################################################################## ################################################################## # Gov 2001/1002/E-2001 Section 12 # 4-23-14 # Soledad Artiz Prillaman ################################################################## ################################################################## # Clear our workspace rm(list=ls()) library(Zelig) library(systemfit) library("ZeligChoice") ################################################### ###### Seemingly Unrelated Regression Models ###### ################################################### # Load the grunfeld data from the Zelig packages data(grunfeld) names(grunfeld) summary(grunfeld) # Run the Seemingly Unrelated Regression for GE investment and # Westinghouse investment # Store our two equations then run model # To run this model we can use the systemfit package formula <- list(mu1 = Ige ~ Fge + Cge, mu2 = Iw ~ Fw + Cw) sur.out <- systemfit(formula = formula, method = "SUR", data = grunfeld) sur.out # To get the variance - covariance matrix vcov.sur<-vcov(sur.out) vcov.sur # We can see that there are dependencies in our parameters across models ####################################### ###### Multinomial Choice Models ###### ####################################### # Load the mexico data from the Zelig package # This has data on vote choice across three parties in Mexico data(mexico) names(mexico) summary(mexico\$vote88) # Let's run a multinomial logit model of voting behavior # Our covariates will be age and gender # Since vote choice is a categorical variable, we need to factorize it # to get a matrix of y's ml.out <- zelig(as.factor(vote88) ~ age + female, model = "mlogit", data = mexico) ml.out\$result #Let's look at the difference in expected vote choice for the youngest and the oldest x.young <- setx(ml.out, age = min(mexico\$age)) x.old<- setx(ml.out, age = max(mexico\$age)) ml.sim <- sim(ml.out, x1 = x.old, x = x.young) summary(ml.sim) # Let's look at our results plot(ml.sim) ########################## ###### Missing Data ###### ########################## # Let's load the Amelia package and the africa data in this package library(Amelia) data(africa) names(africa) summary(africa) # we have a couple of missing values in gdp pc and trade africa[1,] year country gdp_pc infl trade civlib population 1 1972 Burkina Faso 377 -2.92 29.69 0.5 5848380 # Let's run our imputation # x specifies the data to be imputed # all variables in x will be used as predictors in the imputation model # cs and ts specify if cross-sectional or time-series # logs specifies variables which need to be logged # m specifies the number of imputed data sets # the greater m is, the less likely you are to have outliers set.seed(1234) a.out <- amelia(x = africa, cs = "country", ts = "year", logs = "gdp_pc", m = 5) # If you wanted to save these 5 data sets, you could: write.amelia(obj=a.out, file.stem = "a.out") # Now, we just want to estimate a basic regression model with our imputed data # But we have 5 datasets! # Zelif will combine them for us: z.out.imp <- zelig(trade ~ log(population) + log(gdp_pc) + infl + civlib, data = a.out\$imputations, model = "ls") summary(z.out.imp) ### Diagnostics # The missingness map shows us where our missingness is missmap(a.out) # We can plot our amelia output to compare empirical and imputed densities plot(a.out) # There are time-series, cross-sectional plots with imputed values tscsPlot(a.out, var = "trade", cs = "Cameroon") #Overimputation for a specific variable tests the imputation model by # imagining that each observation is missing and generating some imputations to check performance. overimpute(a.out, var = "trade") #The disperse function starts the algorithm at some unlikely values to check # that amelia hasn't found a local rather a global maximum for the likelihood of the complete data. disperse(a.out, dims = 1, m = 5)