r - 在线性混合模型(lme4)中评估ikelihood函数

我目前正在编写一个脚本，以评估用于线性混合模型的(受限)对数似然函数。我需要它来计算将某些参数固定为任意值的模型的可能性。
也许此脚本对您中的某些人也有帮助!

我使用lmer()和lme4中的logLik()来检查我的脚本是否可以正常工作。看起来，事实并非如此!
由于我的教育背景并不真正关心这一级别的统计信息，所以我迷失了方向。

接下来，您将找到一个使用sleepstudy-data的简短示例脚本:

  # * * * * * * * * * * * * * * * * * * * * * * * *
  # * example data

  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 + ( 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 <- nrow(VarCorr(mod0)$group)                  # number of random effects
  n <- nrow(dat)                                  # number of individuals
  m <- length(unique(dat$group))                  # number of goups
  Y <- dat$y                                      # response vector

  X <- cbind(rep(1,n), dat$x)                     # model matrix of fixed effects (n x p)
  beta <- as.numeric(fixef(mod0))                 # fixed effects vector (p x 1)

  Z.sparse <- t(mod0@Zt)                          # model matrix of random effect (sparse format)
  Z <- as.matrix(Z.sparse)                        # model matrix Z (n x q*m)
  b <- as.matrix(ranef(mod0)$group)               # random effects vector (q*m x 1)

  D <- diag(VarCorr(mod0)$group[1:q,1:q], q*m)    # covariance matrix of random effects
  R <- diag(1,nrow(dat))*summary(mod0)@sigma^2    # covariance matrix of residuals
  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)

  # * * * * * * * * * * * * * * * * * * * * * * * *
  # * the likelihood function

  # profile and restricted log-likelihood (Harville, 1997)
  loglik.p <- - (0.5) * (  (log(det(V))) + t((Y - X %*% beta)) %*% solve(V) %*% (Y - X %*% beta)  )
  loglik.r <- 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 <- logLik(mod0, REML=TRUE)
  cbind(loglik.p, loglik.r, loglik.lmer)

也许这里有一些LMM专家可以为您提供帮助？无论如何，您的建议将不胜感激!

编辑:顺便说一句，LMM的似然函数可以在Harville(1977)中找到，(希望如此)可通过以下链接访问:
Maximum likelihood approaches to variance component estimation and to related problems

问候，
西蒙

最佳答案

解决方案(截至2013年3月)是安装lme4的开发版本并使用devFunOnly参数。

自2014年3月14日起，该开发版本以及此功能随lme4 on CRAN一起提供，并且reference guide给出的解释使软件包作者(Ben Bolker)的注释对原始问题的看法更加恭维。