title: Evaluating the ikelihood function in linear mixed models (lme4)
body: <p>I am currently writing a script to evaluate the (restricted) log-likelihood function for use in linear mixed models. I need it to calculate the likelihood of a model with some parameters fixed to arbitrary values.
  Maybe this script is helpful to some of you as well! </p>
  
  <p>I use <code>lmer()</code> from <code>lme4</code> and <code>logLik()</code> to check whether my script works as it should. And as it seems , it does not! 
  As my educational background wasn't really concerned with this level of statistics, I am a bit lost finding the mistake.</p>
  
  <p>Following, you will find a short example  script using the sleepstudy-data:</p>
  
  <pre><code>  # * * * * * * * * * * * * * * * * * * * * * * * * 
    # * example data
  
    library(lme4)
    data(sleepstudy)
    dat &lt;- sleepstudy[ (sleepstudy$Days %in% 0:4) &amp; (sleepstudy$Subject %in% 331:333) ,]
    colnames(dat) &lt;- c("y", "x", "group")
  
    mod0 &lt;- lmer( y ~ 1 + x + ( 1 | group ), data = dat)  
  
  
    # + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + #
    #                                                             #
    #   Evaluating the likelihood-function for a LMM              #
    #   specified as: Y = X*beta + Z*b + e                        #
    #                                                             #
    # + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 
  
    # * * * * * * * * * * * * * * * * * * * * * * * * 
    # * the model parameters
  
    # n is total number of individuals
    # m is total number of groups, indexed by i
    # p is number of fixed effects
    # q is number of random effects
  
    q &lt;- nrow(VarCorr(mod0)$group)                  # number of random effects
    n &lt;- nrow(dat)                                  # number of individuals
    m &lt;- length(unique(dat$group))                  # number of goups
    Y &lt;- dat$y                                      # response vector
  
    X &lt;- cbind(rep(1,n), dat$x)                     # model matrix of fixed effects (n x p)
    beta &lt;- as.numeric(fixef(mod0))                 # fixed effects vector (p x 1)
  
    Z.sparse &lt;- t(mod0@Zt)                          # model matrix of random effect (sparse format)
    Z &lt;- as.matrix(Z.sparse)                        # model matrix Z (n x q*m)
    b &lt;- as.matrix(ranef(mod0)$group)               # random effects vector (q*m x 1)
  
    D &lt;- diag(VarCorr(mod0)$group[1:q,1:q], q*m)    # covariance matrix of random effects
    R &lt;- diag(1,nrow(dat))*summary(mod0)@sigma^2    # covariance matrix of residuals
    V &lt;- Z %*% D %*% t(Z) + R                       # (total) covariance matrix of Y
  
    # check: values in Y can be perfectly matched using lmer's information
    Y.test &lt;- X %*% beta + Z %*% b + resid(mod0)
    cbind(Y, Y.test)
  
    # * * * * * * * * * * * * * * * * * * * * * * * * 
    # * the likelihood function
  
    # profile and restricted log-likelihood (Harville, 1997)
    loglik.p &lt;- - (0.5) * (  (log(det(V))) + t((Y - X %*% beta)) %*% solve(V) %*% (Y - X %*% beta)  )
    loglik.r &lt;- loglik.p - (0.5) * log(det( t(X) %*% solve(V) %*% X ))
  
    #check: value of above function does not match the generic (restricted) log-likelihood of the mer-class object
    loglik.lmer &lt;- logLik(mod0, REML=TRUE)
    cbind(loglik.p, loglik.r, loglik.lmer)
  </code></pre>
  
  <p>Maybe there are some LMM-experts here who can help? In any case your recommendations are greatly appreciated!</p>
  
  <p>edit: BTW, the likelihood function for LMMs can be found in Harville (1977), (hopefully) accessible through this link:
  <a href="http://www.jstor.org/stable/10.2307/2286796">Maximum likelihood approaches to variance component estimation and to related problems</a></p>
  
  <p>Regards, 
  Simon</p>