Gavin Simpson
LightNews LightNews
gsimpson.bsky.social
Gavin Simpson
@gsimpson.bsky.social
2.6K followers 690 following 110 posts
(Palaeo)[ecologist | limnologist] & #fakeStatistican, #rstats user, wielder of #GAMs. He/him/his. Opinions mine…
Posts Replies Media Videos
gsimpson.bsky.social
🌟 plot.gam() also has a new theme: scheme=2

This draws plots with 68 & 95% intervals (by default) and has a ggplot2-like grey plot background
Example plot showing the partial effects of the estimated smooths in a GAM with terms y ~ x0 + s(x1) + s(x2) + x(x3). Three panels are shown, one per smooth. The first panel is for s(x1) and shows an increasing, slightly curved function that increases slightly more as x1 increases. The second panel shows the partial effects for s(x2), which shows much more variation, being initially negative at x2=0, rising to a peak around x2=0.35, then falling to ~-1 at x2=0.5, then remaining roughly flat followed by a gradual decline thereafter. The final panel shows the partial effect of s(x3), which is ~0 for all x3, with uncertainty band covering 0 too. This predictor is known to have no effect in the simulated data. The uncertainty in the partial effects is shown by two credible interval bands; a dark blue central band is a 68% Bayesian credible interval, while the lighter blue outer interval is a 95% Bayesian credible interval. The background of each panel is light grey with white grid lines, in a similar style to ggplot2's default theme.
November 12, 2025 at 11:28 AM
gsimpson.bsky.social
🌟 plot.gam() gains a deriv argument, which if TRUE plots derivatives of univariate smooths instead of the usual partial effect plots

Partial effect plots can be confusing. With deriv=TRUE you see the change in Y (η) for a small change in X, which is comparable with usual interpretations of model 𝛽
Example plot showing the derivative of the estimated smooths in a GAM with terms y ~ x0 + s(x1) + s(x2) + x(x3). Three panels are shown, one per smooth. The first panel is for s(x1) and shows a positive derivative that increases slightly as x1 increases. The second panel shows the derivative for s(x2), which shows much more variation, being initially strongly positive from x2=0, moving gradually to strongly negative by x2=0.35, then back to 0 at x2=0.5, then remaining roughly flat thereafter. This reflects the strongly peaked nature of the fitted smooth s(x2). The final panel shows the derivative of s(x3), which is ~0 for all x3, with uncertainty band covering 0 too. This predictor is known to have no effect in the simulated data. The uncertainty in the derivatives is shown by two credible interval bands; a dark blue central band is a 68% Bayesian credible interval, while the lighter blue outer interval is a 95% Bayesian credible interval. The background of each panel is light grey with white grid lines, in a similar style to ggplot2's default theme.
November 12, 2025 at 11:28 AM
gsimpson.bsky.social
A new release of the mgcv #RStats 📦 is out on CRAN and Simon Wood (U Edinburgh) has added some significant new features despite the small bump in version number:

