评估线性混合模型中的似然函数（lme4）答案

【问题标题】：Evaluating the ikelihood function in linear mixed models (lme4)评估线性混合模型中的似然函数（lme4）
【发布时间】：2013-02-13 15:16:48
【问题描述】：

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

我使用来自lme4 和logLik() 的lmer() 来检查我的脚本是否正常工作。看起来，事实并非如此！由于我的教育背景并不真正关心这个级别的统计数据，所以我有点迷失了。

接下来，您将找到一个使用 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

问候，西蒙

【问题讨论】：

我强烈建议你获取lme4的开发版本（可能来自github，通过devtools），它具有返回的能力（mkDevfunOnly=TRUE）偏差函数
谢谢！我查看了lme4 的github 版本并使用devtools 安装了它。是否有关于 devFunOnly=T 参数及其产生的函数的更多文档？我对必须提供给结果偏差函数的参数特别感兴趣，因为这对我来说又是最重要的一步！
当 \code{devFunOnly} 为 \code{TRUE} 时返回的偏差函数采用单个数字向量参数，表示 \code{theta} 向量。该向量在 Cholesky 参数化中定义了随机效应的方差-协方差函数。对于单个随机效应，这是一个等于随机效应标准差的单个值...
... 对于更复杂或多重随机效果，运行 \code{getME(.,"theta")} 以检索拟合模型的 \code{theta} 向量并检查名称向量可能是确定 \code{theta} 向量的元素与随机效应的 Cholesky 因子的下三角形的元素之间的对应关系的最简单方法。（我刚刚将此添加到文档中。它是否有意义，或者您能提出改进建议吗？）
我忘了说 theta 定义了 scaled 方差-协方差矩阵（即相对于残差方差）。

标签： r statistics multi-level lme4 mixed-models

【解决方案1】：

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

自 2014 年 3 月 14 日起，lme4 on CRAN 提供该开发版本以及此功能，reference guide 提供的解释补充了软件包作者 (Ben Bolker) 对原始问题的解释。

【讨论】：