5
$\begingroup$

Regression Outcome

I had a question regarding negative binomial regression. My data uses COVID-19 deaths (y-intercept), and COVID-19 positive tests in California. I would like to try and use the regression formula to obtain a rate ratio for deaths in each county. I originally intended to use Poisson regression, although there was overdispersion. Both of my predictor variables are count data, However, I am unable to find examples on the internet of others plugging values into the formula to predict an outcome. I am only using two variables, each have 1,417 observation. I understand that the variables I have chosen may not be the best for prediction, although this is just intended to be practice.

How do you calculate the outcome using the formula I created? I would like to plug in the positive tests for each county in California to get an outcome IRR. For example, if San Diego county has 9,200 deaths, I would like to be able to write down my formula and get the IRR for San Diego county by plugging it into my formula. Using STATA I was able to obtain this formula: 38.84652 + 1.000044*(positive tests)

Improve this question
$\endgroup$

1 Answer 1

Reset to default
4
$\begingroup$

Incidence-rate ratios (IRRs) are exponentiated coefficients, so $\exp(b)$ rather than $b$. Standard errors and confidence intervals are similarly transformed.

To predict deaths, you first need to remove the irr option in Stata or take natural logs of the IRRs to recover the coefficients. Then the expected number of deaths is given by $$\begin{align} \mathbf{E}(\mathtt{deaths} \mid \mathtt{tests}) &= \exp (\mathtt{constant + tests} \cdot \mathtt{b_{tests}}) \\ &=\exp (\ln(\mathtt{IRR_{constant}) + tests} \cdot \ln(\mathtt{IRR_{tests}})) \\ &=\mathtt{IRR_{constant} }\cdot \exp(\mathtt{tests} \cdot \ln(\mathtt{IRR_{tests}})). \end{align} $$

Deaths are a rate since they are per length of time that your model uses for estimation. If you are interested in the multiplicative effect of one more test, that would be $\exp(\mathtt{b_{tests}})$. This is one type of IRR and already appears in the table you have. You could also be interested in the effect of $k$ more tests, which would be $\exp(k \cdot \mathtt{b_{tests}})$. This is another type of IRR. This is useful when the effect of a single test is small.

In addition to doing it by hand, you can also use margins, at(positive_tests = (put_number_here)) or nlcom to verify that you did the two steps correctly.

Finally, I would give poisson, vce(robust) a whirl to handle overdispersion. It requires fewer assumptions than the negative binomial, though with overdispersed cross-sectional data this may not matter too much. See more on this here.

Reproducible Stata example:

. sysuse auto, clear
(1978 automobile data)

. nbreg price c.mpg, irr nolog

Negative binomial regression                            Number of obs =     74
                                                        LR chi2(1)    =  24.22
Dispersion: mean                                        Prob > chi2   = 0.0000
Log likelihood = -668.24436                             Pseudo R2     = 0.0178

------------------------------------------------------------------------------
       price |        IRR   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         mpg |   .9651058   .0062012    -5.53   0.000     .9530279    .9773369
       _cons |   12830.24   1831.567    66.26   0.000     9698.898    16972.55
-------------+----------------------------------------------------------------
    /lnalpha |  -2.100561   .1613593                      -2.41682   -1.784303
-------------+----------------------------------------------------------------
       alpha |   .1223877   .0197484                      .0892048     .167914
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to incidence-rate ratios.
Note: _cons estimates baseline incidence rate.
LR test of alpha=0: chibar2(01) = 6.2e+04              Prob >= chibar2 = 0.000

. // (1) Predict prices for 22 mpg cars
. display exp(ln(12830.24)+22*ln(.9651058))
5873.313

. nbreg price c.mpg, nolog

Negative binomial regression                            Number of obs =     74
                                                        LR chi2(1)    =  24.22
Dispersion: mean                                        Prob > chi2   = 0.0000
Log likelihood = -668.24436                             Pseudo R2     = 0.0178

