Selecting GAM model link function and autocorrelation (mgcv)

Question

I'm using a GAM model for continuous dependent and two independent variables recorded at reasonably regular intervals from multiple subjects. The dependent variable Y is strictly >0 and positively skewed.

M1 = bam(log(Y)~s(A)+s(B)+s(id, bs='re'))
M2 = bam(Y~s(A)+s(B)+s(id, bs='re'), link=Gamma(type='log'))

First question is whether it is preferable to use the log transform of the dependent variable (M1) with a gaussian-identity link function in the gam, or use the untransformed data and employ a Gamma-log link function (M2). My reading suggests that while transforming the variable is the traditional method the latter approach of using a generalised link function is these days generally viewed as 'better', and gives more accurate estimates of predicted mean, SD etc. Though the explained deviance is little different between the models (both ~58%), AIC is far lower when using the log-transformed dependent variable (M1:~2000 vs M2:16000). For completeness I did try a gaussian-log link function but it did not fit well.

Second question: there is autocorrelation in the dependent variable so the next model includes an AR(1) rho=0.5 term.

M3 = bam(log(Y)~s(A)+s(B)+s(id, bs='re'), AR.start=ind, rho=0.55)

For a start, the AR function does not work with non-gaussian link function which raises the question of whether it is valid to use an autocorrelation model with log-transformed variables? The AIC from the log-transformed-AR(1) model is -9000. Is this reliable? Again, deviance explained is little different.

The fitted smooth terms are qualitatively similar for the non-autocorrelation models (M1 and M2) so practically it doesn't make much difference to my conclusions. However, the smooths are significantly different for the AR model (M3).

What would be the advice on choosing amongst the models?

Answer 1

Q1: In the $log(y)$ model with Gaussian errors you are modelling the mean of log-transformed $y$, not the mean of $y$, and when back-transformed to the original scale the two won't coincide. That's often why GLM-like models are favoured as you can model on the scale you actually want to model on, not some transformed scale. With the $log(y)$ model there's also nothing forcing model outputs or summaries to be strictly positive, whereas with the Gamma there would be.

Q2: right; in the Gaussian model we can think of the data as being correlated correlated Gaussian random variables and represent the correlations via a correlation matrix. There isn't an easy way to do the same thing in the GLM context.

There is nothing wrong with fitting AR models to log-transformed data or to the residuals of a regression model fitted to such data.

A single AIC is not very informative; you need to compare two or more AIC values.

If A and B are not time variables then you can easily undersmooth (overfit) data that is autocorrelated; when you say the smooths are "significantly different" for M3, I might hazard a guess that they are much smoother/less wiggly than either of the other two models?

If you suspect residual autocorrelation, or you aren't explicitly modelling time in your linear predictor, I'd suggest that fitting with the AR(1) is appropriate.

You can use the gamm() function to fit the M1 type model but estimating $\rho$ rather than fixing it a priori. You'll need to change from the spline-like ranefs to using the random argument to specify the ranefs in the model and use correlation = corAR1(form = ~ time) or correlation = corAR1(form = ~ time | id) to get an AR(1) or an AR(1) nested within each subject (id) respectively.