# GLM tip getting non-linear with splines

September 20, 2019
By

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# GLM tips: get non-linear with splines

This tip is great for a quick non-linear test, before you go all the way
with a GAM or parametric non-linear model.

You’ll need the `splines` library, which comes shipped with R anyway.

First, let’s make up a bit of count data. The underlying ‘true’ model
will be poisson (think count data) with a log link (so slope estimates
are multiplicative of the poisson mean). But we’ll introduce a bit of
non-linearity.

``````n <- 100
set.seed(101)
x <- sort(rnorm(n, sd = 2))
mu <- 2 + 0.1*x - 0.6*x^2 + 0.18*x^3#linear predictor

y <- rpois(n, exp(mu))
plot(x, y)
lines(x, exp(mu))
`````` Now, we could just fit a polynomial, but for real data we wouldn’t know
the mean structure was generated as a polynomial. So we might want to
use something a bit more flexible, like a cubic spline.

So, here’s how to make a cubic spline. We just need to choose the
degrees of freedom. A DF of of 1 will give us a linear fit, higher DFs
allow more bends (‘knots’). We’ll fit a log-linear model, a model with
df = 2 and a model with df= 3 Given we generated our data with a cubic
polynomial, we’d expect the 3 df model will do best

``````library(splines)
#log linear model
m1 <- glm(y ~ x, family = "poisson")
m1pred <- predict(m1, type = "response")

#non-linear models
m2 <- glm(y ~ ns(x,2), family = "poisson")
m2pred <- predict(m2, type = "response")

m3 <- glm(y ~ ns(x,3), family = "poisson")
m3pred <- predict(m3, type = "response")

par(mfrow = c(1,3))
plot(x, y, main = "DF = 1")
lines(x, exp(mu), lwd = 2, col = "grey")
lines(x, m1pred, col = "orange", lwd = 2)

plot(x, y, main = "DF = 2")
lines(x, exp(mu), lwd = 2, col = "grey")
lines(x, m2pred, col = "purple", lwd = 2)

plot(x, y, main = "DF = 3")
lines(x, exp(mu), lwd = 2, col = "grey")
lines(x, m3pred, col = "darkblue", lwd = 2)
`````` The orange line is the naive linear fit, it basically shows no trend.
The purple line (df = 2) does better, but misses the kick up at the end.
The blue line looks closest to the ‘true’ mean function (grey line).

The purple and blue lines are our spline fits with 2 and 3 knots
respectively. Both clearly capture the non-linearity. The grey line is
the ‘true’ mean structure we created above. So our splines just peaks a
bit too hard, but does get the shape right.

We can convince ourselves that the 3 df splines model is better with AIC

``````AIC(m1)

##  492.1094

AIC(m2)

##  434.6807

AIC(m3)

##  362.3653
``````

Yep, the 3 spline model has a much lower AIC despite using more model
D.F., so is better.

## Transfer your spline skills, anywhere

A nice feature of this cubic spline trick is that you can use it
anywhere that takes a model matrix as input. So it will work with `glm`,
`glmer`, `lmer` and any Bayesian GLM method you care to use. Just apply
your normal model selection criteria to find the ‘best’ number of knots.

First, set the knots using `x`:

``````library(splines)
xcs <- ns(x, 3) #3 knots!

##                1         2           3
## [1,]  0.00000000 0.0000000  0.00000000
## [2,] -0.04388799 0.1157184 -0.07117461
## [3,] -0.04775128 0.1262731 -0.07766646

plot(x, xcs[,1], type = 'l', ylim = c(min(xcs), max(xcs)))
lines(x, xcs[,2], col = "salmon")
lines(x, xcs[,3], col = "steelblue")
`````` This splits `x` into three covariates (note the new matrix `xcs` has
three columns), which I’ve plotted above. The cubic spline algorithm
puts bends in the new covariates according to the density of the data.

We can use these new covariates in our model and `glm` will estimate a
coefficient for each one.

Because the `xcs` are a non-linear functions of `x`, fitting a model
against them means we can mix the curves to get a non-linear fit. We do
lose a few degrees of freedom though, because now x is three covariates,

We can just use `xcs` in our model formula as a covariate.

But I prefer just to put the `ns` command directly into the model
formula, that way it is easy to change the knots, as we did above.

If you really want to get into non-linear trend fitting, you should use
a generalized additive model (GAM), such as from package `mgcv`.

In fact, GAM can fit these kind of splines for you too (and more), it
just uses a different method to select the number of knots.

But the splines trick is still handy. For instance, you can quickly
modify an existing linear model to have a non-linear spline. Or use it
in a Bayesian linear model that doesn’t have a GAM equivalent (you can
also spit out model matrices from
mgcv
and use them in
Bayesian models if you want other types of splines).

So I hope you enjoying splining.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.