Regression models on non-linear data

• Regression is a method for predicting a variable, Y, using another, X

Equation of a straight line (1)

$$y= \alpha + \beta x + \epsilon$$

Equation of a straight line (2)

$$y= 2 + 1.5 x + \epsilon$$

What about nonlinear data? (1)

What about nonlinear data? (2)

What about nonlinear data? (3)

What about nonlinear data? (4)

What about nonlinear data? (5)

Degrees of freedom (df)

Within a model, how many 'parts' are free to vary?

E.g. If we have 3 numbers and we know the average is 5:

If we have 2 and 7, our final number is not free to vary. It must be 6: $$\frac{2 + 7 + 6}{3} = 5$$ This means our 'model' is constrained to $n-1$ degrees of freedom

In regression context:

• The number of data points in our model ( $N$ ) limits the df
• Usually the number of predictors ( $k$ ) in our model is considered the df (one of these is the intercept)
• "Residual df" are points left to vary in the model, after considering the df:

$$N-k-1$$ Helpful post on CrossValidated

Overfitting

When our model fits both the underlying relationship and the 'noise' percuiliar to our sample data

• You want to fit the relationship, whilst minimising the noise.

• This helps 'generalizability': meaning it will predict well on new data.

If we allowed total freedom in our model, e.g. a knot at every data point. What would happen?

What if we could do something more suitable?

If we could define a set of functions we are happy with:

• We could use coefficients to 'support' the fit
• Could use penalties to restrict how much they flex

Splines

• How do you (manually) draw a smooth curve?

Mathematical splines

• Smooth, piece-wise polynomials, like a flexible strip for drawing curves.
• 'Knot points' between each section

How smooth?

Can be controlled by number of knots $(k)$, or by a penalty $\gamma$.

Generalized Additive Model

• Regression models where we fit smoothers (like splines) from our data.
• Strictly additive, but smoothers can describe complex relationships.
• In our case:

$$y= \alpha + f(x) + \epsilon$$

$$g(\mu_i) = A_i \theta + f_1(x_1) + f_2(x_{2i}) + f3(x_{3i}, x_{4i}) + ...$$
Where:

What does that mean for me?

• Can build regression models with smoothers, particularly suited to non-linear, or noisy data

• Hastie (1985) used knot every point, Wood (2017) uses reduced-rank version

Issues

• We need to chose the right dimension (degrees of freedom / knots) for our smoothers
• We need to chose the right penalty ( $\lambda$ ) for our smoothers

Consequence

• If you penalise a smooth of $k$ dimensions, it no longer has $k-1$ degrees of freedom as they are reduced
• 'Effective degrees of freedom' - the penalized df of the predictors in the model.

Note:
$df(\lambda) = k$, when $\lambda = 0$
$df(\lambda) \rightarrow 0$, when $\lambda \rightarrow \infty$

mgcv: mixed gam computation vehicle

• Prof. Simon Wood's package, pretty much the standard
• Included in standard R distribution, used in ggplot2 geom_smooth etc.
• Has sensible defaults for dimensions
• Estimates the ideal penalty for smooths by various methods, with REML recommended.

library(mgcv)
my_gam <- gam(Y ~ s(X, bs="cr"), data=dt)

• s() controls smoothers (and other options, t, ti)
• bs="cr" telling it to use cubic regression spline ('basis')
• Default determined from data, but you can alter this e.g. (k=10)
• Penalty (smoothing parameter) estimation method is set to (REML)

Model Output:

summary(my_gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ s(X, bs = "cr")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  43.9659     0.8305   52.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##        edf Ref.df     F p-value    
## s(X) 6.087  7.143 296.3  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.876   Deviance explained = 87.9%
## GCV = 211.94  Scale est. = 206.93    n = 300

Check your model:

gam.check(my_gam)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 4 iterations.
## The RMS GCV score gradient at convergence was 1.107369e-05 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##        k'  edf k-index p-value
## s(X) 9.00 6.09     1.1    0.97

Check your model:

gam.check(my_gam)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 4 iterations.
## The RMS GCV score gradient at convergence was 1.107369e-05 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##        k'  edf k-index p-value
## s(X) 9.00 6.09     1.1    0.93

Is it any better than linear model?

my_lm <- lm(Y ~ X, data=dt)
anova(my_lm, my_gam)

## Analysis of Variance Table
## 
## Model 1: Y ~ X
## Model 2: Y ~ s(X, bs = "cr")
##   Res.Df   RSS     Df Sum of Sq      F    Pr(>F)    
## 1 298.00 88154                                      
## 2 292.91 60613 5.0873     27540 26.161 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(my_lm, my_gam)

##              df      AIC
## my_lm  3.000000 2562.280
## my_gam 8.087281 2460.085

Yes, yes it is!

Advanced elements (1)

The name of the model class is 'Generalized Additive Model' (GAM) - like the 'Generalized Linear Model' (GLM)

• This means we can fit other types of response distribution: e.g. Poisson or Negative Binomial for counts, binomial ("Logistic regression") for binary outcomes etc.
• This is achieved using a 'link' function: log for Poisson, logit for binomial.

The main addition here is the family argument:

my_logistic_gam <- gam(Y ~ s(X, bs="cr"), data=dt, family = "binomial")

As with a GLM, the results are on the scale of the link function, so have to transform predictions back to the response scale:

predict(my_logistic_gam, newdata, type="response")

Advanced elements (2)

Interactions

Sometimes our predictors are not independent of each other, and you may want to isolate the individual and combined effects.

• Smoking and drinking were highly linked: as people would often drink and smoke at pubs.
• A stroke drug trial might have different effects on recover depending on severity of the stroke.

In a GLM, you would add this with a multiplicative interaction term:

my_interaction_model <- glm(Y ~ X * Z , data=dt, family = "binomial")

You can't do the same thing in an additive model, but you can make multidimensional smoothers:

# Within 's' - assuming they are on the same scale
my_interaction_gam1 <- gam(Y ~ s(X,Z), data=dt, family = "binomial")
# As a tensor product, using te, if on different scales
my_interaction_gam2 <- gam(Y ~ te(X,Z), data=dt, family = "binomial")

Summary

• Regression models are concerned with explaining one variable: y, with another: x

• This relationship is assumed to be linear

• If your data are not linear, or noisy, a smoother might be appropriate

• Splines are ideal smoothers, and are polynomials joined at 'knot' points

• GAMs are a framework for regressions using smoothers

• mgcv is a great package for GAMs with various smoothers available

• mgcv estimates the required smoothing penalty for you

• gratia or mgcViz packages are good visualization tool for GAMs

References and Further reading:

GitHub code: https://github.com/chrismainey/GAMworkshop

Simon Wood's comprehensive book:

• WOOD, S. N. 2017. Generalized Additive Models: An Introduction with R, Second Edition, CRC Press.

Noam Ross free online GAM course: https://noamross.github.io/gams-in-r-course/

Kim Larson's great intro article: https://multithreaded.stitchfix.com/blog/2015/07/30/gam/

• HARRELL, F. E., JR. 2001. Regression Modeling Strategies, Springer-Verlag.

• HASTIE, T. & TIBSHIRANI, R. 1986. Generalized Additive Models. Statistical Science, 1, 297-310. 291

• HASTIE, T., TIBSHIRANI, R. & FRIEDMAN, J. 2009. The Elements of Statistical Learning : Data Mining, Inference, and Prediction, Springer.