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

Question

I have checked out the defensive methods as was discussed in this post in order to prevent this error but it still doesn't go away.

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

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

Answer 1

tl;dr you should use glmer instead. Because you haven't named your arguments, R is interpreting them by position (order). lmer's third argument is REML, so R thinks you're specifying REML=binomial, which isn't a legitimate value. family is the third argument to glmer, so this would work (sort of: see below) if you used glmer, but it's usually safer to name the arguments explicitly if there's any possibility of getting confused.

A reproducible example would be nice, but:

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

is a starting point. I foresee a few more problems though:

• if your data are not Bernoulli (0/1) (which I'm guessing not since your response is called Proportion), then you need to include the total number sampled in each trial, e.g. by specifying a weights argument
• you have Plot and Treatment as both fixed and as random-effect grouping variables in your model; that won't work. I see that Crawley really does suggest this in the R book (google books link).

Do not do it the way he suggests, it doesn't make any sense. Replicating:

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)

If we look at the among-group standard deviation, we see that it's zero for both groups; it's exactly zero for block because block is already included in the fixed effects. It need not be for the treatment:block interaction (we have treatment, but not the interaction between block and treatment, in the fixed effects), but is because there's little among-treatment-within-block variation:

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

Conceptually, it makes more sense to treat block as a random effect:

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

Sometimes you might want to treat it as a fixed effect:

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

Now the block-by-treatment interaction variance is effectively zero (because block soaks up more variability if treated as a fixed effect). However, the estimated effect of spraying is very nearly identical.

If you're worried about overdispersion you can add an individual-level random effect (or use MASS::glmmPQL; lme4 no longer fits quasi-likelihood models)

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

The observation-level effect has effectively replaced the treatment-by-block interaction (which is now close to zero). Again, the estimated spraying effect has hardly changed (but its standard error is twice as large ...)