Regression Inference in R
Here we will discuss how to perform standard inference in regression models.
Contents
Setup
We continue the example we started in R_Regression#A first example and which is replicated here:
# This is my first R regression!
setwd("T:/ECLR/R/FirstSteps") # This sets the working directory
mydata <- read.csv("mroz.csv",na.strings = ".") # Opens mroz.csv from working directory
Before we run our initial regression model we shall restrict the dataframe mydata to those data that do not have missing wage information, using the following subset command:
mydata <- subset(mydata, wage!="NA") # select non NA data
Now we can run our initial regression:
# Run a regression
reg_ex1 <- lm(lwage~exper+log(huswage),data=mydata)
reg_ex1_sm <- summary(reg_ex1)
We will introduce inference in this model.
t-tests
We use t-tests to test simple coefficient restrictions on regression coefficients. Let's initially have a look at our regression output
print(reg_ex1_sm)
which delivers the following regression output
Call:
lm(formula = lwage ~ exper + log(huswage), data = mydata)
Residuals:
Min 1Q Median 3Q Max
-3.10089 -0.31219 0.02919 0.37466 2.11402
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.534866 0.139082 3.846 0.000139 ***
exper 0.016684 0.004243 3.933 9.81e-05 ***
log(huswage) 0.236466 0.063684 3.713 0.000232 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7031 on 425 degrees of freedom
Multiple R-squared: 0.05919, Adjusted R-squared: 0.05477
F-statistic: 13.37 on 2 and 425 DF, p-value: 2.338e-06
As you can see, this output contains t-statistics and their associated p-values. These test statistics and their p-values are all associate to the following hypothesis test: [math]H_0: \beta_{i} = 0; H_A: \beta_{i} \neq 0[/math]. Here [math]\beta_{i}[/math] represents the ith unknown population parameter. If you want to test any other hypothesis (rather than the two-sided, equal to 0 hypothesis) you will need to access the regression output in order to calculate
[math]t-stat = \frac{\widehat{\beta}_{i} - \beta_{i}}{se_{\widehat{\beta}_{i}}}[/math]
As was discussed in R_Regression#Accessing Regression Output the easiest way to get to these is to recognise that the coefficients element of reg_ex1_sm contains the parameter estimates in the first column and the standard errors in the second column. So let's say that we wanted to test the following hypothesis:
[math]H_0: \beta_{exper} = 0.01; H_A: \beta_{exper} \gt 0.01[/math]
then we can calculate the relevant test statistic according to:
t_test = (reg_ex1_sm$coefficients[2,1]-0.1)/reg_ex1_sm$coefficients[2,2]
where we recognise that that the experience coefficient is saved in the 2nd row of coefficients. As it turns out the value of this t-test is 1.5755.
F-tests
F-tests are used to test multiple coefficient restrictions on regression coefficients.
Let's say we are interested whether two additional variables age and educ should be included into the model. As a good econometrics student, or even master, you know that to calculate a F-test you need residual sum of squares from a restricted model (that is model reg_ex1) and an unrestricted model. The latter we estimate here:
reg_ex2 <- lm(lwage~exper+log(huswage)+age+educ,data=mydata) reg_ex2_sm <- summary(reg_ex2)
Calculating the F-test is now very easy. What we need is a restricted and unrestricted model. We use the function anova:
print(anova(reg_ex1,reg_ex2))
which delivers the following output:
Analysis of Variance Table
Model 1: lwage ~ exper + log(huswage)
Model 2: lwage ~ exper + log(huswage) + age + educ
Res.Df RSS Df Sum of Sq F Pr(>F)
1 425 210.11
2 423 188.10 2 22.004 24.741 6.895e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The table at the heart of this output delivers the individual residual sum of squares, the F-test statstic and its p-value. The p-value is extremely small which would lead us to reject the null hypothesis, concluding that at least one of age or educ was significant. If you look at the regression output of reg_ex2 you will see that it is the education variable.
Slighly easier, use the lht function
lht1 <- lht(reg_ex2, c("age =0"," educ = 0"))
which requires as input the unrestricted model and the restrictions.
p-values
Often you will want to calculate p-values for test statistics. You need a number of ingredients to do that:
- You need to know the value of a calculated test statistic
- You need to know what the distribution of the test statistic is (assuming that the null hypothesis is true) - this includes knowledge of degrees of freedom parameters if these are required for your distribution.
- You need to know whether you are working with a two-tailed, left-tailed or right-tailed test.
Once you know all these ingredients you can use some R internal functions to get p-values. Let's go to the t-test we calculated earlier (and saved in t_test) to test [math]H_0: \beta_{exper} = 0.01; H_A: \beta_{exper} \gt 0.01[/math]. To get the p-value here (right-tailed area) we can call on the function pt, which calculates probabilities from a t distribution (probability under the t-distribution with 425 degrees of freedom to the left of t_test:
p_t_test = (1-pt(t_test,425)) # right tailed p-value
The result is 0.05794. The (1 - ... ) part ensures that we calculate the right tail. Alternatively (as usual there are several ways to achieve the same result - all equally valied) we could have used an option in the pt function to ask for the right tail:
p_t_test = pt(t_test,425,lower.tail = FALSE)
Of course different test statistics require different distributions. But the principle is the same. The relevant function to get p-values for a F-test is pf. Type ?pf into R to get to the help function where you can see how to use it exactly.
These probability functions exist for a range of common distributions. For a list see [1].
