3
$\begingroup$

I'm using the Conway-Maxwell-Poisson (CMP) distribution to model the amount of nouns in a clause (data is under-dispersed). I've run the model using glmmTMB (family= "compois") but I'm having problems interpreting the coefficients. I understand they are in log scale (just as in any poisson regression), right? Would it be ok to transform them to IRR (=Incidence Rate Ratios)? tab_model() does this automatically for poisson but doesn't do it for CMP, what's the problem here? How should we approach the interpretation of the coefficients in this models?

This is one the models I've been working on:

glmmTMB(nouns ~ lg + ad + sp + lg:ad + sp:ad + (1|sub), family = "compois", data = simulated.data)

Simulated data

I've sampled (with replacement) 200 observations from the original dataset and printed them with dput(simulated.data).

    structure(list(sub = c("lt", "abet", "abet", "alma", "mat", "lisat", 
"nt", "at", "valet", "lt", "abet", "tit", "amt", "fact", "t", 
"lisat", "tit", "abet", "gael", "mat", "jt", "luct", "tit", "at", 
"at", "angt", "angt", "at", "gael", "ct", "mat", "gael", "mat", 
"at", "lisat", "mat", "mat", "angt", "valet", "valet", "gael", 
"tit", "jt", "fact", "valet", "ct", "mat", "amt", "lisat", "lisat", 
"helet", "helet", "ct", "alma", "angt", "alma", "amt", "vera", 
"amt", "angt", "tit", "alma", "nt", "vera", "luct", "t", "luct", 
"luct", "luct", "angt", "fact", "jt", "gael", "mat", "tit", "abet", 
"at", "at", "luct", "tit", "at", "amt", "angt", "angt", "tit", 
"mat", "tit", "at", "lisat", "lt", "tit", "at", "nt", "luct", 
"fact", "gael", "tit", "nt", "at", "at", "amt", "gael", "at", 
"franct", "at", "angt", "valet", "nt", "angt", "angt", "mat", 
"at", "jt", "jt", "angt", "lt", "gael", "at", "at", "amt", "mat", 
"ct", "mat", "angt", "ct", "lt", "t", "at", "t", "luct", "at", 
"ct", "lisat", "at", "angt", "amt", "mat", "fact", "nt", "angt", 
"lt", "t", "fact", "luct", "gael", "angt", "lt", "nt", "t", "helet", 
"jt", "fact", "lt", "tit", "mat", "angt", "franct", "angt", "at", 
"angt", "at", "valet", "gael", "valet", "at", "lt", "mat", "tit", 
"jt", "at", "jt", "valet", "tit", "franct", "abet", "gael", "at", 
"franct", "helet", "mat", "amt", "helet", "fact", "valet", "lt", 
"angt", "gael", "t", "at", "gael", "t", "tit", "at", "amt", "gael", 
"mat", "valet", "tit", "tit", "gael"), nouns = c(0L, 1L, 1L, 2L, 1L, 1L, 0L, 
    0L, 2L, 1L, 0L, 1L, 0L, 4L, 0L, 3L, 0L, 1L, 1L, 0L, 0L, 2L, 1L, 
    2L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 2L, 1L, 2L, 1L, 2L, 1L, 
    0L, 1L, 2L, 1L, 1L, 1L, 0L, 2L, 0L, 0L, 0L, 2L, 3L, 1L, 0L, 1L, 
    3L, 0L, 3L, 1L, 1L, 1L, 1L, 0L, 0L, 2L, 2L, 0L, 1L, 0L, 1L, 0L, 
    0L, 1L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 1L, 1L, 2L, 1L, 0L, 1L, 3L, 
    2L, 1L, 2L, 0L, 3L, 1L, 0L, 2L, 1L, 0L, 0L, 2L, 1L, 0L, 1L, 1L, 
    1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 2L, 1L, 0L, 
    2L, 1L, 1L, 3L, 3L, 1L, 1L, 2L, 0L, 0L, 1L, 2L, 0L, 3L, 0L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 2L, 0L, 4L, 0L, 
    0L, 2L, 0L, 2L, 2L, 0L, 1L, 0L, 1L, 1L, 2L, 0L, 1L, 1L, 0L, 0L, 
    0L, 1L, 2L, 1L, 0L, 2L, 0L, 1L, 1L, 0L, 2L, 0L, 2L, 0L, 1L, 0L, 
    0L, 1L, 2L, 2L, 0L, 0L, 1L, 0L, 1L, 0L, 2L, 2L, 0L, 1L, 1L, 1L, 
    1L), sp = c("C", "A", "A", "C", "A", "A", "A", "A", "C", "C", 
    "C", "A", "C", "C", "A", "A", "C", "A", "C", "C", "A", "C", "C", 
    "C", "C", "C", "C", "A", "A", "C", "C", "C", "C", "C", "C", "A", 
    "C", "C", "A", "C", "C", "C", "C", "C", "C", "A", "C", "C", "C", 
    "A", "A", "A", "C", "C", "C", "C", "A", "A", "A", "A", "A", "A", 
    "C", "C", "A", "C", "A", "A", "A", "C", "A", "A", "C", "A", "A", 
    "A", "A", "C", "C", "A", "C", "A", "A", "A", "A", "A", "A", "A", 
    "A", "C", "A", "C", "C", "C", "A", "A", "A", "A", "C", "A", "A", 
    "C", "A", "C", "A", "C", "A", "A", "A", "C", "A", "A", "C", "A", 
    "A", "C", "C", "C", "A", "A", "C", "A", "C", "C", "C", "C", "A", 
    "C", "A", "C", "A", "A", "A", "C", "A", "C", "C", "A", "C", "C", 
    "C", "A", "C", "C", "C", "C", "C", "C", "A", "C", "A", "C", "A", 
    "A", "A", "A", "A", "A", "A", "C", "C", "A", "A", "C", "C", "C", 
    "C", "C", "A", "A", "A", "A", "C", "A", "C", "A", "A", "C", "C", 
    "A", "C", "A", "A", "A", "C", "C", "A", "A", "C", "C", "C", "A", 
    "A", "A", "A", "C", "C", "A", "A", "C"), ad = c("O", "O", "C", 
    "O", "O", "C", "C", "O", "O", "C", "O", "C", "C", "O", "C", "O", 
    "O", "C", "C", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", 
    "O", "C", "C", "C", "O", "C", "O", "C", "C", "O", "C", "C", "O", 
    "C", "O", "O", "O", "O", "O", "O", "O", "C", "O", "C", "C", "C", 
    "O", "O", "O", "C", "O", "C", "O", "C", "O", "O", "O", "O", "C", 
    "C", "O", "O", "C", "O", "O", "O", "C", "O", "C", "C", "O", "C", 
    "C", "C", "C", "C", "C", "O", "O", "O", "C", "C", "O", "C", "C", 
    "O", "O", "C", "O", "C", "O", "O", "O", "O", "O", "C", "C", "C", 
    "O", "O", "O", "O", "C", "O", "C", "O", "C", "C", "O", "O", "O", 
    "C", "C", "C", "O", "O", "O", "C", "O", "C", "C", "C", "O", "C", 
    "O", "C", "C", "C", "C", "C", "O", "C", "O", "O", "C", "O", "O", 
    "C", "O", "O", "O", "C", "O", "O", "C", "O", "O", "O", "O", "O", 
    "C", "O", "C", "O", "C", "C", "C", "O", "O", "O", "C", "O", "C", 
    "O", "O", "O", "O", "C", "O", "C", "O", "O", "O", "O", "O", "O", 
    "O", "O", "C", "O", "O", "O", "O", "O", "O", "C", "O", "C", "C", 
    "C", "C"), lg = c("SP", "QO", "SP", "SP", "SP", "SP", "QO", "SP", 
    "SP", "SP", "SP", "SP", "CON", "SP", "SP", "SP", "SP", "QO", 
    "CON", "QO", "QO", "SP", "SP", "SP", "SP", "CON", "SP", "SP", 
    "SP", "SP", "SP", "SP", "QO", "SP", "SP", "QO", "SP", "QO", "CON", 
    "CON", "SP", "SP", "SP", "QO", "SP", "SP", "QO", "SP", "SP", 
    "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", 
    "SP", "SP", "SP", "CON", "QO", "SP", "SP", "SP", "SP", "QO", 
    "QO", "SP", "SP", "SP", "SP", "SP", "SP", "QO", "SP", "SP", "SP", 
    "SP", "QO", "SP", "CON", "QO", "CON", "SP", "CON", "SP", "QO", 
    "QO", "SP", "CON", "SP", "QO", "SP", "SP", "SP", "SP", "SP", 
    "SP", "SP", "QO", "SP", "QO", "QO", "SP", "SP", "CON", "QO", 
    "SP", "SP", "SP", "SP", "SP", "SP", "SP", "CON", "QO", "CON", 
    "QO", "SP", "QO", "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", 
    "QO", "SP", "SP", "SP", "CON", "QO", "SP", "SP", "SP", "SP", 
    "SP", "SP", "CON", "SP", "SP", "QO", "SP", "SP", "SP", "SP", 
    "SP", "SP", "SP", "QO", "CON", "SP", "SP", "SP", "SP", "QO", 
    "SP", "SP", "SP", "QO", "CON", "SP", "SP", "SP", "SP", "SP", 
    "SP", "SP", "SP", "SP", "SP", "QO", "SP", "SP", "SP", "QO", "SP", 
    "SP", "SP", "SP", "SP", "SP", "SP", "SP", "CON", "QO", "QO", 
    "QO", "SP", "CON", "SP", "CON", "QO", "SP")), row.names = c(NA, 
    -200L), class = "data.frame")

