“nlme{nlme}”中的嵌套随机效应

Question

我必须使用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相同(至少在结构上)？

Answer 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