---
title: 'Regression III: Lab 3 Answers'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Missing Data and Multiple Imputation

I've got two examples here - one purely cross-sectional and one time-series cross-sectional.  The TSCS one is a study by David Rueda (2008), which was replicated by Lall in his 2016 piece.  The second is from the 2014 American National Election Study. 

### Rudea Replication 

For the sake of this exercise, we're going to replicate one of Rueda's models and assess the quality of the imputations.  You can read in and manage the data like below.  This puts an object named `x` in your workspace that contains the data. 

```{r}
f <- file("https://quantoid.net/files/reg3/rueda.rda")
load(f)
close(f)
tmp <- x[,c("country", "year", "govem", "govpart", "hkcorp", 
"govpxcor", "open", "finop", "dccggx", "unemp", 
"gdpgr", "deca70", "deca80", "deca85", "deca90", "count")]
tmp$count <- as.factor(tmp$count)
```

Next we can make the listwise deleted data.  We can do this as follows: 

```{r}
tmpNA <- na.omit(tmp[,c("country", "year", "govem", "govpart", "hkcorp", 
"govpxcor", "open", "finop", "dccggx", "unemp", 
"gdpgr", "deca70", "deca80", "deca85", "deca90", "count")])
```

  1. Estimate the model of `govem` (government employment) on `govpart` (cabinet partisanship), `hkcorp` (corporatism), `govpart*hkcorp` (the interaction of cabinet partisanship and  corporatism),  `open` (internal openness), `finop` (financial openness), `dccggx` (government debt), `unemp` (unemployment), `gdpgr` (gdp growth), `deca70`, `deca80`, `deca85`, `deca90` and `count` (year and country variables).  

```{r}
mod <- lm(govem ~ govpart*hkcorp  + open + finop + dccggx + 
          unemp + gdpgr + deca70 + deca80 + deca85 + deca90 + count, data=tmpNA)
summary(mod)
```

Rueda used panel corrected standard errors so, I'll show you how you can do that in parallel with the original analysis. 

```{r}
pmod <- pcse(mod, tmpNA$count, tmpNA$year)
summary(pmod)
```


  2. Impute the data three different ways - `mice` with the default settings, `mice` using random forests (make sure that the time and year variables are in the model) as the imputation method and `amelia`.  After doing that, re-estimate the models with the three different imputed datasets.  How do the three models compare? 

```{r}
library(mice)
library(Amelia)
library(mitools)
mi.imp <- mice(tmp[,c("govem", "govpart", "hkcorp", "govpxcor", "open", "finop", "dccggx", "unemp", "gdpgr", "deca70", "deca80", "deca85", "deca90", "count")], m=50, printFlag=FALSE)
```

I talked in class about a way to make sure that transformations are respected in imputations.  Imagine we wanted to use the log of `dccggx` and `unemp` in the imputation models. We also want to make sure that the interaction gets imputed the right way. Imputing these variables is done through what is called 'passive imputation'.  Here, we set the imputation method to the transformation we want to execute.  We could do that with the following: 

```{r}
tmp$logunemp <- log(tmp$unemp)
tmp$logdc <- log(tmp$dccggx)
## run mice one time to get the imputation methods for everything
tmp.imp <- mice(tmp[,c("govem", "govpart", "hkcorp", "govpxcor", 
        "open", "finop", "dccggx", "unemp", "gdpgr", "logunemp",
        "logdc","deca70", "deca80", "deca85", "deca90", "year",
         "count")], m=1, printFlag=FALSE)
pm <- tmp.imp$predictorMatrix
pm[,12:15] <- 0
imp.meth <- tmp.imp$method
imp.meth["logunemp"] <- "~log(unemp)"
imp.meth["logdc"] <- "~log(dccggx)"
imp.meth["govpxcor"] <- "~I(govpart*hkcorp)"
mi.imp2 <- mice(tmp[,c("govem", "govpart", "hkcorp", "govpxcor", 
        "open", "finop", "dccggx", "unemp", "gdpgr", "logunemp",
        "logdc","deca70", "deca80", "deca85", "deca90", 
        "year", "count")], 
        m=50, predictorMatrix=pm, printFlag=FALSE)
```

