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).

Code:
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.

Code:
poisson incidence_per100000 i.country , base irr
Can you help me to spot my mistake?

Thank you in advance!
Martin