A simple linear regression analysis

Fitting a linear regression model in R

If we assume that the trend in attainment is linear, we can represent the relationship between attainment and cohort by a linear regression model of the form:

SCOREi = β0 + β1COHORT90i + ei

where scorei is the attainment score for student i and cohort90i indicates their school year (centred at 1990). The difference between a student’s actual score and that predicted for their cohort is the residual ei.

β0 and β1 are the intercept and slope of the population regression line. Because cohort90 is centred around around 1990, β0 is the attainment score expected for a student taking the exams in 1990. The slope β1 is the expected change in attainment for each 1-year increase in cohort90.

In R, linear regression models are fitted with the lm command. The response variable follows the command which is then followed by a ~ then by a list of explanatory variables (excluding the constant as this is included by default).Below we choose to store the model results in a new object called fit.

 fit <- lm(score ~ cohort90, data = mydata)

We then display the results using the summary command.

 summary(fit)
## 
## Call:
## lm(formula = score ~ cohort90, data = mydata)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -41.31 -11.73   0.56  12.20  46.20 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 30.72873    0.08582  358.04   <2e-16 ***
## cohort90     1.32214    0.01601   82.59   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.8 on 33986 degrees of freedom
## Multiple R-squared:  0.1672, Adjusted R-squared:  0.1671 
## F-statistic:  6822 on 1 and 33986 DF,  p-value: < 2.2e-16

A variety of output is displayed. The top portion of output provides summary statistics for the residuals. The bottom portion of output reports standard model summary statistics. However, of main interest is the table of regression coefficients which is reported in the middle half of the output. The estimated slope coefficient is given next to cohort90 while the estimated intercept is given next to (Intercept). he column labelled Std. Error gives the associated standard errors.

Interpretation

The fitted regression line is:

SCOˆREi = 30.73+1.32COHORT90i

The estimated intercept is 30.73, which is the attainment score predicted for a student with cohort90 = 0, i.e. a student from the 1990 cohort.

The estimated slope is 1.32, which is interpreted as the increase in attainment predicted for each successive cohort. Note that we are assuming a linear trend from 1984 to 1998, i.e. the 1.32 year-on-year increase is assumed to apply across the range.

We can make predictions about attainment for each cohort using the fitted regression line. Note, however, that we should be cautious about making predictions based on a linear model for years where there are no observations, particularly between the 1990 and 1996 cohorts. We have no evidence for a linear trend during this period.

Some examples of predictions are given in below Table, alongside the actual mean attainment scores for these cohorts. The predicted values are close to the observed values, suggesting that the linear assumption is reasonable.

cohort90 scôre Observed mean score -6 (i.e. 1984) 30.73 + (1.32 × -6) = 22.8 23.7 -2 (i.e. 1988) 30.73 + (1.32 × -2) = 28.1 28.5 8 (i.e. 1998) 30.73 + (1.32 × 8) = 41.30 41.3

We can plot the fitted regression line by calculating predicted values and plotting them against cohort90. The predict command is used to calculate the predicted values and to store them in a new object called predscore:

 predscore <- predict(fit)

We use the plot command to plot the predicted values predscore against cohort90. We specify “l” for the type argument to specify that the type of graph should be a line plot:

 plot(mydata$cohort90, predscore, type = "l")

The new object predscore contains the predicted score for every observation in the data. Clearly this variable contains far more information than we need to plot out a straight line, where only two end-points are required. Plotting the predicted score for every observation is computationally inefficient as it slows down the time required to plot the graph. While this is not a particular problem here, it can be a considerable problem for more complex graphs.

A more computationally efficient way to plot this graph is to simply run the command for the subset of data which gives each distinct combinations of the two variables used to plot the graph. To do this, we first need to create a new matrix object which contains the two variables of interest. Here we combine the cohort90 variable with the predscore vector.

 uniquedata <- cbind(cohort90 = mydata$cohort90, predscore = predscore)