Next, we could do the random forests imputation.  For this we would have to make sure to add the year in, and we can put it in as a factor even if we want
```{r}
# find variables with missing data
tmp$year <- as.factor(tmp$year)
tmp.imp <- mice(tmp[,c("govem", "govpart", "hkcorp", "govpxcor", "open",
        "finop", "dccggx", "unemp", "gdpgr", "year", "count", 
        "deca70", "deca80", "deca85", "deca90")], 
        m=1, printFlag=FALSE)
## since we've got all of the variables in propmiss, we simply want to set the 
## imputation method for anything that has missing info to 'rf'
imp.methods <- tmp.imp$method
imp.methods <- gsub("pmm", "rf", imp.methods)
imp.methods["govpxcor"] <- "~I(govpart*hkcorp)"
pm <- tmp.imp$predictorMatrix
pm[,12:15] <- 0

rf.imp <- mice(tmp[,c("govem", "govpart", "hkcorp", "govpxcor", "open",
        "finop", "dccggx", "unemp", "gdpgr", "year", "count", 
        "deca70", "deca80", "deca85", "deca90")], 
        methods=imp.methods, m=50, printFlag=FALSE, 
        predictorMatrix=pm)
```

Finally, we can use `amelia`

```{r}
tmp.am <- x[,c("govem", "govpart", "hkcorp", "open", "finop",         
        "dccggx", "unemp", "gdpgr", "deca70", "deca80", "deca85",
        "deca90", "year", "count")]
am.imp <- amelia(tmp.am, cs = "count", ts="year", 
        logs = c("dccggx", "unemp"), polytime=3, intercs=F, m=50, p2s=0, 
        idvars=c("deca70", "deca80", "deca85", "deca90"))
```

Now, we can re-estimate the models.  I'm going to use `MIcombine` for all of them, but you could use `zelig` or `lm.mids`, too.  

```{r}
# make the default mice datasets
mi.dats <- lapply(1:50, function(i)complete(mi.imp2, i))
rf.dats <- lapply(1:50, function(i)complete(rf.imp, i))
am.dats <- am.imp$imputations

mi.mods <- MIcombine(lapply(mi.dats, function(x)update(mod, data=x)))
rf.mods <- MIcombine(lapply(rf.dats, function(x)update(mod, data=x)))
am.mods <- MIcombine(lapply(am.dats, function(x)update(mod, data=x)))
```

If we wanted to do this with the panel corrected standard errors, it's a touch more complicated
```{r}
## first, we need to get the coefficients from all of the models
mi.b <- lapply(mi.dats, function(x)update(mod, data=x)$coef)
rf.b <- lapply(rf.dats, function(x)update(mod, data=x)$coef)
am.b <- lapply(am.dats, function(x)update(mod, data=x)$coef)

## then we need to get the PC variane-covariande matrices
mi.v <- lapply(1:50, function(x)pcse(update(mod, data=mi.dats[[x]]),
        mi.dats[[x]]$count, mi.dats[[x]]$year)$vcov)
rf.v <- lapply(1:50, function(x)pcse(update(mod, data=rf.dats[[x]]),
        rf.dats[[x]]$count, rf.dats[[x]]$year)$vcov)
am.v <- lapply(1:50, function(x)pcse(update(mod, data=am.dats[[x]]),
        am.dats[[x]]$count, am.dats[[x]]$year)$vcov)

## finally, we can combine the coefficients and variances using Rubin's rules. 
mi.pmods <- MIcombine(mi.b, mi.v)
rf.pmods <- MIcombine(rf.b, rf.v)
am.pmods <- MIcombine(am.b, am.v)
```

Next, we can summarize the models and compare to the original results

```{r}
# summarize the models
summary(mod)
summary(mi.mods)

## make the confidence intervals
ci.lwd <- confint(mod)
ci.mi <- summary(mi.mods)[,3:4]

##identify the significant results 
sig.lwd <- sign(ci.lwd[,1]) == sign(ci.lwd[,2])
sig.mi <- sign(ci.mi[,1]) == sign(ci.mi[,2])

## put the two significant results descriptions together. 
cbind(sig.lwd, sig.mi)
```

If we were going to do this with the panel corrected standard errors:

```{r}
library(lmtest)
ci.lwd <- coefci(mod, vcov.=pmod$vcov)
ci.mi <- summary(mi.pmods)[,3:4]

##identify the significant results 
sig.lwd <- sign(ci.lwd[,1]) == sign(ci.lwd[,2])
sig.mi <- sign(ci.mi[,1]) == sign(ci.mi[,2])

## put the two significant results descriptions together. 
out <- cbind(sig.lwd, sig.mi[-c(1,28)])
rownames(out) <- names(sig.lwd)
out
```

Here, the coefficients of `open` and `finop` both become insignificant in the imputed data model.  We could look at the same thing for the other two models.  Here, I'll just use the PCSE models. 

```{r}
ci.rf <- summary(rf.pmods)[,3:4]
ci.am <- summary(am.pmods)[,3:4]

sig.rf <- sign(ci.rf[,1]) == sign(ci.rf[,2])
sig.am <- sign(ci.am[,1]) == sign(ci.am[,2])

cbind(out, sig.rf[-c(1,28)], sig.am[-c(1,28)])
```

Notice that the largest set of substantive variables is statistically significant using the `amelia` imputations. 

  3. Do the posterior predictive checks to check the quality of the imputations using the original `mice` method (with the default imputation methods) and with `amelia` three methods.  Is there any difference in what you find across the methods?  Pick one missing data point from each variable with missing data and plot out the posterior distribution of the missing observations for the three different methods. 

First, we can do this for `mice`.  There's an easier way of doing this than the way we talked about in class, at least with the `mice` object:

```{r}
## identify the variables with missing observations
missvars <- c("govem", "dccggx", "unemp")

# make the padded out data frame: 
w <- is.na(tmp[,c("govem", "govpart", "hkcorp", 
"govpxcor", "open", "finop", "dccggx", "unemp", "gdpgr", "logunemp", "logdc",
"deca70", "deca80", "deca85", "deca90", "year", "count")])
w[,missvars] <- TRUE

mi.ppd <- update(mi.imp2, where=w, m=500, printFlag=F)
mi.comp <- lapply(1:500, function(x)complete(mi.ppd, i))
```

Now, we can make the ppds for `amelia`

```{r}
newdat <- tmp[,c("govem", "govpart", "hkcorp", 
"govpxcor", "open", "finop", "dccggx", "unemp", "gdpgr", "logunemp", "logdc",
"deca70", "deca80", "deca85", "deca90", "year", "count")]
for(i in 1:length(missvars)){
        tmp.dat <- tmp[,c("govem", "govpart", "hkcorp", 
                "govpxcor", "open", "finop", "dccggx", "unemp", "gdpgr", "logunemp", "logdc",
                "deca70", "deca80", "deca85", "deca90", "year", "count")]
        tmp.dat[[missvars[i]]] <- NA
        newdat <- rbind(newdat, tmp.dat)
}
am.ppd <- amelia(newdat, cs = "count", ts="year", 
        logs = c("dccggx", "unemp"), polytime=3, intercs=F, m=500, p2s=0, 
        idvars=c("deca70", "deca80", "deca85", "deca90"))
d1 <- 1:368
d2 <- 369:736
d3 <- 737:1104
d4 <- 1105:1472

am.comp <- list()
for(i in 1:500){
        tmp.comp <- am.ppd$imputations[[i]][d1, ]
        tmp.comp[[missvars[1]]] <- am.ppd$imputations[[i]][d2, missvars[1]]
        tmp.comp[[missvars[2]]] <- am.ppd$imputations[[i]][d3, missvars[2]]
        tmp.comp[[missvars[3]]] <- am.ppd$imputations[[i]][d4, missvars[3]]
        am.comp[[i]] <- tmp.comp   
}
```

Next, we can estimate the models. 

```{r}
mi.orig <- update(mi.imp2, m=500)
mi.orig.comp <- lapply(1:500, function(i)complete(mi.orig, i))
l1 <- lapply(mi.orig.comp, function(x)update(mod, data=x))
l2 <- lapply(mi.comp, function(x)update(mod, data=x))
## save the coefficients from each model
b1 <- sapply(l1, coef)
b2 <- sapply(l2, coef)
## fin pval as the proportion of times that b1 > b2
diff <- t(b1) - t(b2)
pval <- apply(diff, 2, function(x)mean(x > 0))
## replace one-sided p-value with two-sided p-value
pval <- 2*ifelse(pval > .5, 1-pval, pval)
pval
```

