估计logit模型R中二元和连续系数的平均边际效应

问题描述

我希望使用国家健康访问调查中的数据，通过logistic回归分析各种人口统计学因素对高血压预测概率的平均边际影响。为了澄清，通过平均边际效应，我的意思是我想计算每个X的平均值（如STATA输出）的边际效应。

我的问题是我既有二进制独立变量又有连续独立变量，但据我所读，以二进制平均值为0或1来评估二进制变量是没有意义的。我不知道我不知道如何进行回归分析，在这里我可以求出连续变量的均值，而不是二进制变量的均值。这是我到目前为止的代码。

#Here I create a data frame of the means of the continuous variables 
mean_df=df %>% select(c(AGE,BMICALC,FAMSIZE,YEARSONJOB,HOURSWRK)) %>% summarise_all(mean)


#here is my regression,variables here not in the line of code above are binary 
logit_margin_diabetes <- glm(DIABETES~scale(AGE)+scale(IMMIGRANT)+scale(HOURSWRK)+scale(BELOW_TWICE_POVERTY)
+scale(BMICALC)+scale(FEMALE)+scale(FAMSIZE)+scale(EDUC_1)+scale(EDUC_2)+scale(EDUC_3)+
scale(EDUC_4)+scale(SMOKE)+scale(MARRIED)+scale(HISP)+scale(AFR_AM)+scale(WHITE),data = df,family="binomial")

#This is the stage where I want to apply the logit so it is evaluated at the means of the continuous variables. But I don't know what to do about the binary variables 
marg_mean<-margins(logit_margin_diabetes,data=mean_df)
summary(marg_mean)

抱歉，我很难制作和进行MRE，因为我不知道R中具有此类信息的数据集。但是，如果有人可以提供任何建议，将不胜感激！谢谢。

这是每个注释的修改后输出。但我希望输出也显示SE，AME和p值

margins(logit_margin,at=list(AGE=35.93349,BMICALC=26.90704,FAMSIZE=2.495413,YEARSONJOB=4.538336,HOURSWRK=32.53768,IMMIGRANT=1,BELOW_TWICE_POVERTY=1,FEMALE=1,EDUC_1=1,EDUC_2=1,EDUC_3=1,EDUC_4=1,SMOKE=1,MARRIED=1,HISP=1,AFR_AM=1,WHITE=1))
summary(marg_mean)

这是我运行摘要（marg_mean）后看到的新输出的照片

解决方法

如果您将变量声明为因素，则margins包将自动处理此问题。请参见subsetting section of the vignette或您可以inspect the source code以查看将边际效应计算为因子变量的差异。

请注意，margins的默认设置是计算“平均边际效应”，而不是是“均值的边际效应”。 IMO，默认设置在大多数情况下是最好的，但是如果您坚持要考虑“综合的”平均观察，则可以使用at函数的margins参数很容易。

代码示例。在第一种情况下，vs被视为连续变量。在第二个中，vs被视为二进制变量。

library(margins)

mod1 <- glm(am ~ hp + vs,data=mtcars,family=binomial)
mod2 <- glm(am ~ hp + factor(vs),family=binomial)

margins(mod1)
#> Average marginal effects
#> glm(formula = am ~ hp + vs,family = binomial,data = mtcars)
#>        hp       vs
#>  -0.00203 -0.03193

margins(mod2)
#> Average marginal effects
#> glm(formula = am ~ hp + factor(vs),data = mtcars)
#>        hp      vs1
#>  -0.00203 -0.03154

编辑：这是at参数的示例：

margins(mod1,at=list(hp=200,vs=0))