Converting Scaled Parameters to Original Scale


#1

I have the polynomial regression model

Y_i | \mu, \sigma^2 \sim \mathcal{N}(\mu_i, \sigma^2), i = 1, \dots, n \ \text{independent}

\mu_i = \alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x^2_{i1} + \beta_4 x^2_{i2} + \beta_5 x_{i1} x_{i2}

With appropriate priors on \alpha, \beta_i, and \sigma^2.

My model is as follows:

```{stan output.var="example"}
  data{
    int<lower=1> n;  // number of observations
    int<lower=1> p;  // number of predictors
    matrix[n, p] x;  // design matrix
    vector[n] y;     // responses
  }
  
  parameters{
    real<lower=-10, upper=10> log_sigma;
    real alpha;      // intercept
    vector[p] beta;  // regression coefficients for predictors
  }
  
  transformed parameters{
    real<lower=0> sigma;
    
    sigma = exp(log_sigma);
  }
  
  model{
     y  ~ normal(alpha + x * beta, sigma);  // multiplies matrix x with vector beta 
                                            // to obtain vector of mean values
                                            // the normal() function is vectorised
     alpha ~ normal(0, 1000);
     beta ~ normal(0, 1000);                // vectorised prior
  }
```

And I use the following R code to scale my parameters:

```{r, results='hide'}
library(MASS)
hills[18,3] <- 18.650
X <- with(hills, cbind(dist, climb, climb^2))
x.means <- colMeans(X)
x.sds   <- apply(X, 2, sd)
X.std <- sweep(X, 2, x.means, "-")
X.std <- sweep(X.std, 2, x.sds, "/")
data.in <- list(x=X.std, y=hills$time, n=NROW(X.std), p=NCOL(X.std))
model.fit1 <- sampling(example, data=data.in)
```

Now, the formulae to convert from the scaled parameters back to the original scale is as follows:

Since I’m using Stan, the implementation of this formula to transform the parameters back to the original scale would look like this:

 for(i in 1:n){
        mu[i] = alpha_std + beta1_std*x1_std[i] + beta2_std*x2_std[i] + beta3_std*x1_std[i]^2 + beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
      
        alpha = -beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2 + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);
        beta1 = beta1_std/x1_std[i] - 2*beta3_std*bar_x1/x1_std[i]^2 - beta5*bar_x2/(x1_std[i]*x2_std[i]);
        beta2 = beta2_std/x2_std[i] - 2*beta4_std*bar_x2 - beta5_std*bar_x1/(x1_std[i]*x2_std[i]);
        beta3 = beta3_std/x1_std[i]^2;
        beta4 = beta4_std/x2_std[i]^2;
        beta5 = beta5_std/(x1_std[i]*x2_std[i]);
      }

My original code looked like this, but I was having too many problems with the way Stan deals with vectors and exponents, so I abandoned it:

```{stan output.var="example"}
    data{
      int<lower=1> n;
      vector[n] x1;
      vector[n] x2;
      vector[n] y;
    }
    
    transformed data{
      real bar_x1;
      real x1_sd;
      vector[n] x1_std;
      real bar_x2;
      real x2_sd;
      vector[n] x2_std;
      
      bar_x1 = mean(x1);
      x1_sd = sd(x1);
      x1_std = (x1 - bar_x1)/x1_sd; // centered and scaled
      
      bar_x2 = mean(x2);
      x2_sd = sd(x2);
      x2_std = (x2 - bar_x2)/x2_sd; // centered and scaled
    }
    
    parameters{
      real<lower=-10, upper=10> log_sigma;
      real alpha_std;
      real beta1_std;
      real beta2_std;
      real beta3_std;
      real beta4_std;
      real beta5_std;
    }
    
    transformed parameters {
      real<lower=0> sigma;
      vector[n] mu;
      real alpha;
      real beta1;
      real beta2;
      real beta3;
      real beta4;
      real beta5;

      sigma = exp(log_sigma);
      
      for(i in 1:n){
        mu[i] = alpha_std + beta1_std*x1_std[i] + beta2_std*x2_std[i] + beta3_std*x1_std[i]^2 + beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
      
        alpha = -beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2 + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);
        beta1 = beta1_std/x1_std[i] - 2*beta3_std*bar_x1/x1_std[i]^2 - beta5*bar_x2/(x1_std[i]*x2_std[i]);
        beta2 = beta2_std/x2_std[i] - 2*beta4_std*bar_x2 - beta5_std*bar_x1/(x1_std[i]*x2_std[i]);
        beta3 = beta3_std/x1_std[i]^2;
        beta4 = beta4_std/x2_std[i]^2;
        beta5 = beta5_std/(x1_std[i]*x2_std[i]);
      }
    }
    
    model{
      y ~ normal(mu, sigma);
      alpha ~ normal(0, 1000);
      beta1 ~ normal(0, 1000);
      beta2 ~ normal(0, 1000);
      beta3 ~ normal(0, 1000);
      beta4 ~ normal(0, 1000);
      beta5 ~ normal(0, 1000);
    }
```

