Как мне преодолеть эту ошибку при выполнении функции lmer: Ошибка в if (REML) p else 0L: аргумент не интерпретируется как логический?

Я проверил защитные методы, как обсуждалось в это сообщение, чтобы предотвратить эту ошибку, но она все еще не исчезла.

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

Ошибка в if (REML) p else 0L: аргумент не интерпретируется как логический


person Dominique    schedule 01.08.2016    source источник
comment
Хорошо спасибо. Итак, как мне явно заявить, что весь аргумент является биномиальным? В книге Кроули R он просто использует эту формулу: модель‹-lmer(y~блок+обработка+(1|блок/обработка),биномиальная,данные=данные) сводка(модель)   -  person Dominique    schedule 01.08.2016
comment
Очевидно, это опечатка, и вы действительно хотите использовать функцию glmer? В случае сомнений изучите документацию, например, help("lmer") и help("glmer").   -  person Roland    schedule 01.08.2016


Ответы (1)


tl;dr вместо этого следует использовать glmer. Поскольку вы не назвали свои аргументы, R интерпретирует их по положению (порядку). Третий аргумент lmerREML, поэтому R считает, что вы указываете REML=binomial, что не является допустимым значением. family является третьим аргументом glmer, так что это будет работать (вроде: см. ниже), если вы используете glmer, но обычно безопаснее называть аргументы явно, если есть вероятность запутаться.

Воспроизводимый пример был бы хорош, но:

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

является отправной точкой. Я предвижу еще несколько проблем, хотя:

  • если ваши данные не Бернулли (0/1) (что, как я предполагаю, нет, поскольку ваш ответ называется Proportion), тогда вам нужно включить общее число, отобранное в каждом испытании, например. указав аргумент weights
  • у вас есть Plot и Treatment как фиксированные, так и группирующие переменные со случайным эффектом в вашей модели; это не сработает. Я вижу, что Кроули действительно предлагает это в книге R (ссылка на книги Google).

Не делайте так, как он предлагает, это бессмысленно. Репликация:

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

Эффект уровня наблюдения эффективно заменил взаимодействие «лечение за блоком» (которое сейчас близко к нулю). Опять же, расчетный эффект распыления почти не изменился (но его стандартная ошибка вдвое больше...)

person Ben Bolker    schedule 01.08.2016