4
$\begingroup$

I am undertaking some data analysis in R where I fit a binomial GLM to some data using the glm function. The model is called within a "wrapper" function that manually alters the family and call for the model after it is fit to the data. Here is an example of the kind of command I am using:

CUSTOM.MODEL <- function(formula, options, data) {
  
  [Some commands that create NEWFORMULA from formula and options]
  
  MODEL <- glm(formula = NEWFORMULA, data = data, family = binomial (link = 'cloglog'))
  MODEL$family <- 'custom family'
  MODEL$call   <- sys.call()

MODEL }

Even though this is a binomial GLM, when I call the summary of a model fit with this function I get output where the dispersion parameter is not one --- i.e., the model has adjusted to account for overdispersion. I am not sure why this is happening. I would like to turn this feature off and fit the data to a model with dispersion parameter equal to one ---i.e., force the model to fit without accounting for overdispersion. How do I do this?

$\endgroup$
0
3
$\begingroup$

I don't think this is true. What you get is only a warning, the output is not affeced otherwise, that is, standard errors etc ... are calculated basedc on the binomial likelihood, not a quasibinomial likelihood.

To see for yourself, replace your call with

MODEL <- glm(cbind(Positive, Negative) ~ Var1 + Var2 + Var3, data 
             = DATA, family = quasibinomial(link = 'cloglog'))

and compare outputs.

$\endgroup$
2
  • 1
    $\begingroup$ Thanks for your answer (+1) --- I have revised the question to add some detail that was not previously supplied, and which I think is the source of the problem. Apologies that this information was not there in the first instance. $\endgroup$
    – Ben
    Oct 11 '20 at 5:57
  • 1
    $\begingroup$ I would like to accept an answer for this question, and since yours was accurate for the original question, yours would be my preference. Unfortunately your answer does not answer the updated question, so I do not wish to accept and mislead readers. If you can update this answer to make it accurate for the updated question I will accept it. $\endgroup$
    – Ben
    Oct 22 '20 at 1:17
2
$\begingroup$

The glm function using the binomial family fits a model with no overdispersion (i.e., the dispersion parameter is one). However, if you alter the family element of the resulting model then this can stuff up the summary output so that it gives you something unexpected. Specifically, if you change the family element to anything other than a binomial or Poisson distribution, the summary function will adjust for over-dispersion, even if the original model was fit without it.

For more detail on what is happening here, try fitting the model on its own and you will see that it does not adjust for overdispersion. Then call the function stat:::summary.glm and you will notice that it automatically adjusts for overdispersion unless the object$family$family %in% c("poisson", "binomial") (the object here is the model). So what is happening here is that the summary function is looking at your model and checking its family; then when it sees that it is not a binomial or Poisson family, it computes the summary using a method that adjusts for overdispersion. This means that you will see an unfortunate inconsistency between the outputs when you call MODEL versus when you call summary(MODEL).

One of the takeaway lessons here is that it can be dangerous to monkey with the elements of a model produced by the glm function. If you alter these elements manually, you can get some strange behaviours when you call functions that operate on these models. As to how to "turn this off" you have two options: the simplest thing to do is to revert the wrapper function so that it does not change the family element for the model; another option is to add a new class for your model and then program a custom summary function for that class of object (one that does not adjust for overdispersion).

$\endgroup$
4
  • $\begingroup$ You should ask about this on the R help email list. It mighy be a bug. $\endgroup$ Oct 11 '20 at 9:19
  • $\begingroup$ I'm not sure it would really be classified as a bug per se. The issue is that the stats package was written with a specific set of models in mind, and if it doesn't see the specific models it expects then it adds dispersion in the summary. It would be better if it were programmed in a more general way, but I'm not sure I'd be game to say that its lack of generality is a bug. $\endgroup$
    – Ben
    Oct 11 '20 at 9:28
  • $\begingroup$ Might be. But it merits discussion on that list! $\endgroup$ Oct 11 '20 at 9:32
  • $\begingroup$ Perhaps. I think the response would be that if you go around monkeying with the model family, you can't expect the later functions to work properly on the model. $\endgroup$
    – Ben
    Oct 11 '20 at 12:47
1
$\begingroup$

Another way to replicate this behaviour is by changing the family element before calling the glm function.

### some data for the example
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)

### a custom family 
### made by simply copying the Poisson family
### and changing the family element
###
### this family element is just the name
### the essential parts of the family object 
### are such elements like linkfun, linkinv and variance
f <- poisson()
f$family <- "something_else"

### the two different models
mod <- glm(counts ~ outcome + treatment, family = poisson())
mod2 <- glm(counts ~ outcome + treatment, family = f)

### Results:

### these use dispersion = 1
summary(mod)   ### not defined but selects the default for the family
summary(mod2, dispersion = 1)
### these estimate the dispersion
summary(mod2)

This behavior is not a bug but instead by design. It is described in the R documentation 'Summarizing Generalized Linear Model Fits':

The dispersion of a GLM is not used in the fitting process, but it is needed to find standard errors. If dispersion is not supplied or NULL, the dispersion is taken as 1 for the binomial and Poisson families, and otherwise estimated by the residual Chisquared statistic (calculated from cases with non-zero weights) divided by the residual degrees of freedom.

$\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.