Предельные эффекты от блеска

Question

Вариации по этому вопросу задавались ранее (например, Есть ли способ получить "маргинальные эффекты" от объекта "glmer" ), и большинство из них предлагают использовать ggeffects (или sjPlot). Тем не менее, статистик на нашем факультете испытывает некоторые трудности с согласием с результатами ggeffects.

Я новичок как в статистике, так и в R (так что мне, вероятно, не стоит работать с glmer ...), и у меня возникли проблемы с пониманием того, что я делаю неправильно.

Моя модель в простейшем виде - блеск от lme4: исход (двоичный) ~ FE (двоичный) + (1 | RE) . Фиксированный эффект - это тест, выполненный на некоторых, но не на всех людях из моего случайного эффекта.

Модель выхода

Intercept: 1.2654

FE_2: -0.2305

RE Std. Dev.: 2.896

ggpredict:

ggpredict(model, type = "fe", terms = "FE")

1: 0.780

2: 0.738

Теперь, насколько я понял, я могу получить предельный эффект по вероятности для теста 2 следующим образом:

y <- 1.2654 - 0.2305
prob <- exp(y) / (1 + exp(y))

Что в точности совпадает с выводом ggpredict: 0,738.

Однако он говорит, что это условная вероятность для каждого человека, и что мне нужно «вставить что-то, чего я не понимаю», чтобы получить вероятность, которую я могу обобщить для совокупности. Его быстрый пример решения:

y <- 1.2654 - 0.2305 + rnorm(100000000) * 2.896 
prob_trus <- exp(y) / (1 + exp(y))
mean(prob_trus)
0.62

Использование ggaverage с теми же аргументами, что и в предыдущем ggpredict, дает почти такую ​​же вероятность, 0,639. Итак, наконец-то о моих вопросах:

• Является ли ggaverage тем же, что и его решение, только с некоторой разницей в симуляции и / или ошибкой округления?
• Если его метод - путь, как я могу получить те же результаты от ggeffects?

Answer 1

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

Хедекери другие.(2018) недавно предложили новую идею для получения коэффициентов регрессии с интерпретацией маргинал / популяция.Это реализовано в функции marginal_coefs() пакета R GLMMadaptive , которая подходит для смешанных моделей с использованием адаптивной квадратуры Гаусса.Типичным примером логистической регрессии со смешанными эффектами является:

library("GLMMadaptive")
fm <- mixed_model(fixed = y ~ x1 + x2 + x3, random = ~ x1 | id, data = DF, family = binomial())

marginal_coefs(fm)
marginal_coefs(fm, std_errors = TRUE)

Для дополнительных примеров, например, как получить предельные прогнозы или создать графики эффектов, проверьте виньетку: https://drizopoulos.github.io/GLMMadaptive/articles/Methods_MixMod.html

Answer 2

ggaverage() - это не то же самое, что метод «что-то, чего я не понимаю», подобные результаты обусловлены случайностью.Вы можете прочитать кое-что о различиях между ggpredict() и ggaverage() здесь, в виньетке .

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

Вы увидите разницу при сравнении ggpredict(model, type = "fe", terms = "FE") с ggpredict(model, type = "re", terms = "FE").Обратите внимание на изменение type = "re", которое теперь учитывает случайные эффекты, что приводит к различным прогнозируемым значениям и большим доверительным интервалам.

Вот пример с данными примера из пакета:

data(efc_test)

fit <- glmer(
  negc7d ~ c12hour + e42dep + c161sex + c172code + (1 | grp),
    data = efc_test,
    family = binomial(link = "logit")
  )

ggpredict(fit, "c12hour", type = "fe")

#> # Predicted probabilities for Negative impact with 7 items 
#> # x = average number of hours of care per week 
#> 
#>   x predicted conf.low conf.high
#>   0     0.342    0.249     0.449
#>   5     0.344    0.251     0.450
#>  10     0.346    0.254     0.452
#>  15     0.348    0.256     0.453
#>  20     0.350    0.258     0.455
#>  25     0.352    0.260     0.456
#>  30     0.354    0.261     0.458
#>  35     0.356    0.263     0.460
#>  40     0.357    0.265     0.463
#>  45     0.359    0.266     0.465
#>  ... and 25 more rows.
#> 
#> Adjusted for:
#> *   e42dep = 2.92
#> *  c161sex = 1.76
#> * c172code = 1.97

ggpredict(fit, "c12hour", type = "re")

#> # Predicted probabilities for Negative impact with 7 items 
#> # x = average number of hours of care per week 
#> 
#>   x predicted conf.low conf.high
#>   0     0.472    0.107     0.870
#>   5     0.475    0.108     0.871
#>  10     0.477    0.109     0.872
#>  15     0.479    0.110     0.873
#>  20     0.481    0.111     0.873
#>  25     0.483    0.111     0.874
#>  30     0.485    0.112     0.875
#>  35     0.487    0.113     0.876
#>  40     0.489    0.114     0.877
#>  45     0.491    0.115     0.878
#>  ... and 25 more rows.
#> 
#> Adjusted for:
#> *   e42dep = 2.92
#> *  c161sex = 1.76
#> * c172code = 1.97