Replication of Examples in Chapter 5

We compute the t-statistic based on the following formula,

Upon computing the t-statistic, we compare it with the critical value at the desired significant level, say 5%, which is 1.96 for a two-sided test with the standard normal distribution. Also, we can use the t-statistic to get the p-value.

How can we get all the quantities used in this formula? Of course, you can simply copy them from the output of summary(mod1). But doing so is cumbersome and prone to mistakes. More importantly, we should use the heteroskedasticity-robust standard error of \(\hat{\beta}_1\) instead of the ones reported by default.

Next, I will show you two ways to get the t-statistics using the heteroskedasticity-robust standard errors. The first way is to get all the quantities in Equation (\ref{eq:t-stat-b1}) using R functions, and the second way is to get the appropriate t-statistic through the function coeftest(). Although the second method is much easier than the first one, we can learn how to obtain all the elements in the output of the function lm().

Since hypothesis testing is a very common work in statistics, many R functions have been developed to do it. Here I introduce a function, coeftest(), which is in the package of lmtest loaded automatically through loading the AER package.

coeftest(mod1)

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)
(Intercept) 698.93295    9.46749 73.8245 < 2.2e-16 ***
str          -2.27981    0.47983 -4.7513 2.783e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

By default, it reports the homoskedasticity-only standard errors, the corresponding t-statistics, and the p-values. To get the heteroskedasticity-robust results, we need to add an argument to this function to specify the heteroskedasticity-consistent covariance matrix, which has been defined above as htV <- vcovHC(mod1, type = "HC1").

# coeftest returns a matrix
t_tst <- coeftest(mod1, vcov. = htV); t_tst

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)
(Intercept) 698.93295   10.36436 67.4362 < 2.2e-16 ***
str          -2.27981    0.51949 -4.3886 1.447e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## Get the t-statistic for STR
t_b1 <- t_tst["str", "t value"]; t_b1

[1] -4.388557

Finally, we can construct the 95% confidence interval of \(\beta_1\) using the function of confint(), which uses the homoskedasticity-only standard errors.

# confidence interval with the default homoskedasticity-only SE
confint(mod1, "str")

       2.5 %    97.5 %
str -3.22298 -1.336637

Since there is no existing function to report the confidence interval with heteroskedasticity-robust SE, we can write a user-defined function to do that.

get_confint_rb <- function(lm_obj, param, vcov_ = vcov(lm_obj),
                           level = 0.05){
    ## This function generates a two-sided confidence interval for a
    ## parameter in the linear regression model with a specified
    ## covariance matrix.  The inputs The output

    ## get all the parameters' names and select one based on param
    all_param <- names(coef(lm_obj))
    which_param <- grep(param, all_param)

    ## get the estimated parameter and its standard error
    bhat_param <- coef(lm_obj)[which_param]
    sd_param <- sqrt(vcov_[which_param, which_param])

    ## get the critical value
    cv <- qnorm(1 - level/2)

    ## calculate the confidence interval
    lower <- bhat_param - cv * sd_param
    upper <- bhat_param + cv * sd_param

    conf_interval <- c(lower, upper)
    names(conf_interval) <- c("lower", "upper")
    return(conf_interval)
}

Note that we define the default value of vcov_ to be the homoskedasticity-only covariance matrix in the function get_confint_rb(), and the default significant level is 5%. When computing the confidence interval with the heteroskedasticity-robust standard errors, we need to change vcov_ to a heteroskedasticity-consistent covariance matrix, which is htV.

get_confint_rb(mod1, param = "str", vcov = htV)

    lower     upper
-3.297988 -1.261628

We define a dummy variable representing small classes.

smallclass <- ifelse(df2use$str < 20, 1, 0)

The function ifelse() creates a vector consisting of 1 and 0. The first argument in this function is a condition, df2use$str < 20. If the condition is satisfied for an element in df2use$str, the corresponding element in D is 1, otherwise 0.

Then we can estimate the linear regression of test scores against the dummy variable, and do the zero hypothesis test.

mod2 <- lm(testscr ~ smallclass, data = df2use)
coeftest(mod2, vcov. = vcovHC(mod2, type = "HC1"))

t test of coefficients:

            Estimate Std. Error  t value  Pr(>|t|)
(Intercept) 649.9788     1.3229 491.3317 < 2.2e-16 ***
smallclass    7.3724     1.8236   4.0428 6.288e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We can create a scatterplot to visualize the relationship between test scores and the dummy variable.

plot(smallclass, df2use$testscr,
     xlab = "class size", ylab = "test score")

In Figure 1, since the variable smallclass takes only the value of 1 or 0, many points are overlapped. To make these points more visible, we can create a scatterplot with jittered points. I also adjust the limit of x-axis so that the clusters of points appear close towards the center of the plot.

smallclass_jittered <- jitter(smallclass, amount = 0.05)
plot(smallclass_jittered, df2use$testscr, xlim = c(-0.5, 1.5),
     xlab = "class size", ylab = "test score")
abline(coef(mod2)[1], coef(mod2)[2], col = "blue")