I know similar questions have been asked before, but I've looked through them all and I'm still confused, in part because the other examples seem to be more complicated, so thought I'd throw mine up. I made a very simple GAM and plotted it, but the y axis values are the smoothed values. I'd like to make the y axis have the original units so it's easier to understand.

This is a growth curve, so I'm looking at weight by age separated by sex, with individual ID as a random effect since there are some repeated measures. The plot gives me two plots, one for males and one for females, however I'd love the y axis to be in kilograms like the original dataset is, to better understand and compare the patterns.

gammWeight <- gamm(data = Data, Weight ~ s(Age, by = Sex), 
                   random = list(IndivID=~1))

I know this is likely a simple fix, but I haven't had much luck looking and this is my first time doing GAMs so figured I'd ask the experts.

ETA: I got it! I just added shift = coef(gammWeight$gam)[1] in the plot, so now it's: plot(gammWeight$gam), shift = coef(gammWeight$gam)[1])

I figured I'd keep this up in case anyone is searching and has the same question! I got the answer from Noam Ross's github tutorial.


1 Answer 1


The general solution is to use the predict() for a set of evenly spaced values of Age, repeated for each level of Sex.

But to just get a plot, you could do

plot(gammWeight$gam, shift = coef(gammWeight$gam)[1])

That will only show estimates for the first level of Sex, so the predict() approach is much better as you get actual estimates for all levels of that variable if you create the input data you pass to newdata correctly.

By the way, your model is missing the group means. You want Weight ~ Sex + s(Age, by = Sex)

To use predict() one first needs to create the new data at which the model will be evaluated, then we predict on the link scale, form the required credible interval, then transform to the response scale with the inverse of the link function used in the model. With a family = gaussian() model we can skip the last step as the link in such models is the identity function.

Hence (pseudo-ish code - not tested as there's no reproducible example here):

## new data - sequence of 200 values over range of Age
##          - all levels of Sex
##          - choose a dummy value of the subject
## expand.grid provides all combinations so the Age sequence is repeated for
## each level of Sex
newd <- with(Data, expand.grid(Age = seq(min(Age), max(Age), length = 200),
                               Sex = factor(levels(Sex), levels = levels(Sex)),
                               IndivID = IndivID[1]))

## predict
p <- predict(gammWeight$gam, newdata = newd, se.fit = TRUE)
p <- as.data.frame(p)
pred <- cbind(newd, p)
pred <- transform(pred, lower = fit - (1.96 * se.fit),
                        upper = fit + (1.96 * se.fit))

Now you can plot fit against Age, colouring by Sex say, to get actual model predicted values of the expected weight as a function of Age for each level of Sex. The plot() + shift approach is only adding the constant for the reference level of Sex to each smooth drawn. In other words, it's not showing the differences between sexes in their respective mean weights.

The above assumes you are using the corrected model Weight ~ Sex + s(Age, by = Sex); your model as you have it now is likely wrong because the smooth for each sex is trying to capture both the nonlinear way in which weight varies with age but also to fit the mean weight for that sex.

  • $\begingroup$ thank you for the response! Can you please elaborate on how I should use predict() for plotting? Or do you mean in comparing the curves like in your (super helpful) post? BTW I do have a basic question about that post specifically if you don't mind me asking! I've been using it as my main resource since I'm comparing male and female growth. Super basic q: for the pairwise differences of the curves, what exactly is different when the CI excludes 0? Is it the growth rate? I thought it was the mean at that point but my data suggests otherwise. I've been searching but can't get it. $\endgroup$
    – GRSB.1
    Commented Apr 14, 2022 at 21:05
  • $\begingroup$ Comparing the curves statistically via differencing is a complementary analysis you can do, but it is different from just visualizing the estimated age curves for the sexes. the predict() approach will allow you to do that while the plot() + shift won't. $\endgroup$ Commented Apr 20, 2022 at 9:14
  • $\begingroup$ If you've followed the posts exactly (or used gratia::difference_smooths() to do it for you) &mdash; i.e. you have included the Sex parametric effect for the group means &mdash; then what the difference smooth approach is doing is showing you the difference between the two estimated smooths. As the group means are correctly modelled by parametric factor terms, you could interpret the difference smooth as the estimates difference in growth rates. 1/2 $\endgroup$ Commented Apr 20, 2022 at 9:17
  • $\begingroup$ Some users want differences of predicted values (i.e. including the group means plus the smooths) and this can be done by not zeroing out the relevant columns in the Xp matrix for the relevant parametric terms (i.e. don't zero out the intercept and the parametric term/column for the other sex combination being compared). Then the difference would be in terms of the expected value of the response, which in your case would be the difference in the mean weight as a function of age for the sexes. $\endgroup$ Commented Apr 20, 2022 at 9:20
  • $\begingroup$ Thank you so much! This is all super helpful! I've made a function where I don't zero out the parametric columns (so just not including the X[,!grepl... line). So on those graphs, when the CI doesnt include 0, that means the means are significantly different at those ages, just like the curves are different in the original function? $\endgroup$
    – GRSB.1
    Commented Apr 25, 2022 at 19:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.