Introduction to Generalized Additive Models with R and mgcv

[@thexKCY]

This video is single handedly responsible for me being able to complete one of my dissertation chapters. Thank you for the work you do, it makes such an impact!

[@ftboth]

Thank you! Glad the video was useful.

[@beagle989]

the animation of the functions being added and weighted was amazing, omg

[@ftboth]

Thanks!

[@jennykeeping8918]

I tried learning about GAMS back in 2016 for my Masters and my god, with the skin of my teeth I got by! What an amazing spread of FREE tools now available for learning how to make good models now! Thank you SO much, PhD will now be LEAGUES better.

[@annas-h8x]

Thank you so much for this video. It is such a great and clear explanation of GAMs and is helping me so much in my PhD dissertation!

[@ftboth]

Thanks! I’m glad you found the video useful

[@lauraeli2286]

Gavin you have a brilliant teaching style and I find your explanations are always so clear! Also, the notation you use is I think the best and most straightforward for beginners like myself to understand and I think this is important because I find concepts/aims of the techniques are clear to understand. But unfortunately the notation used in some other videos I've seen detracts from this. So brilliant and much appreciated! Thank you!

[@beginnerpython6183]

Many thanks from South Korea ! One of the best video among others!

[@MichaelTamGP]

Thank you from Australia for this excellent video tutorial!

[@charleydublin7304]

This is great work, thank you for putting together such a brilliant explanation of GAM.

[@HarmonicaTool]

That was a great introduction to an important topic. Great content by a great teacher. Thank you very much. I will go to a big online bookshop and order "the book". Then I will check out your channel. Cheers from Germany.

[@michawalesiak8524]

Big thanks from Poland, best video on GAMs out there

[@doug_sponsler]

@Gavin Simpson You mention in passing (1:49:30) that, for spatial and/or temporal models, the gam.check procedure tends to indicate a need for unreasonably large k values, and that this is due to autocorrelation in the data. Could you elaborate on how you go about selecting an appropriate k value when gam.check is misleading? Perhaps this would be a good topic to address in a blog post, since I'm sure you encounter it often in your climate time series modeling.

[@falco.pellegrino]

Hi Gavin, I have been suggesting this video in my consultations for a while, and I'm finally watching this for the first time. Thank you for taking the time to create this! I noticed a typo at 8:39. Assumption 2 should be Y|X ~ N(\mu, \sgma^2) or \epsilon | X ~ N(0, \sigma^2)

[@ftboth]

Oops. Will have to fix that in my slides. Nice username 😝

[@falco.pellegrino]

@@ftboth I thought I had changed that 😅😂 I guess it shows how little I interact with people on youtube 😅

[@connorrosenberg2736]

Fantastic resource! Highly recommend watching for a wonderful and useful introduction to GAMs. Would love to see a continuation which talks about adding in priors using rstanarm.

[@redshift6668]

Thank you for this excellent video, and for the book recommendation at the end - which is also great.

[@ftboth]

Glad it was helpful!

[@hamishthecat666]

Thank you Gavin, excellent seminar and I learned a lot!

[@laposanti534]

Thanks professor, a very complete tutorial providing a good deal fo intuition.

[@pedro.melo.albuquerque]

Awesome talk! Congratulations.

[@pumpkinbeetle]

This video is a great resource, thank you for creating it. I wanted to confirm whether the standard plot diagnostics for logistic regression are not easy to use and the dharma package is recommended in this case?

[@ftboth]

Thanks; glad you liked it. As to your question: I presume you mean for binary data? If so, yes, many of the typical diagnostic plots are less than useful. The QQ-plot is the most useful of the four produced by mgcv's `gam.check()` or `appraise()` in my gratia 📦. This is because it (with suitable settings for relevant arguments (`method = "simulate"` in `appraise()` for example) it provides reference bands that reflect the kinds of deviations from the 1:1 line you'd expect if the model was correct. It achieves this by simulating residuals from the model, much like DHARMa does (just without the normalization via the ECDF). When I produced this video, DHARMa didn't work with mgcv model fits, but now thanks to Florian, DHARMa's developer and Matteo Fasiolo (developer of the mgcViz 📦) DHARMa can handle many of the families in mgcv, but currently not all of them. I now teach DHARMa as part of my week-long GAM course as a result, and hope to add randomized residuals to gratia at some point in 2024.

