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
как фиксированные, так и группирующие переменные со случайным эффектом в вашей модели; это не сработает. Я вижу, что Кроули действительно предлагает это в книге 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
glmer
? В случае сомнений изучите документацию, например,help("lmer")
иhelp("glmer")
. - person Roland   schedule 01.08.2016