我目前正在编写一个脚本，用于评估线性混合模型中（受限制的）对数似然函数。我需要它来计算一些参数固定为任意值的模型的似然度。
也许这个脚本对你们中的一些人有帮助！
我使用lme4中的lmer()和logLik()来检查我的脚本是否按照预期工作。但似乎并没有！ 由于我的教育背景并不真正关注这个统计水平，所以我有点迷失在找错中。
接下来，您将找到一个使用sleepstudy数据的简短示例脚本：
也许这里有一些LMM专家可以提供帮助？无论如何，您的建议将不胜感激！
编辑：顺便说一下，LMM的似然函数可以在Harville（1977）中找到，（希望）可通过此链接访问： 最大似然方法用于方差分量估计和相关问题 问候， 西蒙
我使用lme4中的lmer()和logLik()来检查我的脚本是否按照预期工作。但似乎并没有！ 由于我的教育背景并不真正关注这个统计水平，所以我有点迷失在找错中。
接下来，您将找到一个使用sleepstudy数据的简短示例脚本：
# * * * * * * * * * * * * * * * * * * * * * * * *
# * 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）中找到，（希望）可通过此链接访问： 最大似然方法用于方差分量估计和相关问题 问候， 西蒙
lme4的开发版本（可能是通过devtools从github获取），该版本具有返回偏差函数的能力（mkDevfunOnly=TRUE）。 - Ben Bolkerlme4的github版本，并使用devtools进行了安装。是否有关于devFunOnly=T参数及其生成的函数的进一步文档？我特别关注必须提供给生成的偏差函数的参数，因为这对我来说是最重要的一步！ - SimonG