# Example 4.9, Parents’ Education in a Birth Weight Equation
# Data set: bwght

# Function for result reporting
source("_report.R")

# Load the data and estimate the model in the background
load("bwght.Rdata")
index=(!is.na(data$motheduc))&(!is.na(data$fatheduc)) # Remove data with mo/fatheduc missing
data=data[index,]
model1=lm(bwght~cigs+parity+faminc+motheduc+fatheduc, data=data)
model2=lm(bwght~cigs+parity+faminc, data=data)
dig1=c(2,3,3,3,3,3,4)
dig2=c(2,3,3,3,4)

# Describe the model and state the hyphothesis
cat("\nModel to estimate: bwght = beta0 + beta1 * cigs + beta2 * parity + beta3 * faminc + beta4 * motheduc + beta5 * fatheduc + u",
    "\nwhere bwght is ", paste(desc[desc[,1]=="bwght",2]),
    "\ncigs is ", paste(desc[desc[,1]=="cigs",2]),
    "\nparity is ", paste(desc[desc[,1]=="parity",2]),
    "\nfaminc is ", paste(desc[desc[,1]=="faminc",2]),
    "\nmotheduc is ", paste(desc[desc[,1]=="motheduc",2]),
    "\nand fatheduc is ", paste(desc[desc[,1]=="fatheduc",2]),
    "\nWe want to test the hypothesis that after controlling for cigs, parity, and faminc, parents’ education has no effect on bwght, i.e. H0: beta4 = 0, beta5 = 0",
    sep="")

# Report results
{
  cat("The estimation result for the unrestricted model is")
  reportreg(model1,dig1)
}
{
  cat("And the estimation result for the restricted model is")
  reportreg(model2,dig2)
}

# Interpretation
cat("The denominator df is n - k - 1 = ", nrow(model1$model)-nrow(summary(model1)$coef),
    "; there are 2 exclusion restrictions, so the numerator df is q = 2. The 5% critical value of this F test is 3.0. Since the F statistic is [(",
    printr(model1,dig1[7]), " - ", printr(model2,dig2[5]), ")/(1 - ", printr(model1,dig1[7]),
    ")] * (", nrow(model1$model)-nrow(summary(model1)$coef), "/2) = ",
    round((as.numeric(printr(model1,dig1[7]))-as.numeric(printr(model2,dig2[5])))/(1-as.numeric(printr(model1,dig1[7])))*(nrow(model1$model)-nrow(summary(model1)$coef))/2,2),
    " < 3.0, we fail to reject H0 at the 5% level, i.e. motheduc and fatheduc are jointly insignificant at the 5% level",
    sep="")