7
$\begingroup$

In "A Practioner's guide to Generalized linear models" in paragraph 1.83 it is stated that:

"In the particular case of a Poisson multiplicative GLM it can be shown that modelling claim counts with an offset term equal to the log of the exposure produced identical results to modelling claim frequencies with prior weights set to be equal to the exposure of each observation."

I am not able to find any further references of this results, so i took upon some empirical testing in which i was not able to find proof that the statement is correct. Can anyone provide some insight in why this results may be right/wrong.

FYI, i used the following R code to test the hypothesis, in which i could not obtain similar results for the two mentioned cases:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y'=y, 'X'=X, 'offset' = offset)
formula = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula, family = "poisson", df, weights = offset) 
#Third model using poisson model on y/offset as reference
dfNew = df
dfNew$y = dfNew$y/offset
fit3 = glm(formula, family = "poisson", dfNew)

#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, fit3$coefficients, c(intercept,coefs))

The coefficient estimates resulting from running this code is given below:

 >  
    (Intercept)       X.1       X.2       X.3        X.4       X.5       X.6
[1,]    1.998277 0.2923091 0.4586666 0.1802960 0.11688860 0.7997154 0.4786655
[2,]    1.588620 0.2708272 0.4540180 0.1901753 0.07284985 0.7928951 0.5100480
[3,]    1.983903 0.2942196 0.4593369 0.1782187 0.11846876 0.8018315 0.4807802
[4,]    2.000000 0.2909240 0.4576965 0.1807591 0.11658183 0.8005451 0.4780123
              X.7       X.8       X.9      X.10
[1,]  0.005772078 0.9154808 0.9078758 0.3512824
[2,] -0.003705015 0.9117014 0.9063845 0.4155601
[3,]  0.007595660 0.9181014 0.9076908 0.3505173
[4,]  0.005881960 0.9150350 0.9084375 0.3511749
> 

and we can observe the coefficients are not identical.

$\endgroup$
  • $\begingroup$ You should't really include rm(list=ls() ) in R code you post here! That could surprise somebody running it, making them angry at you. I removed it. I also edited to include the results from running the code. If you had done this originally, maybe you had gotten som response faster, since few readers will run the code themself. $\endgroup$ – kjetil b halvorsen Aug 24 '17 at 11:39
  • 1
    $\begingroup$ @kjetilbhalvorsen Thanks for the feedback, will do in the future. $\endgroup$ – BDP1 Aug 24 '17 at 11:56
  • $\begingroup$ @BDP1 As it stands now your R code is not testing the claim you are referring to. I suggest you to add the weight term for the fit3 and add an addendum directly to the question. $\endgroup$ – aivanov 1 hour ago
5
$\begingroup$

(with your R code, you could replace "poisson" with "quasipoisson" to avoid all the warnings that get generated. Nothing else of import will change. See (*) below). Your reference use the term "multiplicative glm" which I think just means a glm with log link, since a log link can be thought of as a multiplicative model. Your own example shows that the claim is false, at least as we interpreted it (Since the estimated parameters are not equal). You could write the authors and ask them what they meant. Below I will argue why the claim is false.

Let $\lambda_i$ be the poisson parameter and $\omega_i$ the weights. Let $\eta_i$ be the linear predictor without the offset, and then $\eta_i+\log(\omega_i)$ be the linear predictor with the offset. The poisson probability function is $$ f(y_i) = e^{-\lambda_i} \lambda_i^{y_i}/y_i ! $$ Then the log likelihood function for the model with offset becomes $$ \ell = -\sum_i \omega_i e^{\eta_i} + \sum_i y_i \eta_i +\sum_i y_i\log \omega_i - \sum_i \log y_i! $$ while the log likelihood function for the model with weights becomes $$ \ell^w = -\sum_i \omega_i e^{\eta_i}+\sum_i y_i \omega_i \eta_i -\sum_i \omega_i \log y_i! $$ and this are clearly not the same. So what those authors meant is not clear to me.

(*) Note from the help of R's glm function:

Non-‘NULL’ ‘weights’ can be used to indicate that different observations have different dispersions (with the values in ‘weights’ being inversely proportional to the dispersions); or equivalently, when the elements of ‘weights’ are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations. For a binomial GLM prior weights are used to give the number of trials when the response is the proportion of successes: they would rarely be used for a Poisson GLM.

Looking into the meaning of the weights arguments explains this, it gives little meaning with the poisson family function, which assumes a constant scale parameter $\phi=1$ while the weights arguments modifies $\phi$. This do give more meaning with the quasiposson family function. See the answer to "weight" input in glm and lm functions in R The answer given there also helps in seeing why the likelihood in the weighted case takes the form given above.

The answer given here might be relevant: How is a Poisson rate regression equal to a Poisson regression with corresponding offset term? and is very interesting.

$\endgroup$
  • $\begingroup$ Thanks for the answer. Going trough the likelihood contributions clarifies a lot. Upon searching some more in response to your answer, i found the following question: stats.stackexchange.com/questions/66791/… In which it is shown that by dividing the original response by the offset, and setting the offset variable as weights, the same results can be obtained as the baseline model, where the log(offset) enters the linear predictor with a constant coefficient of 1. $\endgroup$ – BDP1 Aug 24 '17 at 11:58
  • $\begingroup$ I also attempted to derive that the likelihood contributions of this new model are equal to the contributions of the newly mentioned model, but was unable to do so, given the $(y_i/w_i)!$ term that appears for the latter. $\endgroup$ – BDP1 Aug 24 '17 at 12:03
  • $\begingroup$ At least in the experiment in the provided R-script, this claim appears to be true. this can be easilly seen by adding a weights=offset argument in the fit3=glm(...) line $\endgroup$ – BDP1 Aug 24 '17 at 12:07
  • $\begingroup$ I dont understand what you say here, if you think your experiment corroborates, why not edit that change to the expåeriment to the OP and explain why you think it corroborates the claim. I tried to explain why it is not true, what is wrong with my argument? $\endgroup$ – kjetil b halvorsen Aug 24 '17 at 16:26
  • $\begingroup$ Excuses me for the unclarities. I think your answer perfectly answered the initial question. However, given the linked answer, i was curious wheter we can prove that it that claim is valid, the same way you showed that my initial claim is invalid. Since my initial question is answered i Will accept your answer and look for answers on the new questions elsewhere. Thanks a lot for your effort! $\endgroup$ – BDP1 Aug 24 '17 at 16:34
1
$\begingroup$

Sorry not to simply add a comment above, but I do not have enough rep.

The original claim - but modified a little - is infact correct.

The following two models give exactly the same answer in R using a poisson glm with log-link:

  • model y, use an offset being log(offset)
  • model y / offset, use weights equal to offset

adjusting your original code slightly shows identical answers:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y' = y,
                'y_over_offset' = y/offset,
                'X' = X,
                'offset' = offset)

formula_offset = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))
formula_weights = paste("y_over_offset ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula_offset, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula_weights, family = "poisson", df, weights = offset) 


#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, c(intercept,coefs))

Hopefully that should give identical answers.

It is possible to show that the two models are statistically equivalent (there is a CAS paper somewhere which shows this - I will post a link if I have time).

Incidentally, if you are doing penalised regression then the way different packages such as glmnet and H2o measure deviance for the two different ways of defining a model, may lead to different results.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.