Starting with a toy example data (http://www.stata-press.com/data/r14/hbp.dta)

use http://www.stata-press.com/data/r14/hbp, clear

drop if mi(sex)
drop if mi(age_group)
sort city sex age_group year, stable
by city sex age_group: egen obs = total(hbp)
by city sex age_group:  gen pop = _N
by city sex age_group:  keep if _n == 1
drop year race hbp

replace city = 4 if city == 5 // used by the the loop later, so sequential better

bysort city (sex age): egen obs_city = total(obs)
bysort city (sex age): egen pop_city = total(pop)
gen crude = obs_city / pop_city

We have age & sex specific counts of events obs and population denominator pop in 5 different cities. Crude rate in a city is given by crude variable simply dividing observations by population within city and ignoring age structure. Now, I'd like to use information from all cities to create age & sex specific rates and using that - calculated expected number of events given the study population.

I'm trying to achieve that by using city specific standardized rates:

dstdize obs pop age sex, by(city) print saving( "temp", replace )

mat b = r(adj)
matrix list b

gen adjusted = .
gen exp_city_ds = .

forv city = 1/4 {
    local bw = b[1, `city']
    display `bw'
    replace adjusted = `bw' if city == `city'
    replace exp_city_ds = pop_city * `bw' if city == `city'

So far so good, but when I look at results I'd obtain from simple poisson model - the numbers do not add up:

poisson obs i.sex i.age, exp(pop) irr 
predict exp_p, n

bysort city (sex age): egen exp_city_poi = total(exp_p)
drop exp_p

gen dif0 = obs_city - exp_city_ds   
gen dif1 = obs_city - exp_city_poi  
gen dif2 = exp_city_ds - exp_city_poi
su dif? if sex == 1 & age == 2 // first obs by city

saveold "temp_old.dta", v(12) replace // data for R

I also tried to calculate from the same dummy dataset observed number of case using expected function of SpatialEpi package:


data <- read.dta("temp_old.dta") #, convert.factors = FALSE)

expected(data$pop, data$obs, 8)
seq(1, 4, 1)

city <- cbind(seq(1, 4, 1), expected(data$pop, data$obs, 8))
colnames(city) <- c("city", "E")

data <- merge(data, city, by="city")

data$poi_fit <- fitted(glm(obs ~ 1 + offset(log(pop)), data = data, family = "poisson"))

Again - I get different numbers, very close to Poisson fit (both Stata & R) but not the same.

Now, my question: is the reason of such difference lying in my calculations/methods? Or is it expected?

  • $\begingroup$ Both Stata & R code had some errors - they are now fixed. $\endgroup$ – radek Jul 27 '16 at 14:26

You are applying Poisson model to the dataset which does not behave this way -- the assumption of independence in time does not hold.

For the dataset in consideration. We look at the prevalence of hypertension for different cities and age/sex groups.

Try the following R code:

# As example read high blood pressure data
d <- read.dta("hbp.dta")
# Reformat data
d$city <- as.factor(d$city)
d$age_group <- factor(d$age_group) # Remove misleading unused age groups
d$hbp <- as.numeric(d$hbp) - 1
# Fit logistic regression
fm <- glm(hbp ~ city + age_group*sex, data=d)
# Basic model selection
bm <- step(fm)

You will find that the risk age group turns out to be (30-34) which are, indeed, the oldest subjects in the data. Being a woman is slight protective factor and the city3 has the highest prevalence. Therefore to test the Poisson behavior let us subset to, say, city1, men.

ind1 <- which(d$city == "1")
ind2 <- which(d$sex == "Male")
d2 <- d[c(ind1,ind2),]
ag <- aggregate(hbp ~ age_group, d2, sum)

The result is:

age_group hbp
15 - 19   5
20 - 24  12
25 - 29  19
30 - 34  24

Now fit Poisson:

x <- (1:4)
m1 <- glm(ag$hbp ~ x, family="poisson")

The model has $R^2 = 0.96$ so it is good.

  • 1
    $\begingroup$ Thanks. I'm aware that this is perhaps not the perfect dataset to analyse with this method - them main reason of using it was opportunistic availability in Stata. I'm more interested in general concepts and methods of approaching the problem that I could use with data I cannot share openly. $\endgroup$ – radek Jun 23 '16 at 7:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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