0
$\begingroup$

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?

$\endgroup$
5
  • 2
    $\begingroup$ For a start, emulate the examples on the help page. Type ?bs. $\endgroup$
    – whuber
    Aug 5 '19 at 20:51
  • $\begingroup$ @whuber thanks. $\endgroup$
    – Bull
    Aug 5 '19 at 20:52
  • 2
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Aug 5 '19 at 20:52
  • 1
    $\begingroup$ @jbowman Did you mean lm(y~bsx)? If so, this works! I could accept that as an answer. $\endgroup$
    – Bull
    Aug 5 '19 at 20:55
  • $\begingroup$ I actually meant x; the point was to show the collinearity between the variables in des. $\endgroup$
    – jbowman
    Aug 5 '19 at 22:13
2
$\begingroup$

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!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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