# Linear Mixed Model

My model in Stan is as follows:

``````data {
int<lower=0> N;         // number of data
int<lower = 0> p_fix;   // number of covariates, fixed effects
int<lower = 0> ngr;	// number of groups
int<lower = 0> p_ran;   // number of covariates, random effects

real Y[N];  	// response vector
matrix[N, p_fix] X;   	// design matrix (fixed effects)
matrix[N, ngr] G;     	// gropus (dummy) allocation ( e.g. riga i-esima: 1 0 0
// vuol dire che ci sono 3 gruppi e l'individuo i appartiene al gruppo 1)

matrix[N,p_ran ] Z;              // design matrix (random effects)
}

parameters {
vector[p_fix] beta;        	// regression coefficients (fixed effects)
vector[ngr] theta;      	  // (group specific) random effects
matrix[p_ran, ngr] gamma;     // regression coefficients (random)

vector[p_fix] sigma2_beta;	        // variances for the prior on beta
vector[ngr] sigma2_theta;	          // variances for the prior on theta
matrix[p_ran, ngr] sigma2_gamma;    // (group specific) random effects

real<lower=0> sigma2_e;    //error sd

}

transformed parameters
{
vector[N] mu;
for(i in 1:N){
mu[i] =  X[i,] * beta +  Z[i,] * (gamma * G[i,]') + theta' * G[i,]';
//mu[i] =  Z[i,] * (gamma * G[i,]') + theta' * G[i,]';
}
}

model {
// Likelihood
for (s in 1:N)
{
// print(" s = ", s);
Y[s] ~ normal(mu[s],pow(sigma2_e,0.5)) ;
}

for (j in 1:p_fix){
beta[j] ~ normal(0.0, pow(sigma2_beta[j], 0.5));
sigma2_beta[j] ~ inv_gamma(3., 5.);
}

for (j in 1:ngr){
theta[j] ~ normal(0.0, pow(sigma2_theta[j], 0.5));
sigma2_theta[j] ~ inv_gamma(3., 5.);
}

for (j in 1:p_ran) {
for(i in 1:ngr){
gamma[j,i] ~ normal(0.0, pow(sigma2_gamma[j,i], 0.5));
sigma2_gamma[j,i] ~ inv_gamma(3., 5.);
}
}

sigma2_e ~ inv_gamma(2., 10.);
}
``````

I run the model with these values and initial conditions:

``````# initial values
beta0 = rep(0.1, p_fix)
theta0=rep(0.1, ngr)
gamma0=matrix(0.1, p_ran, ngr)

sigma2_beta0=rep(10, p_fix)
sigma2_theta0=rep(10,ngr)
sigma2_gamma0=matrix(10, p_ran,ngr)
sigma2_e0=10

data_lmm <-list(N = N,
p_fix=p_fix,
ngr=ngr,
p_ran=p_ran,
Y=as.vector(Y),
X=as.matrix(X),
G=as.matrix(G),
Z=as.matrix(Z)     # <3
)

#initialization of the parameters

inits <- function()
{
list(
beta =beta0,
theta=theta0,
gamma=gamma0,

sigma2_beta=sigma2_beta0,
sigma2_theta=sigma2_theta0,
sigma2_gamma=sigma2_gamma0,
sigma2_e=sigma2_e0)
}

# run stan model
{tic()
LMM <- stan(file = "LinearMixedModel.stan",
data = data_lmm,
chains = 1,
iter = 15000,
warmup = 5000,
thin= 5,
seed = 42,
init = inits,
algorithm = 'NUTS')
toc()}
``````

but during the sampling phase I obtained these following warnings:

``````Chain 1: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1: Exception: normal_lpdf: Scale parameter is nan, but must be > 0!  (in 'modelfacc061586_LinearMixedModel' at line 63)

Chain 1: If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1: but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1:
Chain 1: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1: Exception: normal_lpdf: Scale parameter is nan, but must be > 0!  (in 'modelfacc061586_LinearMixedModel' at line 63)

Chain 1: If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1: but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1:
Chain 1: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1: Exception: normal_lpdf: Scale parameter is nan, but must be > 0!  (in 'modelfacc061586_LinearMixedModel' at line 51)
``````

and so on…
Can someone explain the warning messages?

Thanks all!

Variance parameters need a positivity constraint

``````vector<lower=0>[p_fix] sigma2_beta;	        // variances for the prior on beta
vector<lower=0>[ngr] sigma2_theta;	          // variances for the prior on theta
matrix<lower=0>[p_ran, ngr] sigma2_gamma;    // (group specific) random effects
``````

Without the constraints, Stan will try to calculate the probability that they are negative but when they are negative `pow(..., 0.5)` does not work and returns `nan` (“not a number”).

1 Like

Thanks very much !