nat_pride_other_data

setwd("~/Downloads/") 
library("foreign")
library("lme4")
library("lattice")
library("sandwich")
library("arm")
library("car")
library("lmtest")

# Load data
base<-read.spss("Inglehart_cut.sav", use.value.labels=TRUE, to.data.frame=TRUE)
# Recode regions
base$RegionsTXTUSA2<-recode(base$Regions, "1='GermanyWest';2='GermanySouth';3='GermanyNorth'; 4='GermanyEast';5='GermanyCenter';6='RussiaPovolzhie';7='RussiaCenter';8='RussiaUral';9='RussiaFarEast';10='RussiaMoskow';11='USANorth';12='USA';13='SpainNonspanish';14='Spainspanish';15='IndiaWest';16='IndiaSouth';17='IndiaCenter';18='IndiaEast';19='IndiaNorth';20='ItalyNorthEast';21='ItalyNorthWest';22='ItalyCenter';23='ItalySouth';24='UkraineEast';25='UkraineSouth';26='UkraineCenter';27='UkraineWest';28='GBEngland';29='GBWales';30='GBScotland'; 31='BrazilNorthwest';32='BrazilNotrheast'; 34='BrazilCenter';35='BrazilSoutheast';36='BrazilSouth'")
# NA for two variables
base$divorced <- ifelse(base$divorced == -99, NA ,base$divorced)
base$children <- ifelse(base$children == -99, NA,base$children)

# Create wave identifier
base$wave<-ifelse(base$wave1==1,1,ifelse(base$wave2==1,2,ifelse(base$wave3==1,3,ifelse(base$wave4==1,4,ifelse(base$wave5==1,5,0)))))

# Subset data, leaving only relevant variables to free memory
basesubset <- base[c("RegionsTXTUSA2","wave","natpride","V235","V238","age","divorced","children")]
base <- NULL

# Fixed effects model (original)
fixed <- lm(natpride ~ age + as.factor(V235) + as.factor(V238) + as.factor(RegionsTXTUSA2) -1, data = basesubset)

# Random effects model (original)
random <- lmer(natpride ~ age + as.factor(V235) + as.factor(V238) + (1|RegionsTXTUSA2), data = basesubset)

# Fixed effects model (with children and divorced)
fixed_family <- lm(natpride ~ age + as.factor(V235) + as.factor(V238) + as.factor(divorced) + as.factor(children) + as.factor(RegionsTXTUSA2) -1, data = basesubset)

# Random effects model (with children and divorced)
random_family <- lmer(natpride ~ age + as.factor(V235) + as.factor(V238) + as.factor(divorced) + as.factor(children) + (1|RegionsTXTUSA2), data = basesubset)

# Huber-Eicker-White Errors 
fixed.se <- coeftest(fixed, vcovHC(fixed,type="HC1"))[,2]
fixed_family.se <- coeftest(fixed_family, vcovHC(fixed_family,type="HC1"))[,2]

# Hausman test function (inspired by p. 8, http://www.kirchkamp.de/oekonometrie/pdf/me-p.pdf)
hausman <- function(fixed,random) {
	rNames <- names(fixef(random))
	fNames <- names(coef(fixed))
	timevarNames <- intersect(rNames,fNames)
	k <- length(timevarNames)
	rV <- vcov(random)
	rownames(rV)=rNames
	colnames(rV)=rNames
	bDiff <- as.matrix((fixef(random))[timevarNames] - coef(fixed)[timevarNames])
	vDiff <- as.matrix(vcov(fixed)[timevarNames,timevarNames] - rV[timevarNames, timevarNames])
	H <- as.numeric(abs(t(bDiff) %*% solve(vDiff) %*% bDiff))
	c(H=H,p.value=pchisq(H,k,lower.tail=FALSE))
}

# Perform Hausman test
hausman(fixed,random)
# Ok in the original specification, can use random effects
hausman(fixed_family,random_family)
# Hausman test fails when children and divorce information is added in the model

# Dot plots
## Random effects
dotplot(ranef(random, condVar=TRUE), binpositions="all",  stackdir = "centerwhole")
dotplot(ranef(random_family, condVar=TRUE), binpositions="all",  stackdir = "centerwhole")

## Fixed effects (with Huber-Eicker-White errors)
end1 <- length(fixed$coefficients)
png("fixed_effects_nat_pride.png",width = 7.5, height = 7.5, units= "in", res = 200)
coefplot(fixed$coefficients[10:end1],fixed.se[10:end1],varnames=fixed$xlevels[[3]],cex.var=0.7,main="Fixed effects estimates on provided data (lines are Huber-Eicker-White 68% and 95% CI)",cex.main=0.9)
dev.off()