用 R 报告调查加权逻辑模型的平均边际效应

Question

我正在使用复杂样本的调查数据来估计二元结果模型。我正在尝试报告 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。

Answer 1

错误信息表明在某些行中有 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 包含 NA。然后我检查了 fit 在你的代码中是如何产生的：

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

NA 应该在 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 列中有多少 NA 以及它们所在的位置。

which(is.na(df46$age))

[1]   17   28  802  888 1045 2401 2898

有 7 个 NA。我怀疑这个数字与错误消息中的数字有关，因为 df46 中有 3034 行。将数字加倍，得到 6068。将 NA 的数目加倍，得到 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

Answer 2

我发现这是怎么回事：原来是包代码中的某种错误。具体可以看一下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 足以避免使用调查加权模型产生此错误。

用 R 报告调查加权逻辑模型的平均边际效应

Reporting average marginal effects of a survey-weighted logit model with R

r

survey

logistic-regression