In order to make better models, sometimes we create new features in our data. This can be as simple as creating a column like BMI or bill_ratio like we’ve done in past classworks, or extracting the day of the week from a date string, but it can also be things like creating polynomial features, step functions, splines, or interactions.
These features can help our model perform better. For polynomial features, step functions, and splines, it also allows us to add non-linearity to our predictions even though we’re using linear regression.
We will use multiple approaches for modelling non-linearity and apply it to the Boston dataset included in the R library MASS. This dataset consists of 506 samples. The response variable is median value of owner-occupied homes in Boston (medv). The dataset has 13 associated predictor variables.
For convenience we can name the response as y and the predictor x. We will also pre-define the labels for the x and y-axes that we will use repeatedly in figures throughout this practical.
y = Boston$medvx = Boston$lstaty.lab ='Median Property Value'x.lab ='Lower Status (%)'
plot( x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, main ="", bty ='l' )
1.1 Polynomial Regression
Start by fitting to the data a degree-2 polynomial using the command lm() and summarizing the results using summary().
poly2 =lm(y ~poly(x, 2, raw =TRUE))summary(poly2)
Call:
lm(formula = y ~ poly(x, 2, raw = TRUE))
Residuals:
Min 1Q Median 3Q Max
-15.2834 -3.8313 -0.5295 2.3095 25.4148
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.862007 0.872084 49.15 <2e-16 ***
poly(x, 2, raw = TRUE)1 -2.332821 0.123803 -18.84 <2e-16 ***
poly(x, 2, raw = TRUE)2 0.043547 0.003745 11.63 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.524 on 503 degrees of freedom
Multiple R-squared: 0.6407, Adjusted R-squared: 0.6393
F-statistic: 448.5 on 2 and 503 DF, p-value: < 2.2e-16
The argument raw = TRUE In terms of fitting the curve poly(x, 2, raw = TRUE)) and poly(x, 2)) will give the same result! They are just based on different (orthogonal) basis but with polynomial regression we are almost never interested in the regression coefficients.
For plotting th results, we need to create an object, which we name sort.x, which has the sorted values of predictor x in a ascending order. Without sort.x we will not be able to produce the plots since in lecture. Then, we need to use predict() with sort.x as input in order to proceed to the next steps.
sort.x =sort(x)sort.x[1:10] # the first 10 sorted values of x
pred2 =predict(poly2, newdata =list(x = sort.x), se =TRUE)names(pred2)
[1] "fit" "se.fit" "df" "residual.scale"
The object pred2 contains fit, which are the fitted values, and se.fit, which are the standard errors of the mean prediction, that we need in order to construct the approximate 95% confidence intervals (of the mean prediction). With this information we can construct the confidence intervals using cbind(). Lets see how the first 10 fitted values and confidence intervals look like.
pred2$fit[1:10] # the first 10 fitted values of the curve
plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="Degree-2 polynomial", bty ='l')lines(sort.x, pred2$fit, lwd =2, col ="firebrick")matlines(sort.x, se.bands2, lwd =1.4, col ="firebrick", lty =3)
Note: We use lines() for pred2$fit because this is a vector, but for se.bands2, which is a matrix, we have to use matlines().
Then we do similar steps to produce a plot of degree-3 up to degree-5 polynomial fits.
poly3 =lm(y ~poly(x, 3))poly4 =lm(y ~poly(x, 4))poly5 =lm(y ~poly(x, 5))pred3 =predict(poly3, newdata =list(x = sort.x), se =TRUE)pred4 =predict(poly4, newdata =list(x = sort.x), se =TRUE)pred5 =predict(poly5, newdata =list(x = sort.x), se =TRUE)se.bands3 =cbind(pred3$fit +2*pred3$se.fit, pred3$fit-2*pred3$se.fit)se.bands4 =cbind(pred4$fit +2*pred4$se.fit, pred4$fit-2*pred4$se.fit)se.bands5 =cbind(pred5$fit +2*pred5$se.fit, pred5$fit-2*pred5$se.fit)par(mfrow =c(2,2))# Degree-2plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="Degree-2 polynomial", bty ='l')lines(sort.x, pred2$fit, lwd =2, col ="firebrick")matlines(sort.x, se.bands2, lwd =2, col ="firebrick", lty =3)# Degree-3plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="Degree-3 polynomial", bty ='l')lines(sort.x, pred3$fit, lwd =2, col ="darkviolet")matlines(sort.x, se.bands3, lwd =2, col ="darkviolet", lty =3)# Degree-4plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="Degree-4 polynomial", bty ='l')lines(sort.x, pred4$fit, lwd =2, col ="royalblue")matlines(sort.x, se.bands4, lwd =2, col ="royalblue", lty =3)# Degree-5plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="Degree-5 polynomial", bty ='l')lines(sort.x, pred5$fit, lwd =2, col ="darkgreen")matlines(sort.x, se.bands5, lwd =2, col ="darkgreen", lty =3)
All four curves look reasonable given the data available. We may choose the degree-2 polynomial since it is simpler and seems to do about as well as the others. However, if we want to base our decision on a more formal procedure, we can use analysis-of-variance (ANOVA). Specifically, we will perform sequential comparisons based on the F-test, comparing first the linear model vs. the quadratic model (degree-2 polynomial), then the quadratic model vs. the cubic model (degree-3 polynomial) and so on. We therefore have to fit the simple linear model, and we also choose to fit the degree-6 polynomial to investigate the effects of an additional predictor as well. We can perform this analysis in RStudio using the command anova() as displayed below.
Analysis of Variance Table
Model 1: y ~ x
Model 2: y ~ poly(x, 2, raw = TRUE)
Model 3: y ~ poly(x, 3)
Model 4: y ~ poly(x, 4)
Model 5: y ~ poly(x, 5)
Model 6: y ~ poly(x, 6)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 504 19472
2 503 15347 1 4125.1 151.8623 < 2.2e-16 ***
3 502 14616 1 731.8 26.9390 3.061e-07 ***
4 501 13968 1 647.8 23.8477 1.406e-06 ***
5 500 13597 1 370.7 13.6453 0.0002452 ***
6 499 13555 1 42.4 1.5596 0.2123125
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Which model would you choose?
1.2 Step Functions
For step function regression we can make use of the command cut(), which automatically assigns samples to intervals given a specific number of intervals. We can check how this works by executing the following syntax:
table(cut(x, 2))
(1.69,19.9] (19.9,38]
430 76
What we see is that cut(x, 2) automatically created a factor with two levels, corresponding to the intervals \((1.69,19.9]\) and \((19.9,38]\) and assigned each entry in x to one of these factors depending on which interval it was in. The command table() tells us that 430 samples of x fall within the first interval and that 76 samples fall within the second interval. Note that cut(x, 2) generated 2 intervals, but this means there is only 1 cutpoint (at 19.9). The number of cutpoints is naturally one less than the number of intervals, but it is important to be aware that cut requires specification of the number of required intervals.
So, we can use cut() within lm() to easily fit regression models with step functions. Below we consider 4 models with 1, 2, 3 and 4 cutpoints (2, 3, 4 and 5 intervals) respectively.
The analysis then is essentially the same as previously. We plot the fitted lines of the four models, along with approximate 95% confidence intervals for the mean predictions.
#| fig-width: 9#| fig-height: 8pred2 =predict(step2, newdata =list(x =sort(x)), se =TRUE)pred3 =predict(step3, newdata =list(x =sort(x)), se =TRUE)pred4 =predict(step4, newdata =list(x =sort(x)), se =TRUE)pred5 =predict(step5, newdata =list(x =sort(x)), se =TRUE)se.bands2 =cbind(pred2$fit +2*pred2$se.fit, pred2$fit-2*pred2$se.fit)se.bands3 =cbind(pred3$fit +2*pred3$se.fit, pred3$fit-2*pred3$se.fit)se.bands4 =cbind(pred4$fit +2*pred4$se.fit, pred4$fit-2*pred4$se.fit)se.bands5 =cbind(pred5$fit +2*pred5$se.fit, pred5$fit-2*pred5$se.fit)par(mfrow =c(2,2))plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="1 cutpoint", bty ='l')lines(sort(x), pred2$fit, lwd =2, col ="firebrick")matlines(sort(x), se.bands2, lwd =1.4, col ="firebrick", lty =3)plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="2 cutpoints", bty ='l')lines(sort(x), pred3$fit, lwd =2, col ="darkviolet")matlines(sort(x), se.bands3, lwd =1.4, col ="darkviolet", lty =3)plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="3 cutpoints", bty ='l')lines(sort(x), pred4$fit, lwd =2, col ="royalblue")matlines(sort(x), se.bands4, lwd =1.4, col ="royalblue", lty =3)plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab,main ="4 cutpoints", bty ='l')lines(sort(x), pred5$fit, lwd =2, col ="darkgreen")matlines(sort(x), se.bands5, lwd =1.4, col ="darkgreen", lty =3)
Note that we do not necessarily need to rely on the automatic selections of cutpoints used by cut(). We can define the intervals if we want to. For instance, if we want cutpoints at 10, 20 and 30 we can do the following
For this analysis we will require package splines.
library(splines)
Initially let’s fit regression splines by specifying knots. From the previous plot it is not clear where exactly we should place knots, so we will make use of the command summary in order to find the 25th, 50th and 75th percentiles ofx, which will be the positions where we will place the knots. We also sort the variable x before fitting the splines.
summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.73 6.95 11.36 12.65 16.95 37.97
cuts =summary(x)[c(2, 3, 5)] cuts
1st Qu. Median 3rd Qu.
6.950 11.360 16.955
sort.x =sort(x)
For a start lets fit a linear spline using our selected placement of knots. For this we can use command lm() and inside it we use the command bs() in which we specify degree = 1 for a linear spline and knots = cuts for the placement of the knots at the three percentiles. We also calculate the corresponding fitted values and confidence intervals exactly in the same way we did in previous practical demonstrations.
plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, main ="Linear Spline", bty ='l')lines(sort.x, pred1$fit, lwd =2, col ="firebrick")matlines(sort.x, se.bands1, lwd =2, col ="firebrick", lty =3)
Using ?bs we see that instead of using the argument knots we can use the argument df, which are the degrees of freedom. Splines have \((d+1)+K\) degrees of freedom, where \(d\) is the degree of the polynomial and \(K\) the number of knots. So in this case we have 1+1+3 = 5 degrees of freedom. Selecting df = 5 in bs() will automatically use 3 knots placed at the 25th, 50th and 75th percentiles. Below we check whether the plot based on df=5 is indeed the same as the previous plot and as we can see it is.
spline1df =lm(y ~ splines::bs(x, degree =1, df =5))pred1df =predict(spline1df, newdata =list(x = sort.x), se =TRUE)se.bands1df =cbind( pred1df$fit +2* pred1df$se.fit, pred1df$fit -2* pred1df$se.fit )par(mfrow =c(1, 2))plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, main ="Linear Spline (with knots)", bty ='l')lines(sort.x, pred1$fit, lwd =2, col ="firebrick")matlines(sort.x, se.bands1, lwd =2, col ="firebrick", lty =3)plot(x, y, cex.lab =1.1, col="darkgrey", xlab = x.lab, ylab = y.lab, main ="Linear Spline (with df)", bty ='l')lines(sort.x, pred1df$fit, lwd =2, col ="firebrick")matlines(sort.x, se.bands1df, lwd =2, col ="firebrick", lty =3)matlines(sort.x, se.bands1, lwd =2, col ="firebrick", lty =3)
