5

我目前正在编写一个脚本来评估用于线性混合模型的（受限）对数似然函数。我需要它来计算模型的可能性，其中一些参数固定为任意值。也许这个脚本对你们中的一些人也有帮助！

我使用lmer()fromlme4和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) 中找到，（希望）可通过此链接访问： 方差分量估计和相关问题的最大似然方法

问候，西蒙

4

1 回答 1

2

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

自 2014 年 3 月 14 日起，CRAN 上的 lme4提供了该开发版本以及此功能，并且参考指南提供了补充包作者 (Ben Bolker) 对原始问题的评论的解释。

于 2014-05-01T16:21:00.483 回答