🌟 scasm() for estimating GAMs with shape constrained smooths. Can be used with any family & smoothness selection is via the EFS method
Model is: b3 <- scasm( y ~ s(x0, bs = "bs", k= k) + s(x1, bs = "sc", xt = "m+", k = k) + s(x2, bs = "bs", k = k) + s(x3, bs = "bs", k = k), family=poisson, bs=200 ) The second smooth `s(x1) is a shape constrained smooth with a positive monotonicity constraint (xt = "m+"). The `bs = 200` arguments uses 200 boostrap samples, which generates bootstrap distributions for each coefficient in the model. These bootstrap samples respect the shape constraints, while the usual +/- 2 SE credible intervals may not. The uncertainty in the partial effects is shown by two credible interval bands; a dark blue central band is a 68% Bayesian credible interval, while the lighter blue outer interval is a 95% Bayesian credible interval. The background of each panel is light grey with white grid lines, in a similar style to ggplot2's default theme.
November 12, 2025 at 11:28 AM
gsimpson.bsky.social
Just updated my manuscript on using #GAMs in #AnimalScience, now on arXiv: doi.org/10.48550/arX...
🐄🐖🪶

Extended examples now show how GAMs go beyond prediction, helping estimate biologically meaningful traits from data.

Code: github.com/gavinsimpson...

🧪 #RStats #mgcv #Statistics #OpenScience
Results of model fitting to the average daily fat content data from @Henderson1990-bd. a) observed average daily fat content (points) and estimated lactation curves from Wood's [-@Wood1967-re] model, a Tweedie GLM, and a Tweedie GAM (lines) with associated 95% confidence (Wood's model) or 95% credible intervals (GLM and GAM). Response residuals for Wood's model (b), Tweedie GLM (c), and Tweedie GAM (d), plus scatter plot smoothers (lines) and 95% credible intervals (shaded ribbons). The fitted lactation curves are like an inverted U, with an extended longer tail to the right (later in lactation). The GAM curve fits the data well, but the fitted curves from Wood's model and the GLM equivalent do not provide good fits to the data, and over predict the amount of fat produced at the peak of lactation, and only grossly capture the decline in fat production later in lactation. The remaining panels show the raw response residuals for the three models, drawing attention to the poor fit; for Wood's model and the GLM there is significant pattern in the residuals, while for the GAM no residual pattern is observed. Quantities of interest derived from Wood's model, a Tweedie GLM, and a Tweedie GAM fitted to the lactation data example: a) the estimated week of peak average daily fat content, b) the estimated average daily fat content at the peak, and c) the rate of change (first derivative) of the lactation curve estimated at a point that is midway between the peak fat content and the end of lacation. The points are the estimated values and the lines are a 95% uncertainty interval. The uncertainty interval is based on the 0.025 and 0.975 percentiles of the bootstrap distribution of model coefficient estimates (Wood's model) or of the posterior distribution (GLM and GAM). Each panel shows three point estimates and an uncertainty range. The three points are the estimates from a GAM, a GLM, and Wood's lactation model. The first panel shows the estimated timing of the peak of lactation, with the GAM capturing the fact that the peak in the data occurs much later in lactation (~ week 11) while the other two models confidently estimate that the peak is in week ~8-9. The GAM estimate has a much wider credible interval, which does include the estimates of Wood's model & the GLM at the extreme end. This reflects the uncertainty in the estimation of the peak timing arising from the data having a wide flat peak. The other panels show the estimates of fat content at the peak, which are broadly similar at ~ 0.7 kg fat per day. The final panel showing the persistency estimate shows the GAM estimate diverging from those of the GLM and Wood's model. Again, the latter two models are overly confident in their estimation of this biologically relevant parameter, despite the fited lactation curve not really following the lactation data. a) Estimated daily growth rate on November 15^th^, 2021 and 95% Bayesian credible interval for the 18 pigs in the pig growth example. b) Posterior distribution of daily growth rate on November 15^th^, 2021, for three pigs (numbers 2, 13, and 17), for whom weight observations ceased before November 1^st^, 2021. In b), the shaded region is the posterior distribution, the point, and thick and thin bars are the posterior median, and 66% and 95% posterior intervals respectively. With the fitted growth curves, we can estimate for any day what the growth rate of each pig was. In this figure I'm showing the estimated growth rate of each pig in the example on November 21st. This growth rate is the first derivative of the fitted growth curve (smooth function). I used posterior sampling to produce the posterior distribution of the growth rate for each pig. These are summarised as a point estimate (median) and ccredible interval in the first panel with most pigs growin at ~1-1.5 kg per day by November 21st, with uncertainties on the order of +/- 0.5 kg per day. The second panel shows the entire posterior distribution of the estimated growth rate for three pigs (2, 13, and 17) for whom there were no weight estimates after November 1st. Here, the model is drawing power from the other pigs to help extrapolate the growth curves for these three pigs, but pig-specific details remain, with the posterior distribution for pig 17 being much more diffuse (wider) than for either pigs 2 or 13, reflecting greater uncertainty for the former animal.
October 29, 2025 at 10:45 AM
gsimpson.bsky.social
A new visitor 🕷️ to our back yard garden 🪴 Not sure what it is yet, but for scale the abdomen is about as large as the end of a human male’s little finger
A spider, hanging from top to bottom of the image. The abdomen is is patterned with a mix of browns, from cream, through tan, to dark brown A spider, hanging from top to bottom of the image. The abdomen is is patterned with a mix of browns, from cream, through tan, to dark brown
August 18, 2025 at 5:57 PM
gsimpson.bsky.social
Packed house at Åben 🍻 for the second night of #PintOfScience in #Aarhus #Denmark 🧬 🧪

Tonight’s talks;

🎙️ Attack of the Clones - Andrew Ferenbach

🎙️ Is there a Serial Killer Gene? Victoria Twiddy
Panoramic photo of people at Åben, Aarhus, 🇩🇰, for the second night of Pint of Science 2025
May 20, 2025 at 5:45 PM
gsimpson.bsky.social
Apricot is getting into the festive spirit 😻🎄💡
Photo of a female Devon Rex cat named Apricot sat on a bundle of Christmas tree lights
December 15, 2024 at 11:45 AM
gsimpson.bsky.social
🚨 New function 🚨

@vincentab.bsky.social's marginaleffects 📦 is great for generating quantities of interest from your model. I've long wanted something similar to plot_predictions() for my gratia 📦

Introducing conditional_values() 🎉

gavinsimpson.github.io/gratia/refer... #RStats #GAMs #mgcv 1/6
The image shows the result of a call to the conditional_values() function where the GAM contains a factor by smooth. In this figure, predicted values from the model conditional on values of x2 are shown for two levels of the factor fac. As we're conditioning on two variables, we include the group means for the 2 levels of fac, plus the smooth for each category. This model also includes a smooth effect of x0. As we don't condition on this covariate, it will be set to a typical representative value (value of observation closest to the median of x0). The code to produce both images is below. This image is created by the last line of the code. # packages library("mgcv") library("gratia") # categorical conditions are also handled df <- data_sim("eg4", seed = 2) m2 <- gam(y ~ fac + s(x2, by = fac) + s(x0), data = df, method = "REML") cv <- conditional_values( m2, condition = list("fac", x2 = "fivenum") ) # plot - we see a discrete x axis with five estimates per level of fac cv |> draw() # in this example we condition on `x2` and `fac %in% c(2,3)` cv <- conditional_values( m2, condition = list("x2", fac = 2:3) ) # plot - smooths of `x2` for `fac == 2` and `fac == 3` cv |> draw() The image shows the result of a call to the conditional_values() function where the GAM contains a factor by smooth. In this figure, predicted values from the model conditional on values of the factor fac are shown for different values of the covariate x2. As we're conditioning on two variables, but focusing on fac, we summarise x2 by default using its Tukey's five number summary, hence there are 5 estimates per category (level) of fac. This model also includes a smooth effect of x0. As we don't condition on this covariate, it will be set to a typical representative value (value of observation closest to the median of x0). The code to produce both images is below. This particular image is created by the first call to the draw(). # packages library("mgcv") library("gratia") # categorical conditions are also handled df <- data_sim("eg4", seed = 2) m2 <- gam(y ~ fac + s(x2, by = fac) + s(x0), data = df, method = "REML") cv <- conditional_values( m2, condition = list("fac", x2 = "fivenum") ) # plot - we see a discrete x axis with five estimates per level of fac cv |> draw() # in this example we condition on `x2` and `fac %in% c(2,3)` cv <- conditional_values( m2, condition = list("x2", fac = 2:3) ) # plot - smooths of `x2` for `fac == 2` and `fac == 3` cv |> draw()
November 30, 2024 at 10:06 AM
gsimpson.bsky.social
Peer review fueled by some Åben SIPA at Aarhus Street Food

Passing time before I pick up my kids off the train from Copenhagen
Photo shows a manuscript, a green fine liner pen, and a pint of session IPA 🍺
November 15, 2024 at 5:06 PM
gsimpson.bsky.social
Inbox achievement unlocked 😎
Screenshot from my iPhone showing the GMail icon with 1337 unread messages
October 28, 2023 at 11:26 AM