Получить p-значения из результатов svyglm при использовании нескольких вменений в R

Question

Я хотел бы получить р-значения из результатов svyglm модель при использовании нескольких вменений. Воспроизводимый пример приведен ниже.

Создать наборы данных

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

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