На что вы устанавливаете коэффициент группировки при использовании glmer / lme4 и ForexInterval? - PullRequest
Новая
       104

На что вы устанавливаете коэффициент группировки при использовании glmer / lme4 и ForexInterval?

1 голос
/ 25 сентября 2019

Проблема: использование многоуровневой модели (смешанных эффектов) и отсутствие уверенности в том, на что установить переменную группировки, чтобы сгенерировать прогнозируемые вероятности для измеренной переменной на уровне группы из модели glmer с использованием функции предсказания merTools.

Цель: Генерация прогнозируемых вероятностей и SE / CI по диапазону значений из переменной уровня группы «второго уровня».

Поиск: советы о том, как правильно выполнить эту или другие рекомендации для генерации прогнозируемых вероятностей и CI.диапазон значений для переменной уровня группы из модели glmer. 

library(lme4)
library(merTools)
library(ggplot2)

hier_data <- data_frame(pass = sample(c(0, 1), size = 1000, replace = T),
                        wt = rnorm(1000),
                        ht = sample(1:5, size = 1000, replace = T, prob = c(.1, .1, .6, .1, .1)),
                        school_funding = rnorm(1000),
                        school = rep(c("A", "B", "C", "D", "E"), each = 200))

mod <- glmer(pass ~ wt + ht + school_funding + (1 | school),
             family = binomial("logit"), data = hier_data)

### Without school: error
ndata <- data.frame(wt = median(hier_data$wt),
                    ht = median(hier_data$ht),
                    school_funding = seq(from = min(hier_data$school_funding), to =max(hier_data$school_funding), length.out = 100))

pp <- cbind(ndata, predictInterval(merMod = mod,
                      newdata = ndata,
                      type = "probability"))

### Problem, when adding school variable: which school?
ndata <- data.frame(wt = median(hier_data$wt),
                    ht = median(hier_data$ht),
                    school_funding = seq(from = min(hier_data$school_funding), to =max(hier_data$school_funding), length.out = 100),
                    school = "A")

pp <- cbind(ndata, predictInterval(merMod = mod,
                                   newdata = ndata,
                                   type = "probability"))

ggplot(pp, aes(x = school_funding, y = fit)) +
    geom_point() +
    geom_errorbar(aes(ymin = lwr, ymax = upr))

1 Ответ

0 голосов
/ 25 сентября 2019

Кажется, что вы пытаетесь достичь effects графиков для ваших переменных с быстрыми интервалами прогнозирования.Прежде всего обратите внимание, что predictInterval не включает неопределенность в оценочные значения параметров дисперсии , theta.Если требуются более точные доверительные интервалы, вы должны использовать функцию bootMer, как описано в ?bootMer, которая использует начальную загрузку для оценки неопределенности.Однако это может быть просто невозможно, поскольку размер и сложность модели возрастают.В качестве альтернативы пакет effects содержит возможность иллюстрировать эффекты объектов merMod (однако документация просто отвратительна).

В целом, когда иллюстрируются эффекты объектов merModвопрос «какие эффекты?».Вас интересуют предельные эффекты или условные эффекты (например, изменчивость случайных эффектов?).Если ваша модель содержит только случайные эффекты первого порядка (без случайных уклонов), и вас интересует неопределенность коэффициента с фиксированным эффектом или влияние на условное среднее, вы можете уйти с помощью любой школы и указать which = "fixed"в predictInterval как 

pp <- cbind(ndata, predictInterval(merMod = mod,
                                   newdata = ndata, #<= any school chosen
                                   type = "probability",
                                   which = "fixed"))

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

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

newData <- expand.grid(wt = median(hier_data$wt), 
                       ht = median(hier_data$ht),
                       school = levels(hier_data$school),
                       school_funding = seq(from = min(hier_data$school_funding), 
                                            to = max(hier_data$school_funding), 
                                            length.out = 100))
pp <- predictInterval(merMod = mod,
                       newdata = newData,
                       type = "probability")
#Split predictions by every column but school
# And calculate estimated means
predictions <- do.call("rbind", lapply(split(as.data.frame(pp), 
                                             newData[, !names(newData) == "school"]), 
                                       colMeans))
rownames(predictions) <- 1:nrow(predictions)
#create a plot
ggplot(as.data.frame(cbind(predictions, funding = newData$school_funding[newData$school == "A"])), 
       aes(x = funding, y = fit, ymax = upr, ymin = lwr)) + 
    geom_point() +
    geom_errorbar()

Для этого примера модель чаще, чем не единственнаяи содержит очень мало групп, и, как таковой, результат вряд ли будет отличной оценкой для предельного эффекта, но за исключением извлечения симуляций из predictInterval этого может быть достаточно.Вероятно, это улучшится с моделями с большим количеством уровней группировки в случайном эффекте.predictInterval, по-видимому, не включает в себя метод для этой ситуации напрямую.

Альтернативой для рассмотрения предельных эффектов было бы предположение о предельном среднем в форме 1/(1+exp(-eta) (что часто предполагается для новых группировокслучайный эффект).Это непосредственно не реализовано в функции predictInterval, но может быть достигнуто путем вычитания случайного эффекта из линейного предиктора и оценки только случайности фиксированных эффектов, как показано ниже: 

pp <- predictInterval(merMod = mod,
                      newdata = ndata, #<= any school chosen
                      type = "linear.prediction",
                      which = "fixed")
#remove random effects
pp <- sweep(pp, 1, predict(mod, newdata = ndata, random.only = TRUE), "-")
pp <- 1/(1+exp(-pp))

, котораязатем можно построить с использованием аналогичных методов.Для меньшего числа групп это могло бы быть лучшим предиктором для предельного среднего (?, Кто-то может исправить меня здесь).

В любом случае добавление небольшого количества x-jitter может улучшить сюжет.

Во всех случаях в ссылках на GLMM FAQ может быть несколько золотых самородков.Болкер и др.

...