I was suggested by StackOverflow that it's probably better to post this question here.

Suppose we have the following linear mixed effects model: enter image description here

How do we fit this nested model in R?

For now, I tried two things:

  1. Rewrite the model as:

enter image description here

Then using lmer function in lme4 package to fit the mixed effect model and put Xi as both random and fixed effect covariate as: lmer(y ~ X-1+(0+X|subject))

But when I pass the result to BIC and do the model selection, it always picks the simplest model, which is not correct.

  1. I tried to regress y_i on X_i first and treat X_i as the fixed effect, then I will get the estimate of the slope, i.e. phi_i vector. Then see phi_i as the new observations and regress on C_i again to get the beta. But it seems not correct since we do not know C_i in the real problem and it looks like C_i and beta jointly decide the coefficients.

So, are there other ways to fit this kind of model in R and where are my mistakes?

Thanks for any help!

Edit: I put the simulation settings here just in case you want to reproduce it.

enter image description here

Edit again: here are some R code related to my first try.


#equally space point not include the start point
time <- function(from, to, length.out) {
length.out <- length.out + 1
result <- seq(from, to, length.out = length.out)
result <- result[-1]

X = array(0, dim = c(N, n, 3))#each X[i,,] is a nx3 matrix
for (i in 1:N){
  for (j in 1:n){
  X[i,j,] <-c(1,time(0,10,n)[j],(time(0,10,n)[j])^2)

y = array(0, dim = c(N, n, 1))

Omega <- matrix(0,nrow=3,ncol=3)
Omega[1,1] = runif(1,0.01,1.01)#under O1, only omega1^2 is not equal to 0

beta <-rep(0,5)
beta[1]= rnorm(1,mean=0.01,sd=1) #mu0
beta[2]= rnorm(1,mean=0.005,sd=1) #mu1
beta[3]= rnorm(1,mean=0.0025,sd=1) #mu2
#under M1, alpha1=alpha2=0

C = array(0, dim = c(N, 3, 5))
for(i in 1:N){

muy = array(0, dim = c(N, n, 1))

sigma2 =1

for (i in 1:N){
  muy[i,,] <- X[i,,]%*%C[i,,]%*%beta
  Cov[i,,] <- X[i,,]%*%Omega%*%t(X[i,,])+ sigma2*diag(n)# 8x8 matrix
  y[i,,] <- mvrnorm(n = 1, muy[i,,], Cov[i,,], tol = 1e-6, empirical = 
FALSE, EISPACK = FALSE)#install package "MASS"

#just change X into X2, which is in a matrix format, easy for computation later

X2 <- rbind(X[1,,],X[2,,])
for(i in 2:(N-1)){
  X2 = rbind(X2,X[i+1,,])

ym <- as.matrix(y)

#add one column for X, i.e. the subject indicator 
X3 = matrix(0,nrow=N*n,ncol=1)
for(i in 1:N){
    for(j in (n*(i-1)+1):(n*i)){


BIC(lmer(ym ~ X2-1+(0+X2|X3)))
BIC(lmer(ym ~ X2[,1:2]-1+(0+X2|X3)))
BIC(lmer(ym ~ X2[,1]-1+(0+X2|X3)))
  • $\begingroup$ Could you provide a reproducible example in R? $\endgroup$ – Dimitris Rizopoulos Feb 6 at 7:14
  • $\begingroup$ @DimitrisRizopoulos Thanks for the response, I put the simulation setting in the question just now. $\endgroup$ – Nan Feb 6 at 15:53
  • $\begingroup$ Thanks but I actually asked if you could put the R code that simulates data and fits the model. $\endgroup$ – Dimitris Rizopoulos Feb 6 at 16:05
  • $\begingroup$ @DimitrisRizopoulos Sure! Thanks! I will re-edit it. $\endgroup$ – Nan Feb 6 at 16:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.