How to Interpret Lm() Coefficient Estimates When Using Bs() Function for Splines

How to interpret lm() coefficient estimates when using bs() function for splines

I would expect a -1 coefficient for the first part and a +1 coefficient for the second part.

I think your question is really about what is a B-spline function. If you want to understand the meaning of coefficients, you need to know what basis functions are for your spline. See the following:

library(splines)
x <- seq(-5, 5, length = 100)
b <- bs(x, degree = 1, knots = 0)  ## returns a basis matrix
str(b)  ## check structure
b1 <- b[, 1]  ## basis 1
b2 <- b[, 2]  ## basis 2
par(mfrow = c(1, 2))
plot(x, b1, type = "l", main = "basis 1: b1")
plot(x, b2, type = "l", main = "basis 2: b2")

basis

Note:

B-splines of degree-1 are tent functions, as you can see from b1;
B-splines of degree-1 are scaled, so that their functional value is between (0, 1);
a knots of a B-spline of degree-1 is where it bends;
B-splines of degree-1 are compact, and are only non-zero over (no more than) three adjacent knots.

You can get the (recursive) expression of B-splines from Definition of B-spline. B-spline of degree 0 is the most basis class, while

B-spline of degree 1 is a linear combination of B-spline of degree 0
B-spline of degree 2 is a linear combination of B-spline of degree 1
B-spline of degree 3 is a linear combination of B-spline of degree 2

(Sorry, I was getting off-topic...)

Your linear regression using B-splines:

y ~ bs(x, degree = 1, knots = 0)

is just doing:

y ~ b1 + b2

Now, you should be able to understand what coefficient you get mean, it means that the spline function is:

-5.12079 * b1 - 0.05545 * b2

In summary table:

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)                       4.93821    0.16117  30.639 1.40e-09 ***
bs(x, degree = 1, knots = c(0))1 -5.12079    0.24026 -21.313 2.47e-08 ***
bs(x, degree = 1, knots = c(0))2 -0.05545    0.21701  -0.256    0.805

You might wonder why the coefficient of b2 is not significant. Well, compare your y and b1: Your y is symmetric V-shape, while b1 is reverse symmetric V-shape. If you first multiply -1 to b1, and rescale it by multiplying 5, (this explains the coefficient -5 for b1), what do you get? Good match, right? So there is no need for b2.

However, if your y is asymmetric, running trough (-5,5) to (0,0), then to (5,10), then you will notice that coefficients for b1 and b2 are both significant. I think the other answer already gave you such example.

Reparametrization of fitted B-spline to piecewise polynomial is demonstrated here: Reparametrize fitted regression spline as piece-wise polynomials and export polynomial coefficients.

Interpreting spline results

You can reverse-engineer the spline formulae without having to go into the R code. It suffices to know that

A spline is a piecewise polynomial function.
Polynomials of degree $d$ are determined by their values at $d+1$ points.
The coefficients of a polynomial can be obtained via linear regression.

Thus, you only have to create $d+1$ points spaced between each pair of successive knots (including the implicit endpoints of the data range), predict the spline values, and regress the prediction against the powers of $x$ up to $x^d$. There will be a separate formula for each spline basis element within each such knot "bin." For instance, in the example below there are three internal knots (for four knot bins) and cubic splines ($d=3$) were used, resulting in $4\times 4=16$ cubic polynomials, each with $d+1=4$ coefficients. Because relatively high powers of $x$ are involved, it is imperative to preserve all the precision in the coefficients. As you might imagine, the full formula for any spline basis element can get pretty long!

As I mentioned quite a while ago, being able to use the output of one program as the input of another (without manual intervention, which can introduce irreproducible errors) is a useful statistical communication skill. This question provides a nice example of how that principle applies: instead of copying those $64$ sixteen-digit coefficients manually, we can hack together a way to convert the splines computed by R into formulas that Excel can understand. All we need do is extract the spline coefficients from R as described above, have it reformat them into Excel-like formulas, and copy and paste those into Excel.

This method will work with any statistical software, even undocumented proprietary software whose source code is unavailable.

