如何克服执行 lmer 函数时出现的错误:Error in if (REML) p else 0L : argument is not interpretable as logical?

How do I overcome this error in execution of lmer function: Error in if (REML) p else 0L : argument is not interpretable as logical?

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

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

Error in if (REML) p else 0L : argument is not interpretable as logical

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 参数
  • 您的模型中有 PlotTreatment 作为固定和随机效应分组变量;那是行不通的。我看到 Crawley really does suggest this in the R book(google 本书 link)。

不要照他说的去做,没有任何意义。复制:

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,但在固定效应中没有 blocktreatment 之间的交互),而是因为它们之间几乎没有-块内处理变化:

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::glmmPQLlme4 不再适合准似然模型)

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 ***

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