Which is why I started using the design matrix methodology.

But now I can’t figure out how to implement this rescaling formula in my new model with the design matrix! In other words, I don’t see how I can include the implementation

mu[i] = alpha_std + beta1_std*x1_std[i] + beta2_std*x2_std[i] + beta3_std*x1_std[i]^2 + beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
          
alpha = -beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2 + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);
beta1 = beta1_std/x1_std[i] - 2*beta3_std*bar_x1/x1_std[i]^2 - beta5*bar_x2/(x1_std[i]*x2_std[i]);
beta2 = beta2_std/x2_std[i] - 2*beta4_std*bar_x2 - beta5_std*bar_x1/(x1_std[i]*x2_std[i]);
beta3 = beta3_std/x1_std[i]^2;
beta4 = beta4_std/x2_std[i]^2;
beta5 = beta5_std/(x1_std[i]*x2_std[i]);

in the new, design matrix model. That is, implement it in the design matrix model in order to convert the scaled values back to the original scale.

I would greatly appreciate it if someone could please take the time to help me with this.


#2

Hi y64c,

Thanks for posting this follow up from your earlier since deleted post.

I haven’t checked your math, as I’d like to address some other points first. Everything below references your original Stan program, which sounds like it should be fine.

  1. I believe Stan is confused about your transformed parameters. Your conversion back to the original scale probably should happen in the generated quantities block, instead of the transformed parameters block since you aren’t using the original scale parameters in the model. I don’t yet have the detailed knowledge to say exactly what Stan is trying to do with the conversion within the transformed parameters block.

  2. You want to put priors on the parameters, not the transformed parameters.

  3. The going recommendation is to put a prior on sigma itself, not a transformed version of it, and not an arbitrarily bounded version (transformed or not) of it.

These recommendations leave you with Stan code like

data{
  int<lower=1> n;
  vector[n] x1;
  vector[n] x2;
  vector[n] y;
}
transformed data{
  real bar_x1;
  real x1_sd;
  vector[n] x1_std;
  real bar_x2;
  real x2_sd;
  vector[n] x2_std;
  real y_sd;

  bar_x1 = mean(x1);
  x1_sd = sd(x1);
  x1_std = (x1 - bar_x1)/x1_sd; // centered and scaled

  bar_x2 = mean(x2);
  x2_sd = sd(x2);
  x2_std = (x2 - bar_x2)/x2_sd; // centered and scaled

  y_sd = sd(y);
}
parameters{
  real<lower=0> sigma;
  real alpha_std;
  real beta1_std;
  real beta2_std;
  real beta3_std;
  real beta4_std;
  real beta5_std;
}
transformed parameters {
  real mu[n];

  for(i in 1:n) {
    mu[i] = alpha_std + beta1_std*x1_std[i]
      + beta2_std*x2_std[i] + beta3_std*x1_std[i]^2
      + beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
  }
}
model{
  alpha_std ~ normal(0, 10);
  beta1_std ~ normal(0, 2.5);
  beta2_std ~ normal(0, 2.5);
  beta3_std ~ normal(0, 2.5);
  beta4_std ~ normal(0, 2.5);
  beta5_std ~ normal(0, 2.5);
  sigma ~ exponential(1 / y_sd);

  y ~ normal(mu, sigma);
}
generated quantities {
  real alpha;
  real beta1;
  real beta2;
  real beta3;
  real beta4;
  real beta5;

  for (i in 1:n) {
    alpha = -beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd
      + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2
      + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);

    beta1 = beta1_std/x1_std[i] - 2*beta3_std*bar_x1/x1_std[i]^2
      - beta5*bar_x2/(x1_std[i]*x2_std[i]);

    beta2 = beta2_std/x2_std[i] - 2*beta4_std*bar_x2
      - beta5_std*bar_x1/(x1_std[i]*x2_std[i]);

    beta3 = beta3_std/x1_std[i]^2;

    beta4 = beta4_std/x2_std[i]^2;

    beta5 = beta5_std/(x1_std[i]*x2_std[i]);
  }
}

