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.
