Dear r-mixed-models experts,
I am looking at an unbalanced two-way random-effects model, a
specialization of the general linear mixed model,
*y = X*Beta + Z*b + eps, b~N_q(0,Sig_theta), eps~N_n(0,sigsq*I).*
My eventual target is to derive ML estimates for a two-way random-effects
model (not necessarily balanced) with a KNOWN value of sigsq, the error
variance (this is for a theoretical problem which assumes sigsq to be
known).
To this end, at the first stage I am trying to compute the ML estimates for
the case of unknown sigsq by myself (rather than using lmer function) using
the computation suggested in section 1.4 of Douglas Bates's book ("lme4:
Mixed-effects modeling with R").
I am working with the Penicillin data (available with the package), and my
code looks like this:
# penicilin data
data(Penicillin)
str(Penicillin)
summary(Penicillin)
xtabs(~ sample + plate, Penicillin)
(fm03 <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin))
# computing the profiled deviance, d.tilde(theta | y) - up to an additive
constant; case of one observation per cell
y <- Penicillin$diameter
#Obtain the appropriate model matrices, X and Z, for a two-way
random-effects (with non-random intercept)
X <- model.matrix(fm03)
Z <- model.matrix(~ plate + sample, data=Penicillin,
contrasts=list(sample=contrasts(Penicillin$sample,contrasts=F),
plate=contrasts(Penicillin$plate,contrasts=F)))[,-1] #
there must be a more elegant way to obtain X
R <- length(unique(Penicillin$plate)) # number of row levels
C <- length(unique(Penicillin$sample)) # number of column levels
q <- R+C
n <- length(y)
profiled.deviance <- function(theta){
theta.A <- theta[1]
theta.B <- theta[2]
Lambda <- diag( c( rep(sqrt(theta.A),R), rep(sqrt(theta.B),C) ) )
L <- t( chol( t(Lambda) %*% t(Z) %*% Z %*% Lambda + diag(q) ) )
# compute min{rsq(theta, beta, u)} over beta and u, where rsq is the
Penalized Residual Sum of Squares
#translate it to a minimization problem for a single vector x=(beta,u) of
dim q+1: min_x{|y - Ax|^2 + |Gamma x|^2}
A <- cbind(X,Z %*% Lambda)
Gamma <- diag( c(rep(0,1),rep(1,q)) )
argmax <- solve( t(A) %*% A + t(Gamma) %*% Gamma ) %*% t(A) %*% y
min.rsq <- sum( (y - A %*% argmax)^2 ) + sum( (Gamma %*% argmax)^2 )
d.tilde <- 2*log( det(L) ) + n*( 1 + log(2*pi*min.rsq/n) )
return(d.tilde)
}
#plug in the ML estimates from the output of
# fm03 <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin)
lambda.A.ml <- 0.71691
lambda.B.ml <- 3.73092
sigsq.ml <- 0.30242
theta.ml <- c(sqrt(lambda.A.ml /sigsq.ml), sqrt(
lambda.B.ml /sigsq.ml))
profiled.deviance(theta.ml)
The problem is that in the last line of the code, I get a different *
deviance* from the number showing up in the output of lmer(...); it's not
off by much - I get 337.7 while the output shows 332.3 - but I'm guessing
it's still large enough a difference to suggest that my computation is
different than what's implemented in the function lmer.
I wonder what I'm doing wrong, and would greatly appreciate any comments...
Thank you so much - I appreciate your time,
Asaf
[[alternative HTML version deleted]]