Expected number of cases, standardizing with age & sex of the whole population

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 )

matrix list b

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:

library(foreign)
library(SpatialEpi)

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? • Both Stata & R code had some errors - they are now fixed. Jul 27 '16 at 14:26 1 Answer 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 library(foreign) 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) summary(bm)  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.

• 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. Jun 23 '16 at 7:52