Regression Inference in R
Here we will discuss how to perform standard inference in regression models.
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 linear hypothesis testing function (lht
) which is also part of the car package:
lht1 <- lht(reg_ex2, c("age = 0"," educ = 0"))
which requires as input the unrestricted model, (here reg_ex2
) and the restrictions. Here we have two restrictions, i.e. that the coefficients to the two variables age
and educ
are 0 and hence both variables irrelevant. You could add more restrictions by adding them into the vector list c( )
.
Look at the output and you will find that it is essentially exactly the same as for the anova
function. The lht
function is more convenient to use for testing purposes, as you will only need the unrestricted model and then add the restriction. Internally the function will then estimate the restricted model, you will not see it.
p-value
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].