Skip to content

LogLikelihood

4 messages · bbonit at tin.it, Andrzej Galecki, Steve Walker

Original
bbonit at tin.it
# 
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) &amp; (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
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:

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Andrzej Galecki
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    
1 day later

    

      
      

        


  

        

      Steve Walker
    

    

      
      #
    

  


  

      
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: