2
$\begingroup$

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

require(stats); require(graphics)
bs(women$height, 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)$?

$\endgroup$

closed as off-topic by Michael Chernick, Peter Flom Aug 14 at 10:40

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Michael Chernick, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Did you try model.matrix() applied to your model? That should give you what you need. $\endgroup$ – Isabella Ghement Aug 13 at 22:53
4
$\begingroup$

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(women$height, 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"))

enter image description here

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.

enter image description here

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.