This document will walk you
through some of the methods you could use to generate pooled model
results that account for both sampling variability and across imputation
variability. The package hot.deck does not come with a set
of functions to do inference, so we will show you how you could use the
data generated by hot.deck in combination with
glm.mids (and similarly lm.mids) from the
mice package, zelig from the
Zelig package and by using MIcombine from the
mitools package on a list of model objects.
The data we will use come from Poe et al.
(1999) dealing with democracy and state repression. First we need
to call the hot.deck routine on the dataset.
library(hot.deck)
data(isq99)
out <- hot.deck(isq99, sdCutoff=3, IDvars = c("IDORIGIN", "YEAR"))
#> Warning in hot.deck(isq99, sdCutoff = 3, IDvars = c("IDORIGIN", "YEAR")): 52 observations with no observed data. These observations were removed
#> Warning in hot.deck(isq99, sdCutoff = 3, IDvars = c("IDORIGIN", "YEAR")): 45 of 4661 imputations with # donors < 5, consider increasing sdCutoff or using method='p.draw'This shows us that there are still 45 observations with fewer than 5
donors. Using a different method or further widening the
sdCutoff parameter may alleviate the problem. If you want
to see the frequency distribution of the number of donors, you could
look at:
numdonors <- sapply(out$donors, length)
numdonors <- sapply(out$donors, length)
numdonors <- ifelse(numdonors > 5, 6, numdonors)
numdonors <- factor(numdonors, levels=1:6, labels=c(1:5, ">5"))
table(numdonors)
#> numdonors
#> 1 2 3 4 5 >5
#> 18 10 11 6 20 4596Before running a model, three variables have to be created from those
existing. Generally, if variables are deterministic functions of other
variables (e.g., transformations, lags, etc…) it is advisable to impute
the constituent variables of the calculations and then do the
calculations after the fact. Here, we need to lag the AI
variable and create percentage change variables for both population and
per-capita GNP. First, to create the lag of AI,
PCGNP and LPOP. To do this, we will make a
little function.
tscslag <- function(dat, x, id, time){
obs <- apply(dat[, c(id, time)], 1, paste, collapse=".")
tm1 <- dat[[time]] - 1
lagobs <- apply(cbind(dat[[id]], tm1), 1, paste, collapse=".")
lagx <- dat[match(lagobs, obs), x]
}
for(i in 1:length(out$data)){
out$data[[i]]$lagAI <- tscslag(out$data[[i]], "AI", "IDORIGIN", "YEAR")
out$data[[i]]$lagPCGNP <- tscslag(out$data[[i]], "PCGNP", "IDORIGIN", "YEAR")
out$data[[i]]$lagLPOP <- tscslag(out$data[[i]], "LPOP", "IDORIGIN", "YEAR")
}Now, we can use the lagged values of PCGNP and
LPOP, to create percentage change variables:
for(i in 1:length(out$data)){
out$data[[i]]$pctchgPCGNP <- with(out$data[[i]], c(PCGNP-lagPCGNP)/lagPCGNP)
out$data[[i]]$pctchgLPOP <- with(out$data[[i]], c(LPOP-lagLPOP)/lagLPOP)
}You can use the MIcombine command from the
mitools package to generate inferences, too. Here, you have
to produce a list of model estimates and the function will combine
across the different results.
# initialize list
out <- hd2amelia(out)
results <- list()
# loop over imputed datasets
for(i in 1:length(out$imputations)){
results[[i]] <- lm(AI ~ lagAI + pctchgPCGNP + PCGNP + pctchgLPOP + LPOP + MIL2 + LEFT +
BRIT + POLRT + CWARCOW + IWARCOW2, data=out$imputations[[i]])
}
summary(mitools::MIcombine(results))
#> Multiple imputation results:
#> MIcombine.default(results)
#> results se (lower upper) missInfo
#> (Intercept) 5.742000e-01 1.392579e-01 2.988367e-01 0.8495632948 18 %
#> lagAI 4.588350e-01 1.942672e-02 4.193226e-01 0.4983474782 38 %
#> pctchgPCGNP 6.443740e-03 3.801600e-03 -1.420511e-03 0.0143079916 46 %
#> PCGNP -2.029774e-05 2.883881e-06 -2.596168e-05 -0.0000146338 9 %
#> pctchgLPOP -3.345152e-01 9.808632e-01 -2.419391e+00 1.7503609836 56 %
#> LPOP 7.244757e-02 8.299065e-03 5.616535e-02 0.0887297991 6 %
#> MIL2 1.224411e-01 4.461532e-02 3.009314e-02 0.2147890798 46 %
#> LEFT -1.590016e-01 4.503555e-02 -2.481225e-01 -0.0698808082 19 %
#> BRIT -1.237613e-01 3.119437e-02 -1.849958e-01 -0.0625267555 7 %
#> POLRT -7.409579e-02 8.942919e-03 -9.181989e-02 -0.0563716977 21 %
#> CWARCOW 6.352371e-01 6.190189e-02 5.090329e-01 0.7614413416 40 %
#> IWARCOW2 2.028580e-01 6.216681e-02 7.683170e-02 0.3288842201 37 %The final method for combining results is to convert the data object
returned by the hot.deck function to an object of class
mids. This can be done with the datalist2mids
function from the miceadds package.
out.mids <- miceadds::datalist2mids(out$imputations)
s <- summary(mice::pool(mice::lm.mids(AI ~ lagAI + pctchgPCGNP + PCGNP + pctchgLPOP + LPOP + MIL2 + LEFT +
BRIT + POLRT + CWARCOW + IWARCOW2, data=out.mids)))
#> Warning in mice::lm.mids(AI ~ lagAI + pctchgPCGNP + PCGNP + pctchgLPOP + : Use with(imp,
#> lm(yourmodel)).
print(s, digits=4)
#> term estimate std.error statistic df p.value
#> 1 (Intercept) 5.184e-01 1.520e-01 3.41164 44.31 1.388e-03
#> 2 lagAI 4.563e-01 1.829e-02 24.94770 54.14 2.305e-31
#> 3 pctchgPCGNP 4.366e-03 5.981e-03 0.73003 9.35 4.833e-01
#> 4 PCGNP -2.133e-05 2.934e-06 -7.27100 177.89 1.090e-11
#> 5 pctchgLPOP -1.325e-02 8.841e-01 -0.01499 17.66 9.882e-01
#> 6 LPOP 7.593e-02 9.271e-03 8.19073 61.53 1.946e-11
#> 7 MIL2 1.189e-01 4.305e-02 2.76173 28.27 9.990e-03
#> 8 LEFT -1.561e-01 4.971e-02 -3.13985 37.96 3.267e-03
#> 9 BRIT -1.191e-01 3.301e-02 -3.60727 138.01 4.316e-04
#> 10 POLRT -7.289e-02 9.105e-03 -8.00628 82.69 6.461e-12
#> 11 CWARCOW 6.349e-01 5.871e-02 10.81350 48.84 1.469e-14
#> 12 IWARCOW2 2.012e-01 5.635e-02 3.57049 114.96 5.213e-04