Estimate incidence rate (looping)

[Lydia Ewurabena Gamey]

Hello all,

I currently have the dataset below.
I want to calculate the incidence rate for each month (pneumonia_age_14days*100000/person_month)

To calculate the incidence rate with 95%CI, I am using the command -
gen rec = _n
global max_num = _N
forvalues num = 1/$max_num {
keep if rec ==`num'
local numerator = pneumonia_age_14days
local denominator = person_month
quietly:cii means `denominator' `numerator', poisson level(95)
gen cimean = r(mean)*100000
gen cise = r(se)*100000
gen cilb = r(lb)*100000
gen ciub = r(ub)*100000
gen diff_lower = cimean - cilb
gen diff_upper = ciub - cimean
}

However, this looping command is unable to loop across all the observations.
The result of the first observation is replicated for the rest of the data.
I believe there is an error with my code and will be glad to get help.

Thanks in advance
Lydia

* Example generated by -dataex-. For more info, type help dataex
clear
input int month long age_cat double person_month int(year pneumonia_age_14days) float rec
17532 2 131352.07236842107 2008  40   1
17532 1          138388.75 2008  69   2
17532 3  317287.0394736842 2008  39   3
17532 4  571625.4934210527 2008  42   4
17563 4  534700.9210526316 2008  32   5
17563 3  298460.6907894737 2008  42   6
17563 1  131521.4802631579 2008  72   7
17563 2 123363.45394736843 2008  38   8
17592 3 321542.00657894736 2008  51   9
17592 4         573293.125 2008  28  10
17592 1 142616.31578947368 2008  63  11
17592 2 133212.99342105264 2008  24  12
17623 4  557192.7631578947 2008  27  13
17623 1 140439.90131578947 2008  57  14
17623 3   313878.552631579 2008  27  15
17623 2 130083.71710526316 2008  27  16
17653 2 135601.90789473685 2008  22  17
17653 1 147662.86184210528 2008  50  18
17653 3  326915.2960526316 2008  35  19
17653 4  577475.1973684211 2008  14  20
17684 1  145085.1644736842 2008  43  21
17684 2 132118.55263157896 2008  24  22
17684 3         318501.875 2008  29  23
17684 4  559900.3289473684 2008  37  24
17714 2 127169.34210526316 2008  24  25
17714 3 308122.89473684214 2008  26  26
17714 1 132643.65131578947 2008  40  27
17714 4  570424.5065789474 2008  25  28
17745 1 134783.61842105264 2008  25  29
17745 3   311575.427631579 2008  20  30
17745 4   569593.552631579 2008  14  31
17745 2 127215.88815789475 2008  15  32
17776 3  303313.9802631579 2008  28  33
17776 4  551564.3421052631 2008  33  34
17776 2 123999.40789473685 2008  29  35
17776 1  132364.6052631579 2008  38  36
17806 3 315241.21710526315 2008  36  37
17806 1 138675.88815789475 2008  76  38
17806 2 128744.27631578948 2008  47  39
17806 4  569557.7631578947 2008  36  40
17837 2 126496.51315789475 2008 106  41
17837 3   306791.677631579 2008  64  42
17837 1         134601.875 2008 130  43
17837 4  551452.5986842106 2008  48  44
17867 3 319386.18421052635 2008  91  45
17867 2 133232.23684210528 2008 100  46
17867 4  571325.1644736843 2008  96  47
17867 1  139460.4605263158 2008 126  48
17898 3  322719.6710526316 2009  49  49
17898 1  140292.0394736842 2009  72  50
17898 4  574609.4407894737 2009  49  51
17898 2 136383.42105263157 2009  52  52
17929 3  293399.4407894737 2009  38  53
17929 2  125371.7105263158 2009  40  54
17929 1 126789.53947368421 2009  60  55
17929 4  520552.1052631579 2009  36  56
17957 3  326640.4605263158 2009  43  57
17957 2 141165.26315789475 2009  58  58
17957 4  576738.8486842106 2009  45  59
17957 1 140047.46710526317 2009  64  60
17988 3 318275.49342105264 2009  26  61
17988 1 135820.98684210528 2009  63  62
17988 4  559917.5328947369 2009  24  63
17988 2 139089.73684210528 2009  37  64
18018 1         140475.625 2009  44  65
18018 4  580253.7828947369 2009  28  66
18018 2  146467.6644736842 2009  35  67
18018 3 331224.57236842107 2009  25  68
18049 4  562847.8947368421 2009  49  69
18049 3  322630.2960526316 2009  42  70
18049 2  144104.2105263158 2009  31  71
18049 1 136232.23684210528 2009  42  72
18079 3  315151.5460526316 2009  31  73
18079 1 140711.34868421053 2009  44  74
18079 2  125820.6907894737 2009  23  75
18079 4  566282.9605263158 2009  21  76
18110 2 127644.57236842105 2009  14  77
18110 3  316422.0394736842 2009  15  78
18110 1  140222.3355263158 2009  35  79
18110 4  566871.4473684211 2009  14  80
18141 3  307825.9868421053 2009  46  81
18141 4  549018.8157894737 2009  55  82
18141 2 125837.13815789475 2009  37  83
18141 1  135167.1052631579 2009  47  84
18171 2  131499.6052631579 2009  61  85
18171 4  564211.3486842106 2009  79  86
18171 3  317774.6381578948 2009  71  87
18171 1 138931.21710526317 2009  84  88
18202 2  128681.6447368421 2009  82  89
18202 4  543685.6907894737 2009  54  90
18202 1 133923.45394736843 2009 127  91
18202 3 307627.59868421056 2009  81  92
18232 1 138648.45394736843 2009 175  93
18232 2  135146.3815789474 2009 132  94
18232 3  319515.8552631579 2009  96  95
18232 4  562873.5855263158 2009  50  96
18263 3 321961.90789473685 2010  54  97
18263 2  138073.5855263158 2010  43  98
18263 4  565481.0855263158 2010  31  99
18263 1 138984.24342105264 2010  98 100
end
format %tdNN-CCYY month
label values age_cat age

