# 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="")