Find difference in means at particular points in a by-factor GAMM

data visualizationgeneralized-additive-modelmgcvr

This is a bit long but I wanted to give enough photos and code to explain what I'm doing and the issues I'm having. The main question is in the last paragraph.

I'm examining differences in male and female growth in a wildlife species using GAMMs, and while the GAMM output compares male and female growth across the curve (via the parametric coefficient telling me the significance of the partial effect of sex), I am wondering how I can compare the difference at different ages/x values. It was suggested I use this method because I have some repeat measures of individuals but not all repeats (we captured and measured whoever we saw on measurement day each year, so some animals are measured once in their life and some are represented in >5 age categories).

I've tried a few methods but all are giving me a different answer than my GAM is, and I'm wondering which is correct.

My code is as follows:

gamm(data = data, Weight ~ Sex + s(Age, by = Sex), random = list(SubjectID = ~1), method = "REML")

This is the output:

And here are the effects plots:

(Note that some of the graphs and outputs say "AgeReset10Up", this is just because I have a few age categories, I changed it to "Age" in a few graphs and in this description for ease of reading but they're the same thing).

The parametric coefficient and partial effects plot tells me that there is a significant effect of sex, saying that over the growth curve females weigh more than males, but I want to know at what ages this is actually true, so that I can say for example females are significantly bigger than males at age x but not age y. I tried to find this out in three ways, but all say that there are no differences.

First, when I make a ggplot of the raw data and the GAM line overlapping, the CIs overlap at all points, which seems to indicate there is no significant difference in weight at any age.

Second, when I run a linear regression
lmer(Weight ~ Sex*Age + (1|SubjectID), data = data), there is no point at which the age*sex interaction that is significant, see below. Running the LMER was suggested to me as a way to find difference in means at each age, kind of like running individual t-tests (the reference age is age 10 btw).

And third, I had used Gavin Simpson's awesome function and tutorial to find the differences in the smooths, and when I posted a question about it he kindly suggested that I modify that code by not zeroing out the intercept and parametric lines and that that would give me a difference in means when the CI excludes 0, like in the post. I ran that, and there is no point at which the CI excludes 0, suggesting the means are the same at every point.

How is it that there's an overall significant effect of sex, with females bigger than males in the GAM, but when examining the actual mean at each age using several different methods there's no difference? A friend who specializes in GAMs said that the GAM output saying there's a significant difference is right, and that the other methods must be off. So am I just misunderstanding something? Is there a better way to find means?

Best Answer

I believe this

... and there is no point at which the CI excludes 0, suggesting the means are the same at every point

is one source of confusion.

On the face of it, the data and smooths suggest that the difference of means is not exactly zero. Put another way, whatever difference in means there is between the two groups is small relative to the uncertainty in the estimates of those mean. You observed a difference but cannot say with any level of confidence that this observed difference is a real difference in the population of subjects you might have sampled.

Why are each of the model terms significant but the difference not significant? This comes from propagating the uncertainty in all of the things you estimated into the computation of the difference of means as a function of age.

Controlling for the smooth effect of Age by sex, you detect a small difference between the two groups in the mean response. The two smooth terms were assessed to be unlikely to have arise if the true functions were flat constant functions.

Computing a difference of means as a function of age requires you to bring together all these estimated effects and propagate their uncertainties into your estimates of the difference in the response between the two groups as a function of age.

You could try a more direct estimation of the thing you are interested in by using the ordered-factor method.

data <- transform(data, oSex = ordered(Sex))
m <- gam(Weight ~ oSex + s(Age) + s(Age, by = oSex) + s(Subject, bs = "re"),
         method = "REML"

where s(Age) is the smooth effect of Age for the reference level of Sex, while s(Age, by = oSex) is the smooth difference between the smooth effect of Age in the reference level and the other level of Sex.

This shouldn't change the outcome much though; I suspect you'll see a small effect of Sex overall once you condition on the smooth effects of Age in both groups, but you won't find large ("significant") differences when you ask the more specific question of "At what ages are the groups different?".

Related Solutions

Solved – GAM factor smooth interaction–include main effect smooth

You need to be careful with ordered factors here in mgcv as they aren't doing what I think you want to be fitting.

If you pass an ordered factor to by, then gam() etc set up a smooth for all the levels except the reference level, and further more they are set up as smooth differences between the reference level and the level for a specific smooth. What is happening in your first model is that the reference level of age is modelled as a constant term (it is the intercept), with the effect of age for the other levels being smooth differences from this constant.

In the second model, you add s(age), which then models the smooth effect of age in the reference level. Now, the by smooths model smooth differences from this no-longer-constant reference smooth.

I suspect that in the second model, all the levels of sex respond similarly to age hence there are no large deviations from the smooth for the reference level of sex and hence the terms are not significant. In the first model, the effect of age for the reference level was constant, so the difference smooths picked up the actual non-linear effect of age and hence were significantly different from zero.

If you just want to estimate a model with separate smooth function of age for each level of sex I would use an unordered factor (factor(..., ordered = FALSE), not ordered() or factor(..., ordered = TRUE). The the model would be:

y ~ fsex + s(age, by = fsex)

where fsex is the unordered factor.

If you want the model to be explicitly set up like ANOVA contrasts (estimate an effect for the reference level then have differences between individual levels and the reference), then you need to fit the model as per your second example with and ordered factor

y ~ osex + s(age) + s(age, by = osex)

where osex is the ordered factor. But note that in this model, s(age) is not the main smooth effect of age. It is the smooth effect of age in the reference level of osex.

Change y-axis value in R GAM plot from smoothed to actual response variable

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.

Best Answer

Related Solutions

Solved – GAM factor smooth interaction–include main effect smooth

Change y-axis value in R GAM plot from smoothed to actual response variable

Related Question