# Example 4.10, Salary-Pension Tradeoff for Teachers
# Data set: meap93

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

# Load the data and estimate the model in the background
load("meap93.Rdata")
data$bs=data$benefits/data$salary    # Create the b/s ratio
model1=lm(lsalary~bs, data=data)
model2=lm(lsalary~bs+lenroll+lstaff, data=data)
model3=lm(lsalary~bs+lenroll+lstaff+droprate+gradrate, data=data)
dig1=c(3,3)
dig2=c(3,3,4,3)
dig3=c(3,3,4,3,5,5)

# Prepare for result reporting

# Create row names
temp0=c(names(model3$coef)[-1],"intercept")
temp1=vector()
for (n in c(1:6)) {
  temp1[2*n-1]=temp0[n]
  temp1[2*n]=""
}
temp1[13:14]=c("Observations","R-squared")
# Append columns of results to the table
for (n in c(1:3)) {
  temp2=vector()
  # Fill in coefficients and SEs
  nvar=nrow(summary(eval(parse(text=paste("model",n,sep=""))))$coef)-1
  for (m in c(1:nvar)) {
    if (as.numeric(printcoef(eval(parse(text=paste("model",n,sep=""))),m+1,eval(parse(text=paste("dig",n,sep="")))[m+1]))>0)
      temp2[2*m-1]=paste("  ",printcoef(eval(parse(text=paste("model",n,sep=""))),m+1,eval(parse(text=paste("dig",n,sep="")))[m+1])," ",sep="")
    else
      temp2[2*m-1]=paste(" ",printcoef(eval(parse(text=paste("model",n,sep=""))),m+1,eval(parse(text=paste("dig",n,sep="")))[m+1])," ",sep="")
    temp2[2*m]=paste(" (",printse(eval(parse(text=paste("model",n,sep=""))),m+1,eval(parse(text=paste("dig",n,sep="")))[m+1]),")",sep="")
}
  nemp=10-2*nvar
  # Fill the empty cells
  if (nemp>0) for (p in c(1:nemp)) temp2[2*nvar+p]="   -"
  # Fill the intercept coefficients/SEs, no. of obs. and R^2
  temp2[11]=paste(" ",printcoef(eval(parse(text=paste("model",n,sep=""))),1,eval(parse(text=paste("dig",n,sep="")))[1])," ",sep="")
  temp2[12]=paste(" (",printse(eval(parse(text=paste("model",n,sep=""))),1,eval(parse(text=paste("dig",n,sep="")))[1]),")",sep="")
  temp2[13]=nrow(eval(parse(text=paste("model",n,sep="")))$model)
  temp2[14]=printr(eval(parse(text=paste("model",n,sep=""))),3)
  temp1=data.frame(temp1,temp2)
}
table=temp1
names(table)=c("Independent Variables","(1)","(2)","(3)")

# Describe the model and summarize bs
cat("In this example, we want to test the salary-benefits tradeoff. To do this, we estimate three models, in all of which lsalary is the dependent variable and the benefits to salary ratio, b/s, is an independent variable, but different controlling variables are used. (Please refer to the textbook for the derivation of the relationship.)",
    "\nModel 1: lsalary = beta0 + beta1 * b/s + u",
    "\nModel 2: lsalary = beta0 + beta1 * b/s + beta2 * lenroll + beta3 * lstaff + u",
    "\nModel 2: lsalary = beta0 + beta1 * b/s + beta2 * lenroll + beta3 * lstaff + beta4 * droprate + beta5 * gradrate + u",
    "\nwhere lenroll is ", paste(desc[desc[,1]=="lenroll",2]), " (enroll: ", paste(desc[desc[,1]=="enroll",2]), ")",
    "\nlstaff is ", paste(desc[desc[,1]=="lstaff",2]), " (staff: ", paste(desc[desc[,1]=="staff",2]), ")",
    "\ndroprate is ", paste(desc[desc[,1]=="droprate",2]),
    "\nand gradrate is ", paste(desc[desc[,1]=="gradrate",2]),
    sep="")
summary(data$bs)

# Report results
{
  cat("The estimation results are summarized in the table below:\n\n",
      "\tDependent Variable: lsalary\n")
  print(table,row.names=F,right=F)
}

# Interpretation
cat("To test the benefits-salary tradeoff, we test H0: beta1 = -1 against H1: beta1 ≠ -1",
    "\nColumn 1:\nWhen controlling for no other variables, the coefficient estimate on bs is ",
    printcoef(model1,2,dig1[2]), ", and the t statistic for testing H0 is (",
    printcoef(model1,2,dig1[2]), " + 1)/", printse(model1,2,dig1[2]),
    " = ", round((as.numeric(printcoef(model1,2,dig1[2]))+1)/as.numeric(printse(model1,2,dig1[2])),3),
    ", so we fail to reject H0",
    "\nColumn 2:\nWhen controlling for lenroll and lstaff, the coefficient estimate on bs is ",
    printcoef(model2,2,dig2[2]), ", and the t statistic is ",
    round((as.numeric(printcoef(model2,2,dig2[2]))+1)/as.numeric(printse(model2,2,dig2[2])),2),
    ", and we reject H0 in favor of H1 at the 5% level",
    "\nColumn 3:\nWhen also controlling for droprate and gradrate, the coefficient estimate on bs is ",
    printcoef(model3,2,dig3[2]), ", and the t statistic is ",
    round((as.numeric(printcoef(model3,2,dig3[2]))+1)/as.numeric(printse(model3,2,dig3[2])),2),
    ", and we again reject H0 in favor of H1 at the 5% level. Besides, the new controlling variables, droprate and gradrate, are very statistically insignificant",
    sep="")