[R-sig-ME] LogLikelihood
Steve Walker
steve.walker at utoronto.ca
Tue Jan 27 04:26:06 CET 2015
Hope this helps clear things up:
with(getME(mod0, c("n", "L", "X", "beta", "Z", "Lambda", "u", "y")), {
mu <- as.numeric((X %*% beta) + (Z %*% Lambda %*% u))
r2 <- sum((y-mu)^2) + sum(u^2)
ldL2 <- 2*determinant(L, logarithm = TRUE)$modulus
-0.5*(ldL2 + n*(1 + log((2*pi*r2)/n)))
})
logLik(mod0)
Other useful references include the lme4pureR package and the lmer paper:
https://github.com/lme4/lme4pureR/blob/master/R/JSS.R
http://arxiv.org/pdf/1406.5823v1.pdf
Equation 34 of the paper is minus twice the log-likelihood.
Cheers,
Steve
On 2015-01-25 12:27 PM, bbonit at tin.it wrote:
>
>
> Dear list, my name is Gianluca Bonitta
> I'm trying to build up the Loglikelihood of the following model.
> For check it I had used logLik(mod0,REML=F) like "gold standard"
> Like You see there is a difference # diff logLik(mod0,REML=F) - mylog = 0.6339805
> Can somebody help to resolve my mistake ?
> Maybe professor Bolker or professor Bates that are the "fathers" of lme4 pack
> thank You in advance
> Best
> Gianluca
>
> ########################################################################################
> library(lme4)
> data(sleepstudy)
> dat <- sleepstudy[ (sleepstudy$Days %in% 0:4) & (sleepstudy$Subject %in% 331:333) ,]
> colnames(dat) <- c("y", "x", "group")
> mod0 <- lmer( y ~ 1 + x +( x | group ), data = dat,REML="F")
> ########################################################################################
>
> q <- 2 # number of random effects
> n <- nrow(dat) # number of individuals
> m <- length(unique(dat$group)) # number of groups
> Y <- dat$y # response vector
> R <- diag(1,nrow(dat))*summary(mod0)$sigma^2 # covariance matrix of residuals
> beta <- as.numeric(fixef(mod0)) # fixed effects vector (p x 1)
> a<-rep(c(597.1903,60.05023),m) # variance rand effects
> ranef(mod0)$group
> b <-c(17.94432, -3.753130,-33.31148, 10.294328,15.36716, -6.541198) # random effect estimated
> D <-matrix(-0.97,6,6) # random effect estimated correlation
> diag(D) <-a
> X <- cbind(rep(1,n), dat$x) # model matrix of fixed effects (n x p)
> Z.sparse<- getME(mod0,"Z") # model matrix of random effect (sparse format)
> Z <- as.matrix(Z.sparse)
> V <-Z%*% D %*% t(Z) + R # (total) covariance matrix of Y
> # check: values in Y can be perfectly matched using lmer's information
> Y.test <- X %*% beta + Z %*% b + resid(mod0)
> cbind(Y, Y.test)
> mu = X %*% beta + Z %*% b
> ###############################################################################################
> ll = -n/2*log(2*pi) - sum(log(diag(chol(V)))) - .5 * t(Y- mu) %*% chol2inv(chol(V)) %*% (Y-mu);
> logLik(mod0,REML=F)
> ll
> ####################################à
> # diff 'log Lik.' 0.6339805 (df=6)
>
> logLik(mod0,REML=F) -ll
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list