了解混合效應模型的邊際可能性

我正在閱讀一些關於混合模型的信息，但我不確定所使用的術語以及它們如何組合在一起。Pinheiro 在他的書“S 和 S-Plus 中的混合效應模型”的第 62 頁上描述了似然函數。

第二個方程的第一項被描述為條件密度，第二個是邊際密度.

我一直在嘗試為簡單的隨機效應模型生成這些對數似然 (ll)，因為我認為這將有助於我的理解，但我一定是誤解了推導。

---

我嘗試計算 ll 的示例。

示例模型

library(lme4)
model <- lmer(angle  ~ temp + (1|replicate), data=cake, REML=FALSE)

條件對數似然

我試圖計算這個模型的條件對數似然：它似乎來自 我應該能夠通過在單元/複製級別的預測中找到數據的密度來計算這一點。

sum(dnorm(cake$angle,
         predict(model), # predictions at replicate unit (XB + Zb) level
         sd=sigma(model),
         log=TRUE))
#[1] -801.6044

# Which seems to agree with
cAIC4::cAIC(model)$loglikelihood
# [1] -801.6044

# or should we really be using a multivariate normal density
# but doesn't make a difference as variance is \sigma^2 I_n
dmvnorm(cake$angle, predict(model), diag(sigma(model)^2,270, 270), log=TRUE)
#[1] -801.6044

邊際對數似然

我嘗試計算邊際對數似然，lme4給出為

logLik(model)
#'log Lik.' -834.1132 (df=4)

採取與以前類似的方法 我應該能夠通過在人口水平的預測中找到數據的密度來計算這一點，但它並不接近。

sum(dnorm(cake$angle,
         predict(model, re.form=NA), # predictions at population (XB) level
         sd=sigma(model),
         log=TRUE))
# [1] -1019.202

所以也許我需要使用第二個等式並且需要使用y的條件模型和****b的邊際，但仍然沒有接近。

sum(
 dnorm(cake$angle, predict(model), sd=sigma(model), log=TRUE) , # conditional
 dnorm(0,  ranef(model)$replicate[[1]], # RE predictions
           sd=sqrt(VarCorr(model)$replicate), log=TRUE) 
 ) 
# [1] -849.6086

編輯：接下來去…

對於線性混合模型,和我認為似然計算的方差應該估計為 var(Y) =，但又錯了！

所以在代碼中

z = getME(model, "Z")
zt =  getME(model, "Zt")
psi = bdiag(replicate(15, VarCorr(model)$replicate, simplify=FALSE))

betw = z%*%psi%*%zt
err = Diagonal(270, sigma(model)^2)
v = betw + err

sum(dnorm(cake$angle,
     predict(model, re.form=NA), 
     sd=sqrt(diag(v)),
     log=TRUE))
# [1] -935.652

我的問題：

• 你能告訴我在計算邊際可能性時哪裡出錯了。我真的不需要代碼來重現 ll，而是更多地描述為什麼我嘗試的方法不起作用。

非常感謝。

PS。我確實查看了生成這些值的函數，lme4:::logLik.merMod並lme4:::devCrit看到作者使用了一些困難/技術方法，這再次導致我需要幫助，為什麼我的方法是錯誤的。

lme4:::logLik.merMod我可以通過意識到Y的邊際分佈是多元正態 (MVN)來重現返回的邊際對數似然（因為b的邊際分佈和****Y的條件是 MVN）。

所以這段代碼會重現

library(mvtnorm)
dmvnorm(cake$angle, predict(model, re.form=NA), as.matrix(v), log=TRUE)
#[1] -834.1132

其中cake$angle是觀察值，predict(model, re.form=NA)是總體預測（使用固定效應係數計算），並且v是邊際模型的方差（如問題所示）。

---

關於我的問題中失敗的努力的一些評論。

在計算條件對數似然時，我使用了單變量正態密度函數，而我可以/應該使用多變量。在這種情況下，它沒有區別，因為單位方差是

mvtnorm::dmvnorm(cake$angle, predict(model), diag(sigma(model)^2,270, 270), log=TRUE)
#[1] -801.6044

試圖計算邊際分佈

• 首先嘗試使用總體（XB）級別的預測，但方差不正確（它忽略了單位間誤差）

sum(dnorm(cake$angle, predict(model, re.form=NA), sd=sigma(model), log=TRUE))

• 我認為第二次嘗試只是胡說八道
• 第三次嘗試使用單變量而不是多變量正態分佈。

引用自：https://stats.stackexchange.com/questions/271903