Bayes HLM

library(lme4)
library(blme)
library(arm)


cps <- read.csv('cps.csv')
cps <- subset(cps,wageworkerlastyear==1 & !is.na(lhourlywage))
cps$state<-as.integer(cps$statefip)
cps$nchildnum<-as.integer(cps$nchild)-1
cps$statefip<-as.factor(cps$statefip)

states <- (sort(unique(as.character(cps$statefip))))

#regular OLS model
modelOLS <- lm(lannualhoursworked ~ lhourlywage + lhourlywage*statefip + statefip + lnonlaborincome + sex + race + nchildnum + age + ageage , data = cps)
summary(modelOLS)

modelOLS$coefficients[79:128]
stateCoefficients <- tail(modelOLS$coefficients,50)
stateCoefficients<-setNames(stateCoefficients,states[-1])
wald.test(b=coef(modelOLS),Sigma=vcov(modelOLS),Terms= 79:128)
dotplot(rev(stateCoefficients), xlab = 'State Effects on lhourlywage')


#varying intercepts, no predictor
model0 <- lmer(lannualhoursworked ~ 1 + (1 | statefip ), data = cps)
summary(model0)
coef(model0)
ranef(model0)
confint(model0, level = 0.95, method="Wald")



#varying intercepts, predictor
model1 <- lmer(lannualhoursworked ~  lhourlywage  + lnonlaborincome + sex + race + nchildnum + age + ageage + (1 | statefip ), data = cps)
summary(model1)
coef(model1)
ranef(model1)
confint(model1, level = 0.95, method="Wald")



#same intercepts, varying slopes
model2 <- lmer(lannualhoursworked ~  lhourlywage  + lnonlaborincome + sex + race + nchildnum + age + ageage + (lhourlywage - 1 | statefip ), data = cps)
summary(model2)
coef(model2)
ranef(model2)
confint(model2, level = 0.95, method="Wald")
dotplot(rev(setNames(ranef(model2)$statefip$lhourlywage,states)), xlab = 'State Effects on lhourlywage')


#varying intercepts, varying slopes
model3 <- lmer(lannualhoursworked ~ lhourlywage  + lnonlaborincome + sex + race + nchildnum + age + ageage + (1 + lhourlywage | statefip ), data = cps)
summary(model3)
coef(model3)
ranef(model3)
confint(model3, level = 0.95, method="Wald")
dotplot(ranef(model3))
dotplot(rev(setNames(ranef(model3)$statefip$lhourlywage,states)), xlab = 'State Effects on lhourlywage')


#uniform variance prior, normal mean
model4 <- blmer(lannualhoursworked ~ lhourlywage  + lnonlaborincome + sex + race + nchildnum + age + ageage + (1 + lhourlywage | statefip ),
                data = cps, cov.prior = NULL, fixef.prior = normal)
summary(model4)
coef(model4)
ranef(model4)
confint(model4, level = 0.95, method="Wald")
dotplot(ranef(model4))
dotplot(rev(setNames(ranef(model4)$statefip$lhourlywage,states)), xlab = 'State Effects on lhourlywage')



#uniform variance prior, normal mean
model5 <- blmer(lannualhoursworked ~ lhourlywage  + lnonlaborincome + sex + race + nchildnum + age + ageage + (1 + lhourlywage | statefip ),
                data = cps, cov.prior = invgamma, fixef.prior = normal)

summary(model5)
coef(model5)
ranef(model5)
dotplot(ranef(model5))
dotplot(rev(setNames(ranef(model5)$statefip$lhourlywage,states)), xlab = 'State Effects on lhourlywage')