如何在执行 lmer 函数时克服此错误： if (REML) p else 0L 中的错误：参数不可解释为逻辑？答案

【问题标题】：How do I overcome this error in execution of lmer function: Error in if (REML) p else 0L : argument is not interpretable as logical?如何在执行 lmer 函数时克服此错误： if (REML) p else 0L 中的错误：参数不可解释为逻辑？
【发布时间】：2016-12-06 09:03:53
【问题描述】：

我已经检查了this post 中讨论的防御方法，以防止此错误，但它仍然没有消失。

model<-lmer(Proportion~Plot+Treatment+(1|Plot/Treatment),binomial,data=data)

if (REML) p else 0L 中的错误：参数不能解释为逻辑

【问题讨论】：

好的，谢谢。那么我如何明确声明整个论点是二项式的呢？在 Crawley 的 R 书中，他简单地使用了这个公式：model
显然这是一个错字，您真的想使用glmer 函数吗？如有疑问，请研究文档，例如 help("lmer") 和 help("glmer")。

标签： r lme4

【解决方案1】：

tl;dr 你应该改用glmer。因为你还没有命名你的参数，R 是按位置（顺序）解释它们。 lmer 的第三个参数是 REML，所以 R 认为你指定了 REML=binomial，这不是一个合法的值。 family 是 glmer 的第三个参数，所以如果你使用glmer，这将起作用（有点：见下文），但如果有任何可能，明确命名参数通常更安全变得混乱。

一个可重现的例子会很好，但是：

model <- glmer(Proportion~Plot+Treatment+(1|Plot/Treatment),
    family=binomial,data=data)

是一个起点。不过，我预见到还有一些问题：

如果您的数据不是伯努利 (0/1)（我猜不是因为您的回复称为 Proportion），那么您需要包括每次试验中采样的总数，例如通过指定 weights 参数
您的模型中有Plot 和Treatment 作为固定变量和随机效应分组变量；那是行不通的。我看到 Crawley really does suggest this in the R book（谷歌图书链接）。

不要按照他的建议去做，这没有任何意义。复制：

library(RCurl)
url <- "https://raw.githubusercontent.com/jejoenje/Crawley/master/Data/insects.txt"
dd <- read.delim(text=getURL(url),header=TRUE)
## fix typo because I'm obsessive:
levels(dd$treatment) <- c("control","sprayed") 
library(lme4)
model <- glmer(cbind(dead,alive)~block+treatment+(1|block/treatment),
               data=dd,family=binomial)

如果我们查看组间标准差，我们会发现两组的标准差都为零；对于block，它完全为零，因为block 已经包含在固定效果中。不一定是treatment:block 交互作用（我们有treatment，但不是block 和treatment 之间的交互作用，在固定效应中），而是因为块间处理内的变化很小：

VarCorr(model)
##  Groups          Name        Std.Dev.  
##  treatment:block (Intercept) 2.8736e-09
##  block           (Intercept) 0.0000e+00

从概念上讲，将块视为随机效应更有意义：

dd <- transform(dd,prop=dead/(alive+dead),ntot=alive+dead)
model1 <- glmer(prop~treatment+(1|block/treatment),
               weights=ntot,
               data=dd,family=binomial)
summary(model)
## ...
## Formula: prop ~ treatment + (1 | block/treatment)
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  treatment:block (Intercept) 0.02421  0.1556  
##  block           (Intercept) 0.18769  0.4332  
## Number of obs: 48, groups:  treatment:block, 12; block, 6
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -1.1640     0.2042  -5.701 1.19e-08 ***
## treatmentsprayed   3.2434     0.1528  21.230  < 2e-16 ***

有时您可能想将其视为固定效果：

model2 <- update(model1,.~treatment+block+(1|block:treatment))
summary(model2)
## Random effects:
##  Groups          Name        Variance  Std.Dev. 
##  block:treatment (Intercept) 5.216e-18 2.284e-09
## Number of obs: 48, groups:  block:treatment, 12
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.5076     0.0739  -6.868 6.50e-12 ***
## treatmentsprayed   3.2676     0.1182  27.642  < 2e-16 ***

现在，逐个处理的交互方差实际上为零（因为如果将其视为固定效应，块会吸收更多的可变性）。然而，喷洒的估计效果几乎相同。

如果您担心过度分散，您可以添加个体级别的随机效应（或使用MASS::glmmPQL；lme4 不再适合准似然模型）

dd <- transform(dd,obs=factor(seq(1:nrow(dd))))
model3 <- update(model1,.~.+(1|obs))

## Random effects:
##  Groups          Name        Variance  Std.Dev. 
##  obs             (Intercept) 4.647e-01 6.817e-01
##  treatment:block (Intercept) 1.138e-09 3.373e-05
##  block           (Intercept) 1.813e-01 4.258e-01
## Number of obs: 48, groups:  obs, 48; treatment:block, 12; block, 6
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -1.1807     0.2411  -4.897 9.74e-07 ***
## treatmentsprayed   3.3481     0.2457  13.626  < 2e-16 ***

观察级效应有效地取代了逐块治疗的相互作用（现在接近于零）。同样，估计的喷洒效果几乎没有变化（但它的标准误差是两倍大......）

【讨论】：