# fit the linear mixed effects model

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

Suppose we have the following linear mixed effects model:

How do we fit this nested model in R?

For now, I tried two things:

1. Rewrite the model as:

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.

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

N=20
n=5

#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]
return(result)
}

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){
C[1,1,1]=C[i,2,2]=C[i,3,3]=1
}

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

sigma2 =1

for (i in 1:N){
C[i,2,4]=C[i,3,5]=rnorm(1,mean=0,sd=1)
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)){
X3[j]=i
}
}

#library(lme4)

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)))

• Could you provide a reproducible example in R? – Dimitris Rizopoulos Feb 6 at 7:14
• @DimitrisRizopoulos Thanks for the response, I put the simulation setting in the question just now. – Nan Feb 6 at 15:53
• Thanks but I actually asked if you could put the R code that simulates data and fits the model. – Dimitris Rizopoulos Feb 6 at 16:05
• @DimitrisRizopoulos Sure! Thanks! I will re-edit it. – Nan Feb 6 at 16:18