# Calculates incidence rate in China
rm(list=ls())
setwd('~/Documents/Research/3Paper2/4ERL/Code/Data_and_Analysis/Health/WHO')

##== Deaths
dth0 = read.table('Morticd9', header=TRUE, sep=',')
# China mortality in the latest, and A35 (Selected Urban and Rural Areas)
dth1 = dth0[dth0$Country==3068 & dth0$SubDiv=='A35',]
yr = max(unique(dth1$Year))
print(yr)
dth2 = dth1[dth1$Year==yr,]

# Disease matching
dss = c('IHD', 'CEV', 'COPD', 'LC', 'ALRI')
nds = length(dss)
ids1 = is.element(dth2[,'Cause'],'C046')|is.element(dth2[,'Cause'],'C047')
ids2 = is.element(dth2[,'Cause'],'C051')
ids3 = is.element(dth2[,'Cause'],'C054')
ids4 = is.element(dth2[,'Cause'],'C028')
ids5 = is.element(dth2[,'Cause'],'C053')

# Deaths by 5-yr age group from 30+ to 80+ 
#               30-34      75-79     90+
idt  = match(c('Deaths12','Deaths21','Deaths25'),colnames(dth2))

# IHD
# 30-79 by 5-yr 
dth3 = as.matrix(dth2[ids1,idt[1]:idt[2]])
# both-sex and both-disease deaths 
dth4 = colSums(dth3,na.rm=TRUE)
# 80+ both-sex and both-disease deaths 
dth5 = sum(as.matrix(dth2[ids1,(idt[2]+1):idt[3]]),na.rm=TRUE)
dth.1 = c(dth4,dth5) 

# CEV
# 30-79 by 5-yr 
dth3 = as.matrix(dth2[ids2,idt[1]:idt[2]])
# both-sex and both-disease deaths 
dth4 = colSums(dth3,na.rm=TRUE)
# 80+ both-sex and both-disease deaths 
dth5 = sum(as.matrix(dth2[ids2,(idt[2]+1):idt[3]]),na.rm=TRUE)
dth.2 = c(dth4,dth5)

# 30+ deaths for COPD and LC
dth.3 = sum(as.matrix(dth2[ids3,idt[1]:idt[3]]),na.rm=TRUE)
dth.4 = sum(as.matrix(dth2[ids4,idt[1]:idt[3]]),na.rm=TRUE)

# 5- deaths for ALRI
#              0 year    4 year
idt = match(c("Deaths2","Deaths6"),colnames(dth2))
dth.5 = sum(as.matrix(dth2[ids5,idt[1]:idt[2]]),na.rm=TRUE)

dth = list(dth.1,dth.2,dth.3,dth.4,dth.5)

##== Population
pop0 = read.table('pop', header=TRUE, sep=',') 
pop1 = pop0[pop0$Country==3068 & pop0$Year==yr & pop0$SubDiv=='A35',]

#               30-34   75-79   90+
ipp  = match(c('Pop12','Pop21','Pop25'),colnames(pop1))
# pop by age group 30-79 & 80+
pop2 = as.matrix(pop1[,ipp[1]:ipp[2]])
pop3 = colSums(pop2,na.rm=TRUE)
pop4 = sum(as.matrix(pop1[,(ipp[2]+1):ipp[3]]),na.rm=TRUE)
pop.1 = c(pop3,pop4) 
# pop 30+
pop.3 = sum(as.matrix(pop1[,ipp[1]:ipp[3]]),na.rm=TRUE)
# pop 5-
#              0 year 4 year
ipp = match(c("Pop2","Pop6"),colnames(pop1))
pop.5 = sum(as.matrix(pop1[,ipp[1]:ipp[2]]),na.rm=TRUE)

pop = list(pop.1,pop.1,pop.3,pop.3,pop.5)

##== Incidence rate
pop.total = array() 
dth.total = array()
y0.ch.total = array() 
y0.ch = list() #by age group

for (k in 1:nds){
    pop.total[k] = sum(pop[[k]])
    dth.total[k] = sum(dth[[k]])
    y0.ch.total[k] = dth.total[k]/pop.total[k]
    y0.ch[[k]] = dth[[k]]/pop[[k]]
}
dth.total
# [1]  62403 142851 102523  38131   3060
pop.total
# [1] 61155452 61155452 61155452 61155452  7306977

names(y0.ch.total) = dss
names(y0.ch) = dss
print(round(y0.ch.total*1000,3))
  # IHD   CEV  COPD    LC  ALRI 
# 1.020 2.336 1.676 0.624 0.419 

save(y0.ch,y0.ch.total,file="y0_pm_China.Rdata")