# How to predict by hand in R using splines regression? [closed]

The R package splines allows one to fit a non linear model using splines. For instance,

require(stats); require(graphics)
bs(womenheight, df = 5) summary(fm1 <- lm(weight ~ bs(height, df = 5), data = women)) ## example of safe prediction plot(women, xlab = "Height (in)", ylab = "Weight (lb)") ht <- seq(57, 73, length.out = 200) lines(ht, predict(fm1, data.frame(height = ht)))  This produces the following estimates Call: lm(formula = weight ~ bs(height, df = 5), data = women) Residuals: Min 1Q Median 3Q Max -0.31764 -0.13441 0.03922 0.11096 0.35086 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 114.8799 0.2167 530.146 < 2e-16 *** bs(height, df = 5)1 3.4657 0.4595 7.543 3.53e-05 *** bs(height, df = 5)2 13.0300 0.3965 32.860 1.10e-10 *** bs(height, df = 5)3 27.6161 0.4571 60.415 4.70e-13 *** bs(height, df = 5)4 40.8481 0.3866 105.669 3.09e-15 *** bs(height, df = 5)5 49.1296 0.3090 158.979 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.2276 on 9 degrees of freedom Multiple R-squared: 0.9999, Adjusted R-squared: 0.9998 F-statistic: 1.298e+04 on 5 and 9 DF, p-value: < 2.2e-16  If I want to predict using these estimates, what should I put in the predictive model? $$\hat{y} = \hat{\beta}^T \textbf{(?)}.$$ I know that I can obtain the predictions using the command predict, but I want to understand what is this command doing. Is it $$\textbf{(?)} = bs(x, df = 5)$$? • Did you try model.matrix() applied to your model? That should give you what you need. – Isabella Ghement Aug 13 '19 at 22:53 ## 1 Answer As per my comment, once you fit your model, you can extract the values of the predictors included in the model using the model.matrix() function: require(stats) require(splines) require(graphics) fm1 <- lm(weight ~ bs(height, df = 5), data = women) summary(fm1) model.matrix(fm1)  The R output produced by model.matrix() for your model is as follows: > model.matrix(fm1) (Intercept) bs(height, df = 5)1 bs(height, df = 5)2 bs(height, df = 5)3 bs(height, df = 5)4 bs(height, df = 5)5 1 1 0.000000e+00 0.000000000 0.000000000 0.000000e+00 0.000000000 2 1 4.534439e-01 0.059857872 0.001639942 0.000000e+00 0.000000000 3 1 5.969388e-01 0.203352770 0.013119534 0.000000e+00 0.000000000 4 1 5.338010e-01 0.376366618 0.044278426 0.000000e+00 0.000000000 5 1 3.673469e-01 0.524781341 0.104956268 0.000000e+00 0.000000000 6 1 2.001640e-01 0.595025510 0.204719388 9.110787e-05 0.000000000 7 1 9.110787e-02 0.566326531 0.336734694 5.830904e-03 0.000000000 8 1 3.125000e-02 0.468750000 0.468750000 3.125000e-02 0.000000000 9 1 5.830904e-03 0.336734694 0.566326531 9.110787e-02 0.000000000 10 1 9.110787e-05 0.204719388 0.595025510 2.001640e-01 0.000000000 11 1 0.000000e+00 0.104956268 0.524781341 3.673469e-01 0.002915452 12 1 0.000000e+00 0.044278426 0.376366618 5.338010e-01 0.045553936 13 1 0.000000e+00 0.013119534 0.203352770 5.969388e-01 0.186588921 14 1 0.000000e+00 0.001639942 0.059857872 4.534439e-01 0.485058309 15 1 0.000000e+00 0.000000000 0.000000000 0.000000e+00 1.000000000 attr(,"assign") [1] 0 1 1 1 1 1  According to this output, the weight can be predicted from height via a linear combination of basis functions as follows: \begin{align*} weight = \beta_0 + \beta_1*bs(height, df = 5)1 + \beta_2*bs(height, df = 5)2 + \\ \beta_3*bs(height, df = 5)3 + \beta_4* bs(height, df = 5)4 + \\ \beta_5 * bs(height, df = 5)5 \end{align*} However, since the regression coefficients $$\beta_0$$ through $$\beta_5$$ are unknown, you will need to replace them with their estimated values $$b0$$ through $$b5$$ reported in the model summary when computing the predicted weight. For example, from summary(fm1), you will find that $$b0 = 114.8799$$, $$b1 = 3.4657$$, etc. You can therefore predict weight as a function of height - for all the heights observed in the data - using the R command, which is R's translation of the above equation): pred <- model.matrix(fm1) %*% coef(fm1)  The predicted weights are as follows: > pred [,1] 1 114.8799 2 117.2767 3 119.9608 4 122.8568 5 125.8895 6 128.9841 7 132.1124 8 135.3176 9 138.6491 10 142.1564 11 145.8886 12 149.8935 13 154.2176 14 158.9074 15 164.0095  Finally, you can plot these predicted weights and compare them against your "safe" predictions to see they look identical: plot(women, xlab = "Height (in)", ylab = "Weight (lb)") ht <- seq(57, 73, length.out = 200) lines(ht, predict(fm1, data.frame(height = ht)), lwd=3, col="grey") lines(womenheight, pred, col = "magenta", lty=2, lwd = 2)

legend("topleft", c("Safe prediction","model.matrix() prediction"),
lty=c(1,2), lwd=c(3,2), col=c("grey","magenta"))


You can plot the basis functions used in the approximation of the effect of height on weight with these commands:

matplot(women\$height, model.matrix(fm1)[,-1], type='l',
col=rainbow(5), lwd=2)
legend("top", legend = dimnames(model.matrix(fm1))[[2]][-1],
col=rainbow(5), lty=1:5, bty="n", lwd=2)



The resulting plot is shown below. From this plot, if you want to predict weight for a height of 70 inches, say, you would evaluate each of the 5 basis functions at 70 inches and then multiply the results with the corresponding basis function coefficients reported by R and add the estimated intercept on top of it:

\begin{align*} weight = b_0 + b_1*bs(height = 70, df = 5)1 + b_2*bs(height = 70, df = 5)2 + \\ b_3*bs(height = 70, df = 5)3 + b_4* bs(height = 70, df = 5)4 + \\ b_5 * bs(height = 70, df = 5)5 \end{align*}

Here, the $$b$$'s are used to denote the estimated values of the $$\beta$$'s.

• Great answer. Thanks. – Splinter Aug 14 '19 at 7:41
• You're welcome, @Splinter! 😊 – Isabella Ghement Aug 14 '19 at 22:23