Binomial model with GAM (mgcv) using weights

generalized-additive-modelmgcvr

I have been having an issue when trying to fit a binomial GAM to data. There are two ways these models can be coded, (i) providing a proportion as the response variable, and the number of trials as weights; and (ii) providing two columns, with successes and failures. I have reason to want to weight my data points (independently of the number of samples). However, I have noticed that if I use approach (ii) and add weights (using the weights argument), I get very odd results indeed. Furthermore, if I supply the same weights in relative terms (but different absolute magnitudes), I get very different output. This does not happen when using an equivalent GLM model (or, indeed, when using the gam package). How can I provide a set of weights for the data points?

Here is a MRE:

library('mgcv')

# Random data.
x = 1:100
y_binom = cbind(rpois(100, 5 + x/2), rpois(100, 100))
w = sample(seq_len(100), 100, replace = TRUE)

# GAM models.
m1 = gam(y_binom ~ s(x), family = 'binomial')
m2 = gam(y_binom ~ s(x), weights = w / mean(w), family = 'binomial')
m3 = gam(y_binom ~ s(x), weights = w / sum(w), family = 'binomial')
m4 = gam(y_binom ~ s(x), weights = w * 100, family = 'binomial')

ms = list(m1, m2, m3, m4)

# Different RMSEs.
lapply(X = ms, FUN = function(x) return(sqrt(mean(x$residuals^2))))

# Different predictions, e.g.
plot(predict(m2), predict(m3))


# This does not happen with GLMs.
m1 = glm(y_binom ~ x, family = 'binomial')
m2 = glm(y_binom ~ x, weights = w / mean(w), family = 'binomial')
m3 = glm(y_binom ~ x, weights = w / sum(w), family = 'binomial')
m4 = glm(y_binom ~ x, weights = w * 100, family = 'binomial')

ms = list(m1, m2, m3, m4)

# Same RMSEs (for m2-m4).
lapply(X = ms, FUN = function(x) return(sqrt(mean(x$residuals^2))))

# Same predictions, e.g.
plot(predict(m2), predict(m3))

Best Answer

I think you are seeing a difference because of an issue where smooths have difficulty and not any inherent problem in the GLM part of the model; your choice of weights is changing the magnitude of the log-likelihood which is resulting in slightly different models being returned.

I'll get back to that shortly. First, the "problem" goes away if you just fit a common or garden GLM with gam():

library('mgcv')

# Random data
set.seed(1)
x <- 1:100
y_binom <- cbind(rpois(100, 5 + x/2), rpois(100, 100))
w <- sample(seq_len(100), 100, replace = TRUE)

gam_m <- gam(y_binom ~ x, weights = w / mean(w), family = 'binomial')
glm_m <- glm(y_binom ~ x, weights = w / mean(w), family = 'binomial')

Exactly the same model is fitted

> logLik(gam_m)
'log Lik.' -295.6122 (df=2)
> logLik(glm_m)
'log Lik.' -295.6122 (df=2)
> coef(gam_m)
(Intercept)           x 
 -2.1698127   0.0174864 
> coef(glm_m)
(Intercept)           x 
 -2.1698127   0.0174864

and even if you change the magnitude of the log-likelihood by using a different normalization of the weights, you get the same fitted model even though the log+likelihood is different:

gam_other <- gam(y_binom ~ x, weights = w / sum(w), family = 'binomial')

> logLik(gam_other)
'log Lik.' -2.956122 (df=2)
> coef(gam_other)
(Intercept)           x 
 -2.1698127   0.0174864

The behaviour of glm() is that same in this regard:

> logLik(glm(y_binom ~ x, weights = w / sum(w), family = 'binomial'))
'log Lik.' -2.956122 (df=2)

# compare with logLik(gam_other)

This might break down in cases where the optimisation is more marginal, and this is what's happening with gam(). Using my gratia package we can easily compare the two GAMs fitted above:

# using your GAM m2 and m3 as examples
library(gratia)
comp <- compare_smooths(m2, m3)
draw(comp)

which produces

enter image description here

Note that by default, that smooths in those plots include a correction related to bias introduced when the smooth is estimated to be linear.

As you can see, the two fits are different; with one optimization penalising the smooth all the way back to a linear function and the other not quite penalizing as far. With more data, the extra complexity involved in fitting this model over a GLM (where in the GAM we're having to select smoothness parameters), would be overcome and I would expect the change to the log-likelihood to not have such a dramatic effect.

This situation is one where a some of the theory about GAMs starts to get a little looser there's work to try to correct or account for these issues, but often it can be difficult to tell the difference between something that is linear or slightly non-linear on the scale of the link function. Here the true function is slightly non-linear on the scale of the link function but m3 isn't able to identify this, in part I think because the weights are dominating the likelihood calculation.

Related Question