Why is my manual calculation of the log-likelihood for a 3-level model different than what nlme provides?

Question

In short: I want to manually calculate the log-likelihood of a 3-level multilevel/mixed/hierarchical model, but my result is different from what nlme gives. I don't understand why. Examples of the code and formulae below

---

Background

I am trying to identify influential cases in a mixed-model. I came across a method proposed by Coffman & Millsap (2006) who suggest looking into the individual contribution to the overall log-likelihood. For various reasons I find this method preferable to other methods that I’ve found and after some searching I came across a function that implements this method in R, thelogLik1 function in the nlmeU package.

The function works great for a two-level random effect model. However, it has not yet been extended to accommodate further levels of nesting. Now, my dataset contains longitudinal measures of students within schools, so I would love to apply this method to a three-level model. I’ve looked into the function and the book by Galecki & Burzykowski (2013) that describes how logLik1 utilizes the likelihood function.

The paper by Longford (1987) describes an extension of the log-likelihood function for $h = 2, ..., m$ levels. The likelihood function is described as follows:

$$\lambda = -\frac{N}{2}\log (2\pi) - \frac{1}{2} \sum_{i=1}^{N} \log [\det(\mathbf{V})] - \frac{1}{2} (\mathbf{y} - \mathbf{X \beta})'\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X \beta}) $$

Were variance of $\mathbf{y}$ is denoted:

$$\mathbf{V} = \sigma_e^2I_N + \sum_{l=2}^{h}\sum_{i=1}^{n_i}\mathbf{Z}_{li}\Theta_l \mathbf{Z}_{li}' $$

and $\Theta_l$ is the variance-covariance matrix at level l and $\mathbf{Z}_{li}$ the matrix of random-effect regressors. An adaption to the logLik1 function should be as simple as adding the matrix Z and D to the sum that calculates V. I've tried to make these adaptions but the sum of the individual log-likelihoods does not correspond to the overall log-likelihood given by nlme.

Function

The function is changed by adding D2 and Z2 to the original code:

# Modified logLik1() function (simplified for readability)
logLik2.lme <- function(modfit, dt1){ 
  m            <- modfit$modelStruct                 # Model structure, this gets an additional parameter for the top level variance

  sigma        <- modfit$sigma                       # sigma, residual standard deviation
  D            <- as.matrix(m$reStruct[[1]])         # "subject" ;    as.matrix(m$reStruct[[2]]) will get you the top level intercept variance
  D            <- D  * sigma^2                       # Matrix D, design matrix of random effects; VarCorr(modfit)[,c(1,3)] * sigma^2

  D2  <- as.matrix(m$reStruct[[2]]) * sigma^2        # Matrix D2, variance-covariance matrices for random effects at level 2 

  vecR         <- rep(sigma, nrow(dt1))
  vecR2        <- vecR^2
  R            <- diag(vecR2, nrow=length(vecR))        # R_i matrix   Residual variance matrix

  # For now I constructed Z and Z2 to handle a two-level random intercept model, this is easily extended.
  n            <- nrow(dt1)                 # No. of obs for subject dt1

  Z            <- model.matrix(m$reStruc, data=dt1)[,1]     # Z_i matrix gets an additional column of 1's to indicate top level intercept, take first column
  Z2           <- model.matrix(m$reStruc, data=dt1)[,2]     # Z2 second column of Z_i matrix to indicate top level intercept

  V            <- R + Z %*% D %*% t(Z) +  Z2 %*% D2 %*% t(Z2)   # V_i matrix; changes with an addition level, because m changes, D changes and Z
  #J <- matrix(1,ncol = n, nrow = n); V <- D[1]*J+ D2[1]*J  + R   # simplified Longford (1987) 2.4 

  predY        <- predict(modfit, dt1, level=0)         # Predict fixed level = 0 (Y - XB); These are an individual's residuals without random effects
  #predict(modfit, dt1, level=0)  + ranef(modfit)[1,1] + ranef(modfit)[1,2] * c(0,2:7) == predict(modfit, dt1, level=1)

  dvName       <- as.character(formula(modfit)[[2]]) # Dependent variable name
  r            <- dt1[[dvName]] - predY              # Residuals; observed - predicted


  lLik         <- n*log(2*pi) + log(det(V)) +  t(r) %*% solve(V) %*% r # page 256, formula 13.27 Galecki & Burzykowski (2013)
  return(-0.5 * as.numeric(lLik)) # Divide everything by 2
}

Example

An example shows that the likelihood remains unchanged (even though nlme says the likelihood is increased in the second model):

require(nlmeU)

# Construct a dataset; simple no regressors
set.seed(1010101)
# id's
school.id = rep(1:10, each = 70)
stud.id = rep(1:100, each = 7)
time <- rep(1:7,100)
# random effects
school <- rnorm(10,0,.5)
stud <- rnorm(100,0,1)
e <- rnorm(700,0,.2)
# dependend variable
math <- school[school.id] + stud[stud.id] + e
df1<-as.data.frame(cbind(school.id,stud.id,time,math))

# Create 2 level model
mod2 <- lme(math~ 1,random= list(stud.id = ~1),data=df1,method="ML",na.action=na.exclude)
summary(mod2)

# Create 3 level model
mod3 <- lme(math~ 1,random= list(school.id = ~ 1, stud.id = ~1),data=df1,method="ML",na.action=na.exclude)
summary(mod3)

# Apply the old function to the two-level model and the new function to the three-level model
lLik.i0 <- by(df1,df1$stud.id, function(dfi) logLik1(mod2,dfi))
lLik.i0 <- as.vector(lLik.i0)

lLik.i <- by(df1,df1$stud.id, function(dfi) logLik2.lme(mod3,dfi))
lLik.i <- as.vector(lLik.i)

# This should be the same?
sum(lLik.i0);logLik(mod2)
sum(lLik.i);logLik(mod3)

Please help me, what am I missing?

Answer 1

The problem stems from the fact that the Z matrix should be reconstructed to account for the nested random effects when computing the V matrix. In this particular example, it should be

Z <- kronecker(diag(10), rep(1, 7))
DD <- matrix(c(D2), 10, 10)
diag(DD) <- diag(DD) + c(D)
V <- Z %*% DD %*% t(Z) + R

Also you will need

lLik.i <- by(df1, df1$school.id, function(dfi) logLik2.lme(mod3, dfi))

Here I've hardcoded it but you could program it to work in the general case.