Dear R,

 

Using Henderson's equation, I can estimate the variance components for
models with a random intercept (see below), but I have not been able to
figure out how to calculate the variance estimates for a model with a random
slope.  Using a very small toy data set, I have the following code for one
iteration only for an intercept only model.  Can someone can tell me how to
modify it so that it yields a random slopes model?

 

 

 

 

#Data and Packages

library(lme4)

library(MASS)  #for generalized inverse

z=c(1,1,1,1,2,2,2,2,3,3,3,3)

x=c(-.5,-.5,.5,.5,-.5,-.5,.5,.5,-.5,-.5,.5,.5)

y=c(1, 2, 3.5, 4, 3.5, 4.5, 5.5, 6, 6.5, 7.5, 8.5, 9)

 

#Fixed Predictor with Random Intercept

mix.1 <-lmer(y ~ 1 + (1|z))

summary(mix.1)

 

#Get matrices

X= getME(mix.1,"X")

Z = as.matrix(getME(mix.1,"Z"))

 

#Calculate Preliminary Matrices

yy=t(y)%*%y

Xy=t(X)%*%y

Zy=t(Z)%*%y

XX <-t(X)%*%X

XZ <-t(X)%*%Z

ZX <-t(Z)%*%X

ZZ <-t(Z)%*%Z

D  <-diag(1,ncol(Z))                                                

N <-length(y)

k=ncol(X)

 

#initial values

sigma <-c(7,.5)                                           

 

#Begin iteration here  -- #D needs to become a block diagonal matrix for
random coefficients

D <-sigma[1]*D                                                      

R <-sigma[2]*diag(1,N)                                          

Im=diag(1,ncol(Z)) 

V=sigma[2]*Im+ZZ%*%D                                       

W=(solve(sigma[2]*(Im+ZZ%*%D))*sigma[2] )


XZD <-t(X) %*% Z %*% D

A= rbind(cbind(XX, XZD),cbind(ZX, V));A              

ab=rbind(Xy,Zy)  

BB=ginv(A)%*%ab  

r=nrow(BB)                                              

B=BB[1:k];B


v=cbind(BB[(k+1):r])

u=D%*%v;u                                                                 

 

 

#REML FOR NEXT ITERATION

rx=ncol(XX)

ZMZ <-ZZ-t(XZ)%*%ginv(XX)%*%XZ                            

T <-solve(Im+(ZMZ%*%D/sigma[2]))                          

trT <-sum(diag(T))

q=ncol(Z)

sig.1 <-(t(u)%*%u)/(q-trT)

sig.2 <-(yy-(t(B)%*%Xy)-(t(u)%*%Zy))/(N-rx)

sig.REML <-c(sig.1,sig.2);sig.REML

 

#Substitute returned values for sigma and perform second iteration

 

 

 


	[[alternative HTML version deleted]]