[@pumpkinbeetle]

@ftboth - thank you for the detailed response. I am finding analysis of binary data is quite painful. Will keep an eye out for updates to gratia and use DHARMa in the meantime. Thanks!

[@charliezhu163]

Great video! One question: Is it possible to output the explicit coefficients of the piecewise fitted spline for each factor? (as not being to see the spline functions makes it more of a black box compared to the traditional GLM)

[@ftboth]

Yes, just use `coef(my_model)` and there's a `gratia::smooth_coefs()` (in the GitHub version of gratia, soon to be on CRAN) if you want to extract the coefficients for specific smooths

[@abalfazlzareei]

Thanks so much! This was great!

[@nicholastrandberg6346]

Thanks for this excellent video!

[@teamgoulag4381]

Hello, thank you very much for this course and the fluid and very clear explanation. I wanted to know if we could apply the GAM to spatial data (altitude, longitude)? I thank you for your answers.

[@ftboth]

Yes, GAMs can be applied to such data. If you want interactions between such variables then tensor product smooths would be appropriate as altitude and longitude are in different units

[@ssperandei]

Simply fantastic! Thank you!

[@JoyLetters]

Hi! thank you for sharing this excellent video. is there any difference between generalised additive models and generalised additive mixed model?

[@ftboth]

Loosely GAMM is reserved for GAMs with random effects, but seeing as you can add a random effect as a smooth to a GAM, these distinctions get are very blurred

[@cristianrperdomogarcia5745]

Saludos, me has ayudado en algunos foros y lo agradezco. Aunque no tengo muy buen inglés lo veré muchas veces si es necesario para entender aún más los gam. Éxitos !!!

[@cristianrperdomogarcia5745]

