Воспроизведите «поля» из Stata и рассчитайте разницу между двумя прогнозируемыми средними значениями с 95% ДИ
Я надеюсь повторить то, что создал здесь потрясающий Чак Хубер (поля). В этом примере он смог рассчитать прогнозируемый вес при рождении для курящих и некурящих матерей. Кодекс штата звучал примерно так:
regress bwt age i.smoke
margins smoke
Где bwt = вес при рождении, возраст = возраст и i.smoke = курильщик против некурящего в качестве факторной переменной.
я зашел так далеко
library(tidyverse)
library(margins)
set.seed(42)
n <- 1000
dat <- data.frame(id=1:n,
treat = factor(sample(c('Treat','Control'), n, rep=TRUE, prob=c(.5, .5))),
age=sample(18:80, n, replace=TRUE),
sex = factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4))),
smoke=factor(sample(c("Never", 'Former', 'Current'), n, rep=TRUE, prob=c(.25, .6, .15))),
score=runif (n, min=16, max=45))
mod <- lm(score~treat+age+sex+smoke, data=dat)
summary(margins(mod, variables = "treat"))
Выходные данные предоставляют средний предельный эффект (AME) для treat, но только для одного уровня этой переменной.
Мои вопросы:
- Есть ли способ, чтобы выходные данные отображали AME и 95% CI дляtreatTreat иtreatControl (а именно, два уровня этого фактора)?
- Как я могу вычислить среднюю разницу и 95% ДИ между этими двумя AME дляtreatTreat иtreatControl?
Большое спасибо,
Сандро
Ответы 1
Пакет margins больше не разрабатывается активно и по состоянию на 23 апреля 2024 г. не был доступен в CRAN. Это архивирование может быть временным.
Вместо этого вы можете использовать пакет marginaleffects (отказ от ответственности: я автор). На веб-сайте содержится более 30 глав с обучающими материалами и практическими примерами, включая вводную страницу «Начало работы»: https://marginaleffects.com
Прежде чем показывать код, нам нужно договориться о терминологии. Команда margins smoke в Stata вычисляет «предельные прогнозы», то есть средние прогнозы по подгруппам. Это не то же самое, что средний предельный эффект, который рассчитывает пакет margins в R. AME в этом пакете лучше называть «Средним наклоном» и определяется как частная производная итогового уравнения по интересующей переменной.
Вы можете легко вычислить все эти величины с помощью marginaleffects. Сначала подгоните модель:
library(dplyr)
library(marginaleffects)
set.seed(42)
n <- 1000
dat <- data.frame(id=1:n,
treat = factor(sample(c('Treat','Control'), n, rep=TRUE, prob=c(.5, .5))),
age=sample(18:80, n, replace=TRUE),
sex = factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4))),
smoke=factor(sample(c("Never", 'Former', 'Current'), n, rep=TRUE, prob=c(.25, .6, .15))),
score=runif (n, min=16, max=45))
mod <- lm(score~treat+age+sex+smoke, data=dat)
Теперь вычислите средние прогнозы для каждой подгруппы переменной treat. Я показываю код dplyr, чтобы вы точно понимали, что происходит:
avg_predictions(mod, by = "treat")
treat Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
Treat 30.9 0.38 81.4 <0.001 Inf 30.2 31.7
Control 30.9 0.36 86.0 <0.001 Inf 30.2 31.6
Columns: treat, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response
dat |> mutate(prediction = predict(mod, newdata = dat)) |>
summarize(prediction = mean(prediction), .by = treat)
treat prediction
1 Treat 30.90587
2 Control 30.92291
Вы можете использовать аргумент hypothesis для проверки гипотезы, сравнивающей две величины:
avg_predictions(mod, by = "treat", hypothesis = "sequential")
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
Control - Treat 0.017 0.523 0.0326 0.974 0.0 -1.01 1.04
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response
dat |> mutate(prediction = predict(mod, newdata = dat)) |>
summarize(prediction = mean(prediction), .by = treat) |>
summarize(diff = diff(prediction))
diff
1 0.01703766
Альтернативно, вы можете сравнить средние контрфактические прогнозы, используя функцию avg_comparisons():
avg_comparisons(mod, variables = "treat")
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
treat Treat - Control -0.0108 0.523 -0.0207 0.983 0.0 -1.04 1.01
Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type: response
d0 <- mutate(dat, treat = "Control")
d1 <- mutate(dat, treat = "Treat")
mean(predict(mod, newdata = d1) - predict(mod, newdata = d0))
[1] -0.01083374
Спасибо @ Винсент. Всегда приятно услышать от разработчика пакета объяснение функций его пакета.