# make a model of (neg) in Philadelphia
# Kyle Gorman <kgorman.ling.upenn.edu>

# libraries
library(arm)
library(Design)
library(languageR)

# read in data
dat = read.csv('lcv_neg.csv', header=1)

# look at collinearities
cor(dat$Occ, dat$Res)
cor(dat$Occ, dat$Sc1)
cor(dat$Occ, dat$Sc2)
cor(dat$Occ, dat$Mo)
cor(dat$Res, dat$Sc1)
cor(dat$Res, dat$Sc2)
cor(dat$Res, dat$Mo)
cor(dat$Sc1, dat$Sc2) # huge
cor(dat$Sc1, dat$Mo) 
cor(dat$Sc2, dat$Mo) 
collin.fnc(dat, c('Occ', 'Res'))$cnumber # just an example
vif(glm(Neg ~ Occ + Res, data=dat, family=binomial)) # another one

# plot 'em
pdf('pairs.pdf')
pairs( ~ Occ + Res + Sc1 + Sc2 + Mo, data=dat)
lower.panel=panel.smooth, upper.panel=panel.cor
dev.off()

# do some residualization
dat$rRes = residuals(lm(Res ~ Occ, data=dat))
dat$rSc2 = residuals(lm(Sc2 ~ Occ + rRes, data=dat))
dat$rSc1 = residuals(lm(Sc1 ~ Occ + rRes + rSc2, data=dat))
dat$rMo = residuals(lm(Mo ~ Occ + rRes + rSc2 + rSc1, data=dat))

# standardize
dat$sOcc = (dat$Occ - mean(dat$Occ)) / (2 * sd(dat$Occ))
dat$srRes = (dat$rRes - mean(dat$rRes)) / (2 * sd(dat$rRes))
dat$srSc2 = (dat$rSc2 - mean(dat$rSc2)) / (2 * sd(dat$rSc2))
dat$srSc1 = (dat$rSc1 - mean(dat$rSc1)) / (2 * sd(dat$rSc1))
dat$srMo = (dat$rMo - mean(dat$rMo)) / (2 * sd(dat$rMo))

# fit a flat model
m0 = glm(Neg ~ sOcc + srRes + srSc2 + srSc1 + srMo + Sex + Style + rcs(Age, 4), family=binomial, data=dat)
# bootstrap validation analysis based on Baayen, 2008 p. 212, not shown here.

# fit a hierarchical model
m1 = lmer(Neg ~ sOcc + srRes + srSc2 + srSc1 + srMo + Sex + Style + rcs(Age, 4) + (1 | Name), family=binomial, data=dat)
source('http://ling.upenn.edu/~kgorman/R/ranOutliers.R') # get this if necessary
ranOutliers(m1)
pvals.fnc(m1)

# plot making is easy, see the longer draft