```
Improve this question
$\endgroup$
2
  • $\begingroup$ There! I can't share the original data because it's not mine and it's not public but I've built a similar df with 200 obs. $\endgroup$
    – Leandra
    Commented Feb 2, 2021 at 15:01
  • $\begingroup$ Added too, thanks! $\endgroup$
    – Leandra
    Commented Feb 2, 2021 at 15:08

1 Answer 1

Reset to default
2
$\begingroup$

First off, you are right that family="compois" will by default use a log link as is standard in Poisson regression. You can learn this from the output of ?family_glmmTMB:

compois(link = "log")

which tells you the default link is "log". Thus, you can do the standard exponential transformation of parameter estimates to obtain the multiplicative effect of a unit increment in the corresponding predictor on the mean of the response (which is what is modeled).

And yes, you can calculate IRRs between two predictor value vectors based on the fitted means. For instance, let's look at the first two rows of your simulated.data:

> simulated.data[1:2,]
  sub nouns sp ad lg
1  c*     0  C  O SP
2  l*     1  A  O QO

The fitted means for these two predictor value vectors would be (note the type="response"!):

> (responses <- predict(model,newdata=simulated.data[1:2,],type="response"))
[1] 1.041901 1.246722

so the IRR between the two would be

> responses[2]/responses[1]
[1] 1.196584

You can also obtain this ratio by writing out the model, multiplying predictor values by parameter estimates, exponentiating and then taking the ratio... but it's likely easier to just call predict.glmmTMB() as I did.

As to why tab_model() doesn't do this: that is just a case of not every add-on package playing nice with every other one, in which case you just have to do your own calculations.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for your explanation! So just by running exp(fixef(model)) I could get IRRs, right? Let me see if I'm interpreting your example correctly, then a clause of subject l* in conditions A (sp), O (ad), and QO (lg) has 1.196584 times the 'risk' of including a noun than a clause of subject c* in conditions C (sp), O (ad) and SP (lg). Is this right? $\endgroup$
    – Leandra
    Commented Feb 2, 2021 at 15:56
  • $\begingroup$ Taking exponentials of the parameter estimates amounts to comparing the situation where exactly one predictor is set to 1, and the others are 0. Especially in the presence of interactions, that may not be what you want - in such a situation, you have to check carefully which situations you want to compare in your ratio. ... $\endgroup$
    – Stephan Kolassa
    Commented Feb 2, 2021 at 16:04
  • 2
    $\begingroup$ ... Your interpretation is almost correct. Remember that we are not modeling risks, but counts. So it's not a relative risk, but the expected number of nouns that is increased by a factor of 1.196584. $\endgroup$
    – Stephan Kolassa
    Commented Feb 2, 2021 at 16:05
  • $\begingroup$ "Taking exponentials of the parameter estimates amounts to comparing the situation where exactly one predictor is set to 1, and the others are 0." I thought I'd be comparing one predictor against the reference level. So for instance I could compare the effect of lg=SP against lg=QO (in a previous model, before including the interactions). $\endgroup$
    – Leandra
    Commented Feb 2, 2021 at 16:11
  • 2
    $\begingroup$ Yes, that is correct, I was sloppy in not distinguishing between categorical predictors and their dummy coding. You still need to be careful about your interactions. For instance, since you have a lg:ad interaction, the exponential of the parameter estimates for the lg factor levels gives the IRR when ad is at the reference level only. For the IRR of a lg factor level at other settings of ad, you need to either collect the parameter estimates yourself, or use predict(). $\endgroup$
    – Stephan Kolassa
    Commented Feb 2, 2021 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.