# How to simulate non-linear regression model with splines?

I want to simulate a non linear regression model of the type

$$y_i = 1 + x_i + bs(x_i)+ e_i,$$

$$i=1,...,1000$$, where $$bs(x_i)$$ is a b-spline. I have tried to use the splines R package but when I try to fit the non linear regression model, the estimate of one of the entries of the spline is NA, and the rest of the estimates are far from the true value:

library(splines)

n = 1000
set.seed(1234)
x <- rnorm(n)

bsx <- bs(x)
eps <- rnorm(n,0,0.25)

y = 1 + x + bsx[,1] + bsx[,2] + bsx[,3] + eps

des <- cbind(x,bsx)

lm(y~ des)

Call:
lm(formula = y ~ des)

Coefficients:
(Intercept)         desx         des1         des2         des3
1.6017       1.1633       0.5384       0.2188           NA


What am I doing wrong?

• For a start, emulate the examples on the help page. Type ?bs.
– whuber
Aug 5 '19 at 20:51
• @whuber thanks.
– Bull
Aug 5 '19 at 20:52
• Try regressing lm(x~bsx) and see what happens... x is already captured as a linear combination of the columns of bsx, so you don't need to have it in the regression. Aug 5 '19 at 20:52
• @jbowman Did you mean lm(y~bsx)? If so, this works! I could accept that as an answer.
– Bull
Aug 5 '19 at 20:55
• I actually meant x; the point was to show the collinearity between the variables in des. Aug 5 '19 at 22:13

The problem is caused by the fact that x is implicitly captured in the variables returned by the bs function, so adding it as a separate regressor results in perfect multicollinearity.

To see this, we can calculate the eigenvalues of the covariance matrix of des"

> eigen(cov(des))$values [1] 1.024286e+00 7.888118e-03 1.515617e-03 1.098538e-18  where the smallest eigenvalue is effectively zero. Leaving out x results in: > eigen(cov(bsx))$values
[1] 0.0320042724 0.0060866427 0.0009164426


so the issue is evidently the addition of the column x to the columns of bsx.

Running the regression without x results in:

> summary(lm(y~bsx))

Call:
lm(formula = y ~ bsx)

Residuals:
Min       1Q   Median       3Q      Max
-0.78588 -0.15857  0.00394  0.16535  0.71671

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.34905    0.09557  -24.58   <2e-16 ***
bsx1         3.09459    0.20696   14.95   <2e-16 ***
bsx2         5.33126    0.09630   55.36   <2e-16 ***
bsx3         7.66870    0.15867   48.33   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.245 on 996 degrees of freedom
Multiple R-squared:  0.9554,    Adjusted R-squared:  0.9552
F-statistic:  7109 on 3 and 996 DF,  p-value: < 2.2e-16


... a much more useful result!