Here is an example taken from the question, but modified to have knots at three internal points ($200, 500, 800$) as well as at the endpoints $(1, 1000)$. The plots show R's version followed by Excel's rendering. Very little customization was performed in either environment (apart from specifying colors in R to match Excel's default colors approximately).

R plots

Excel plots

(The vertical gray gridlines in the R version show where the internal knots are.)

Here is the full R code. It's an unsophisticated hack, relying entirely on the paste function to accomplish the string manipulation. (A better way would be to create a formula template and fill it in using string matching and substitution commands.)

#
# Create and display a spline basis.
#
x <- 1:1000
n <- ns(x, knots=c(200, 500, 800))

colors <- c("Orange", "Gray", "tomato2", "deepskyblue3")
plot(range(x), range(n), type="n", main="R Version",
     xlab="x", ylab="Spline value")
for (k in attr(n, "knots")) abline(v=k, col="Gray", lty=2)
for (j in 1:ncol(n)) {
  lines(x, n[,j], col=colors[j], lwd=2)
}
#
# Export this basis in Excel-readable format.
#
ns.formula <- function(n, ref="A1") {
  ref.p <- paste("I(", ref, sep="")
  knots <- sort(c(attr(n, "Boundary.knots"), attr(n, "knots")))
  d <- attr(n, "degree")
  f <- sapply(2:length(knots), function(i) {
    s.pre <- paste("IF(AND(", knots[i-1], "<=", ref, ", ", ref, "<", knots[i], "), ", 
                   sep="")
    x <- seq(knots[i-1], knots[i], length.out=d+1)
    y <- predict(n, x)
    apply(y, 2, function(z) {
      s.f <- paste("z ~ x+", paste("I(x", 2:d, sep="^", collapse=")+"), ")", sep="")
      f <- as.formula(s.f)
      b.hat <- coef(lm(f))
      s <- paste(c(b.hat[1], 
            sapply(1:d, function(j) paste(b.hat[j+1], "*", ref, "^", j, sep=""))), 
            collapse=" + ")
      paste(s.pre, s, ", 0)", sep="")
    })
  })
  apply(f, 1, function(s) paste(s, collapse=" + "))
}
ns.formula(n) # Each line of this output is one basis formula: paste into Excel

The first spline output formula (out of the four produced here) is

"IF(AND(1<=A1, A1<200), -1.26037447288906e-08 + 3.78112341937071e-08*A1^1 + -3.78112341940948e-08*A1^2 + 1.26037447313669e-08*A1^3, 0) + IF(AND(200<=A1, A1<500), 0.278894459758071 + -0.00418337927419299*A1^1 + 2.08792741929417e-05*A1^2 + -2.22580643138594e-08*A1^3, 0) + IF(AND(500<=A1, A1<800), -5.28222778473101 + 0.0291833541927414*A1^1 + -4.58541927409268e-05*A1^2 + 2.22309136420529e-08*A1^3, 0) + IF(AND(800<=A1, A1<1000), 12.500000000002 + -0.0375000000000067*A1^1 + 3.75000000000076e-05*A1^2 + -1.25000000000028e-08*A1^3, 0)"

For this to work in Excel, all you need do is remove the surrounding quotation marks and prefix it with an "=" sign. (With a bit more effort you could have R write a file which, when imported by Excel, contains copies of these formulas in all the right places.) Paste it into a formula box and then drag that cell around until "A1" references the first $x$ value where the spline is to be computed. Copy and paste (or drag and drop) that cell to compute values for other cells. I filled cells B2:E:102 with these formulas, referencing $x$ values in cells A2:A102.

Excel snippet

interpretation of the output of R function bs() (B-spline basis matrix)

The matrix b

#              1         2         3
# [1,] 0.0000000 0.0000000 0.0000000    
# [2,] 0.8270677 0.0000000 0.0000000    
# [3,] 0.8198433 0.1801567 0.0000000    
# [4,] 0.0000000 0.7286085 0.2713915    
# [5,] 0.0000000 0.0000000 1.0000000

is actually just the matrix of the values of the three basis functions in each point of x, which should have been obvious to me since it's exactly the same interpretation as for a polynomial linear model. As a matter of fact, since the boundary knots are

bknots <- attr(b,"Boundary.knots")
# [1]  0.0 77.4

and the internal knots are

iknots <- attr(b,"knots")
# 33.33333% 66.66667% 
#  13.30000  38.83333

then the three basis functions, as shown here, are:

knots <- c(bknots[1],iknots,bknots[2])
y1 <- c(0,1,0,0)
y2 <- c(0,0,1,0)
y3 <- c(0,0,0,1)
par(mfrow = c(1, 3))
plot(knots, y1, type = "l", main = "basis 1: b1")
plot(knots, y2, type = "l", main = "basis 2: b2")
plot(knots, b3, type = "l", main = "basis 3: b3")

Sample Image

Now, consider b[,1]

#              1
# [1,] 0.0000000
# [2,] 0.8270677
# [3,] 0.8198433
# [4,] 0.0000000
# [5,] 0.0000000

These must be the values of b1 in x <- c(0.0, 11.0, 17.9, 49.3, 77.4). As a matter of fact, b1 is 0 in knots[1] = 0 and 1 in knots[2] = 13.3000, meaning that in x[2] (11.0) the value must be 11/13.3 = 0.8270677, as expected. Similarly, since b1 is 0 for knots[3] = 38.83333, the value in x[3] (17.9) must be (38.83333-13.3)/17.9 = 0.8198433. Since x[4], x[5] > knots[3] = 38.83333, b1 is 0 there. A similar interpretation can be given for the other two columns.

How to extract the underlying coefficients from fitting a linear b spline regression in R?

You can extract the coefficients manually from fit.spline like this

summary(fit.spline)

Call:
lm(formula = wage ~ bs(age, knots = 30, degree = 1), data = Wage)
Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                         54.19       4.05    13.4   <2e-16 ***
bs(age, knots = 30, degree = 1)1    58.43       4.61    12.7   <2e-16 ***
bs(age, knots = 30, degree = 1)2    68.73       4.54    15.1   <2e-16 ***
---

range(Wage$age)
## [1] 18 80
## coefficients of the first model
a1 <- seq(18, 30, length.out = 10)
b1 <- seq(54.19, 58.43+54.19, length.out = 10)
## coefficients of the second model
a2 <- seq(30, 80, length.out = 10)
b2 <- seq(54.19 + 58.43, 54.19 + 68.73, length.out = 10)
plot(Wage$age, Wage$wage, col="gray", xlim = c(0, 90))
lines(x = a1, y = b1, col = "blue" )
lines(x = a2, y = b2, col = "red")

If you want the coefficients of the slope like in a linear model then you can simply use

b1 <- (58.43)/(30 - 18)
b2 <- (68.73 - 58.43)/(80 - 30)

Note that in fit.spline the intercept means the value of wage when age = 18 whereas in a linear model the intercept means the value wage when age = 0.

How to interpret lm() coefficients when formula=y~exp(x)?

This code demonstrates that your approach should work:

    set.seed( 100 )

    x <- rnorm(10)
    y <- runif(10)

    m <- lm( y~exp(x) )

    cf <- coef(m)

    yp1 <- predict( m, newdata=data.frame(x=x) )
    yf <- fitted.values(m)

    stopifnot( max( abs(yp1 - yf) ) < 1e-10 )

    yp2 <- cf[1] + cf[2] * exp(x)

    stopifnot( max( abs(yp1 - yp2) ) < 1e-10 )

    cat( "Hi-Ho Silver - Away!\n" )

Export fitted regression splines (constructed by 'bs' or 'ns') as piecewise polynomials

I was constantly asked to wrap up the idea in my original answer into a user-friendly function, able to reparametrize a fitted linear or generalized linear model with a bs or ns term. Eventually I rolled out a small R package SplinesUtils at https://github.com/ZheyuanLi/SplinesUtils (with a PDF version package manual). You can install it via

## make sure you have the `devtools` package avaiable
devtools::install_github("ZheyuanLi/SplinesUtils")

The function to be used here is RegBsplineAsPiecePoly.

library(SplinesUtils)

library(splines)
library(ISLR)
fit.spline <- lm(wage ~ bs(age, knots=c(42), degree=2), data = Wage)

ans1 <- RegBsplineAsPiecePoly(fit.spline, "bs(age, knots = c(42), degree = 2)")
ans1
#2 piecewise polynomials of degree 2 are constructed!
#Use 'summary' to export all of them.
#The first 2 are printed below.
#8.2e-15 + 4.96 * (x - 18) + 0.0991 * (x - 18) ^ 2
#61.9 + 0.2 * (x - 42) + 0.0224 * (x - 42) ^ 2

## coefficients as a matrix
ans1$PiecePoly$coef
#              [,1]        [,2]
#[1,]  8.204641e-15 61.91542748
#[2,]  4.959286e+00  0.20033307
#[3,] -9.914485e-02 -0.02240887

## knots
ans1$knots
#[1] 18 42 80

The function defaults to parametrize piecewise polynomials in shifted form (see ?PiecePoly). You can set shift = FALSE for a non-shifted version.

ans2 <- RegBsplineAsPiecePoly(fit.spline, "bs(age, knots = c(42), degree = 2)",
                              shift = FALSE)
ans2
#2 piecewise polynomials of degree 2 are constructed!
#Use 'summary' to export all of them.
#The first 2 are printed below.
#-121 + 8.53 * x + 0.0991 * x ^ 2
#14 + 2.08 * x + 0.0224 * x ^ 2

## coefficients as a matrix
ans2$PiecePoly$coef
#              [,1]        [,2]
#[1,] -121.39007747 13.97219046
#[2,]    8.52850050  2.08267822
#[3,]   -0.09914485 -0.02240887

You can predict the splines with predict.

xg <- 18:80
yg1 <- predict(ans1, xg)  ## use shifted form
yg2 <- predict(ans2, xg)  ## use non-shifted form
all.equal(yg1, yg2)
#[1] TRUE

But since there is an intercept in the model, the predicted values would differ from model prediction by the intercept.

yh <- predict(fit.spline, data.frame(age = xg))
intercept <- coef(fit.spline)[[1]]
all.equal(yh, yg1 + intercept, check.attributes = FALSE)
#[1] TRUE

The package has summary, print, plot, predict and solve methods for a "PiecePoly" class. Explore the package for more.

Multivariate regression splines in R

You can't use splines::bs in this case, as it is strictly for construction of a univariate spline. If you do bs(mat) where mat is a matrix, it is just doing bs(c(mat)). For example,

mat <- matrix(runif(8), 4, 2)
identical(bs(mat), bs(c(mat)))
# [1] TRUE

This explains why you get double number of rows, when doing bs(cbind(mtcars$wt,mtcars$hp).

To create a 2D spline, the simplest way is to create additive spline:

lm(mpg ~ factor(gear) + bs(wt, knots = 5) + bs(hp, knots = 4), mtcars)

but this may not be what you want. Then consider interaction:

model <- lm(mpg ~ factor(gear) + bs(wt, knots = 5):bs(hp, knots = 4), mtcars)

The bs(wt, knots = 5):bs(hp, knots = 4) is forming row-wise Kronecker product between the two design matrices. Since bs(wt, knots = 5) is a matrix of 4 columns, and bs(hp, knots = 4) is a matrix of 3 columns, the interaction has 4 * 3 = 12 columns.

Alternatively, consider using mgcv package. In mgcv, multivariate splines can be constructed in two ways:

isotropic thin-plate splines;
scale invariant tensor product splines.

Clearly you want the second here, as wt and hp have different units. To construct tensor product splines, we can use:

library(mgcv)
fit <- gam(mpg ~ factor(gear)
                 + s(wt, bs = 'cr', k = 4, fx = TRUE)
                 + s(hp, bs = 'cr', k = 4, fx = TRUE)
                 + ti(wt, hp, bs = 'cr', k = c(4, 4), d = c(1, 1), fx = TRUE),
                 data = mtcars)

Here, I specially set fx = TRUE to disable penalized regression.

I don't want to write an extensive answer to introduce mgcv. For how s, ti and gam work, just read documentation. If you need to bridge the gap in theory, grab Simon Wood's book published in 2006: Generalized Additive Models: an introduction with R.

A practical example of mgcv usage?

I had an answer Cubic spline method for longitudinal series data which might help you get familiar with mgcv. But as an introductory example, it only shows how to work with univariate spline. Fortunately, that is also the key. Tensor product spline is constructed from univariate spline.

My other answers related to mgcv is more of theoretical aspect; while not all my answers related to spline is making reference to mgcv. So that question & answer is the best I could give you at this stage.

Would the scale invariant tensor product splines be equivalent to radial smoothing or would that be the isotropic thin-place splines?

Radial smoothing is equivalent to thin-plate spline, as the basis function for a thin-plate spline is radial. That is why it is isotropic and can be used in spatial regression.

Tensor product spline is scale invariant, as it is constructed as (pairwise) multiplication of univariate spline basis.

How to Interpret Lm() Coefficient Estimates When Using Bs() Function for Splines