@@ftboth I have seen the video and I understood that when we put in s (x, f, bs = 'fs') this is understood as slope and random intercept. that is, it is the same as placing s (f, bs = "re") and s (x, by = f, bs = "re." Or is there a difference?

[@TimSter15]

Without this, I wouldn't be able to do my undergraduate dissertation; you have NO idea how valuable these free seminars and videos are to students. Thank you so much Gavin for this detailed exploration! One question, around the 1:15:00 mark you've put that the mgcv package contains other non-exponential families but the Gamma distribution IS an example of an exponential... or am I getting confused with distributions here?

[@TimSter15]

@@ftboth Thanks for replying so quickly and confirming that!

[@TimSter15]

@@ftboth one other quick thing I need to confirm; how do you put a fixed factor in a GAM? If you don't specify otherwise, does it just assume that a variable you plug in is fixed?

[@TimSter15]

@@ftboth ahhh I see. Thanks

[@TimSter15]

@@ftboth for the summary, what do the parametric coefficients tell us? Is it the same as for the smoothed covariates?

[@constantin1481]

Thank you. Great (and complete) tutorial. Am I correct to say that when you use information criterion to choose between two models with a different number of knots or spline types you are basically treating them as random effects and that's why you are able to use REML=T? Also how "correct" is it to choose between fixed effects using ML and then refitting the final model with REML? If it is sound, can you point me toward research that says as much, please?

[@ftboth]

We normally recommend `method = "REML"` - it isn't fitting the model as a random effect model, but it is minimising a REML criterion that arises from viewing penalised splines as random effects. You shouldn't then compare models using their likelihoods if they differ in their fixed effects. If I did that, it was an oversight and I should have fitted the models with `method = "ML"` before doing those comparisons

[@sitendurocks]

Hey gavin! just downloaded gratia, looking forward to this project!

[@satyamgopal]

Thanks for wonderful video. Looks like selecting right basis function is key here. Any guidance for beginners.. on how to select..?

[@satyamgopal]

@@ftboth Thanks for quick response. As per direction on stack exchange, let me try with TP and CR respectively. To see how model works. I want to apply this to forecast sales based on historical data along with promotion calendar..let's see how it goes

[@ftboth]

The actual choice of basis (the type of spline) is not usually important in standard settings. The default thin plate spline is a good choice except in big data settings and removed the need for worrying about knot placement. More important is the size of the basis - how many functions to use.

[@satyamgopal]

Make sense. Thanks

[@jimdurant1192]

Thanks for great talk, it was extremely helpful. Is there any thoughts on the cenGAM package for censored data analysis with mgcv?

[@ftboth]

Great suggestion! I might do something on this; I'm not as familiar with that package but have fitted some censored GAMs using {brms} 📦 for a paper that came out a year or so ago

[@jimdurant1192]

@@ftboth I checked out the brms for my problem. Really helps. Do you ever teach in person or virtual classes on gams?

[@hakimb3855]

Wonderful presentation. Thanks a lot. @Gavin would you generously give your permission to reproduce some of your figures for a presentation?

[@ftboth]

Thanks; all the materials should be under a Creative Commons Attribution (CC-BY) license, so you are free to use whatever you want from the slides so long as you cite the source.

[@yadavghimirey]

Thank you Gavin! How about wrapping a categorical variable using bs = 'fs' argument?

[@ftboth]

Not sure of the context; an `fs` smooth, say `s(x, f, bs = "fs")` is a combination of a continuous variable (the thing you smooth) and a categorical (factor) variable so that you end up with a smooth of `x` for each level of factor `f`. If you are just talking about a term involving just a factor, then `s(f, bs = "re")` would give you a random intercept per level of `f` just like you would with `(1 | f)` in `lmer()` or `glmer()`. In general you can't smooth a categorical variable, though there is some work on smoothing ordinal factors.

[@yadavghimirey]

@@ftboth Thank you for the reply.. If I understood correctly, I can't just do s(f,bs='fs') because f is a factor variable, right? Also, is there a support group for mgcv? I didn't find any. Would be greatly helpful to ask questions.

[@yadavghimirey]

@@AnivetGavin Thank you!

[@ftboth]

@yadavghimirey (Sorry; originally replied using my teaching KZbin account) You can actually do `s(f, bs = "fs")` because Sinmon Wood has a bit of code in the constructor for the `fs` basis to catch this usage; if you do this the code just drops out to the constructor you would get if you had use `s(f, bs = "re")`, for random intercepts. No, there isn't an mgcv support group. Instead I would suggest you ask questions on StackOverflow (tag them mgcv or gam; use this for questions about the software) or CrossValidated (again tag the questions mgcv or generalized-additive-model; use this if the question is statistical in nature and not just limited to using mgcv). There are already plenty of questions under those tags at both locations so lots of existing help is out there, and I routinely answer questions on mgcv/GAMs in both places.

[@srnlo3937]

thanks for this very usefull webinar ! However I didn't found how to set the lambda coeffeicient .Is that the same than the gamma in mgcv package ? thanks !

[@srnlo3937]

@@ftboth Humm in the webinar I understood that "lambda" refers to the whole model and k to the individuals termes. In the mgcv sp should be furnished as a vector with one value for each term, so how to provide parameter taking into account the wiggliness of whole model ?

[@siosaiavaihola3252]

Hi Prof. Gavin, thank you so much for your work/presentation here, exactly what I need for my phd. Just a question Sir, can you please guide to the animation code in time 2:44:60 Galveston Bay slide 127. With appreciation.

[@ftboth]

The code is on line 1875 of the .Rmd file for the slides: github.com/gavinsimpson/intro-gam-webinar-2020/blob/master/gam-intro.Rmd

[@siosaiavaihola3252]

@@ftboth appreciate this very much

[@marianagelambidesiato7495]

thank you! Learned a lot :)

