Блог

Моя glmer модель, которая содержит два предиктора и термин взаимодействия, страдает от полного разделения. Следуя рекомендациям Бена Болкера здесь и здесь, я затем подгоняю модель bglmer , накладывая нулевые средние значения нормы на фиксированные эффекты. Мой код выглядит следующим образом:

 bglmer(Binary_outcome ~ (1|Subject)   Factor1   Factor2   Factor1:Factor2, 
       mydata, 
       control=glmerControl(optimizer="bobyqa"), 
       family = binomial, 
       fixef.prior = normal(sd = c(3, 3, 3)))
  

И Factor1, и Factor2 являются факторными переменными с четырьмя уровнями каждая. Для моего кода я следовал примеру здесь. Насколько я понимаю, теперь я устанавливаю нормальные значения нулевого значения с SD, равным 3, на все элементы моей структуры фиксированных эффектов.

Код, похоже, сработал, но я совершенно не уверен, действительно ли то, что я сделал, правильно. Является ли 3 SD общей рекомендацией для полного разделения? И как бы я указал fixef.priors , что это относится только к термину взаимодействия? (Полное разделение относится к конкретной комбинации Factor1 и Factor2 , а не Factor1 или Factor2 в целом). Или я должен в любом случае устанавливать значения фиксированного эффекта для всех трех элементов, если речь идет о взаимодействии?

Если сомневаетесь, экспериментируйте. tl;dr вероятно, можно штрафовать все, но вы можете, если хотите, штрафовать взаимодействия и не штрафовать основные эффекты (или почти; вы должны указать конечный sd, но он может быть очень большим). sd = 3 является разумным; sd = 2.5 по умолчанию для параметров, отличных от перехвата. Если вы не подбираете прогностическую модель и не хотите утруждать себя перекрестной проверкой, чтобы выбрать оптимальную силу наказания, на самом деле нет автоматического способа принятия решения; вам просто нужно выбрать sd, который делает все ваши параметры «разумными», не слишком приближая ваши четко определенные параметры к нулю.

Я смоделировал данные, немного похожие на ваши (2 фактора по 4 уровня каждый, двоичный ответ, один случайный термин перехвата) и попробовал разные сценарии наказания. Я настроил истинные параметры так, чтобы большинство параметров были разумными (бета = 1), но был один большой основной эффект и один большой параметр взаимодействия (бета = 12).

• нет наказания ( glmer )
• наказание с sd=3
• штраф за все термины с sd=1
• очень слабое наказание за основные условия эффекта ( sd=1e4 ) и среднее наказание за взаимодействия ( sd=2 ) [«смешанный»]

От help("bmerDist-class") :

При указании стандартных отклонений вектор длиной меньше числа фиксированных эффектов будет иметь повторяющийся хвост, в то время как предполагается, что первый элемент применяется только к термину перехвата. Таким образом, по умолчанию ‘c (10, 2.5)’, перехват получает стандартное отклонение 10, а различным наклонам присваивается стандартное отклонение 2.5

 library(lme4)
library(blme)
library(broom.mixed)
library(dotwhisker)
library(colorspace)

dd <- expand.grid(Subject=factor(1:10),
                  Factor1=letters[1:4],
                  Factor2=LETTERS[1:4],
                  rep=1:20
                  )
bvec <- rep(1,16)
bvec[c(2,10)] <- 12
form <- response ~ (1|Subject)   Factor1   Factor2   Factor1:Factor2
set.seed(101)
dd$response <- simulate(form[-2],
         newdata=dd,
         newparams=list(beta=bvec, theta=1),
         family=binomial,
         weights=rep(1,nrow(dd)))[[1]]

b0 <- glmer(form, 
            dd,
            family = binomial)
bfun <- function(sd) {
   bglmer(form, dd, family = binomial, fixef.prior = normal(sd = sd))
}
b1 <- bfun(3)
b2 <- bfun(1)
## eight intercept   main effects first, then eight interaction parameters
b3 <- bfun(rep(c(1e4,2),c(8,8))

theme_set(theme_bw())
dwplot(list(unpenalized=b0,sd3=b1,sd1=b2,mixed=b3),effect="fixed")  
    coord_cartesian(xlim=c(-5,5)) 
    geom_vline(xintercept=c(0,1),lty=2,colour="darkgray")  
    scale_colour_discrete_qualitative(guide=guide_legend(reverse=TRUE))