Next, we could do the omnibus test

```{r}
coef <- colMeans(diff)[c(1:8,28)]
v <- var(diff[,c(1:8,28)])
linearHypothesis(model=list(df.residual=NULL),
  hypothesis.matrix=diag(length(coef)),
  vcov.=v, coef.=coef, rhs=rep(0, length(coef)))
```

Now, we could do the same thing for `amelia`

```{r}
am.orig.comp <- lapply(1:500, function(i)am.ppd$imputations[[i]][d1, ])
l1 <- lapply(am.orig.comp, function(x)update(mod, data=x))
l2 <- lapply(am.comp, function(x)update(mod, data=x))
## save the coefficients from each model
b1 <- sapply(l1, coef)
b2 <- sapply(l2, coef)
## fin pval as the proportion of times that b1 > b2
diff <- t(b1) - t(b2)
pval <- apply(diff, 2, function(x)mean(x > 0))
## replace one-sided p-value with two-sided p-value
pval <- 2*ifelse(pval > .5, 1-pval, pval)
pval
```

Notice, we do quite a bit better here - the quality of these imputations is higher. 

```{r}
coef <- colMeans(diff)[c(1:8,28)]
v <- var(diff[,c(1:8,28)])
linearHypothesis(model=list(df.residual=NULL),
  hypothesis.matrix=diag(length(coef)),
  vcov.=v, coef.=coef, rhs=rep(0, length(coef)))
```

The omnibus result, though is the same.  This would probably suggest trying to fit a better imputation model.  The `rf` model might work, but it will take a pretty long time to run.  If you want, try this at home when you've got some time to let your computer run.  

### ANES 2014

Using `lab3b_data.dta`, do multiple imputation (with 5 imputations) on all of the variables in the dataset.  You can download the data with: 

```{r}
library(haven)
dat <- read_dta('http://www.quantoid.net/files/reg3/lab3b_data.dta')
searchVarLabels(dat, "")
dat$dem_edugroup <- as_factor(dat$dem_edugroup)
levels(dat$dem_edugroup) <- c("ltHS", "HS", "Some Coll", "BA", "Grad Deg")
dat$dem_agegrp <- as_factor(dat$dem_agegrp)
dat$dem_agegrp <- car:::recode(dat$dem_agegrp, 
        'c("17-20", "21-24", "25-29") = "17-30"; c("30-34", "35-39", "40-44", "45-49") = "30-49"; 
         c("50-54", "55-59", "60-64", "65-69") = "50-69"; c("70-74", "Over 75") = "70+"', 
        as.factor = TRUE)
dat$gender_respondent <- as_factor(dat$gender_respondent)
dat <- as.data.frame(dat)
```


  1. See how the coefficients change in a model of \verb"libcpre_self" on all the other variables in the dataset from listwise deletion to multiple imputation.  Try `mice`, `mice` with random forest imputation and `amelia` each with 5 imputations.  Are there differences between the models?   

```{r}
mod <- lm(libcpre_self ~ ., data=dat)
summary(mod)
```

Now, we can do the imputations
```{r}
nes.mi <- mice(dat, m=5)
meths <- nes.mi$method
meths <- gsub("[a-z]*", "rf", meths)
nes.rf <- mice(dat, m=5, method=meths)
nes.am <- amelia(dat, m=5, noms=c("dem_edugroup", "dem_agegrp", "gender_respondent", "relig_chmember"))


nes.mi.dats <- lapply(1:5, function(i)complete(nes.mi, i))
nes.rf.dats <- lapply(1:5, function(i)complete(nes.rf, i))
nes.am.dats <- nes.am$imputations

nes.mi.mods <- MIcombine(lapply(nes.mi.dats, function(x)update(mod, data=x)))
nes.rf.mods <- MIcombine(lapply(nes.rf.dats, function(x)update(mod, data=x)))
nes.am.mods <- MIcombine(lapply(nes.am.dats, function(x)update(mod, data=x)))

summary(nes.mi.mods)
summary(nes.rf.mods)
summary(nes.am.mods)

mi.mi <- as.numeric(gsub(" %", "", summary(nes.mi.mods)[,5]))
rf.mi <- as.numeric(gsub(" %", "", summary(nes.rf.mods)[,5]))
am.mi <- as.numeric(gsub(" %", "", summary(nes.am.mods)[,5]))

# show the missing info for the mi and rf models relative to amelia
mi.mi/am.mi
rf.mi/am.mi

mean(mi.mi/am.mi)
mean(rf.mi/am.mi)
```