[@karunamayiholisticinc]

Can I clarify if I understand it then number of "k" means the number of basis functions between each knot? so greater the "k" more basis functions we have, greater the wiggliness?

[@ftboth]

`k` is the number of basis functions to use (before any identifiability constraints). It is not the number of basis functions between each knot; it is the total number of basis functions used to represent the function. It's best not to think of knots as the default thin plate regression spline isn't really knot-based. As for the other point, the larger `k` the larger the upper limit on possible wiggliness. It doesn't hold that increasing `k` leads to increased wiggliness because the penalty will penalize excess wiggliness. Hence my saying above that `k` controls the maximum amount of wiggliness allowed. Hope that helps?

[@karunamayiholisticinc]

@@ftboth Thanks. yes I understood it. I watched videos about basis functions and spline models from Stanford Playlist. It represents one function over entire data with different basis functions part of it. Natural spline has a knot at every data point which smoothens out due to penalty of reducing wiggliness represented by the change in slope or double derivative. In continuous cubic spline models there is discontinuity in the third derivative at knots and there are extra basis functions that gets added based on knot position but that is for continuous cubic spline model. The natural spine works differently with knots at every data point and smoothens the wiggliness to fit an overall function with effective degrees of freedom that accounts for the amount of smoothness. I hope I got I right.

[@ftboth]

What you are describing is a smoothing spline. Smoothing splines require N coefficients for N data. In the regression modelling context, using as many coefficients as there are data is i) unnecessary, and ii) would exclude models with more than 1 smooth function. Smoothing splines can use many kinds of spline. You mention a natural spline; this is formed from piecewise cubic polynomials (so just like the cubic spline you mention) but a natural spline is distinguished by have 0 second derivatives at the boundary knots (hence any extrapolation will be linear beyond the limits of the data). The above doesn't really apply to mgcv as we are never fitting smooth functions with as many coefficients as data. Instead we use fewer basis functions in GAMs. Another difference is that several of the smooths in mgcv (including the default thin plate regression spline, `bs = "tp"`) don't really work this way; these bases typically do have a knot at every (unique) data point, but as discussed, this is wasteful and not useful in a modelling context. Instead, the full basis (with a knot at each (unique) data) is subject to an eigendecomposition and the K eigenvectors associated with the K largest eigenvalues are retained as the basis. Hence even with these types of smooths we only use K basis functions. Hence my point about it best being to think of `k` as being the number of basis functions in the basis (before any identifiability constraints). These splines that use fewer basis functions are typically termed "regression splines". We borrow the idea of penalizing wigliness from smoothing splines, so the smoothing parameter(s) control the resulting wiggliness of the regression spline. Hence we use the term "penalised regression spline" to describe the general approach taken by mgcv for GAMs (and this is typical for the use of smooth functions in regression models.)

[@ramansingh-cg7rs]

Hi Gavin, many thanks for posting this online, do you know if this can be implemented in C++

[@ftboth]

I'm not aware of a C++ application that can fit these models, but much of the GAM code is C and or uses standard matrix algebra libraries.

[@dr.fragrance3185]

I had a question about random smooths. I have x- and y-coordinates. The simplified model for the smoothing terms is: Y ~ s(X) but for a random smooth, should I do: s(X, Participant) or s(Participant) or s(X, by = Participant) and what is the difference? Is one better? They all seem to make different plots!

[@ftboth]

