Как сделать пошаговую модель со случайным эффектом (lme4 + lmerTest?) - PullRequest
Новая
       58

Как сделать пошаговую модель со случайным эффектом (lme4 + lmerTest?)

1 голос
/ 11 апреля 2019

Я пытаюсь выполнить пошаговую модель со случайным эффектом, из которого я могу получить значение BIC.

Пакет lmerTest сказал, что он работает с lme4, но я могу заставить его работать, только если яудалить одну из моих независимых переменных из модели (это фактор с двумя опциями (TM))

Код ошибки:

Ошибка в $<- (*tmp*, формула, значение = Термины): нет метода для присвоения подмножеств этого класса S4

или

Ошибка в as_lmerModLmerTest (модель): модель не класса 'lmerMod': не может привести к классу 'lmerModLmerTest

Я где-то читал, что это может быть связано с drop1, но я все еще не понял этого.Я также открыт для предложений других пакетов и функций.

Раньше, когда пробовал full.model <- lm (... все работало. После перехода на lmer больше не было. </p>

Код, который я сейчас использую: 

full.model <- lme4::lmer(dep ~ TM + ind + (1 | dorp),  data=test)  #lmerTest:: give same outcome

step.model<- lmerTest::step(full.model, direction="both",k=log(16))   # n=16

summary(step.model)

BIC(step.model)

#Example dataset

test <- data.frame(TM = as.factor(c(rep("org", 3), rep("min", 3),rep("org", 3), rep("min", 3),rep("org", 3), rep("min", 3))),
                   dep = runif(18,0,20),
                   ind = runif(18,0,7),
                   dorp = as.factor(c(rep(1,6),rep(2,6),rep(3,6))))

1 Ответ

0 голосов
/ 12 апреля 2019

Проблема в том, что lmerTest::step.lmerModLmerTest прерывается, когда все случайные эффекты исключаются из модели на этапе выбора случайных эффектов. Это, вероятно, не должно (я думаю, что более ранние версии пакета не могут), но это не слишком сложно обойти. Вы можете указать, что модель случайных эффектов не должна быть упрощена (step(full.model, reduce.random=FALSE)), или , когда вы столкнетесь с этой ошибкой, отбросьте компоненты случайных эффектов модели и затем используйте step() на Полученная линейная модель: 

fixmodel <- lm(formula(full.model,fixed.only=TRUE),
               data=eval(getCall(full.model)$data))
step(fixmodel)

(поскольку он включает eval(), это будет работать только в среде, где R может найти фрейм данных, на который ссылается аргумент data=).

Я отправил вопрос об этой проблеме.

Кроме того (смущает), stats::step имеет другие аргументы / делает другие предположения, чем step.lmerModLmerTest в пакете lmerTest. stats::step определяется как 

step(object, scope, scale = 0,
     direction = c("both", "backward", "forward"),
     trace = 1, keep = NULL, steps = 1000, k = 2, ...)

, в то время как step.lmerModLmerTest использует 

step(object, ddf = c("Satterthwaite",
  "Kenward-Roger"), alpha.random = 0.1, alpha.fixed = 0.05,
  reduce.fixed = TRUE, reduce.random = TRUE, keep, ...)

В частности, аргумент direction не применяется (step.lmerModLmerTest делает только обратное исключение); не делает k (я думаю, step.lmerModLmerTest использует AIC, но мне придется перепроверить). 

set.seed(1001)
dd <- data.frame(x1=rnorm(500),x2=rnorm(500),
                 x3=rnorm(500),f=factor(rep(1:50,each=10)))
library(lme4)
dd$y <- simulate(~x1+x2+x3+(1|f),
                 newdata=dd,
                 newparams=list(theta=1,beta=c(1,2,0,0),
                                sigma=1),
                 family=gaussian)[[1]]
library(lmerTest)
full.model <- lmer(y~x1+x2+x3+(1|f), data=dd)
step.model<- step(full.model)

step.model имеет класс step_list; есть метод печати, но нет итогового метода.

...