###########################################################################
# PUBLG100: Introduction to Quantitative Methods
#
# Week 8 Seminar: Panel Data, Time-Series Cross-Section Models
#
#

# Change your working directory
setwd("N:/PUBLG100")

# Check your working directory
getwd()


## ------------------------------------------------------------------------
install.packages("plm")
install.packages("tidyr")

## ----message = FALSE-----------------------------------------------------
library(plm)
library(lmtest)
library(texreg)
library(tidyr)
library(dplyr)

## ------------------------------------------------------------------------
# clear environment
rm(list = ls())

## ------------------------------------------------------------------------
wdi_original <- read.csv("http://uclspp.github.io/PUBLG100/data/WDI_Data.csv", na.strings = "..")

## ----eval = FALSE--------------------------------------------------------
## View(wdi_original)

## ------------------------------------------------------------------------
wdi <- wdi_original %>%
  select(-1) %>% 
  filter(Country.Code != "") %>%
  gather(Year, Series.Value, X1960..YR1960.:X2015..YR2015.) %>%
  spread(Series.Code, Series.Value) %>%
  mutate(Year = as.numeric(substring(Year, 2, 5)))

## ------------------------------------------------------------------------
wdi <- rename(wdi,
              CountryName = Country.Name,
              CountryCode = Country.Code,
              SafeWaterAccess = SH.H2O.SAFE.ZS,
              MaternalMortality = SH.STA.MMRT,
              PregnantWomenWithAnemia = SH.PRG.ANEM,
              HealthExpenditure = SH.XPD.TOTL.ZS)

## ----eval = FALSE--------------------------------------------------------
## View(wdi)

## ------------------------------------------------------------------------
fixed_effects <- plm(MaternalMortality ~ SafeWaterAccess + HealthExpenditure + PregnantWomenWithAnemia, 
                     data = wdi, 
                     index = c("CountryCode", "Year"), 
                     model = "within", 
                     effect = "individual")
summary(fixed_effects)

## ------------------------------------------------------------------------
plmtest(fixed_effects, effect="individual")

## ------------------------------------------------------------------------
time_effects <- plm(MaternalMortality ~ SafeWaterAccess + HealthExpenditure + PregnantWomenWithAnemia, 
                    data = wdi, 
                    index = c("CountryCode", "Year"), 
                    model = "within", 
                    effect = "time")
summary(time_effects)

## ------------------------------------------------------------------------
plmtest(time_effects, effect="time")

## ------------------------------------------------------------------------
twoway_effects <- plm(MaternalMortality ~ SafeWaterAccess + HealthExpenditure + PregnantWomenWithAnemia, 
                      data = wdi, 
                      index = c("CountryCode", "Year"), 
                      model = "within", 
                      effect = "twoways")
summary(twoway_effects)

## ------------------------------------------------------------------------
screenreg(list(fixed_effects, time_effects, twoway_effects), 
          custom.model.names = c("Country Fixed Effects", "Time Fixed Effects", "Twoway Fixed Effects"))

## ------------------------------------------------------------------------
# clear the environment first
rm(list = ls())

guns_data <- read.csv("http://uclspp.github.io/PUBLG100/data/guns.csv")

## ------------------------------------------------------------------------
guns_data$shall_law <- factor(guns_data$shall, levels = c(0, 1), labels =c("NO", "YES"))

## ------------------------------------------------------------------------
fixed_effects <- 
  plm(mur ~ shall_law + incarc_rate + pm1029, 
      data = guns_data, 
      index = c("stateid", "year"), 
      model = "within", 
      effect = "individual")

summary(fixed_effects)

## ------------------------------------------------------------------------
plmtest(fixed_effects, effect="individual")

## ------------------------------------------------------------------------
twoway_effects <- 
  plm(mur ~ shall_law + incarc_rate + pm1029, 
      data = guns_data, 
      index = c("stateid", "year"), 
      model = "within", 
      effect = "twoways")

summary(twoway_effects)

## ------------------------------------------------------------------------
pbgtest(twoway_effects)

## ------------------------------------------------------------------------
twoway_effects_hac <- coeftest(twoway_effects, 
                               vcov = vcovHC(twoway_effects, method = "arellano", type = "HC3"))

screenreg(list(twoway_effects, twoway_effects_hac),
          custom.model.names = c("Twoway Fixed Effects", "Twoway Fixed Effects (HAC)"))

## ------------------------------------------------------------------------
ldv_model <- 
  plm(mur ~ lag(mur) + shall_law + incarc_rate + pm1029, 
      data = guns_data, 
      index = c("stateid", "year"), 
      model = "within", 
      effect = "twoways")

summary(ldv_model)

## ------------------------------------------------------------------------
pcdtest(twoway_effects)

## ------------------------------------------------------------------------
twoway_effects_pcse <- coeftest(twoway_effects, 
                                vcov = vcovBK(twoway_effects, type="HC3", cluster = "group")) 

twoway_effects_pcse

## ------------------------------------------------------------------------
twoway_effects_scc <- coeftest(twoway_effects,
                               vcov = vcovSCC(twoway_effects, type="HC3", cluster = "group"))

twoway_effects_scc

## ------------------------------------------------------------------------
screenreg(list(fixed_effects, twoway_effects, ldv_model, twoway_effects_pcse, twoway_effects_scc), 
          custom.model.names = c("Country Effects", "Twoway Fixed Effects", "LDV", "PCSE", "SCC"))