## This script caculates the incidence rates for oz-related diseases according to GBD in China
rm(list=ls())
library(abind)
setwd('~/Documents/Research/3Paper2/4ERL/Code/Data_and_Analysis/Health/WHO')

##== Read disease-specific deaths
dth0 = read.table('Morticd9', header=TRUE, sep=',')
# China
sort(unique(dth0$Country))
dth1 = dth0[dth0$Country==3068,]
# Latest year
yr = max(unique(dth1$Year))
print(yr)
dth2 = dth1[dth1$Year==yr,]
# Subdiv
unique(dth2$SubDiv)
# A35: Selected Urban and Rural Areas
# A41: Selected Rural Areas
# A51: Selected Urban Areas
# deaths of all disease
dth.all = dth2[dth2$SubDiv=='A35',]

# Disease matching
dth.all[,'Cause']
dss = c('Circ', 'Resp')
nds = length(dss)
ids1 = is.element(dth.all[,'Cause'],'C041')|is.element(dth.all[,'Cause'],'C035')
ids2 = is.element(dth.all[,'Cause'],'C052')

# Deaths of 30+
# Index of Deaths12 and Deaths25
# Deaths12: Deaths at age 30-34 yearshea
# Deaths25: Deaths at age 95 years and above 
idt  = match(c('Deaths12','Deaths25'),colnames(dth.all))
dth.1 = sum(as.matrix(dth.all[ids1,idt[1]:idt[2]]),na.rm=TRUE)
dth.2 = sum(as.matrix(dth.all[ids2,idt[1]:idt[2]]),na.rm=TRUE)
dth = abind(dth.1,dth.2)
dth
# [1] 276555 122255

##== Population
pop0 = read.table('pop', header=TRUE, sep=',') 
# 3068: China country code
# A35: Selected Urban and Rural Areas
pop1 = pop0[pop0$Country==3068 & pop0$Year==yr & pop0$SubDiv=='A35',]
head(pop1)

# Pop from 30+
ipp = match(c('Pop12','Pop25'),colnames(pop1))
pop2 = pop1[,ipp[1]:ipp[2]]
pop = sum(as.matrix(pop2),na.rm=TRUE)
pop
# [1] 61155452

# Incidence rate
y0.ch = dth/pop
names(y0.ch) = dss
print(round(y0.ch*1000,3))
# circulatory respiratory 
      # 4.522       1.999
save(y0.ch, file="y0_oz_China.Rdata")