[R-sig-ME] LogLikelihood

[Andrzej Galecki]

Hello Gianluca,

There are two random effects (q=2).

Matrix D should be 2 by 2, not 6 by 6.

Did not check the rest of your code, but this is an obvious mistake/error.

Best wishes

Andrzej Galecki


On Sun, Jan 25, 2015 at 12:27 PM, bbonit at tin.it <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
>

	[[alternative HTML version deleted]]