使用 R 报告调查加权 logit 模型的平均边际效应答案

【问题标题】：Reporting average marginal effects of a survey-weighted logit model with R使用 R 报告调查加权 logit 模型的平均边际效应
【发布时间】：2021-10-22 11:11:16
【问题描述】：

我正在使用复杂样本的调查数据来估计二元结果模型。我正在尝试报告 logit 模型的平均边际效应，我通过 R 中 survey 包的svyglm 估计。但是，当我使用包中的margins 时出现以下错误同名：

margins(fit, design = lapop) %>% summary()

Error in h(simpleError(msg, call)) : error in evaluating the argument 'object' in selecting a method for function 'summary': arguments imply differing number of rows: 6068, 6054

似乎它不是summary 函数，因为在执行带有参数的边距命令时会弹出错误。我试图完全忽略调查权重，并向我展示相等的系数和 AME，但不是标准误差。显然，我不能通过忽略调查权重来展示这项工作。所以我想我真正需要的是标准错误。

我一直在阅读该主题并没有找到明确的解决方案，我怀疑这可能与模型中 X 的缺失值有关，但与任何其他线性模型一样，R 应该只是使用完整的案例。

我不确定是否有人对此有所了解，或者我是否应该只报告没有标准错误（因此没有 p 值）的 AME。我已经上传了一个MWE，如果有人感兴趣，可以找到here。

【问题讨论】：

标签： r logistic-regression survey

【解决方案1】：

我发现这是怎么回事：原来这是包代码中的某种错误。你可以看看具体的here。今天的解决方案是安装 Tomasz Żółtak 的 prediction 和 margins 包的分支，直到他在 Github 上的拉取请求被合并。

devtools::install_github("tzoltak/prediction")
devtools::install_github("tzoltak/margins")

如果您还没有安装 devtools 软件包，则必须在安装后完成此操作。

install.packages('devtools')

执行此操作后，如果模型的数据框在模型的部分或全部协变量上缺少值，则在模型上运行 margins() 应该不再产生错误。因此，它将呈现平均部分效应及其相应的调查加权标准误差。查看 MWE here。

希望将来直接从 CRAN 调用 margins 将足以避免调查加权模型出现此错误。

【讨论】：

请添加更多详细信息以扩展您的答案，例如工作代码或文档引用。

【解决方案2】：

错误消息表明在某些行中存在NAs，R 不会自动排除。首先，我尝试使用fit 和lapop 变量来重现错误消息，并且确实弹出了错误。

margins(fit, design = lapop)

#Error in data.frame(..., check.rows = FALSE, check.names = FALSE, fix.empty.names = FALSE,  : 
# arguments imply differing number of rows: 6068, 6054

然后，我尝试确认哪个变量有问题NAs。

margins(fit)

#Note: Estimating marginal effects without survey weights. Specify 'design' to adjust for weighting.
#Error in data.frame(..., check.rows = FALSE, check.names = FALSE, fix.empty.names = FALSE,  : 
# arguments imply differing number of rows: 6068, 6054

弹出相同的错误消息，所以我相信fit 包含NAs。然后我检查了fit 在你的代码中是如何产生的：

fit<-svyglm(ctol ~ y16 + age,
            design = lapop,
            family = quasibinomial(link = 'logit'))

NAs 应该位于 ctol、y16 或 age 列中的任何一个中。然后，我在age中找到了NAs

> str(df46$age)

 dbl+lbl [1:3034] 30, 62, 25, 38, 24, 76, 39, 16, 71, 62, 29, 27, 60, 41, 22, 20, NA, 5...
 @ labels: Named num [1:4] NA 888 988 0
  ..- attr(*, "names")= chr [1:4] "Don't Know" "ns" "nr" "No sabe/No responde"
 @ label : chr "Age"

然后，我检查了age 列中有多少个NAs 以及它们的位置。

which(is.na(df46$age))

[1]   17   28  802  888 1045 2401 2898

有 7 个NAs。我怀疑这个数字与错误消息中的数字有关，因为df46 中有 3034 行。将数字加倍，得到 6068。将 NAs 的数量加倍，得到 14，6068-14 = 6054，即错误消息中显示的确切数字。

然后，我尝试排除df46中的七行以获得完整的案例，并使用完整的案例创建lapop和fit。

ind = which(is.na(df46$age))
df46_complete = df46[-ind,]
lapop<-svydesign(ids = ~ upm, 
                 strata = ~ estratopri, 
                 weights = ~ weight1500, 
                 nest = T,
                 data = df46_complete)
fit<-svyglm(ctol ~ y16 + age,
            design = lapop,
            family = quasibinomial(link = 'logit'))

最后，当我运行margins()时没有弹出错误：

margins(fit, design = lapop) %>% summary()

# factor     AME     SE       z      p   lower   upper
#    age -0.0026 0.0004 -6.0633 0.0000 -0.0035 -0.0018
#    y16  0.1323 0.0187  7.0638 0.0000  0.0962  0.1696

【讨论】：

非常感谢您的回答！我想投票，但它不会让我，我是 Stack Overflow 的新手。还有一个问题，这是否意味着对于我运行的每个模型，我都必须删除所有缺失值？这似乎也只是调查权重的问题，并非没有。再次感谢您！
您可以签入选项（'na.action'）。如果它显示“na.omit”，那么您很可能需要在运行“glm”建模之前忽略缺失值的行。我试图将“na.omit”选项更改为“na.pass”并运行您的模型，但此选项导致了另一种类型的错误。因此，对于您的模型，省略 NA 似乎是唯一可行的方法。
我为你投了赞成票，@DanielSanchéz。我确实想指出答案通常是不正确的，尽管在这种情况下可以。问题是缺少数据，但要获得正确的标准错误，您应该在声明调查设计之后，而不是之前，删除缺少的观察结果。
谢谢@ThomasLumley，但是我发现使用边距包的修改版本（即我的回答）要实用得多。到目前为止，我一直在关注另一个答案，但是为每个模型创建新的数据框是非常不切实际的。我什至没有想到这个问题的解决方案应该反过来。谢谢！