[R-sig-ME] LogLikelihood
Andrzej Galecki
agalecki at umich.edu
Sun Jan 25 18:58:13 CET 2015
More precisely D should be positive definite and have 2 by 2 blocks on the
diagonal.
AG
On Sun, Jan 25, 2015 at 12:51 PM, Andrzej Galecki <agalecki at umich.edu>
wrote:
> 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]]
More information about the R-sig-mixed-models
mailing list