在使用R进行多重插补时，从svyglm结果获取p值

Question

我希望从使用多重插补的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

- scottsmith

Answer 1

如果你需要 $p$ 值，那么你需要计算它们。所有必要的信息都包含在 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 

- Thomas Lumley