###################################################################################
# Carl 25 Feb 2016
# estimate the exposure over calendar years 2009-2011, per person who was integer
# age A on the 2010 census date (1 Aug 2010)
#
# results are 100x100 matrices Em and Ef for males and females, respectively
#
# For example,
#
# Em %*% [2010 male census pop ages 0..99] = [approx. 2009-2010 male exposure at 0..99]
#
###################################################################################
rm(list=ls())
library(dplyr)
graphics.off()
if (.Platform$OS.type == 'windows') windows(record=TRUE)
## construct survival functions l(x) for males and females, based on HMD standard
## the calculated multipliers won't be very sensitive to the standard chosen
## the key point is to account for non-survival at high ages: a 97.58-year old
## observed in the Census (t=2010.58) represents *more* than one expected person-year
# of exposure at age 96 over [2009.0, 2010.0), for example, because we only see the survivors
mort = read.csv('HMDstd.csv') %>%
mutate( sf = exp(-exp(f)), sm = exp(-exp(m)),
lf = c(1, head(cumprod(sf),-1)),
lm = c(1, head(cumprod(sm),-1))
)
lx.male = approxfun( x=mort$age, y=mort$lm, yleft=1, yright=0)
lx.female = approxfun( x=mort$age, y=mort$lf, yleft=1, yright=0)
## calculate the person-years of exposure at ages [A,A+h) x times [0,bigT]
## per exactly x-year-old at census date bigC
## if 2009.0 is t=0, 1 Aug 2010 Census was at t=1.58
PY = function(x, A , h=1, bigC=1.58, bigT=3, this.sex='m', dt=.05) {
tA = bigC - x + A # time of A-th birthday
tgrid = seq(dt/2, bigT-dt/2, dt) # grid of times at which exposure is possible
this.lx <- lx.male
if (this.sex == 'f') this.lx <- lx.female
# per x year old at time bigC, how many p-y of exposure in ages [A,A+n) x times [0,T)
sum( (tgrid >= tA) * (tgrid < (tA+h)) * this.lx(x-bigC+tgrid) * dt) / this.lx(x)
}
## build a (long) data frame with all combinations of exact census ages (in 1/5ths of a year)
## and one-year age groups. The expos variables will contain the expected person-years
## lived at ages [A,A+1) over years [0,3], per x-year-old at the Census date of 1.58
census.age = seq(-1.9, 99.9,.20) # possible ages at time C
age.group = 0:99
df = expand.grid( x=census.age, A=age.group)
for (i in 1:nrow(df)) {
df[i,'m.expos'] = PY(x=df$x[i], A=df$A[i], this.sex='m')
df[i,'f.expos'] = PY(x=df$x[i], A=df$A[i], this.sex='f')
}
## construct a multiplier matrix to convert Census populations by single years of age
## to period exposure. Aggregate exact census ages into integer age groups and then\
## calculate the average exposure. For example
## Em [100x100] %*% male.census.pop [100x1] = male period exposure [100x1]
tmp = df %>%
mutate(intx = floor(x)) %>%
group_by(intx,A) %>%
summarize(f.expos=mean(f.expos), m.expos=mean(m.expos))
Em = matrix(round(tmp$m.expos,2), nrow=100,
dimnames=list(paste0('A',0:99),paste0('x',-2:99)))
Ef = matrix(round(tmp$f.expos,2), nrow=100,
dimnames=list(paste0('A',0:99),paste0('x',-2:99) ))
## collapse the first 3 columns (floor(x) = -2, -1, 0) by summing. This is equivalent to assuming
## that the cohorts born over the two years after the Census will be identical in size to current 0-yr-olds
W = rbind( cbind(1, matrix(0, 2, 99)) ,
diag(100))
Em = Em %*% W
Ef = Ef %*% W
save(Em,Ef, file='EmEf.RData')
########### examples
big.df <- read.csv("../big2010.df.csv",header=T, skip=10, sep=",",
dec=".", stringsAsFactors = FALSE)
muni.df <- read.csv("../muni.df.csv", header=T, sep=",", dec=".", stringsAsFactors = FALSE)
BR.df <- left_join(big.df, muni.df, by='munid')
for (this.muniname in c('Fernando de Noronha','Batatais', 'Ribeirão Preto', 'São Paulo')) {
tmp = BR.df %>% filter(sex=='m', muniname==this.muniname)
N = tmp$pop[-101]
X = Em %*% N
plot( 0:99, N*3, type='h', lwd=2)
lines(0:99, X, lwd=3, col='seagreen')
title(paste(this.muniname,'Males'))
abline(v= seq(0,100,10), lty=2)
plot(0:99, log( head(tmp$death,-1)*3 / X), ylim=c(-9,0), pch='X', col='red')
points(0:99, log( head(tmp$death,-1) / N))
rug((0:99)[head(tmp$death, -1) == 0])
title(paste(this.muniname,'Males'))
}