Gamma family puts plot_nonlinear upside down

Interfaces
,
1

Hi,

I’m fitting a smooth, using stan_gamm4, to all positive data and have played around a fair bit with different data transformations and families. From my trials it seems like the Gamma family is best suited for my data unless I want to log transform the data and run it with the gaussian family (inverse.gaussian is a decent fit to untransformed data as well, but worse than Gamma and seems to take longer to fit). If I can avoid transformations then I would prefer that.
However, when I use plot_nonlinear function, the smooth from the Gamma model comes out upside down relative to the plot of the log transformed model smooth. Why is that and can it somehow be corrected? Does it have to do with the link function (1/mu^2) for Gamma?

I’ve uploaded some dummy/example data that works for demonstrating the issue
example.data.Rdata (8.0 KB)
.

example.gamma <- stan_gamm4(y ~ s(x, bs = 'ts'), 
                            data = example.data, family = Gamma(), 
                            warmup = 1000, iter = 3000, thin = 3, chains = 4, cores = 4, 
                            adapt_delta = 0.99)

example.data$ylog <- log(example.data$y)
example.gaussian <- stan_gamm4(ylog ~ s(x, bs = 'ts'), 
                               data = example.data, family = gaussian(), 
                               warmup = 1000, iter = 3000, thin = 3, chains = 4, cores = 4, 
                               adapt_delta = 0.99)

plot_nonlinear(example.gamma, smooths = 's(x)')
plot_nonlinear(example.gaussian, smooths = 's(x)')

Windows 10
rstanarm 2.32.1

2

Both conditional smooth plots are showing what you’re asking for, so there’s nothing you need to correct. You’re right that the “flipped” direction of the conditional smooth function is because of the differences in the link functions used in your two models.

Your first model uses Gamma(link = "inverse") (the default). Remember that the inverse link is a decreasing function of the (smooth) predictor—bigger values on the link scale imply smaller values on the natural scale of the response. In contrast, the identity link in your second model (i.e., gaussian(link = "identity")) is an increasing function of the predictor. The same increasing property would hold if you had used, for example, Gamma(link = "log").

Like your last question, I definitely recommend incorporating posterior expectation/prediction plots into your post-fitting routine. Dissecting out individual terms is an important part of model interrogation, but their net impact on the natural response scale gets very difficult to infer in isolation from a complex model (smooths, non-identity link functions, etc.).

Close, but not quite. Drop the square on the expected value in the denominator. 1/E[Y]=X\beta in typical regression notation.

1 Like