The above code runs but, I think the conversion back to the original scale idea might be different from what you are hoping for. It’s my understanding that Bayesian linear regression will not recover the same values under these two parameterizations, even after attempting to convert things back to the original scale. If this idea is correct, it would be best then to use auto-scaled default priors, as per rstanarm: https://cran.r-project.org/web/packages/rstanarm/vignettes/priors.html

Hope that helps. Cheers.


#3

If you are happy to transform the parameters in R, one could always extract() the samples and transform them back in R as well. Otherwise, returning the coefficients of the design matrix to their original scale is probably best done with a loop (or even two) in the generated quantities block. (Edit: I see @roualdes beat me to the punch here).

This is entirely correct, and the \text{log}(\sigma) \sim \text{U}[-10, 10] is very much a hangover from developing the Stan code alongside the BUGS code.

N.B. for anyone else answering this question, This is an assessment for a unit I TA for.


#4

@roualdes
@hhau

Thank you both for the great advice.

@hhau
I didn’t even notice that the BUGS document demonstrated the transformation in R. Instead, I was looking at the way it was done in the Stan code and trying to replicate that.

It’s great that you’re on here. Any chance I can email you in the future if anything is unclear?


#5

You should be using the QR reparameterization anytime you are tempted to use a polynomial regression.

Normal with a standard deviation of 1000 is unappropriate.


#6

Yes, but reply speed may vary considerably!

Upon rereading your assignment I see that you have to do the (un)scaling inside Stan; so yes, best to do so in the generated quantities block


#7

@hhau What are you reading that indicates that? I can’t see anything that explicitly says it must be in Stan?

@bgoodri As I understand it, QR reparameterization increases the efficiency of posterior sampling, which is also the purpose of parameter standardisation. So if we use QR reparameterization, then we do not get scaled parameter values. For this task, however, I need the scaled parameter values, so It seems that I have to go without the QR reparameterization.


#8

In the the QR reparameterization, the columns of Q all have the same scale, in addition to being orthogonal to each other. Both of these are critical to the ability to sample efficiently, even if the columns of X are highly correlated as they often are when doing polynomial regression. If you feel that you have to have coefficients that correspond to scaled columns of X, do the QR reparameterization, invert the transformation at the end to get unscaled coefficients, and then scale them however you want.


#9

@roualdes Ok, so I found numerous errors in my previous implementation of the transformation, and after fixing it, I got this:

alpha = alpha_std - beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd
      + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2
      + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);

    beta1 = beta1_std/x1_sd - 2*beta3_std*bar_x1/x1_sd^2
      - beta5*bar_x2/(x1_sd*x2_sd);

    beta2 = beta2_std/x2_sd - 2*beta4_std*bar_x2/x2_sd^2
      - beta5_std*bar_x1/(x1_sd*x2_sd);

    beta3 = beta3_std/x1_sd^2;

    beta4 = beta4_std/x2_sd^2;

    beta5 = beta5_std/(x1_sd*x2_sd); 

Notice that the indices are gone. This is because, likely due to lack of sleep, I mixed up the x1_sd and x2_sd variables with the vectors x1_std and x2_std. So I’ve gone ahead and replaced all of the x1_std and x2_std with x1_sd and x2_sd, which also eliminates the need for the for() loop, since we no longer have indices.

However, I’m unsure if I should be getting rid of the for() loop in the generated quantities block now?

The issue is, if I get rid of the for loop, then my code is as follows:

