##################################
# Import data and only keep 1995 #
##################################

library(rms)
library(sandwich)
library(lmtest)
library(car)

data <- stata.get("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/CourseBios312/salary.dta")
data95 <- subset(data, year=="95")



**************************************************
* Convert string variables to numeric indicators *
**************************************************


data95$male <- as.numeric(data95$sex=="M")
data95$arts <- as.numeric(data95$field=="Arts")
data95$other <- as.numeric(data95$field=="Other")

with(data95, table(male, sex))
with(data95, table(arts, field))
with(data95, table(other, field))

data95$malearts <- data95$male*data95$arts
data95$maleother <- data95$male*data95$other

#### The regression model

m1 <- lm(salary ~ male + arts + other + malearts + maleother, data=data95)
coeftest(m1, vcov=sandwich)

# Compare to using robcov() with Deisgn() package to generate
#   robust standard errors
m2 <- ols(salary ~ male + arts + other + malearts + maleother, data=data95, x=T)
robcov(m2)

# In females, arts compared to professional
waldtest(m1, "arts", vcov=sandwich)

* In females, other compared to professional
waldtest(m1, "other", vcov=sandwich)

* In females, arts compared to other
linearHypothesis(m1, "arts-other=0", vcov=sandwich)

* Field effect in females
waldtest(m1, c("arts","other"), vcov=sandwich)

* In males, arts compared to professional
linearHypothesis(m1, "arts+malearts=0", vcov=sandwich)

* In males, other compared to professional
linearHypothesis(m1, "other+maleother=0", vcov=sandwich)

* In males, arts compared to other
linearHypothesis(m1, "arts+malearts=other+maleother", vcov=sandwich)


* Field effect in males
linearHypothesis(m1, c("arts+malearts=0", "other+maleother=0"), vcov=sandwich)

* Salary by gender, contolling for field
waldtest(m1, c("male","malearts","maleother"), vcov=sandwich)

* Salary by field, controlling for gender
waldtest(m1, c("arts","other","malearts","maleother"), vcov=sandwich)

* Effect modifcation (both questions)
waldtest(m1, c("malearts","maleother"), vcov=sandwich)

linearHypothesis(m1, c("malearts=maleother"), vcov=sandwich)