##== This script calculates the baseline mortality in China from WHO
rm(list=ls())
setwd('~/Documents/Research/3Chn/NCC/Final/Li&Zhang_2018/Data_and_Scripts')

##== Incidence from WHO
# Deaths
dth0 = read.table('Health/GBD/y0_WHO/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))
dth2 = dth1[dth1$Year==yr,]

# Disease matching
dss = c('IHD', 'STROKE', '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+ 
idt  = match(c('Deaths12','Deaths21','Deaths25'),colnames(dth2))
dth3 = as.matrix(dth2[ids1,idt[1]:idt[2]])
dth4 = apply(dth3,2,sum,na.rm=TRUE)
dth5 = sum(as.matrix(dth2[ids1,(idt[2]+1):idt[3]]),na.rm=TRUE)
dth.1 = c(dth4,dth5) 

dth3 = as.matrix(dth2[ids2,idt[1]:idt[2]])
dth4 = apply(dth3,2,sum,na.rm=TRUE)
dth5 = sum(as.matrix(dth2[ids2,(idt[2]+1):idt[3]]),na.rm=TRUE)
dth.2 = c(dth4,dth5)

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)

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('Health/GBD/y0_WHO/pop', header=TRUE, sep=',') 
pop1 = pop0[pop0$Country==3068 & pop0$Year==yr & pop0$SubDiv=='A35',]

ipp  = match(c('Pop12','Pop21','Pop25'),colnames(pop1))
# pop by age group 30-75 & 80+
pop2 = as.matrix(pop1[,ipp[1]:ipp[2]])
pop3 = apply(pop2,2,sum,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-
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
y0 = array()   #total
y0.ch = list() #by age group

for(k in 1:nds){
   y0[k] = sum(dth[[k]])/sum(pop[[k]])
   y0.ch[[k]] = dth[[k]]/pop[[k]]
}

names(y0) = dss
names(y0.ch) = dss
print(round(y0*1000,3))
   # IHD STROKE   COPD     LC   ALRI 
 # 1.020  2.336  1.676  0.624  0.419
save(y0.ch,file = "Health/GBD/y0.Rdata")