在 R 中使用多重插补时从 svyglm 的结果中获取 p 值答案

【问题标题】：Get p-values from results of svyglm when using multiple imputations in R在 R 中使用多重插补时从 svyglm 的结果中获取 p 值
【发布时间】：2022-03-18 08:10:31
【问题描述】：

我想在使用多重插补时从svyglm 模型的结果中获取 p 值。下面是一个可重现的示例。

创建数据集

library(tibble)
library(survey)
library(mitools)

# Data set 1
# Note that I am excluding the "income" variable from the "df"s and creating  
# it separately so that it varies between the data sets. This simulates the 
# variation with multiple imputation. Since I am using the same seed
# (i.e., 123), all the other variables will be the same, the only one that 
# will vary will be "income."

set.seed(123)

df1 <- tibble(id      = seq(1, 100, by = 1),
              gender  = as.factor(rbinom(n = 100, size = 1, prob = 0.50)),
              working = as.factor(rbinom(n = 100, size = 1, prob = 0.40)),
              pweight = sample(50:500, 100,  replace   = TRUE))


# Create random income variable.

set.seed(456)

income <- tibble(income = sample(0:100000, 100))

# Bind it to df1

df1 <- cbind(df1, income)


# Data set 2

set.seed(123)

df2 <- tibble(id      = seq(1, 100, by = 1),
              gender  = as.factor(rbinom(n = 100, size = 1, prob = 0.50)),
              working = as.factor(rbinom(n = 100, size = 1, prob = 0.40)),
              pweight = sample(50:500, 100,  replace   = TRUE))

set.seed(789)

income <- tibble(income = sample(0:100000, 100))

df2 <- cbind(df2, income)


# Data set 3

set.seed(123)

df3 <- tibble(id      = seq(1, 100, by = 1),
              gender  = as.factor(rbinom(n = 100, size = 1, prob = 0.50)),
              working = as.factor(rbinom(n = 100, size = 1, prob = 0.40)),
              pweight = sample(50:500, 100,  replace   = TRUE))

set.seed(101)

income <- tibble(income = sample(0:100000, 100))

df3 <- cbind(df3, income)

应用权重并运行模型

# Apply weights via svydesign

imputation <- svydesign(id      = ~id,
                        weights = ~pweight,
                        data    = imputationList(list(df1, 
                                                      df2, 
                                                      df3)))


# Logit model with weights and imputations

logitImp <- with(imputation, svyglm(gender ~ working + income,
                                    family = binomial()))


# Combine results across MI datasets

summary(MIcombine(logitImp))

结果：

                 results           se        (lower       upper) missInfo
(Intercept)  6.824145e-02 9.549646e-01 -2.573937e+00 2.710420e+00     79 %
working1    -5.468836e-02 4.721469e-01 -9.800804e-01 8.707037e-01      0 %
income      -5.776083e-06 1.764326e-05 -5.984829e-05 4.829612e-05     86 %

有没有办法添加 p 值或从这个输出中计算它们？我不太熟悉使用汇总数据计算 p 值。

我想要你通常只使用svyglm 得到的输出。例如，如果我只使用上面的df1，我会得到：

df1Design <- svydesign(id      = ~id,
                       weights = ~pweight,
                       data    = df1)


df1Logit <- svyglm(gender ~ working + income,
                   family = binomial(),
                   data = df1,
                   design = df1Design)

summary(df1Logit)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.226e-01  5.178e-01  -1.395    0.166
working1    -4.428e-02  4.561e-01  -0.097    0.923
income       9.834e-06  8.079e-06   1.217    0.226

【问题讨论】：

标签： r survey imputation

【解决方案1】：

如果需要，您需要计算 $p$-values。您需要的所有部分都在MIcombine返回的对象中

比如运行example(MIcombine)之后

> a<-MIcombine(models)
> coef(a)
(Intercept)        wave         sex    wave:sex 
-2.25974358  0.24055250  0.64905222 -0.03725422 
> SE(a)
(Intercept)        wave         sex    wave:sex 
 0.26830731  0.06587423  0.34919264  0.08609199 
> coef(a)/SE(a)
(Intercept)        wave         sex    wave:sex 
 -8.4222216   3.6516935   1.8587225  -0.4327257 
> pt( abs(coef(a)/SE(a)),df=a$df,lower.tail=FALSE)*2
 (Intercept)         wave          sex     wave:sex 
5.874387e-17 3.067194e-04 6.307325e-02 6.653235e-01

【讨论】：