[Jared Greathouse]

I'm confused. Why are we looping over observations? Doesn't this do the solution you seek?

[Bjarte Aagnes]

First, scale the exposure, not the events (used for calculation of CI).

You have aggregated data by month × age_cat, so to "to calculate the incidence rate for each month" you must aggregate.

replace person_month = person_month * 10^-5
bys month : ci means pneumonia_age_14days, exposure(person_month)

To collect the results from the ci command:

Alternative 1: frame post

frame create results long(month) double(mean lb ub)

levelsof month 

foreach month in `r(levels)'  {

    ci means pneumonia_age_14days ///
        if month == `month' , exposure(person_month)
    
    frame post results ( `month' ) ( r(mean) ) ( r(lb)  ) (  r(ub) )
}

frame results : save results, replace

Alternative 2: collect

levelsof month 

foreach level in `r(levels)' {
    
    lab define month `level' "`: display %tdNN-CCYY  `level' '" , modify
}

lab val month month 
 
qbys month : ///
    collect : ///
        ci means pneumonia_age_14days, exposure(person_month)
               
collect style cell result[mean lb ub ] , nformat(%4.1f)

collect label levels result "mean" "r", modify
collect label levels result "lb" "lb", modify
collect label levels result "ub" "ub", modify

collect layout (month)(result[mean] result[lb] result[ub])

. frame results: list

     +-------------------------------------------+
     | month        mean          lb          ub |
     |-------------------------------------------|
  1. | 17532   16.398347   14.149445   18.903118 |
  2. | 17563   16.911041   14.555738    19.53892 |
  3. | 17592   14.179981   12.104895   16.508725 |
  4. | 17623   12.088351   10.155696    14.28175 |
  5. | 17653   10.188142   8.4538498   12.173529 |
     |-------------------------------------------|
  6. | 17684   11.509114   9.6363404   13.639601 |
  7. | 17714   10.102249   8.3404425    12.12622 |
  8. | 17745   6.4732365    5.082883   8.1265557 |
  9. | 17776    11.51864   9.6097251     13.6957 |
 10. | 17806   16.923864   14.631739   19.473208 |
     |-------------------------------------------|
 11. | 17837   31.089675   27.908658   34.533921 |
 12. | 17867   35.499275   32.157725   39.093719 |
 13. | 17898   18.909637   16.503817   21.567535 |
 14. | 17929   16.320975   13.985957   18.934363 |
 15. | 17957   17.727622   15.410884   20.294339 |
     |-------------------------------------------|
 16. | 17988    13.00837   11.009951   15.264643 |
 17. | 18018   11.014487   9.2157304   13.061777 |
 18. | 18049   14.067416   11.996792   16.392791 |
 19. | 18079   10.366156   8.5875005   12.404651 |
 20. | 18110   6.7757717   5.3559663   8.4564697 |
     |-------------------------------------------|
 21. | 18141   16.549641   14.250674   19.113899 |
 22. | 18171   25.598377   22.760293   28.692521 |
 23. | 18202   30.881975   27.704391   34.324096 |
 24. | 18232   39.180605   35.655182   42.960317 |
 25. | 18263   19.407457   16.959504   22.109498 |
     +-------------------------------------------+

. collect layout

Collection: default
      Rows: month
   Columns: result[mean] result[lb] result[ub]
   Table 1: 25 x 3

------------------------
        |    r   lb   ub
--------+---------------
01-2008 | 16.4 14.1 18.9
02-2008 | 16.9 14.6 19.5
03-2008 | 14.2 12.1 16.5
04-2008 | 12.1 10.2 14.3
05-2008 | 10.2  8.5 12.2
06-2008 | 11.5  9.6 13.6
07-2008 | 10.1  8.3 12.1
08-2008 |  6.5  5.1  8.1
09-2008 | 11.5  9.6 13.7
10-2008 | 16.9 14.6 19.5
11-2008 | 31.1 27.9 34.5
12-2008 | 35.5 32.2 39.1
01-2009 | 18.9 16.5 21.6
02-2009 | 16.3 14.0 18.9
03-2009 | 17.7 15.4 20.3
04-2009 | 13.0 11.0 15.3
05-2009 | 11.0  9.2 13.1
06-2009 | 14.1 12.0 16.4
07-2009 | 10.4  8.6 12.4
08-2009 |  6.8  5.4  8.5
09-2009 | 16.5 14.3 19.1
10-2009 | 25.6 22.8 28.7
11-2009 | 30.9 27.7 34.3
12-2009 | 39.2 35.7 43.0
01-2010 | 19.4 17.0 22.1
------------------------

[Lydia Ewurabena Gamey]

Thanks for the feedback.

Jared, the link doesn't provide the solution.
The same command needs to be executed for each month(over 500 times)- looping allows me to run the command once.

Bjarte, the code provided does not calculate the incidence rate with 95%CI.

[Jared Greathouse]

Lydia Ewurabena Gamey why didn't Bjarte Aagnes ' command do what you ask? I'm not at my computer at the moment, so post the exact results you got with your dataex after using the syntax and say why it didn't give the results you expected.

[Lydia Ewurabena Gamey]

Sure!

Executing my command-

local max = _N

forvalues num = 1 / `max' {
local numerator = pneumonia_age_14days
local denominator = person_month
quietly: cii means `denominator' `numerator', poisson level(95)
gen cimean = r(mean)*100000
gen cise = r(se)*100000
gen cilb = r(lb)*100000
gen ciub = r(ub)*100000
gen diff_lower = cimean - cilb
gen diff_upper = ciub - cimean if rec==`num'
}

gives the results below

* Example generated by -dataex-. For more info, type help dataex
clear
input int month long age_cat double person_month int(year pneumonia_age_14days) float(rec cimean cise cilb ciub diff_lower diff_upper)
17532 2 131352.07236842107 2008 40  1 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17532 1          138388.75 2008 69  2 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17532 3  317287.0394736842 2008 39  3 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17532 4  571625.4934210527 2008 42  4 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17563 4  534700.9210526316 2008 32  5 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17563 3  298460.6907894737 2008 42  6 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17563 1  131521.4802631579 2008 72  7 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17563 2 123363.45394736843 2008 38  8 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17592 3 321542.00657894736 2008 51  9 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17592 4         573293.125 2008 28 10 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17592 1 142616.31578947368 2008 63 11 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17592 2 133212.99342105264 2008 24 12 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17623 4  557192.7631578947 2008 27 13 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17623 1 140439.90131578947 2008 57 14 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17623 3   313878.552631579 2008 27 15 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17623 2 130083.71710526316 2008 27 16 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17653 2 135601.90789473685 2008 22 17 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17653 1 147662.86184210528 2008 50 18 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17653 3  326915.2960526316 2008 35 19 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
17653 4  577475.1973684211 2008 14 20 30.452507 4.814964 21.75572 41.46767 8.696789 11.015163
end
format %tdNN-CCYY month
label values age_cat age

The cimean estimated here is the same when one manually calculates the incidence rate using the formula (pneumonia_age_14days*100000/person_month) for the first cell.
But as you can see, the same cimean is generated for all months.

[Bjarte Aagnes]

Seems you want age-group-specific rates per month. Example code in #3 could be modified to give rates by month × age_cat. But, the following code is closer to your example code:

sort month age_cat

gen double r  = .
gen double lb = .
gen double ub = . 
    
local events pneumonia_age_14days
local exposure person_month  
local scale = 10^-5
    
qui forvalues n = 1/`=_N' {
       
    cii means `=`exposure'[`n']*`scale'' `=`events'[`n']' , ///
        poisson level(95)
    
    replace r = r(mean) if _n==`n'
    replace lb = r(lb)  if _n==`n'
    replace ub = r(ub)  if _n==`n'
}

list month age_cat pneumonia_age_14days person_month r lb ub ///
    in 1/8 , sepby(month)

     +------------------------------------------------------------------------------+
     |   month   age_cat   pneumo~s   person_~h           r          lb          ub |
     |------------------------------------------------------------------------------|
  1. | 01-2008         1         69   138388.75   49.859544   38.793701   63.100444 |
  2. | 01-2008         2         40   131352.07   30.452508   21.755718   41.467672 |
  3. | 01-2008         3         39   317287.04   12.291709   8.7406068   16.803171 |
  4. | 01-2008         4         42   571625.49    7.347468   5.2954086   9.9316422 |
     |------------------------------------------------------------------------------|
  5. | 02-2008         1         72   131521.48   54.743909   42.833738   68.940894 |
  6. | 02-2008         2         38   123363.45   30.803288    21.79824   42.279919 |
  7. | 02-2008         3         42   298460.69   14.072205   10.142008   19.021533 |
  8. | 02-2008         4         32   534700.92    5.984654   4.0934989   8.4485458 |
     +------------------------------------------------------------------------------+

test:

. cii means 138388.75 69 , poisson

                                                            Poisson exact    
    Variable |   Exposure        Mean    Std. err.       [95% conf. interval]
-------------+---------------------------------------------------------------
             |   138388.8    .0004986      .00006        .0003879     .000631

. di %3.1f r(mean) * 10^5  _skip(5)  %3.1f r(lb) * 10^5 _skip(5)  %3.1f r(ub) * 10^5
49.9     38.8     63.1

[Lydia Ewurabena Gamey]

Dear Bjarte,

This is great!! Exactly what I needed!!
Thanks for your patience, time taken to work through this with me and sharing knowledge.
Wish you the best!!