------------------------------------------------------------------------------
       price | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         mpg |  -.0355175   .0064254    -5.53   0.000    -.0481111   -.0229239
       _cons |    9.45956   .1427539    66.26   0.000     9.179768    9.739353
-------------+----------------------------------------------------------------
    /lnalpha |  -2.100561   .1613593                      -2.41682   -1.784303
-------------+----------------------------------------------------------------
       alpha |   .1223877   .0197484                      .0892048     .167914
------------------------------------------------------------------------------
LR test of alpha=0: chibar2(01) = 6.2e+04              Prob >= chibar2 = 0.000

. display exp(_b[_cons] + _b[mpg]*22)
5873.3192

. margins, at(mpg == 22)

Adjusted predictions                                        Number of obs = 74
Model VCE: OIM

Expression: Predicted number of events, predict()
At: mpg = 22

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   5873.319   240.4986    24.42   0.000     5401.951    6344.688
------------------------------------------------------------------------------

. nlcom expct_price:exp(_b[_cons] + _b[mpg]*22)

 expct_price: exp(_b[_cons] + _b[mpg]*22)

------------------------------------------------------------------------------
       price | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
 expct_price |   5873.319   240.4986    24.42   0.000     5401.951    6344.688
------------------------------------------------------------------------------

. // (2) Some IRRs
. display exp(_b[mpg]*1) // 1 more mpg
.96510585

. nlcom irr_1:exp(_b[mpg]*1) // now with more useful info

       irr_1: exp(_b[mpg]*1)

------------------------------------------------------------------------------
       price | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       irr_1 |   .9651058   .0062012   155.63   0.000     .9529517      .97726
------------------------------------------------------------------------------

. nlcom irr_10:exp(_b[mpg]*10) // now with more info

      irr_10: exp(_b[mpg]*10)

------------------------------------------------------------------------------
       price | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      irr_10 |   .7010508   .0450456    15.56   0.000     .6127631    .7893384
------------------------------------------------------------------------------
Improve this answer
$\endgroup$
12
  • $\begingroup$ I will use the margins command, thank you. I was hoping to get an outcome that express the IRR. Using my formula, would this be accurate: 9562.26 =38.84652 + 1.000044*(9523) ? This assumes that a county in California has 9523 positive tests. The interpretation is what is confusing to me. It looks like those in this county are 9562 times more likely to get covid? Is this a correct interpretation or am I missing a step? $\endgroup$
    – Karl Bailey Williams
    Commented Aug 11 at 19:03
  • $\begingroup$ You forgot to log the exponentiated IIR coefficients, so your number is too high. This quantity is the expected number of deaths given that number of positive tests. If you want deaths per additional positive test, that’s a multiplicative factor of 1.000044, so slightly more since it's greater than 1. You can also use margins, dydx(tests) at(tests=k) to get an additive effect on deaths from the (k+1)st test. It should also be a small number. $\endgroup$
    – dimitriy
    Commented Aug 11 at 21:42
  • $\begingroup$ You can learn more about IRR here: stats.stackexchange.com/a/481541/7071 $\endgroup$
    – dimitriy
    Commented Aug 11 at 21:44
  • $\begingroup$ I see what you are saying. So, by your instruction the formula would be ln(38.84652 + 1.000044*(9523)) with the outcome being 9.16558. $\endgroup$
    – Karl Bailey Williams
    Commented Aug 11 at 22:57
  • 1
    $\begingroup$ No, that's the expected number of deaths, not an IRR. Let's say there were zero positive tests. Then you would expect to see $\exp(\ln(38.84652))=38.84652$ baseline deaths. At 9623 tests, you would expect to see $\exp(\ln(38.84652))\cdot\exp(\ln(1.000044)*9523)=38.84652\cdot 1.5204446 = 59.063982$. The tests "scale up" the number of baselines deaths by 1.52 or 52%. This is because $\exp(a+b)=\exp(a) \cdot \exp(b)$. If there was one additional positive test, then you would expect to see $59.063982 \cdot 1.000044 = 59.066581$ deaths. $\endgroup$
    – dimitriy
    Commented Aug 12 at 3:18

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.