I would like to get p-values from the results of a svyglm model when using multiple imputations. A reproducible example is below.
Create data sets
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 and run model
# 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:
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 %
Is there a way to add p-values or a way to calculate them from this output? I am less familiar calculating p-values with pooled data.
I would like output that you normally get just using svyglm. For example, if I just use df1 from above, I get:
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