Comparing incidence counts between countries

[Martin Mueller]

Dear all

I have a very simple dataset with three variables.
Each observation is the report of one of approx. 1000 ambulance accidents in one of three countries.
The variables are (accident_number, country and inhabitants_of_country).

So far, I collapsed the data, used the - cii means, poisson - command to get confidence intervals for the total number of accidents in a country, and divided it by the number of inhabitants of the country.
So I get an incidence rate, e.g. per 100'000 inhabitants with a 95% CI, that I can compare (see below).

gen accident_counter = 1
collapse (firstnm) inhabitants_of_country (sum)  accident_counter, by(country )

// Get the incidence + 95% CI
gen ci_poisson_lower = .
gen ci_poisson_upper = .
su country
forval i = 1/`r(N)' {
  su accident_counter if _n == `i'
  cii means 1 `r(mean)' , poisson    
  replace ci_poisson_upper = `r(ub)' in `i'
  replace ci_poisson_lower = `r(lb)' in `i'
}

replace ci_poisson_upper = round(ci_poisson_upper / inhabitants_of_country * 100000,.001)
replace ci_poisson_lower = round(ci_poisson_lower / inhabitants_of_country  * 100000, .001)
gen incidence_per100000 = round(accident_counter / inhabitants_of_country  * 100000,  .001)

I wonder if I can calculate a p-value for the obtained incidences and 95% CI.

My first approach was to use the - poisson - command, but the results of the p-value are not compatible with the 95% CI.

poisson incidence_per100000 i.country , base irr

Can you help me to spot my mistake?

Thank you in advance!
Martin

[Clyde Schechter]

I would be more inclined to do it this way:

gen accident_counter = 1
collapse (firstnm) inhabitants_of_country (sum) accident_counter, by(country )
poisson accident_counter i.country, exposure(inhabitants_of_country)
margins country

This will give you estimates of the accident rate per population in each country with a confidence interval. I don't know what kind of p-values you are looking for here. Do you wish to test a null hypothesis that these rates are zero--that seems like an absurd straw man. Perhaps you want to test whether they are equal across the three countries? You could do that with -testparm i.country, equal-.

[Martin Mueller]

This is exactly, what I was looking for, including the -testparm i.country, equal-.


One remaining question. I wonder, how do I interpret the margin command output, e.g. the - poisson accident_counter i.country, exposure(total_inh) base irr - command shows an IRR between country 1 and 2 of 0.97 (95% CI: 0.75-1.26), p= 0.819. Country 1 has 10 times for inhabitants.

The output of the margin command is:

Adjusted predictions                            Number of obs     =          3
Model VCE    : OIM

Expression   : Predicted number of events, predict()

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+------------------------------------------------------------------
    country    |
    country_1  |   240.7882   9.854807    24.43   0.000     221.4732    260.1033
    country_2  |   233.5551   29.66153     7.87   0.000     175.4196    291.6906
    country_3  |   97.96965   19.59393     5.00   0.000     59.56625     136.373
--------------------------------------------------------------------------------

Is there a way I can interpret the predicted number of events per inhabitants or is it more a constructed (adjusted) population size because otherwise the IRR would be 10 times higher in country 2 compared to country 1, no?

[Clyde Schechter]

The default output of -margins- after -poisson- is n, the number of incidents in each country. I forgot about that when responding in #2. To get the incidence rates in each country it should be:

margins country, predict(ir)

Sorry about that!