Linear Mixed Model

Modeling
1

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!

2

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
3

Thanks very much !
Your advice results useful to solve problems.