Assuming Participant is a factor, `s(Participant, bs = "re")` isn't a random smooth. `Participant + s(X, by = Participant)` would be more like a fixed effect smooth; you get one smooth per level of `Participant` and absent any other terms in the model you wouldn't be able to predict for new levels of the factor. `s(X, Participant, bs = "fs")` would be a definition of a random smooth, but again, absent any other mode terms you wouldn't be able to predict for an unobserved level of the factor. Of your models, option 1, the fs smooth, would fit one smooth per level of the factor but all the smooths share the same smoothing parameter, so they have the same wiggliness, option 2, isn't a random smooth as mentioned, it's just a random intercept that describes the Participant mean response, and option 3, the factor by smooth, is like option 1, but each smooth gets it's own smoothing parameter so the wiggliness can vary between smooths for different levels of the factor. They plot differently because if the different assumptions and usages behind each one. As for which is better, that depends on what you are trying to achieve by fitting one of these models.

[@DYT32]

Professor Gavin Simpson, thank you for this great lecture. I'll try my best to explain myself in English: I've been trying to use mgcv::gam to describe some data that have a stepper pattern at low values of the continuous predictor, but with a smoother pattern at the high values of the continuous predictor. There is a lot of variance between levels of random effect. What kind of models do you recommend? I've explored some basis functions (tp, cc, cr, ps, gp; ad wasn't possible - the only likely is 'cr'), but I always get too much 'wiggliness' in the estimated pattern at the stepper segment - similar to your example at slide 57 figure "wrong functional form". I am very grateful for your generosity in sharing your knowledge with us and for your time in reading my question. Muito Obrigado Danilo Pereira Mori, a biologist at the theoretical ecology lab - Institute of Biosciences, University of São Paulo, Brazil.

[@caviper1]

Gracias. Muy claro.

[@danielarobles7927]

Hello professor, I have one question, in the complex example you gave at the end, why you didn't scale the predictors? It is not necessary for gams?

[@ftboth]

It's not necessary to scale the predictors in any model like this; what you might gain in comparison of the estimated effect sizes is lost elsewhere because your model is fitted to scaled covariates. As all the covariates in the Galveston Bay example were reasonably scaled already (none of the covariates was in units orders or magnitude larger/smaller than the others) there was little point in scaling them. HTH.

[@sorayaalvarez6280]

Awesome!

[@TrangVo-pz5xt]

Hi, I have a question related to the tensor smooth term. So, it is true to say that te(x,y) ~ s(x,y) + ti(x,y) if x,y are in different units (such as, x: temperature and y: precipitation)? Because I want to account the main effect of the relationship between P & T without their interactive effect, so I got confused how could I construct the model. Any help is much appreciated! Thanks

[@ftboth]

No, that's not right; te(x,z) ~= s(x) + s(z) + ti(x,z). You don't want s(x,z) if x and z are in different units. The reason it's not exactly equal is due to the extra smoothing parameters in the s(x) + s(z) + ti(x,z) form. If you just want additive main effects of P and T use y ~ s(P) + s(T)

[@haraldurkarlsson1147]

The code for reading in the data did not work for me. Possibly due to a change in the data since this video was recorded. I did manage to download some HadCRUT5.0.1.0 in, however, as .csv files and work with those.

[@ftboth]

Yeah; they've changed the format *again*

[@Alex-ll7qc]

Hi, professor - thank you so much for your great lesson! - just one quick question, how can we incorporate complex survey design features within this analysis? :)

[@ftboth]

I don't know; those designs usually require specialised software, like the survey package in R, so you'd need to look at those methods and how you might include penalised splines there.

[@franklingarciafernandez5387]

Thanks

[@bezaeshetu5454]

Dear Prof. Simson, Thank you for a nice and clear explanation. I have a few questions, 1. there are different types of spline bases(tinplate spline, cubic spline, p-spline, ...), when and with what condition we use them. 2. I have a data set and the response variable is binary (yes, no) and am trying to analyze using mgcv package in r. my question here is one variable is significant when I do gam.check and when I increase the number of knots "K", I got this error message "Error in smooth.construct.cr.smooth.spec(object, dk$data, dk$knots) : dd has insufficient unique values to support 40 knots: reduce k." I also use the tinplate spline but still the same. Could you assist please. Thank You.

[@ftboth]

