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

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? 1 Answer 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.