Предельные эффекты logit в данных взвешенного опроса с использованием R

Question

Я пытаюсь оценить маргинальное влияние модели логита, в которой у меня есть несколько дихотомических объясняющих переменных.

Скажем, модель оценивается

logit<- svyglm (if_member ~ if_female + dummy_agegroup_2 + dummy_agegroup_3 + dummy_education_2 + dummy_education_3 + dummy_education_4, family = quasibinomial(link = "logit"), design = survey_design) 

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

Answer 1

Вы ищете маргинальный номер effects или marginal прогнозов?

Как следует из названия, функция marginpred() возвращает прогнозов. Аргумент для predictat представляет собой кадр данных с контрольными переменными и переменными, которые находятся в модели. Подчеркнем, что: контрольные переменные должны быть исключены из модели.

library("survey") 

odds2prob <- function(x) x/(x + 1) 
prob2odds <- function(x) x/(1 - x) 
expit <- function(x) odds2prob(exp(x)) 
logit <- function(x) log(prob2odds(x)) 

set.seed(1) 

survey_data <- data.frame(
    if_female = rbinom(n = 100, size = 1, prob = 0.5), 
    agegroup = factor(sample(x = 1:3, size = 100, replace = TRUE)), 
    education = NA_integer_, 
    if_member = NA_integer_) 
survey_data["agegroup"] <- relevel(survey_data$agegroup, ref = 3) 
# Different probabilities between female and male persons 
survey_data[survey_data$if_female == 0, "education"] <- sample(
    x = 1:4, 
    size = sum(survey_data$if_female == 0), 
    replace = TRUE, 
    prob = c(0.1, 0.1, 0.5, 0.3)) 
survey_data[survey_data$if_female == 1, "education"] <-sample(
    x = 1:4, 
    size = sum(survey_data$if_female == 1), 
    replace = TRUE, 
    prob = c(0.1, 0.1, 0.3, 0.5)) 
survey_data["if_member"] <- rbinom(n = 100, size = 1, prob = 
            expit((survey_data$education - 3)/2)) 
survey_data["education"] <- factor(survey_data$education) 
survey_data["education"] <- relevel(survey_data$education, ref = 3) 
survey_design <- svydesign(ids = ~ 1, data = survey_data) 

logit <- svyglm(if_member ~ if_female + agegroup + education, 
       family = quasibinomial(link = "logit"), 
       design = survey_design) 
exp(cbind(`odds ratio` = coef(logit), confint(logit))) 
newdf <- data.frame(if_female = 0:1, education = c(3, 3), agegroup = = c(3, 3)) 
# Fails 
mp <- marginpred(model = logit, adjustfor = ~ agegroup + education, 
       predictat = newdf, se = TRUE, type = "response") 
logit2 <- svyglm(if_member ~ if_female, 
       family = quasibinomial(link = "logit"), 
       design = survey_design) 
mp <- marginpred(model = logit2, adjustfor = ~ agegroup + education, 
       predictat = newdf, se = TRUE, type = "response") 
# Probability for male and for female persons controlling for agegroup and education 
cbind(prob = mp, confint(mp)) 

Вот как я оцениваю маргинальным эффекты с survey пакета:

# Probability difference between female and male persons 
# when agegroup and education are set to 3 
svycontrast(full_model, quote(
    (exp(`(Intercept)` + if_female)/(exp(`(Intercept)` + if_female) + 1)) - 
    (exp(`(Intercept)`)/(exp(`(Intercept)`) + 1)))) 
# Can't use custom functions like expit :_(

Есть, вероятно, более умные способы, но я надеюсь, что это помогает.

Обратите внимание, что разница между вероятностями, предсказанными marginpred(), отличается от разницы, оцененной svycontrast(). Вероятности, предсказанные marginpred(), по-видимому, не зависят от изменения значения управляющих переменных (например, education = c(4, 4) вместо education = c(3, 3)), но оценки от svycontrast() зависят от модели регрессии.