3
$\begingroup$

I`m using a GLM to determine the influence of climate variables in the incidence of a disease in 6 different cities across time (2007-2020).

I'm using a negative binomial regression, since the dependent variable is a count.

Dependent variable: monthly number of cases (casos)

Offset: monthly population in each city (populacao)

Year (ano) is a independent continuous variable (0, 1, 2 ...)

Each city (d1, d2 ...) is a independent dummy variable one-hot coded.

dataset

I want to determine the temporal trend of a disease in each city without adjust for climate variables and after adjusting with climate variables.

I'm removing the intercept because I don't want that one city becomes the reference level and I need coefficients for each city.

I have 2 questions:

  1. Why I'm getting NAs in the last coefficient interaction? Are there any errors in my model?

  2. The year (ano) exp(coef) is 1,07. Is the interpretation that there is a upward trend in the incidence of disease across cities correct?

summary(m1 <- glm.nb(casos ~ 0 + ano + ano*(d1 + d2 + d3 + d4 + 
              d5 + d6) + offset(log(populacao)), data = projeto10))

Call:
glm.nb(formula = casos ~ 0 + ano + ano * (d1 + d2 + d3 + d4 + 
    d5 + d6) + offset(log(populacao)), data = projeto10, 
    init.theta = 1.202601993, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6965  -1.0652  -0.7363   0.3766   4.0045  

Coefficients: (1 not defined because of singularities)
        Estimate Std. Error  z value Pr(>|z|)    
ano      0.06878    0.03153    2.181  0.02916 *  
d1     -10.22251    0.08561 -119.414  < 2e-16 ***
d2     -10.41664    0.08864 -117.510  < 2e-16 ***
d3     -11.44485    0.10843 -105.548  < 2e-16 ***
d4     -11.59067    0.12339  -93.938  < 2e-16 ***
d5     -11.60676    0.12727  -91.198  < 2e-16 ***
d6     -11.77710    0.12766  -92.250  < 2e-16 ***
ano:d1  -0.09080    0.03802   -2.389  0.01692 *  
ano:d2  -0.10842    0.03845   -2.820  0.00480 ** 
ano:d3  -0.06582    0.04144   -1.588  0.11223    
ano:d4  -0.14021    0.04388   -3.196  0.00139 ** 
ano:d5   0.01989    0.04447    0.447  0.65471    
ano:d6        NA         NA       NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.2026) family taken to be 1)

    Null deviance: 24450.50  on 1008  degrees of freedom
Residual deviance:   989.11  on  996  degrees of freedom
AIC: 2863

Number of Fisher Scoring iterations: 1


              Theta:  1.203 
          Std. Err.:  0.134 

 2 x log-likelihood:  -2837.036 
> exp(coef(m1))
         ano           d1           d2           d3           d4           d5           d6       ano:d1       ano:d2       ano:d3 
1.071203e+00 3.634290e-05 2.993030e-05 1.070447e-05 9.252024e-06 9.104293e-06 7.678393e-06 9.131966e-01 8.972502e-01 9.362954e-01 
      ano:d4       ano:d5       ano:d6 
8.691727e-01 1.020087e+00           NA 
$\endgroup$
1
  • $\begingroup$ For the future, please note that code fences must be at beginning of line (with nothing more on the same line) $\endgroup$ Dec 10, 2022 at 19:36

2 Answers 2

6
$\begingroup$

You are presumably trying to fit a slope and a trend for each city. You have six cities, so that's six slopes and six trends. But you have attempted to fit a linear model with 13 coefficients, so obviously one of the coefficients must be redundant.

The redundancy is not caused by removing the intercept but rather by trying to force R to estimate 6 interactions terms with your six one-hot coded dummy variables when there are actually only 5 degrees of freedom available for interaction.

Your problem would be easier if you created a factor indicating the city (lets call it City) instead of doing your own coding of dummy variables. The model you presumably want is:

m1 <- glm.nb(casos ~ 0 + City + City:ano + offset(log(populacao)))

which will estimate a separate intercept and slope for each city. An added advantage of this model is that the coefficients for the different cities are statistically independent, so forming confidence intervals for the slopes or for pairwise comparisons between them is easy.

Note the use of : instead of *, which tells R to include only the simple nested interaction without trying to expand out a main effect for ano, which is what has caused all your problems.

I can already tell from the model summary posted in your question what the trend will be for each city. The trend will be positive (upward) for cities 5 and 6, there's no trend for city 3, and the trend is negative (downward) for cities 1, 2 and 4. The slope for city 6 is 0.06878 (P = 0.029). The slope for city 5 is 0.06878 + 0.01989 = 0.08867, which will presumably also be statistically significant. The slope for city 4 is 0.06878 - 0.14021 = -0.07143, which will almost certainly be statistically significant in the opposite direction (you need to fit my model to get the exact p-value). The trends for cities 4 and 6 are significantly different (P = 0.00139), so pooling them would be questionable.

It appears that your cities fall into three groups. If I were you, I would examine the characteristics of the cities, especially their geographical location, to understand why the trends might appear different.

The easiest way to estimate an overall trend averaged across all cities is to fit the simpler model

m1 <- glm.nb(casos ~ 0 + City + ano + offset(log(populacao)))

This model allows a different baseline for each city but assumes a common trend. Almost certainly, the overall trend will not be significant because the conflicting results for the individual cities will cancel out. Concluding no trend may be scientifically questionable however given the significant trends for individual cities.

$\endgroup$
6
  • $\begingroup$ Thanks for you answer. According to your model, how would I get the main effect of ano across all cities? $\endgroup$ Dec 10, 2022 at 22:27
  • $\begingroup$ When you say "main effect", do you perhaps mean the overall regression slope for all cities combined? I can show you how to get that. As a contrast, it would just be a weighted average of the individual slopes. But I question the value of it. What utility is the overall slope when the individual cities so obviously have different slopes, which implies that the overall slope would be misleading for individual cities? Your question said that you needed coefficients for each city. $\endgroup$ Dec 10, 2022 at 22:44
  • $\begingroup$ In my model I need to know what is the temporal trend of the incidence of disease. So I need that the coefficient for ano represent if there is a up/downward trend. The average trend, considering that there is different trends across the cities. $\endgroup$ Dec 10, 2022 at 22:54
  • $\begingroup$ The trend isn't even in the same direction for each city. You have positive trends for cities 5 and 6, no trend for city 3, and negative trends for cities 1, 2 and 4. Hence your cities appear to fall into three groups. $\endgroup$ Dec 10, 2022 at 22:54
  • $\begingroup$ The slope parameter is the temporal trend. My model tells you whether there is a upward or downward trend for each city. As I noted in my previous comment, the trend is up for cities 5 and 6 and down for the cities 1, 2 and 4. If you try fitting my model, you will see that. Almost certainly the overall trend will not be statistically significant if you try to average these conflicting results. $\endgroup$ Dec 10, 2022 at 22:57
8
$\begingroup$

I'm removing the intercept because I don't want that one city becomes the reference level and I need coefficients for each city.

That is why you are getting an NA estimate for a coefficient. In some simple cases you can use that trick to get coefficients for all levels of a categorical predictor. A single categorical predictor with 6 levels can be represented either as an intercept plus 5 coefficients, or as 6 coefficients without an intercept. Either can work.

Once you add interactions you get into trouble like this. Your model and data only allow for 12 linearly independent coefficients. If you include an intercept, you get the intercept, a coefficient for age, 5 for cities, and 5 interaction coefficients. That's 12 total. When you try to remove the intercept you get a coefficient for age, 6 for cities, and you ask for 6 interaction coefficients. That would be 13, so one is omitted and listed as NA.

Your interpretation of the ano coefficient is not correct. It doesn't represent an overall association of outcome with time. Under the usual coding of predictors in R, it's the coefficient at some value of d, presumably for d6. The interaction coefficients are the differences from that ano coefficient for the other cities. d1, d2 and d4 all have significantly lower ano coefficients, numerically lower than 0.

There is no need to remove the intercept to accomplish what you need. It's better to fit the model with an intercept, then use post-modeling tools like those in the emmeans package to get predictions for whatever combinations of predictor values and comparisons among cities, etc., that you want. You can also use such tools to get an estimate of an overall trend in time across all cities, although that is tricky and potentially misleading when there are interactions.

$\endgroup$
3
  • $\begingroup$ Thanks for answering me. But if I maintain the intercept, I will be setting a reference level. Is it possible to not have a reference level? Won't this affect the results when doing the post-analysis with emmeans? $\endgroup$ Dec 10, 2022 at 21:45
  • 1
    $\begingroup$ @ÁlvaroSouza it isn't possible to omit a reference level when you have an interaction like this. I'd recommend letting R code the cities as a single 6-level categorical predictor instead of dummy coding the way that you have. That might make it easier for downstream software to understand the results. The emmeans functions will know how to work with the intercept and provide the estimates for all cities and any comparisons that you want. $\endgroup$
    – EdM
    Dec 10, 2022 at 21:53
  • 2
    $\begingroup$ Actually OP would still get a NA coefficient even if they remove the 0+ from the model. The redundancy isn't caused by the intercept (or not) but by the fact that they've tried to hardwire 6 terms for interaction. You're right of course that coding cities as a 6-level factor would have been better -- that would avoid the redundancy. $\endgroup$ Dec 11, 2022 at 3:40

Your Answer

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

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