This suggests that with either of the `mice` impuation methods, the `amelia` imputations generate lower amounts of missing information.  If you looked at `mi.mi/rf.mi` you would find that the average is in favor of the default `mice` model.  

  2. How do the results change if you have use 25 imputations instead of 5? 

Repeat above with `m=25`. 

### Finite Mixture Models

Using the data above, estimate a finite mixture model where you assume that `indsocial` and `indspend` operationalize the two different theories (don't use `relig_chmember` and `inc_incgroup_pre` in the models, yet.)  Evaluate the resulting model against the linear model where there is an interaction between `indsocial` and `indspend` including the other controls. 


  1. Run an OLS regression of `libcpre_self_num` on `indsocial` (a composite of social policy attitudes) and `indspend` (a composite of spending policy attitudes), their interaction and the controls `gender_respondent`, `dem_edugroup` and `dem_agegrp_num`.  

```{r, fig.align='center', fig.height=6, fig.width=9, out.width='100%'}
library(flexmix)
dat <- read.dta('http://www.quantoid.net/files/reg3/lab3_data.dta')
dat$libcpre_self_num <- as.numeric(dat$libcpre_self)
dat$dem_agegrp_num <- as.numeric(dat$dem_agegrp)
mod <- lm(libcpre_self_num ~ indsocial*indspend + gender_respondent + dem_edugroup + dem_agegrp_num, data=dat)
DAintfun2(mod, c("indsocial", "indspend"), hist=T, scale.hist=.3)
```


  2. Estimate a finite mixture model where `indsocial` is the variable of interest in one component and `indspend` is the variable of interest in the other.  That is set the `indspend` coefficient to zero in the model with `indsocial` and `indsocial`'s coefficient to zero in the model with `indspend`.  Fix the coefficients on the other regressors to be constant across the two components. 

```{r}
model <- FLXMRglmfix(family = "gaussian", nested=list(k=c(1,1), 
	formula = c(~indsocial,
		 ~indspend)), fixed= ~ dem_agegrp_num + gender_respondent + dem_edugroup)
out.a <- stepFlexmix(libcpre_self_num ~ 1 , k=2, model=model, data=dat, nrep=20)
mod.refit <- refit(out.a)
```

	
  3. How do the two different models fit?  Which do you think it better?  

```{r}
fit1 <- predict(out.a)$Comp.1[[1]]
fit2 <- predict(out.a)$Comp.2[[1]]
post <- out.a@posterior$scaled
predfix1 <- rowSums(cbind(fit1, fit2)*post)
## Correlations of fitted values and observed values from the mixture model
cor(predfix1, dat$libcpre_self)^2
## Correlation using only the fits from best predicted theory
predfix2 <- fit1
predfix2[which(post[,2] > .5)] <- fit2[which(post[,2] > .5)]
cor(predfix2, dat$libcpre_self)^2
## correlations of fitted values and observed from the linear model 
cor(mod$fitted, dat$libcpre_self)^2
```

  4. Consider that income (an economic predictor operationalized by `inc_incgroup_pre`) and church membership (a social predictor operationalized by `relig_chmember`) might tell us something about the probabilities of being in one or the other group.  Incorporate that information and see whether, in fact, they do provide information about group membership. 

```{r}
out.b <- stepFlexmix(libcpre_self ~ 1 , k=2, model=model, data=dat, nrep=20, 
	concomitant = FLXPmultinom( ~ inc_incgroup_pre + relig_chmember))
out.b.refit <- refit(out.b)
out.b.refit@concomitant
```