data{
  int<lower=1> n;
  vector[n] x1;
  vector[n] x2;
  vector[n] y;
}
transformed data{
  real bar_x1;
  real x1_sd;
  vector[n] x1_std;
  real bar_x2;
  real x2_sd;
  vector[n] x2_std;
  real y_sd;

  bar_x1 = mean(x1);
  x1_sd = sd(x1);
  x1_std = (x1 - bar_x1)/x1_sd; // centered and scaled

  bar_x2 = mean(x2);
  x2_sd = sd(x2);
  x2_std = (x2 - bar_x2)/x2_sd; // centered and scaled

  y_sd = sd(y);
}
parameters{
  real<lower=0> sigma;
  real alpha_std;
  real beta1_std;
  real beta2_std;
  real beta3_std;
  real beta4_std;
  real beta5_std;
}
transformed parameters {
  real mu[n];

  for(i in 1:n) {
    mu[i] = alpha_std + beta1_std*x1_std[i]
      + beta2_std*x2_std[i] + beta3_std*x1_std[i]^2
      + beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
  }
}
model{
  alpha_std ~ normal(0, 10);
  beta1_std ~ normal(0, 2.5);
  beta2_std ~ normal(0, 2.5);
  beta3_std ~ normal(0, 2.5);
  beta4_std ~ normal(0, 2.5);
  beta5_std ~ normal(0, 2.5);
  sigma ~ exponential(1 / y_sd);

  y ~ normal(mu, sigma);
}
generated quantities {
  real alpha;
  real beta1;
  real beta2;
  real beta3;
  real beta4;
  real beta5;
  
    alpha = alpha_std - beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd
      + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2
      + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);

    beta1 = beta1_std/x1_sd - 2*beta3_std*bar_x1/x1_sd^2
      - beta5*bar_x2/(x1_sd*x2_sd);

    beta2 = beta2_std/x2_sd - 2*beta4_std*bar_x2/x2_sd^2
      - beta5_std*bar_x1/(x1_sd*x2_sd);

    beta3 = beta3_std/x1_sd^2;

    beta4 = beta4_std/x2_sd^2;

    beta5 = beta5_std/(x1_sd*x2_sd);
}

And, when I run the model with the data

x1 <- hills$dist
x2 <- hills$climb
y <- hills$time
n <- length(x1)
data.in <- list(x1 = x1, x2 = x2, y = y, n = n)
model.fit <- sampling(MLRQR, data.in)

, I get the following error:

starting worker pid=85884 on localhost:11196 at 22:00:33.126
starting worker pid=85894 on localhost:11196 at 22:00:33.444
starting worker pid=85904 on localhost:11196 at 22:00:33.734
starting worker pid=85914 on localhost:11196 at 22:00:34.012

SAMPLING FOR MODEL 'stan-eae9207ce7d3' NOW (CHAIN 1).

Gradient evaluation took 4.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.48 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.327414 seconds (Warm-up)
               0.235446 seconds (Sampling)
               0.56286 seconds (Total)


SAMPLING FOR MODEL 'stan-eae9207ce7d3' NOW (CHAIN 2).

Gradient evaluation took 7.1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.71 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.30663 seconds (Warm-up)
               0.220162 seconds (Sampling)
               0.526792 seconds (Total)


SAMPLING FOR MODEL 'stan-eae9207ce7d3' NOW (CHAIN 3).

Gradient evaluation took 2.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.231989 seconds (Warm-up)
               0.180229 seconds (Sampling)
               0.412218 seconds (Total)


SAMPLING FOR MODEL 'stan-eae9207ce7d3' NOW (CHAIN 4).

Gradient evaluation took 6.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.225621 seconds (Warm-up)
               0.153597 seconds (Sampling)
               0.379218 seconds (Total)

 Show Traceback
Error in checkForRemoteErrors(val) : 4 nodes produced errors; first error: invalid class “stanfit” object: The following variables have undefined values: beta1

However, If I keep the for() loop in the generated quantities block, despite the fact that there no longer are any indices in the generated quantities block, the code seems to work. However, I have no idea whether this is the correct way to go about things; in other words, I have no idea whether running the for() loop still, despite the fact that there are no longer any indices, is the correct way to go about this?

The code with the for() loop in the generated quantities block, despite the fact that there are no longer any indices in the transformation formula, is as follows:

