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


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.


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


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


# Construct a dataset; simple no regressors
# 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

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

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

# 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?

Please help me, what am I missing?


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.

  • $\begingroup$ Thank you for your answer. Can you extend your answer please and explain the R Code? $\endgroup$ – Ferdi Sep 20 '18 at 12:22
  • $\begingroup$ Thank you for your answer, it helps a lot! Your suggestion works for calculating the likelihood for each level 2 unit. How would you adapt this to calculate the likelihood for each individual (level 1 unit)? I've tried some things but keeping running into incompatible matrices. $\endgroup$ – Niek Sep 21 '18 at 9:50
  • 1
    $\begingroup$ The independent units are the schools (level 2 units), and therefore the likelihood of the whole model is written as the product of the likelihood contributions for each school. AFAIC see this cannot be written as the product of the likelihood contributions of the subjects (level 1 units) because the V matrix is a full covariance matrix (i.e., it is not block-diagonal). $\endgroup$ – Dimitris Rizopoulos Sep 21 '18 at 10:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.