We then reduce this matrix to keep only unique records using the unique command

 uniquedata <- unique(uniquedata)

We now repeat the plot command, but this time using the unique records of the two variables which are found in the uniquedata object. Note how the [,1] and [,2] syntax is used to refer to column one (cohort90) and column two (predscore) of the uniquedata matrix.

 plot(uniquedata[, 1], uniquedata[, 2], xlab = "cohort", ylab = "predscore", type = "l")

The following graph will appear which is the same as the computationally less efficient one plotted above. However, this graph is drawn considerably quicker as it is based on 6 observations rather than 33988. Note that in the plot command we specified the xlab and ylab arguments in order to give labels to the x-axis and the y-axis.

The residual (unexplained) variance parameter σ2e is not directly estimated by lm, rather it is calculated post-estimation and is reported, as a standard deviation, in the lower portion of the lm output next to Residual standard error:. The variance can be calculated straightforwardly by squaring Residual standard error: to give σˆ2e = 15.82 = 249.64.

R-squared, the proportion of the total variance in score that is explained by cohort90, is also reported in the lower portion of the lm output as Multiple R-squared: 0.1672. Thus 16.7% of the variance in attainment scores is explained by differences between cohorts.

Note that R2 is equal to the squared correlation coefficient: (0.409)^2 = 0.167.

Hypothesis testing

We can test the statistical significance of the linear trend using a t-test. The null hypothesis for the test is

H0: β1 = 0, i.e. the slope of the population regression line is zero.

The test statistic is calculated by dividing the estimated slope by its standard error:

t = 1.322/0.016 = 82.6

which is clearly strongly significant (much larger than 1.96, the critical value at the 5% level of significance). We therefore conclude that there is strong evidence of a linear trend.

Fortunately, lm automatically calculates t-test statistics for us and so we do not need to calculate these by hand. t-test statistics for each explanatory variable are reported in the table of regression coefficients in the column labelled t value. We also do not need to compare manually these t-test statistics to tables of critical values as the output additionally reports p-values for the significance of the coefficient of each explanatory variable. These are reported in the column labelled Pr(>|t|). Variables with p-values smaller than 0.05 are significant at the 5% level.

Model checking

Before carrying out further analysis, we should check the model assumptions. In particular, given the non-normality of score (due to a high proportion of zero scores), we will check the distribution of the residuals.

We will look at two plots:

-standardised residuals vs. normal scores (a normal plot), and -standardised residuals vs. fixed part prediction (scôrei).

To calculate standardised residuals we use the rstandard command and we store these in a new object named estd:

 estd <- rstandard(fit)

To produce the standardised residuals vs. normal scores plot we use the qqnorm command. We follow this command with the qqline command in order to superimpose a 45 degree line. The qqnorm command produces a normal QQ plot of the values in the variable (here estd). qqline adds a line to a normal quantile-quantile plot which passes through the first and third quartiles.

 qqnorm(estd)

 qqline(estd)

To produce the standardised residuals vs. fixed part prediction (scôrei) plot:

plot(predscore, estd)

If the residuals are normally distributed, then the first plot should show them lying on the 45 degree line. There seems to be some departure from linearity at the top and bottom ends; we will explore this further using a histogram.

The second plot is used to check for heteroskedasticity. If the assumption of constant residual variance (homoskedasticity) holds, the vertical spread of points should be the same for each predicted value. This assumption seems reasonable here.

To obtain a histogram of the standardised residuals:

hist(estd)

The histogram shows an approximately normal distribution and so we do not consider a transformation of score.

None of the above plots show any obvious outliers.

References

Crawley, MJ (2007) The R Book. John Wiley & Sons, Ltd, Chichester.

Croxford, L. and Raffe, D. (2006) “Education Markets and Social Class Inequality: A Comparison of Trends in England, Scotland and Wales”. In R. Teese (Ed.) Inequality Revisited. Berlin: Springer.