Thanks Beza. 1. I have largely covered this in an Answer on CrossValidated: stats.stackexchange.com/a/305446/1390 2. You can only fit as many basis functions as you have unique values of the variable you are smoothing (minus some amount due to constraints on the spline). There is nothing you can do about this; this is an upper limit, an it may vary by 1 or 2 depending on which spline type you choose and what order of penalty is used, but that isn't going to substantially change the maximum value of `k` that you can use. This test is just a heuristic so you don't have to treat it as some hard and fast rule that you must get the p value above 0.05. There isn't anything you can do substantially to fix the issue in your case as you can't simply make `k` larger. What I would do is go back to the original `k` you used and then fit another model where you have increased `k` by some amount (doubled it for example) and then I would plot the smooths for the variable of interest from the two models. Are they similar or different? If they are similar then perhaps the test is showing a false positive and you can safely move on. I would also look to see if the other smooths change substantially or not; basically are my results sensitive to the choice of `k for this one variable? You could do more; you could take the residuals from the first model (where you used the original value of `k`) and then model those (with `family = "gaussian" now as the deviance residuals will be kind of spread about 0) and fit with `residual ~ s(x)` replacing `x` with the variable flagged by the `gam.check()` test. Is that smooth "significant" when you look at the summary? Plot that smooth and see what it looks like. Does it deviate from the 0 line markedly? Finally, at least in terms of what comes to mind right now, consider missing covariates; you could have ignored an important covariate that may be related to the variable giving you trouble in `gam.check()`.

[@bezaeshetu5454]

@@ftboth Thank You prof. 1. what is the criteria to use GAM when the family is not Gaussian (for example, family = binomial). when it is gaussian we plot the response versus the predictor then if the pattern is unknown then we use GAM. but when the response is categorical variable (0, 1) and if we do not have time variable for example to plot with, then what do we do. 2. is the plots in gam.check can be used to check the diagnostic of the model in this case. I am asking this because I saw the plots are not quite interpretable since the response is categorical.

[@fabiofavoretto4153]

Completely unrelated question... which software did you used for this webinar?

[@WWEAnimations619]

Hi, when I try to load in the data (using read_delim as in the script), I get a warning that there was one or more parsing issues and then the temperature column is just full of NAs. Do you know how to resolve this? Thanks.

[@ftboth]

They've changed the data structure. You can find updated code to read these data in the materials for my GAM course: github.com/day-5/time-series-gam-example.R

[@an1k199]

Great video, thank you! But i have a question: if i have hours after sunset (from 0 to 7, not countinous, just one value for an hour) on the x and individuals counted on the y-axis, what model should i use? I tried negative binomial Distribution and k is set on maximum (with bs = "cr" or "tp"), but the Output of gam.check always says: k-index: 0.59 and p is significant, but i cant increase k anymore. What to do next?

[@ftboth]

That typically indicates unmodelled structure in the response, so you may be missing one or more covariates that have effects on the response. Try looking at your residuals to see if that sheds some light on what might be missing

[@satyamgopal]

Gavin , Any special care need to be taken in feature engineering in GAM context..?

[@ftboth]

Beyond the usual things you need to worry about, I guess the main GAM-related thing you need to be careful about are concurved features, which is a bit like collinearity in LMs/GLMs etc but now you need to worry about features that could be combinations of smooth functions of one another, or "concurvity".

[@celestecanazza5505]

Hi, if I want to compare gam models using cross validation (mgcv package) how can I do ? Can I write the full model in cv with select=T and method=reml? I have to compare some models in term of rmse (gam, glmer,glm)

[@hansmeiser6078]

I understood this for a while, but nowadays I always get only a straight, linear line... strange.

[@yaladandanikhila2944]

Is any fuction for gam model to run all possible models as olsr in lm.. Please help me out..

[@ftboth]

I don't think so and I would not recommend this as a way of doing statistics in general

[@AnaheimElectronicsInc.]

1:42:30