data{
  int<lower=1> n;
  vector[n] x1;
  vector[n] x2;
  vector[n] y;
}
transformed data{
  real bar_x1;
  real x1_sd;
  vector[n] x1_std;
  real bar_x2;
  real x2_sd;
  vector[n] x2_std;
  real y_sd;

  bar_x1 = mean(x1);
  x1_sd = sd(x1);
  x1_std = (x1 - bar_x1)/x1_sd; // centered and scaled

  bar_x2 = mean(x2);
  x2_sd = sd(x2);
  x2_std = (x2 - bar_x2)/x2_sd; // centered and scaled

  y_sd = sd(y);
}
parameters{
  real<lower=0> sigma;
  real alpha_std;
  real beta1_std;
  real beta2_std;
  real beta3_std;
  real beta4_std;
  real beta5_std;
}
transformed parameters {
  real mu[n];

  for(i in 1:n) {
    mu[i] = alpha_std + beta1_std*x1_std[i]
      + beta2_std*x2_std[i] + beta3_std*x1_std[i]^2
      + beta4_std*x2_std[i]^2 + beta5_std*x1_std[i]*x2_std[i];
  }
}
model{
  alpha_std ~ normal(0, 10);
  beta1_std ~ normal(0, 2.5);
  beta2_std ~ normal(0, 2.5);
  beta3_std ~ normal(0, 2.5);
  beta4_std ~ normal(0, 2.5);
  beta5_std ~ normal(0, 2.5);
  sigma ~ exponential(1 / y_sd);

  y ~ normal(mu, sigma);
}
generated quantities {
  real alpha;
  real beta1;
  real beta2;
  real beta3;
  real beta4;
  real beta5;
  
  for (i in 1:n) {
    alpha = alpha_std - beta1_std*bar_x1/x1_sd - beta2_std*bar_x2/x2_sd
      + (beta3_std*bar_x1^2)/x1_sd^2 + (beta4_std*bar_x2^2)/x2_sd^2
      + (beta5_std*bar_x2*bar_x1)/(x1_sd*x2_sd);

    beta1 = beta1_std/x1_sd - 2*beta3_std*bar_x1/x1_sd^2
      - beta5*bar_x2/(x1_sd*x2_sd);

    beta2 = beta2_std/x2_sd - 2*beta4_std*bar_x2/x2_sd^2
      - beta5_std*bar_x1/(x1_sd*x2_sd);

    beta3 = beta3_std/x1_sd^2;

    beta4 = beta4_std/x2_sd^2;

    beta5 = beta5_std/(x1_sd*x2_sd);
  }
}

And it seems to run fine:

library(MASS)
hills[18, 3] <- 18.65
x1 <- hills$dist
x2 <- hills$climb
y <- hills$time
n <- length(x1)
data.in <- list(x1 = x1, x2 = x2, y = y, n = n)
model.fit <- sampling(MLRQR, data.in)

print(model.fit, pars = c("alpha_std", "alpha", "beta1_std", "beta2_std", "beta3_std", "beta4_std", "beta5_std", "beta1", "beta2", "beta3", "beta4", "beta5", "sigma"), probs = c(0.05, 0.5, 0.95), digits = 5)

So

(1) either I should be leaving the for() loop in the generated quantities block, despite the fact that I have have gotten rid of the vector indexation, or

(2) I should be now be getting rid of the for() loop in the generated quantities block, since I have now gotten rid of the vector indexation, but, in which case, there seems to be an error.

Personally, despite the fact the code with the for() loop remaining in the generated quantities block works, it’s making me nervous, since I don’t understand why it would have to remain after we’ve gotten rid of any need for vector operations. But this could very well just be because I’m misunderstanding how the generated quantities block works …

I would greatly appreciate it if people could please clarify this.


#10

It tells you:

which is because

uses beta5 before it is defined. You may have meant beta5_std.


#11

Thanks for the clarification. I’ll look into the QR reparameterization after I’ve finished doing it this (easier) way. I don’t have a lot of experience with Stan nor Bayesian statistics in general, so I’ll start with the easier (and less correct) way and, after I understand that, try it with the QR reparameterization.


#12

Oh my goodness, I apologise for making such a stupid mistake. You are, of course, correct, and It seems that the code is working fine now. Thank you for the assistance.

What do you think of the output? The fact that beta4 is 0.00000 for everything makes me suspect that I’ve made an error. Or is this actually fine?


#13

I would suspect error


#14

Hmm, ok. I’ll leave this thread as is and start a new one for a review of the mathematics and outputs. Thank you all for the input.