Можно ли построить прогнозы модели R glmer с использованием Python? - PullRequest
Новая
       68

Можно ли построить прогнозы модели R glmer с использованием Python?

2 голосов
/ 04 мая 2019

У меня есть модель glmer в R, для которой я хочу построить прогнозы.Я нашел функцию plot_model из библиотеки sjPlot, и она отлично работает.

Вот MWE: 

library(lattice)

cbpp$response <- sample(c(0,1), replace=TRUE, size=nrow(cbpp))
gm1 <- glmer(response ~ size + incidence + (1 | herd),
              data = cbpp, family = binomial)

Например, вызов plot_model(gm1, type = "pred", show.data = TRUE) дает следующий рисунок:

enter image description here

Однако я не знаком с R, и мне трудно пытаться контролировать эстетику сюжета и вычерчивать несколько моделей на одной фигуре (уже задавался вопросвопрос по этому вопросу здесь ).Я знаком с Python и matplotlib, и заставить эти цифры работать в среде Python было бы намного проще для меня.

Я предполагаю, что одним из способов достижения этого было бы взять значения y (прогнозируемые вероятности возникновения пожара) из R и экспортировать их, чтобы я мог прочитать их на Python, чтобы построить их для каждого ковариата (evi prev) в этом примере.Однако я не уверен, как это сделать.Кроме того, я попытался прочитать исходный код sjPlot, чтобы выяснить, как он строит прогнозы, но также не смог понять его.

Ответы [ 2 ]

2 голосов
/ 06 мая 2019

ggpredict() фактически возвращает больше значений (и вдоль оси x, т. Е. Для рассматриваемого термина - size в вашем примере - они с четным интервалом), но только печатает меньше значений . 

library(lme4)
#> Loading required package: Matrix
library(ggeffects)

cbpp$response <- sample(c(0,1), replace=TRUE, size=nrow(cbpp))
gm1 <- glmer(response ~ size + incidence + (1 | herd), data = cbpp, family = binomial)

pr1 <- ggpredict(gm1, term = "size")

pr1
#> 
#> # Predicted probabilities of response
#> # x = size
#> 
#>   x predicted std.error conf.low conf.high
#>   2     0.632     0.717    0.297     0.875
#>   6     0.610     0.550    0.347     0.821
#>  10     0.587     0.407    0.390     0.759
#>  14     0.563     0.321    0.407     0.708
#>  18     0.539     0.339    0.376     0.695
#>  22     0.515     0.448    0.306     0.719
#>  26     0.491     0.601    0.229     0.758
#>  34     0.444     0.951    0.110     0.837
#> 
#> Adjusted for:
#> * incidence = 1.77
#> *      herd = 0 (population-level)
#> Standard errors are on link-scale (untransformed).

as.data.frame(pr1)
#>     x predicted std.error  conf.low conf.high group
#> 1   2 0.6323758 0.7168742 0.2967912 0.8751705     1
#> 2   4 0.6211339 0.6316777 0.3221952 0.8497229     1
#> 3   6 0.6097603 0.5501862 0.3470481 0.8212222     1
#> 4   8 0.5982662 0.4743133 0.3701925 0.7904902     1
#> 5  10 0.5866630 0.4072118 0.3898523 0.7592017     1
#> 6  12 0.5749627 0.3539066 0.4033525 0.7302266     1
#> 7  14 0.5631779 0.3213384 0.4071542 0.7076259     1
#> 8  16 0.5513213 0.3159857 0.3981187 0.6953669     1
#> 9  18 0.5394060 0.3391396 0.3759558 0.6947993     1
#> 10 20 0.5274456 0.3857000 0.3438768 0.7038817     1
#> 11 22 0.5154536 0.4484344 0.3063836 0.7192510     1
#> 12 24 0.5034437 0.5215385 0.2672889 0.7380720     1
#> 13 26 0.4914299 0.6012416 0.2292244 0.7584368     1
#> 14 28 0.4794260 0.6852450 0.1938167 0.7791488     1
#> 15 30 0.4674458 0.7721464 0.1619513 0.7994688     1
#> 16 32 0.4555030 0.8610687 0.1339908 0.8189431     1
#> 17 34 0.4436111 0.9514457 0.1099435 0.8373008     1

Создано в 2019-05-06 пакетом Представления (v0.2.1)

Существуют некоторые виньетки, которые показывают различные особенности пакета, этот здесь демонстрирует, как вычислить предельные эффекты при определенных значениях / уровнях фокусных терминов.

Рецепт, опубликованный Беном, который показывает, как рассчитать доверительные интервалы (обусловленные или не обусловленные случайными эффектами), реализован в ggpredict(), короткая виньетка, объясняющая различия , находится здесь .

2 голосов
/ 04 мая 2019

Самый простой способ сделать это, вероятно, с ggeffects::ggpredict().

Что-то вроде 

library(ggeffects)
pred_frame <- ggpredict(myModel, term="evi_prev")

должно создать кадр данных с предсказаниями, нижним и верхним доверительными уровнями.Я не уверен, будет ли он делать прогнозы для равномерно распределенных значений вдоль оси x (что было бы неплохо), или как его обмануть.(Если вы предоставите воспроизводимый пример, я мог бы дать ему шанс.)

Игра с MWE, которую вы разместили, предполагает, что трудно получить прогнозы для равномерно распределенных значений (или, в более общем случае, для значений, которые не 'т в исходных данных);Я пробовал что-то вроде terms="size [1:35]", но это ограничивает диапазон прогнозируемых значений, а не заполняет их.

По сути, встроенный метод predict() для merModобъекты могут использоваться (возможно, с newdata для указания, например, равномерно распределенных значений) для получения прогнозов [используйте type="response" для получения прогнозов по вероятности, а не по шкале логарифмов];доверительные интервалы сложнее, но могут быть получены с помощью показанного рецепта здесь

...