R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(wooldridge)
library(sandwich)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
data("meap00_01")
head(meap00_01) 
##   dcode bcode math4 read4 lunch enroll    exppp  lenroll   lexppp
## 1  1010  4937  83.3  77.8 40.60    468 5870.673 6.148468 8.677725
## 2  3010   790  85.7  60.0 12.75    251 4824.836 5.525453 8.481532
## 3  3010  1403  77.3  59.1 17.08    439 4358.772 6.084499 8.379946
## 4  3010  4056  85.2  67.0 23.17    561 4701.396 6.329721 8.455615
## 5  3020   922  91.8  67.7 25.57    442 4311.984 6.091310 8.369153
## 6  3020  2864  87.7  80.0 27.30    381 4507.331 5.942800 8.413461
colnames(meap00_01) 
## [1] "dcode"   "bcode"   "math4"   "read4"   "lunch"   "enroll"  "exppp"  
## [8] "lenroll" "lexppp"
# Estimate OLS Model (Part i)
model_ols <- lm(math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01)
summary(model_ols) # View OLS results
## 
## Call:
## lm(formula = math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.326  -8.758   0.914   9.164  51.416 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 91.93240   19.96170   4.605 4.42e-06 ***
## lunch       -0.44874    0.01464 -30.648  < 2e-16 ***
## log(enroll) -5.39915    0.94041  -5.741 1.11e-08 ***
## log(exppp)   3.52475    2.09785   1.680   0.0931 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.3 on 1688 degrees of freedom
## Multiple R-squared:  0.3729, Adjusted R-squared:  0.3718 
## F-statistic: 334.6 on 3 and 1688 DF,  p-value: < 2.2e-16
# Calculate Robust Standard Errors
robust_se <- coeftest(model_ols, vcov = vcovHC(model_ols, type = "HC1"))
print(robust_se) # Compare robust and usual standard errors
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 91.932404  23.087103   3.9820 7.125e-05 ***
## lunch       -0.448743   0.016584 -27.0583 < 2.2e-16 ***
## log(enroll) -5.399152   1.131175  -4.7730 1.971e-06 ***
## log(exppp)   3.524751   2.353737   1.4975    0.1344    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Perform White Test for Heteroskedasticity (Part ii)
bptest(model_ols, ~ fitted(model_ols) + I(fitted(model_ols)^2), data = meap00_01)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 229.78, df = 2, p-value < 2.2e-16
# Weighted Least Squares (WLS) Estimation (Part iii)
# Regress log(residuals^2) on math4, math4^2
meap00_01$residuals_sq <- residuals(model_ols)^2
meap00_01$log_residuals_sq <- log(meap00_01$residuals_sq)

aux_model <- lm(log_residuals_sq ~ math4 + I(math4^2), data = meap00_01)
summary(aux_model) # View auxiliary regression results
## 
## Call:
## lm(formula = log_residuals_sq ~ math4 + I(math4^2), data = meap00_01)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.1803  -0.8389   0.3460   1.4690   3.7304 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.7704029  0.4412731   26.67   <2e-16 ***
## math4       -0.2361655  0.0142496  -16.57   <2e-16 ***
## I(math4^2)   0.0016576  0.0001105   15.00   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.135 on 1689 degrees of freedom
## Multiple R-squared:  0.1594, Adjusted R-squared:  0.1584 
## F-statistic: 160.1 on 2 and 1689 DF,  p-value: < 2.2e-16
# Calculate Fitted Values and Weights
meap00_01$weights <- exp(fitted(aux_model))

# Perform WLS
model_wls <- lm(math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01, weights = 1 / meap00_01$weights)
summary(model_wls) # View WLS results
## 
## Call:
## lm(formula = math4 ~ lunch + log(enroll) + log(exppp), data = meap00_01, 
##     weights = 1/meap00_01$weights)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1285 -1.4897  0.1499  1.3713  4.6178 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.34591   15.05791   3.211 0.001349 ** 
## lunch       -0.25518    0.01182 -21.596  < 2e-16 ***
## log(enroll) -1.65287    0.71742  -2.304 0.021348 *  
## log(exppp)   5.36162    1.56806   3.419 0.000643 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.706 on 1688 degrees of freedom
## Multiple R-squared:  0.2179, Adjusted R-squared:  0.2165 
## F-statistic: 156.8 on 3 and 1688 DF,  p-value: < 2.2e-16
# Compare WLS and OLS Coefficients (Part iii continued)
ols_coef <- summary(model_ols)$coefficients
wls_coef <- summary(model_wls)$coefficients

print("OLS Coefficients:")
## [1] "OLS Coefficients:"
print(ols_coef)
##               Estimate  Std. Error    t value      Pr(>|t|)
## (Intercept) 91.9324043 19.96169801   4.605440  4.424155e-06
## lunch       -0.4487434  0.01464203 -30.647632 2.240954e-164
## log(enroll) -5.3991517  0.94041163  -5.741264  1.112079e-08
## log(exppp)   3.5247506  2.09784606   1.680176  9.310811e-02
print("WLS Coefficients:")
## [1] "WLS Coefficients:"
print(wls_coef)
##               Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept) 48.3459141 15.05790988   3.210666 1.349192e-03
## lunch       -0.2551755  0.01181573 -21.596250 1.566017e-91
## log(enroll) -1.6528723  0.71741850  -2.303917 2.134828e-02
## log(exppp)   5.3616196  1.56806479   3.419259 6.428607e-04
# Compute Robust Standard Errors for WLS (Part iv)
robust_se_wls <- coeftest(model_wls, vcov = vcovHC(model_wls, type = "HC1"))
print("WLS Robust Standard Errors:")
## [1] "WLS Robust Standard Errors:"
print(robust_se_wls)
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 48.345914  15.198596   3.1809 0.0014946 ** 
## lunch       -0.255175   0.011958 -21.3396 < 2.2e-16 ***
## log(enroll) -1.652872   0.702104  -2.3542 0.0186781 *  
## log(exppp)   5.361620   1.588233   3.3758 0.0007526 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Compare Precision (OLS vs. WLS for log(exppp)) (Part v)
print("Effect of Spending (log(exppp)):")
## [1] "Effect of Spending (log(exppp)):"
print("OLS Estimate:")
## [1] "OLS Estimate:"
print(ols_coef["log(exppp)", ])
##   Estimate Std. Error    t value   Pr(>|t|) 
## 3.52475064 2.09784606 1.68017602 0.09310811
print("WLS Estimate:")
## [1] "WLS Estimate:"
print(wls_coef["log(exppp)", ])
##     Estimate   Std. Error      t value     Pr(>|t|) 
## 5.3616195672 1.5680647921 3.4192589452 0.0006428607

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.