My question has to do with how to manually compute the likelihood of a mixed-effect model.
I understand how to determine the likelihood of a fixed-effect model manually.
For example, if I make up some x and y data:
N=100
x = runif(N,0,1); x
y = 5*x + 3 + rnorm(N,0,1); y
Then, I run a simple linear regression:
lm.0 = lm(y~x); summary(lm.0)
I can use dnorm() to compute the probability of each residual:
dnorm(lm.0$residuals,mean(lm.0$residuals), sd(lm.0$residuals))
I can determine the likelihood of the regression by taking the product of all those probabilities:
prod(dnorm(lm.0$residuals,mean(lm.0$residuals), sd(lm.0$residuals)))
I can then determine the negative log-likelihood by taking the negative log of the above:
-log(prod(dnorm(lm.0$residuals,mean(lm.0$residuals), sd(lm.0$residuals))))
Alternatively, I can find the negative sum of the log of the probabilities to determine the negative log likelihood:
-sum(dnorm(lm.0$residuals,mean(lm.0$residuals),sd = sd(lm.0$residuals),log=T))
And this negative log-likelihood that I've now computed in two different ways can be verified using the 'canned' logLik() command:
as.numeric(-logLik(lm.0))
But, my question is, "How can I do the same thing for a mixed-effect model? Namely, how can I determine the likelihood of a mixed-effect model "manually"?
Perhaps, a fake dataset would be useful. Here, I create a fake dataset with 10 groups of 10 values for a variable called, "fish" within each group:
set.seed(45)
fish = runif(100,0,50); fish; hist(fish,col=4)
group = rep(seq(1,10),10); group
fake10 = as.data.frame(cbind(fish,group)); fake10; str(fake10)
fake10$group = as.factor(fake10$group); str(fake10)
fake10 = fake10[order(fake10[,2]),]; head(fake10)
rownames(fake10) = NULL; head(fake10,25)
We could summarize the mean, sd and count for each group if we wanted to:
mean.per.group = tapply(fake10$fish,fake10$group,mean); mean.per.group
sd.per.group = tapply(fake10$fish,fake10$group,sd); sd.per.group
length.per.group = tapply(fake10$fish,fake10$group,length); length.per.group
fake10.pivot = rbind(mean.per.group,sd.per.group,length.per.group); fake10.pivot
I could run a simple fixed-effect linear regression without regard for the grouping variable:
lm.fake10.000= lm(fish~1,fake10); summary(lm.fake10.000)
And compute the negative log-likelihood as:
-sum(dnorm(lm.fake10.000$residuals,mean(lm.fake10.000$residuals),sd(lm.fake10.000$residuals),log=T))
And I could verify that using the 'canned' logLik command:
-logLik(lm.fake10.000)
But, how would one manually compute the negative log likelihood if this data-set were modeled as a mixed-effect model using lmer()?
library(lme4)
lmer.fake10 = lmer(fish ~ 1 + (1|group),fake10,REML=T)
summary(lmer.fake10)
I understand the negative log-likelihood could be computed using the canned logLik() command:
-logLik(lmer.fake10)
But I'd like to know how to compute the negative log-likelihood manually (going back to dnorm() at least).
In researching this question, I've found two sources (stats books) that seem to provide the answer, but alas, I'm unable to decipher the notation in each of them sufficiently to answer my question. Here are the two sources:
Verbeke, G and G. Molenberghs. 2000. Linear Mixed Models for Longitudinal Data. Springer Series in Statistics. Springer. page 45 (equation 5.7)
Zuur, A.F. et al. 2009. Mixed Effects Models and Extensions in Ecology with R. Springer. page 121 (equation at top of the page).
