`lme {nlme}` 中的嵌套随机效果答案

【问题标题】：Nested random effects in `lme {nlme}``lme {nlme}` 中的嵌套随机效果
【发布时间】：2017-04-28 10:30:54
【问题描述】：

我必须使用lme 命令来拟合具有交互随机效应但没有边际随机效应的 LMM。也就是说，我想在 oats.lmer（见下文）中拟合模型，但使用 nlme 包中的函数 lme。

代码是

require("nlme")
require("lme4")

oats.lmer <- lmer(yield~nitro + (1|Block:Variety), data = Oats)

summary(oats.lmer)
#Linear mixed model fit by REML ['lmerMod']
#Formula: yield ~ nitro + (1 | Block:Variety)
#   Data: Oats
# 
#REML criterion at convergence: 598.1
#
#Scaled residuals: 
#     Min       1Q   Median       3Q      Max 
#-1.66482 -0.72807 -0.00079  0.56416  1.85467 
# 
#Random effects:
# Groups        Name        Variance Std.Dev.
# Block:Variety (Intercept) 306.8    17.51   
# Residual                  165.6    12.87   
#Number of obs: 72, groups:  Block:Variety, 18
#
#Fixed effects:
#            Estimate Std. Error t value
#(Intercept)   81.872      4.846   16.90
#nitro         73.667      6.781   10.86
#
#Correlation of Fixed Effects:
#      (Intr)
#nitro -0.420

我开始玩这个

oats.lme <- lme(yield~nitro, data = Oats, random = (~1|Block/Variety))
summary(oats.lme)
#Linear mixed-effects model fit by REML
# Data: Oats 
#       AIC      BIC    logLik
#  603.0418 614.2842 -296.5209
#
#Random effects:
# Formula: ~1 | Block
#        (Intercept)
#StdDev:    14.50596
#
# Formula: ~1 | Variety %in% Block
#        (Intercept) Residual
#StdDev:    11.00468 12.86696
#
#Fixed effects: yield ~ nitro 
#               Value Std.Error DF  t-value p-value
#(Intercept) 81.87222  6.945273 53 11.78819       0
#nitro       73.66667  6.781483 53 10.86291       0
# Correlation: 
#      (Intr)
#nitro -0.293
#
#Standardized Within-Group Residuals:
#        Min          Q1         Med          Q3         Max 
#-1.74380770 -0.66475227  0.01710423  0.54298809  1.80298890 
#
#Number of Observations: 72
#Number of Groups: 
#             Block Variety %in% Block 
#                 6                 18

但问题是它还对Variety 产生了边际随机效应，我想省略它。

问题是：如何指定oats.lme 中的随机效应，使oats.lme 与oats.lmer 相同（至少在结构上）？

【问题讨论】：

您的意思是“Block 的边际随机效应”吗？

标签： r regression mixed-models nlme

【解决方案1】：

可以这么简单：

library(nlme)
data(Oats)

## construct an auxiliary factor `f` for interaction / nesting effect
Oats$f <- with(Oats, Block:Variety)
## use `random = ~ 1 | f`
lme(yield ~ nitro, data = Oats, random = ~ 1 | f)

#Linear mixed-effects model fit by REML
#  Data: Oats 
#  Log-restricted-likelihood: -299.0328
#  Fixed: yield ~ nitro 
#(Intercept)       nitro 
#   81.87222    73.66667 
#
#Random effects:
# Formula: ~1 | f
#        (Intercept) Residual
#StdDev:    17.51489 12.86695
#
#Number of Observations: 72
#Number of Groups: 18

【讨论】：