# Interpreting tweedie-location-scale result from R mgcv package

I came across the Tweedie-location-scale family function in R (twlss function from mgcv). It looks appealing to model my data. I understand the general concepts is that "twlss" family allows modelling scale and power parameter to the Tweedie distribution. However, after I tried it on data, I am puzzled by the returned values. For example, if I fit data with the command:

gam(list(y ~ s(x), ~1, ~1), family = twlss(), ...)

I expected the fitted power parameter and scale to be a constant "~1". My understanding is that the power must be between [1, 2] as required by Tweedie distribution. But sometimes, I get result saying fitted intercept.1 field to be a negative number. The documentation didn't mention whether intercept.1 is indeed the power parameter, but I guess they are related (for example, some kind of transformation, perhaps?)

I wonder if anyone ha used this package and knows how to interpret the fitted parameters. In particular, how to extract the fitted power and scale parameters from the fitted model.

Many thanks.

New contributor
xing_yu is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.

These location-scale-shape models are often tricky to work with in {mgcv} because the parameterisations chosen by Simon Wood, while numerically stable (as possible), don't often fit into the standard link/inverse link framework common to other GLMs. This coupled with the fact that the documentation for these general families often is terse or just lacking detail makes them difficult to work with unless you are an expert.

Doing a bit of digging (because I am not an expert), it seems that what is returned in the second column of fitted() or predict() for a model fitted using twlss() is not the power parameter itself but rather $$\theta$$ (which is on an unbounded scale) which needs to have an unstated inverse function applied to it to transform it to the bounded [a,b] scale.

I've written a few functions to help you convert from what is returned with the model fit to the parameters of the Tweedie distribution, $$\mu$$, $$\text{p}$$, and $$\phi$$, but for generality, it works on the fitted values of the model (or the predicted values if given some new data).

The first function just recovers the bounds on the power parameter $$\text{a}$$ and $$\text{b}$$ that were used when the model was fitted:

get_tw_ab <- function(model) {
fam <- family(model)
if (fam$family != "twlss") { stop("'model' wasn't fitted with 'twlss()' family.", call. = FALSE) } rfun <- fam$residuals
a <- get(".a", envir = environment(rfun))
b <- get(".b", envir = environment(rfun))
c(a, b)
}


Next we have a helper function that does the conversion from the $$\theta$$ scale to the $$\text{p}$$ (power) scale, given values for $$\text{a}$$ and $$\text{b}$$:

theta_2_power <- function(p, a, b) {
i <- p > 0
exp_theta_pos <- exp(-p[i])
exp_theta_neg <- exp(p[!i])
p[i] <- (b + a * exp_theta_pos) / (1 + exp_theta_pos)
p[!i] <- (b * exp_theta_neg + a) / (1 + exp_theta_neg)
p
}


These are used internally in the final function which, returns the parameters of the fitted Tweedie distribution for each observation or new data combination on their real scales

tw_params <- function(model, data = NULL) {
data <- if (is.null(data)) {
} else {
predict(model, newdata = data, type = "link")
}
ab <- get_tw_ab(model)
data[, 2] <- theta_2_power(data[, 2], a = ab[1], b = ab[2])
data[, 3] <- exp(data[, 3])
colnames(data) <- c("mu", "power", "scale")
data
}


Here's an example building on the one in ?twlss

set.seed(3)
n <- 400
## Simulate data...
dat <- gamSim(1, n = n, dist = "poisson", scale = .2)
dat$$y <- rTweedie(exp(dat$$f), p = 1.3, phi = 0.5) ## Tweedie response

## Fit a fixed p Tweedie, with wrong link ...
b <- gam(list(y ~ s(x0) + s(x1) + s(x2) + s(x3), ~1, ~1),
family=twlss(), data=dat)


What you get back for the intercepts for the power and scale parameters in the summary() output are the second and third columns of the matrix returned by fitted() or predict()

head(predict(b, type = "link"))

 head(predict(b, type = "link"))
[,1]       [,2]       [,3]
1 0.3770615 -0.4588997 -0.9378858
2 1.8091046 -0.4588997 -0.9378858
3 1.1973448 -0.4588997 -0.9378858
4 2.5295548 -0.4588997 -0.9378858
5 2.1003765 -0.4588997 -0.9378858
6 0.6554773 -0.4588997 -0.9378858


here the first column is on the link scale for the linear predictor of $$\mu$$, which by default is the log scale. If you asked for fitted(b) you'd get it on the natural non-log scale.

Notice that columns 2 and 3 of this matrix are constants in this case because we fitted a LSS model with constants in the linear predictors for the power and scale parameters.

The full suite of parameters on their natural/true scales are given y tw_params():

head(tw_params(b))

> head(tw_params(b))
mu    power     scale
1  1.457994 1.389502 0.3914546
2  6.104979 1.389502 0.3914546
3  3.311313 1.389502 0.3914546
4 12.547918 1.389502 0.3914546
5  8.169245 1.389502 0.3914546
6  1.926062 1.389502 0.3914546


These are pretty close to the values used to simulate the data (power = 1.3, and phi = 0.5)