Homework 2
Fikrewold Bitew
Nov 13, 2017
Used data from Ethiopan Demographic and Health Surveys (EDHS 2000 - 2016)
dhs_2016<-read.dta("C:/Users/atq765/Google Drive/Dem7283/DHS Resources/ET_2016_DHS_08152017_2025_103714/etir70dt/ETIR70FL.dta", convert.factors = F)
dhs_2011<-read.dta("C:/Users/atq765/Google Drive/Dem7283/DHS Resources/ET_2011_DHS_01242017_168_103714/etir61dt/ETIR61FL.dta", convert.factors = F)
dhs_2005<-read.dta("C:/Users/atq765/Google Drive/Dem7283/DHS Resources/ET_2005_DHS_02252017_1140_103714/etir51dt/ETIR51FL.dta", convert.factors = F)
dhs_2000<-read.dta("C:/Users/atq765/Google Drive/Dem7283/DHS Resources/ET_2000_DHS_02252017_1142_103714/etir41dt/ETIR41FL.dta", convert.factors = F)## identify the working variables and the all the edhs datasets
edhs_16<-subset(dhs_2016,select=c("caseid","v445","b0_01", "b1_01", "b2_01", "b3_01", "b4_01", "b5_01","b6_01","b7_01","b11_01","m1_1","m10_1", "m14_1", "m15_1","m18_1","mm2_01", "m70_1","h1_1","h2_1","h3_1", "h4_1", "h5_1", "h6_1", "h7_1", "h8_1","h9_1", "h10_1", "h33_1","v002", "v003", "v005", "v007", "v008","v011","v012", "v013", "v021", "v022", "v024", "v025","v102", "v106", "v212","v113", "v115","v116", "v119", "v120", "bidx_01","bidx_02","bidx_03", "bord_01",
"v121", "v122", "v123", "v124", "v125", "v127", "v129", "v153", "v190" , "v201","v211", "v212", "v213", "v217", "v302", "v364", "v302a", "v501","v511","v525", "v623", "v701","v705", "v717"))
## "hw70_1", "hw71_1", "hw72_1", "hw73_1" child height weight BMI found in 2011 and 2016 EDHS
edhs_11<-subset(dhs_2011,select=c("caseid","v445","b0_01", "b1_01", "b2_01", "b3_01", "b4_01", "b5_01","b6_01","b7_01","b11_01","m1_1", "m10_1", "m14_1", "m15_1","m18_1", "mm2_01", "m70_1","h1_1","h2_1","h3_1", "h4_1", "h5_1", "h6_1", "h7_1", "h8_1","h9_1", "h10_1", "h33_1","v002", "v003", "v005", "v007", "v008","v011","v012", "v013", "v021", "v022", "v024", "v025","v102", "v106", "v212","v113", "v115","v116", "v119", "v120", "bidx_01","bidx_02","bidx_03","bord_01",
"v121", "v122", "v123", "v124", "v125", "v127", "v129", "v153", "v190" , "v201","v211", "v212", "v213", "v217", "v302","v302a", "v364","v501","v511","v525", "v623", "v701","v705", "v717"))
edhs_05<-subset(dhs_2005,select=c("caseid","v445","b0_01", "b1_01", "b2_01", "b3_01", "b4_01", "b5_01","b6_01","b7_01", "b11_01","m1_1", "m10_1","m14_1", "m15_1","m18_1","mm2_01", "h1_1","h2_1","h3_1", "h4_1", "h5_1", "h6_1", "h7_1", "h8_1","h9_1", "h10_1", "h33_1", "v002", "v003", "v005", "v007", "v008","v011","v012", "v013", "v021", "v022", "v024", "v025","v102", "v106", "v212","v113", "v115","v116", "v119", "v120", "bidx_01","bidx_02","bidx_03","bord_01","v121", "v122", "v123", "v124", "v125", "v127", "v129", "v153", "v190" , "v201","v211", "v212", "v213", "v217", "v302", "v364","v501","v511","v525", "v623", "v701","v705", "v717"))
edhs_00<-subset(dhs_2000,select=c("caseid","v445","b0_01", "b1_01", "b2_01", "b3_01", "b4_01", "b5_01","b6_01","b7_01","b11_01","m1_1","m10_1","m14_1", "m15_1","m18_1","mm2_01","h1_1","h2_1","h3_1", "h4_1", "h5_1", "h6_1", "h7_1", "h8_1","h9_1", "h10_1", "h33_1","v002", "v003", "v005", "v007", "v008","v011","v012", "v013", "v021", "v022", "v024", "v025","v102", "v106", "v212","v113", "v115","v116", "v119", "v120", "bidx_01","bidx_02","bidx_03","bord_01","v121", "v122", "v123", "v124", "v125", "v127", "v129", "v153", "v201","v211", "v212", "v213", "v217", "v302", "v364", "v501","v511","v525", "v623", "v701","v705", "v717"))
###Add the variable "Year"
edhs_16$year <- edhs_16$v007
edhs_11$year <- edhs_11$v007
edhs_05$year <- edhs_05$v007
edhs_00$year <- edhs_00$v007
edhs_00$v190 <- NA
edhs_00$h33_1<- NA
edhs_00$v302a<-NA
edhs_05$v302a<-NA
edhs_00$m70_1<-NA
edhs_05$m70_1<-NA
edhs_16$yr <-'2016'
edhs_11$yr <-'2011'
edhs_05$yr <-'2005'
edhs_00$yr <-'2000'
edhs_all<-rbind(edhs_00, edhs_05, edhs_11, edhs_16)edhs_all$ever_use<-ifelse(edhs_all$yr==2016|edhs_all$yr==2011, recode(edhs_all$v302a, recodes='0="No"; 1:2="Yes"; else=NA', as.factor.result = F),recode(edhs_all$v302, recodes='0="No"; 1:3="Yes";else=NA', as.factor.result = F ))
edhs_all$region<- ifelse(edhs_all$yr==2016, recode(edhs_all$v024, recodes='1="Tigray"; 2="Afar"; 3= "Amhara"; 4= "Oromia"; 5="Somali"; 6="Ben-Gumuz"; 7="SNNP"; 8="Gambella"; 9="Harari"; 10="Addis Ababa"; 11="Dire Dawa";else=NA', as.factor.result = F),recode(edhs_all$v024, recodes='1="Tigray"; 2="Afar"; 3= "Amhara"; 4= "Oromia"; 5="Somali"; 6="Ben-Gumuz"; 7="SNNP"; 12="Gambella"; 13="Harari"; 14="Addis Ababa"; 15="Dire Dawa";else=NA', as.factor.result = F ))
# Marital Status
edhs_all$marst_new<-recode(edhs_all$v501, recodes='0="Notin Union"; 1:2="In Union"; 3:5="Notin Union"; else=NA', as.factor.result = T)
edhs_all$marst<-recode(edhs_all$v501, recodes='0="1nevermarried"; 1:2="0in_Union"; 3:5="3widowed/divorced/separated"; else=NA', as.factor.result = T)
edhs_all$child_sex<-recode(edhs_all$b4_01, recodes='1="Male"; 2="Female"; else=NA', as.factor.result = T)
edhs_all$sibling_alive<-recode(edhs_all$mm2_01, recodes='0="Dead"; 1="Alive"; else=NA', as.factor.result = T)
edhs_all$TT_pregnancy<-recode(edhs_all$m1_1, recodes='0="0Never"; 1="1TT"; 2:8="2+TT"; else=NA', as.factor.result = T)
edhs_all$prev_b_interval<-recode(edhs_all$b11_01, recodes='0:23="<2 years"; 24:47="2-3 years"; 48:252="4+ years"; else=NA', as.factor.result = T)
edhs_all$b_order<-recode(edhs_all$bord_01, recodes='1="1"; 2:3="2-3"; 4:8="4+"; else=NA', as.factor.result = T)
edhs_all$rank_order<- ifelse(edhs_all$bord_01== 1,"1st birth rank",ifelse((edhs_all$bord_01==2|edhs_all$bord_01==3) & edhs_all$b11_01>24 ,"2nd/3rd birth rank, >2 years",ifelse((edhs_all$bord_01==2|edhs_all$bord_01==3) & edhs_all$b11_01<=24,"2nd/3rd birth rank, <=2 years",
ifelse((edhs_all$bord_01==2|edhs_all$bord_01==3) & edhs_all$b11_01<=24,"2nd/3rd birth rank, <=2 years",
ifelse(edhs_all$bord_01==4 & edhs_all$b11_01>24 ,"4th birth rank, >2 years",ifelse(edhs_all$bord_01==4 & edhs_all$b11_01<=24 ,"4th birth rank, <=2 years",
"5+ birth order"))))))
edhs_all$bmi<-ifelse(edhs_all$v445==9999 | edhs_all$v445==9998 , NA, edhs_all$v445/100)
edhs_all$binbmi<-recode(edhs_all$bmi, recodes='0:18.4999="<18.5"; 18.5:60="18.5+"; else=NA', as.factor.result = F)
#Education
edhs_all$educ<-recode(edhs_all$v106, recodes="0='0Noeducation'; 1='1Primary'; 2:3='2Secodary and Higher'; else=NA", as.factor.result=F)
#Vaccination
edhs_all$ever_vaccin<-recode(edhs_all$h10_1, recodes="0='No'; 1='Yes'; else=NA", as.factor.result=T)
edhs_all$bcg_vaccin<-recode(edhs_all$h2_1, recodes="0='No'; 1:3='Yes'; else=NA", as.factor.result=T)
edhs_all$vitA_vaccin<-recode(edhs_all$h33_1, recodes="0='No'; 1:3='Yes'; else=NA", as.factor.result=T)
#Antenatal Care Visit
edhs_all$post_natal<-recode(edhs_all$m70_1, recodes="0='1No'; 1='0Yes'; else=NA", as.factor.result=T)
edhs_all$Post_Natal_Care<-recode(edhs_all$m70_1, recodes="0='No'; 1='Yes'; else=NA", as.factor.result=T)
edhs_all$antenatal_visits<-recode(edhs_all$m14_1, recodes="0='0_No Visit'; 1:4='1-4 visits'; 5:98='5+ Visits'; else=NA", as.factor.result=T)
edhs_all$ever_antenatal<-recode(edhs_all$m14_1, recodes="0='0No'; 1:98='Yes'; else=NA", as.factor.result=T)
edhs_all$contraceptive_use<-recode(edhs_all$v364, recodes="1='Yes'; 2:5='No'; else=NA", as.factor.result=T)
edhs_all$home_delivery<-recode(edhs_all$m15_1, recodes="11:12='1Home'; 13:96='Health Sector'; else=NA", as.factor.result=T)
edhs_all$fac_CS_delivery<-recode(edhs_all$m15_1, recodes="11:12='Home'; 22:25='Fac with CS delivery';
32:35='Fac with CS delivery'; 21='Fac without CS delivery';
31='Fac without CS delivery';26='Fac without CS delivery';
36='Fac without CS delivery';41:46='Fac without CS delivery';
96='Fac without CS delivery';else=NA",
as.factor.result=T)
##edhs_all$contraceptive_use<-recode(edhs_all$v302, recodes="0:2='1No'; 3='2Yes'; else=NA", as.factor.result=T)
edhs_all$child_wanted<-recode(edhs_all$m10_1, recodes="1='1Then'; 2='2Later'; 3='3Not at all'; else=NA", as.factor.result=T)
#Partner's/Husband's education level
edhs_all$partner_educ<-recode(edhs_all$v701, recodes="0='0Noeducation'; 1='1Primary'; 2:3='Secodary and Higher'; else=NA", as.factor.result=T)
#Occupation
edhs_all$occup<-recode(edhs_all$v717, recodes="0='1Unemployed'; 1:3='0Non-manual/Professional'; 4:9='2Agricultural/Manual(skilled/unskilled)'; else=NA", as.factor.result=T)
edhs_all$occup_new<-recode(edhs_all$v717, recodes="0='1Unemployed'; 1:9='Employed'; else=NA", as.factor.result=T)
#Partner's/Husband's occupation
edhs_all$partner_occup<-recode(edhs_all$v705, recodes="0 ='1Unemployed'; 1:3='2Non-manual/Professional'; 4:9='3Agricultural/Manual(skilled/unskilled)'; else=NA", as.factor.result=T)
#Age
edhs_all$mother_age<-recode(edhs_all$v013, recodes="1='<20'; 2:3='20-29'; 4:5='30-39';6:7='40-49';else='NA'", as.factor.result=T)
#Age
edhs_all$agefb<-recode(edhs_all$v212, recodes="0:19='<20'; 20:49='20+'; else='NA'", as.factor.result=T)
#Children Ever Born
edhs_all$parity<-recode(edhs_all$v201, recodes="1:2='1-2'; 3:4='3-4'; 5:18='5+';else=NA", as.factor.result=T)
#wealth Index
edhs_all$wealthindex<-recode(edhs_all$v190, recodes="1:2='1Poor'; 3='2Middle'; 4:5='3Rich'; else=NA", as.factor.result=T)
#Residence
edhs_all$residence<-recode(edhs_all$v025, recodes='1="Urban"; 2="Rural"; else=NA', as.factor.result = T)
edhs_all$urban<-recode(edhs_all$v102, recodes='1=1; else=0')
edhs_all$rural<-recode(edhs_all$v102, recodes='2=1; else=0')
edhs_all$post_natal_yes<-recode(edhs_all$m70_1, recodes='1=1; else=0')
edhs_all$kidid=paste(edhs_all$caseid, edhs_all$bidx, sep="-")
edhs_all$rural=ifelse(edhs_all$v025==2,1,0)
edhs_all$poor=ifelse(edhs_all$v190<=2,1,0)
edhs_all$Sec_educ=ifelse(edhs_all$v106>1,1,0)
edhs_all$prim_above=ifelse(edhs_all$v106>0,1,0)
weight=edhs_all$v005/1000000
psu=edhs_all$v021
strata=edhs_all$v022
#Yearedhs_all<-edhs_all[edhs_all$v201!=0,]
edhs_all$death.age<-ifelse(edhs_all$b5_01==1,((((edhs_all$v008))+1900)-(((edhs_all$b3_01))+1900)),edhs_all$b7_01)
#censoring indicator for death by age 5, in months (60 months)
edhs_all$d.event<-ifelse(is.na(edhs_all$b7_01)==T|edhs_all$b7_01>60,0,1)
edhs_all$d.eventfac<-factor(edhs_all$d.event); levels(edhs_all$d.eventfac)<-c("Alive at 5", "Dead by 5")
table(edhs_all$d.eventfac)##
## Alive at 5 Dead by 5
## 37704 2948table(edhs_all$d.eventfac)##
## Alive at 5 Dead by 5
## 37704 2948Kaplan-Meier function of under-five mortality of child by Survey Year EDHS (2000 - 2016)
edhs_all<-edhs_all[edhs_all$v201!=0,]
fit8<-survfit(Surv(death.age, d.event)~yr, data=edhs_all)
ggsurvplot(fit8,xlab = "Time in Months", ylim=c(.85, 1), conf.int = F, risk.table = F, title = "Survivorship Function for Under-5 Mortality by Survey Year
EDHS 2000- 2016")
#two group compairison
survdiff(Surv(death.age, d.event)~yr, data=edhs_all)## Call:
## survdiff(formula = Surv(death.age, d.event) ~ yr, data = edhs_all)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## yr=2000 10143 1060 727 152.620 205.566
## yr=2005 9339 692 674 0.466 0.613
## yr=2011 10896 663 793 21.266 29.514
## yr=2016 10274 533 754 64.756 88.280
##
## Chisq= 243 on 3 degrees of freedom, p= 0Kaplan-Meier function of under-five mortality of child by Place of Residence
edhs_all<-edhs_all[edhs_all$v201!=0,]
fit8<-survfit(Surv(death.age, d.event)~residence, data=edhs_all)
ggsurvplot(fit8,xlab = "Time in Months", ylim=c(.85, 1), conf.int = F, risk.table = F, title = "Survivorship Function for Under-5 Mortality by Survey Year
EDHS 2000- 2016")
PART A
Event Variable-ESTIMATING NEONATAL MORTALITY
edhs_all<-edhs_all[edhs_all$v201!=0&edhs_all$yr>2005,]
edhs_all$death.age.neo<-ifelse(edhs_all$b5_01==1,((((edhs_all$v008))+1900)-(((edhs_all$b3_01))+1900)),edhs_all$b7_01)
#censoring indicator for death by age 5, in months (60 months)
edhs_all$d.event.neo<-ifelse(is.na(edhs_all$b7_01)==T|edhs_all$b7_01>1,0,1)
edhs_all$d.eventfac.neo<-factor(edhs_all$d.event.neo); levels(edhs_all$d.eventfac.neo)<-c("Alive at 30 days", "Dead by 30 days")
table(edhs_all$d.eventfac.neo)##
## Alive at 30 days Dead by 30 days
## 20553 617table(edhs_all$d.eventfac.neo)##
## Alive at 30 days Dead by 30 days
## 20553 617PART B
Event Variable-ESTIMATING INFANT MORTALITY
edhs_all<-edhs_all[edhs_all$v201!=0&edhs_all$yr>2005,]
edhs_all$death.age.inf<-ifelse(edhs_all$b5_01==1,((((edhs_all$v008))+1900)-(((edhs_all$b3_01))+1900)),edhs_all$b7_01)
#censoring indicator for death by age 5, in months (60 months)
edhs_all$d.event.inf<-ifelse(is.na(edhs_all$b7_01)==T|edhs_all$b7_01>12,0,1)
edhs_all$d.eventfac.inf<-factor(edhs_all$d.event.inf); levels(edhs_all$d.eventfac.inf)<-c("Alive at 1 years", "Dead by 1 year")
table(edhs_all$d.eventfac.inf)##
## Alive at 1 years Dead by 1 year
## 20196 974table(edhs_all$d.eventfac.inf)##
## Alive at 1 years Dead by 1 year
## 20196 974PART C
Event Variable - Under 5 Mortality 1) Under 5 Mortality: Probability of dying between birth and exactly 5 years of age, expressed per 1,000 live births
If the child is dead at the time of interview, then the variable B5_01!=1, then the age at death in months is the variable B7_01.
If the child is alive at the time of interview, then the variable B5_01==1, and the age at death is censored. If the age at death is censored, then the age at the date of interview (censored age at death) is the date of the interview - date of birth (in months).
Since some health related indicators of under-five mortality like (postnatal care visits in the past two months) that are found to be valuable risk factors for under 5 mortality in other developing countries are incorporated only in the 2011 and 2016 EDHS datasets. In this analysis, we therefore employed only the 2011 and 2016 EDHS datasets.
edhs_all<-edhs_all[edhs_all$v201!=0&edhs_all$yr>2005,]
edhs_all$death.age<-ifelse(edhs_all$b5_01==1,((((edhs_all$v008))+1900)-(((edhs_all$b3_01))+1900)),edhs_all$b7_01)
#censoring indicator for death by age 5, in months (60 months)
edhs_all$d.event<-ifelse(is.na(edhs_all$b7_01)==T|edhs_all$b7_01>60,0,1)
edhs_all$d.eventfac<-factor(edhs_all$d.event); levels(edhs_all$d.eventfac)<-c("Alive at 5", "Dead by 5")
table(edhs_all$d.eventfac)##
## Alive at 5 Dead by 5
## 19974 1196table(edhs_all$d.eventfac)##
## Alive at 5 Dead by 5
## 19974 1196a) Estimating Survival Time Functions from Ethiopian EDHS 2011& 2016
#Here we see the data
head(edhs_all[,c("death.age","d.event")], n=20)## death.age d.event
## 16520 75 0
## 41002 67 0
## 51002 43 0
## 61002 42 0
## 71002 34 0
## 81002 123 0
## 91002 0 1
## 101002 97 0
## 121002 3 0
## 131002 110 0
## 141001 9 0
## 151001 42 0
## 16522 201 0
## 18102 0 1
## 19102 35 0
## 21102 51 0
## 24102 38 0
## 26102 0 0
## 27102 0 1
## 29102 48 0#The Surv() object
head(Surv(edhs_all$death.age, edhs_all$d.event), n=20)## [1] 75+ 67+ 43+ 42+ 34+ 123+ 0 97+ 3+ 110+ 9+ 42+ 201+ 0
## [15] 35+ 51+ 38+ 0+ 0 48+mort<-survfit(Surv(death.age, d.event)~1, data=edhs_all,conf.type="none")
plot(mort, ylim=c(.9,1), xlim=c(0,60), main="Survival Function for Under 5 Child Mortality", xlab="Child Age in Months")
summary(mort) ## Call: survfit(formula = Surv(death.age, d.event) ~ 1, data = edhs_all,
## conf.type = "none")
##
## time n.risk n.event survival std.err
## 0 21170 554 0.974 0.00110
## 1 20402 63 0.971 0.00116
## 2 19927 41 0.969 0.00120
## 3 19494 44 0.967 0.00124
## 4 19032 18 0.966 0.00126
## 5 18612 25 0.964 0.00128
## 6 18223 53 0.962 0.00133
## 7 17757 23 0.960 0.00136
## 8 17362 28 0.959 0.00139
## 9 16968 23 0.958 0.00141
## 10 16609 8 0.957 0.00142
## 11 16319 19 0.956 0.00144
## 12 16025 75 0.951 0.00152
## 13 15546 4 0.951 0.00153
## 14 15161 6 0.951 0.00153
## 15 14802 4 0.951 0.00154
## 16 14477 4 0.950 0.00154
## 17 14099 4 0.950 0.00155
## 18 13757 10 0.949 0.00156
## 19 13448 4 0.949 0.00157
## 20 13156 3 0.949 0.00158
## 21 12901 2 0.949 0.00158
## 22 12657 1 0.949 0.00158
## 23 12417 4 0.948 0.00159
## 24 12188 87 0.942 0.00173
## 33 9696 1 0.941 0.00174
## 36 9112 48 0.937 0.00187
## 48 7078 24 0.933 0.00197
## 54 6227 1 0.933 0.00198
## 60 5672 15 0.931 0.002071000*(1-summary(mort)$surv[36])## [1] NAThe Overall under-five mortality rate is 60.6 deaths per 1000 live births s
b) Estimating Hazard Functions from Ethiopian DHS 2016
library(muhaz)
haz<-kphaz.fit(time=edhs_all$death.age, status=edhs_all$d.event, method = "product-limit")
kphaz.plot(haz, main="Hazard function plot")
data.frame(haz)## time haz var
## 1 0.5 3.122226e-03 1.547402e-07
## 2 1.5 2.080708e-03 1.055971e-07
## 3 2.5 2.283937e-03 1.185597e-07
## 4 3.5 9.580398e-04 5.099315e-08
## 5 4.5 1.356998e-03 7.366020e-08
## 6 5.5 2.943701e-03 1.635058e-07
## 7 6.5 1.311146e-03 7.474594e-08
## 8 7.5 1.633702e-03 9.532519e-08
## 9 8.5 1.369585e-03 8.155728e-08
## 10 9.5 4.870003e-04 2.964645e-08
## 11 10.5 1.174372e-03 7.258848e-08
## 12 11.5 4.735599e-03 2.990339e-07
## 13 12.5 2.597993e-04 1.687528e-08
## 14 13.5 3.999526e-04 2.666134e-08
## 15 14.5 2.726624e-04 1.858662e-08
## 16 15.5 2.817236e-04 1.984267e-08
## 17 16.5 2.863696e-04 2.050231e-08
## 18 17.5 7.365089e-04 5.424664e-08
## 19 18.5 2.996825e-04 2.245331e-08
## 20 19.5 2.308575e-04 1.776526e-08
## 21 20.5 1.562134e-04 1.220131e-08
## 22 21.5 7.910454e-05 6.257528e-09
## 23 22.5 3.257326e-04 2.652594e-08
## 24 23.5 7.244146e-03 6.032363e-07
## 25 28.5 1.152545e-05 1.328360e-10
## 26 34.5 1.777095e-03 6.579678e-08
## 27 42.0 2.845552e-04 3.373913e-09
## 28 51.0 2.689690e-05 7.234435e-10
## 29 57.0 4.443896e-04 1.316580e-08This illustrates, that while the largest drop in survivorship occurred between 0 and 5 years, the hazard is actually higher in the 1-12 month and 25-30 months range, illustrating the conditionality of that probability.
c) Estimating Cummulative Hazard Functions from Ethiopian DHS 2016
#cumulative hazard
plot(cumsum(haz$haz)~haz$time, main = "Cumulative Hazard function", ylab="H(t)",xlab="Time in Months", type="l",xlim=c(0,60), lwd=2,col=3)
#Survival function, I just store this in an object so I can use it
surv<-mortd) Estimating the PDF from Ethiopian DHS 2016
#here is a cheap version of the pdf
ft<- -diff(mort$surv)
plot(ft, xlim=c(.5,59.5), type="s", ylab="f(t)",xlab="Time in Months",main="Probability Density Function")
#here is the cumulative distribution function
Ft<-cumsum(ft)
plot(Ft, xlim=c(0.5,60), type="s", ylab="F(t)",xlab="Time in Months", main="Cumulative Distribution Function")
A-1 Fit the Cox Model for NEONATAL MORTALITY
edhs_all<-edhs_all[edhs_all$v201!=0&edhs_all$yr>2005,]
options(survey.lonely.psu = "adjust")
edhs_all$weight<-edhs_all$v005/1000000
options(survey.lonely.psu = "adjust")
edhs_all$yrstratum<-paste(edhs_all$yr, edhs_all$v022, sep="-")
edhs_all$yrpsu<-paste(edhs_all$yr, edhs_all$v021, sep="-")
des1<-svydesign(ids=~yrpsu, strata=~yrstratum, weights=~weight, data = edhs_all[is.na(edhs_all$weight)==F,], nest=T)
cox.neo<-svycoxph(Surv(death.age.neo,d.event.neo)~yr+mother_age+educ+occup+marst+child_sex+child_wanted+TT_pregnancy+educ+contraceptive_use+antenatal_visits+post_natal+wealthindex+rank_order+child_wanted+fac_CS_delivery+sibling_alive+residence+v201, design=des1)B-1 Fit the Cox Model for INFANT MORTALITY
edhs_all<-edhs_all[edhs_all$v201!=0&edhs_all$yr>2005,]
options(survey.lonely.psu = "adjust")
edhs_all$weight<-edhs_all$v005/1000000
options(survey.lonely.psu = "adjust")
edhs_all$yrstratum<-paste(edhs_all$yr, edhs_all$v022, sep="-")
edhs_all$yrpsu<-paste(edhs_all$yr, edhs_all$v021, sep="-")
des2<-svydesign(ids=~yrpsu, strata=~yrstratum, weights=~weight, data = edhs_all[is.na(edhs_all$weight)==F,], nest=T)
cox.inf<-svycoxph(Surv(death.age.inf,d.event.inf)~yr+mother_age+educ+occup+marst+child_sex+child_wanted+TT_pregnancy+educ+contraceptive_use+antenatal_visits+post_natal+wealthindex+rank_order+child_wanted+fac_CS_delivery+sibling_alive+residence+region+v201, design=des2)C-1) Fit the Cox Model for CHILD MORTALITY
edhs_all<-edhs_all[edhs_all$v201!=0&edhs_all$yr!=2000,]
options(survey.lonely.psu = "adjust")
edhs_all$weight<-edhs_all$v005/1000000
options(survey.lonely.psu = "adjust")
edhs_all$yrstratum<-paste(edhs_all$yr, edhs_all$v022, sep="-")
edhs_all$yrpsu<-paste(edhs_all$yr, edhs_all$v021, sep="-")
des3<-svydesign(ids=~yrpsu, strata=~yrstratum, weights=~weight, data = edhs_all[is.na(edhs_all$weight)==F,], nest=T)
cox.child<-svycoxph(Surv(death.age,d.event)~yr+mother_age+educ+occup+marst+child_sex+child_wanted+TT_pregnancy+educ+contraceptive_use+antenatal_visits+post_natal+wealthindex+rank_order+child_wanted+fac_CS_delivery+sibling_alive+residence+region+v201, design=des3)edhs_all<-edhs_all[edhs_all$v201!=0,]
fit8<-survfit(Surv(death.age, d.event)~mother_age, data=edhs_all)
ggsurvplot(fit8,xlab = "Time in Months", xlim=c(0,60), ylim=c(.85, 1), conf.int = F, risk.table = F, title = "Survivorship Function for Under-5 Mortality by Post-natal Care visits within two months EDHS 2000- 2016")
#two group compairison
survdiff(Surv(death.age, d.event)~mother_age, data=edhs_all)## Call:
## survdiff(formula = Surv(death.age, d.event) ~ mother_age, data = edhs_all)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## mother_age=<20 780 48 35.1 4.71 4.96
## mother_age=20-29 8213 335 429.1 20.65 33.07
## mother_age=30-39 7658 382 440.4 7.76 12.45
## mother_age=40-49 4519 431 291.3 67.00 91.32
##
## Chisq= 103 on 3 degrees of freedom, p= 0edhs_all<-edhs_all[edhs_all$v201!=0,]
fit8<-survfit(Surv(death.age, d.event)~Post_Natal_Care, data=edhs_all)
ggsurvplot(fit8,xlab = "Time in Months", ylim=c(.9, 1), conf.int = F, risk.table = F, title = "Survivorship Function for Under-5 Mortality by Post-natal Care visits within two months EDHS 2000- 2016")
Summary results from COX Proportional Hazard Model to find out risk factors for Neonatal Mortality, Infant Mortality and Child Mortality saparately
summary(cox.neo)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (1239) clusters.
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = edhs_all[is.na(edhs_all$weight) == F, ], nest = T)## Call:
## svycoxph(formula = Surv(death.age.neo, d.event.neo) ~ yr + mother_age +
## educ + occup + marst + child_sex + child_wanted + TT_pregnancy +
## educ + contraceptive_use + antenatal_visits + post_natal +
## wealthindex + rank_order + child_wanted + fac_CS_delivery +
## sibling_alive + residence + v201, design = des1)
##
## n= 12991, number of events= 302
## (8179 observations deleted due to missingness)
##
## coef exp(coef) se(coef)
## yr2016 0.029500 1.029939 0.196886
## mother_age20-29 -0.833551 0.434504 0.356820
## mother_age30-39 -1.301045 0.272247 0.451493
## mother_age40-49 -1.313466 0.268887 0.549603
## educ1Primary 0.160050 1.173570 0.211557
## educ2Secodary and Higher 0.009963 1.010013 0.633023
## occup1Unemployed 0.121712 1.129429 0.271864
## occup2Agricultural/Manual(skilled/unskilled) 0.482625 1.620323 0.281479
## marst1nevermarried -0.604161 0.546533 0.785658
## marst3widowed/divorced/separated 0.694552 2.002811 0.269031
## child_sexMale 0.667848 1.950036 0.168384
## child_wanted2Later -0.349118 0.705310 0.249695
## child_wanted3Not at all -0.033522 0.967033 0.292269
## TT_pregnancy1TT -0.024392 0.975903 0.294098
## TT_pregnancy2+TT -0.532753 0.586987 0.233034
## contraceptive_useYes -0.249778 0.778974 0.225649
## antenatal_visits1-4 visits 0.014140 1.014240 0.204273
## antenatal_visits5+ Visits -0.179984 0.835284 0.376351
## post_natal1No 0.744886 2.106201 0.426789
## wealthindex2Middle -0.120453 0.886519 0.248245
## wealthindex3Rich 0.062292 1.064273 0.223287
## rank_order2nd/3rd birth rank, <=2 years 0.261529 1.298914 0.426918
## rank_order2nd/3rd birth rank, >2 years 0.336207 1.399628 0.328716
## rank_order4th birth rank, <=2 years -0.269790 0.763540 0.632838
## rank_order4th birth rank, >2 years -0.001822 0.998180 0.465026
## rank_order5+ birth order 0.068437 1.070833 0.533237
## fac_CS_deliveryFac without CS delivery 0.819637 2.269676 0.449290
## fac_CS_deliveryHome 0.027568 1.027952 0.347954
## sibling_aliveDead 0.268595 1.308125 0.187120
## residenceUrban -0.095466 0.908949 0.359997
## v201 0.145309 1.156397 0.074352
## z Pr(>|z|)
## yr2016 0.150 0.88090
## mother_age20-29 -2.336 0.01949 *
## mother_age30-39 -2.882 0.00396 **
## mother_age40-49 -2.390 0.01686 *
## educ1Primary 0.757 0.44933
## educ2Secodary and Higher 0.016 0.98744
## occup1Unemployed 0.448 0.65437
## occup2Agricultural/Manual(skilled/unskilled) 1.715 0.08642 .
## marst1nevermarried -0.769 0.44190
## marst3widowed/divorced/separated 2.582 0.00983 **
## child_sexMale 3.966 7.3e-05 ***
## child_wanted2Later -1.398 0.16206
## child_wanted3Not at all -0.115 0.90869
## TT_pregnancy1TT -0.083 0.93390
## TT_pregnancy2+TT -2.286 0.02224 *
## contraceptive_useYes -1.107 0.26832
## antenatal_visits1-4 visits 0.069 0.94481
## antenatal_visits5+ Visits -0.478 0.63248
## post_natal1No 1.745 0.08093 .
## wealthindex2Middle -0.485 0.62752
## wealthindex3Rich 0.279 0.78026
## rank_order2nd/3rd birth rank, <=2 years 0.613 0.54014
## rank_order2nd/3rd birth rank, >2 years 1.023 0.30641
## rank_order4th birth rank, <=2 years -0.426 0.66988
## rank_order4th birth rank, >2 years -0.004 0.99687
## rank_order5+ birth order 0.128 0.89788
## fac_CS_deliveryFac without CS delivery 1.824 0.06811 .
## fac_CS_deliveryHome 0.079 0.93685
## sibling_aliveDead 1.435 0.15117
## residenceUrban -0.265 0.79087
## v201 1.954 0.05066 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef)
## yr2016 1.0299 0.9709
## mother_age20-29 0.4345 2.3015
## mother_age30-39 0.2722 3.6731
## mother_age40-49 0.2689 3.7190
## educ1Primary 1.1736 0.8521
## educ2Secodary and Higher 1.0100 0.9901
## occup1Unemployed 1.1294 0.8854
## occup2Agricultural/Manual(skilled/unskilled) 1.6203 0.6172
## marst1nevermarried 0.5465 1.8297
## marst3widowed/divorced/separated 2.0028 0.4993
## child_sexMale 1.9500 0.5128
## child_wanted2Later 0.7053 1.4178
## child_wanted3Not at all 0.9670 1.0341
## TT_pregnancy1TT 0.9759 1.0247
## TT_pregnancy2+TT 0.5870 1.7036
## contraceptive_useYes 0.7790 1.2837
## antenatal_visits1-4 visits 1.0142 0.9860
## antenatal_visits5+ Visits 0.8353 1.1972
## post_natal1No 2.1062 0.4748
## wealthindex2Middle 0.8865 1.1280
## wealthindex3Rich 1.0643 0.9396
## rank_order2nd/3rd birth rank, <=2 years 1.2989 0.7699
## rank_order2nd/3rd birth rank, >2 years 1.3996 0.7145
## rank_order4th birth rank, <=2 years 0.7635 1.3097
## rank_order4th birth rank, >2 years 0.9982 1.0018
## rank_order5+ birth order 1.0708 0.9339
## fac_CS_deliveryFac without CS delivery 2.2697 0.4406
## fac_CS_deliveryHome 1.0280 0.9728
## sibling_aliveDead 1.3081 0.7645
## residenceUrban 0.9089 1.1002
## v201 1.1564 0.8648
## lower .95 upper .95
## yr2016 0.70020 1.5150
## mother_age20-29 0.21591 0.8744
## mother_age30-39 0.11237 0.6596
## mother_age40-49 0.09157 0.7896
## educ1Primary 0.77523 1.7766
## educ2Secodary and Higher 0.29208 3.4927
## occup1Unemployed 0.66290 1.9243
## occup2Agricultural/Manual(skilled/unskilled) 0.93327 2.8132
## marst1nevermarried 0.11718 2.5490
## marst3widowed/divorced/separated 1.18206 3.3934
## child_sexMale 1.40189 2.7125
## child_wanted2Later 0.43235 1.1506
## child_wanted3Not at all 0.54533 1.7148
## TT_pregnancy1TT 0.54837 1.7368
## TT_pregnancy2+TT 0.37177 0.9268
## contraceptive_useYes 0.50055 1.2123
## antenatal_visits1-4 visits 0.67962 1.5136
## antenatal_visits5+ Visits 0.39947 1.7466
## post_natal1No 0.91247 4.8616
## wealthindex2Middle 0.54498 1.4421
## wealthindex3Rich 0.68705 1.6486
## rank_order2nd/3rd birth rank, <=2 years 0.56258 2.9990
## rank_order2nd/3rd birth rank, >2 years 0.73487 2.6657
## rank_order4th birth rank, <=2 years 0.22088 2.6394
## rank_order4th birth rank, >2 years 0.40122 2.4834
## rank_order5+ birth order 0.37656 3.0452
## fac_CS_deliveryFac without CS delivery 0.94086 5.4752
## fac_CS_deliveryHome 0.51975 2.0331
## sibling_aliveDead 0.90651 1.8877
## residenceUrban 0.44886 1.8406
## v201 0.99958 1.3378
##
## Concordance= 0.689 (se = 0.016 )
## Rsquare= NA (max possible= NA )
## Likelihood ratio test= NA on 31 df, p=NA
## Wald test = 85.78 on 31 df, p=4.887e-07
## Score (logrank) test = NA on 31 df, p=NAsummary(cox.inf)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (1239) clusters.
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = edhs_all[is.na(edhs_all$weight) == F, ], nest = T)## Call:
## svycoxph(formula = Surv(death.age.inf, d.event.inf) ~ yr + mother_age +
## educ + occup + marst + child_sex + child_wanted + TT_pregnancy +
## educ + contraceptive_use + antenatal_visits + post_natal +
## wealthindex + rank_order + child_wanted + fac_CS_delivery +
## sibling_alive + residence + region + v201, design = des2)
##
## n= 12991, number of events= 465
## (8179 observations deleted due to missingness)
##
## coef exp(coef) se(coef)
## yr2016 0.136779 1.146575 0.162491
## mother_age20-29 -0.949289 0.387016 0.281488
## mother_age30-39 -1.392908 0.248352 0.355389
## mother_age40-49 -1.415394 0.242830 0.412551
## educ1Primary 0.111646 1.118117 0.154420
## educ2Secodary and Higher -0.307330 0.735408 0.565540
## occup1Unemployed -0.146494 0.863731 0.199852
## occup2Agricultural/Manual(skilled/unskilled) 0.120976 1.128598 0.228640
## marst1nevermarried -0.829463 0.436283 0.688954
## marst3widowed/divorced/separated 0.488135 1.629274 0.231073
## child_sexMale 0.562810 1.755599 0.134403
## child_wanted2Later -0.248097 0.780284 0.190657
## child_wanted3Not at all -0.050881 0.950392 0.238108
## TT_pregnancy1TT 0.132449 1.141621 0.241889
## TT_pregnancy2+TT -0.380023 0.683845 0.188010
## contraceptive_useYes -0.402439 0.668687 0.188374
## antenatal_visits1-4 visits 0.014300 1.014402 0.178872
## antenatal_visits5+ Visits -0.221893 0.801001 0.314970
## post_natal1No 0.782296 2.186487 0.347766
## wealthindex2Middle -0.180405 0.834932 0.206743
## wealthindex3Rich 0.106547 1.112430 0.166282
## rank_order2nd/3rd birth rank, <=2 years 0.266288 1.305111 0.330522
## rank_order2nd/3rd birth rank, >2 years 0.273746 1.314880 0.273632
## rank_order4th birth rank, <=2 years -0.212528 0.808538 0.529259
## rank_order4th birth rank, >2 years -0.106760 0.898741 0.384411
## rank_order5+ birth order 0.157154 1.170176 0.436060
## fac_CS_deliveryFac without CS delivery 0.909270 2.482510 0.400438
## fac_CS_deliveryHome 0.310828 1.364554 0.296250
## sibling_aliveDead 0.198718 1.219838 0.151909
## residenceUrban -0.014634 0.985473 0.291026
## regionAfar 0.419680 1.521474 0.504367
## regionAmhara 0.247630 1.280986 0.520659
## regionBen-Gumuz -0.102493 0.902584 0.514454
## regionDire Dawa 0.361389 1.435322 0.532780
## regionGambella 0.252920 1.287780 0.495955
## regionHarari 0.126354 1.134684 0.507226
## regionOromia 0.414595 1.513758 0.502422
## regionSNNP 0.384055 1.468226 0.500539
## regionSomali 0.206186 1.228982 0.494105
## regionTigray 0.006832 1.006855 0.498660
## v201 0.129508 1.138269 0.056868
## z Pr(>|z|)
## yr2016 0.842 0.399920
## mother_age20-29 -3.372 0.000745 ***
## mother_age30-39 -3.919 8.88e-05 ***
## mother_age40-49 -3.431 0.000602 ***
## educ1Primary 0.723 0.469678
## educ2Secodary and Higher -0.543 0.586835
## occup1Unemployed -0.733 0.463549
## occup2Agricultural/Manual(skilled/unskilled) 0.529 0.596727
## marst1nevermarried -1.204 0.228610
## marst3widowed/divorced/separated 2.112 0.034646 *
## child_sexMale 4.187 2.82e-05 ***
## child_wanted2Later -1.301 0.193165
## child_wanted3Not at all -0.214 0.830791
## TT_pregnancy1TT 0.548 0.583994
## TT_pregnancy2+TT -2.021 0.043250 *
## contraceptive_useYes -2.136 0.032648 *
## antenatal_visits1-4 visits 0.080 0.936283
## antenatal_visits5+ Visits -0.704 0.481128
## post_natal1No 2.249 0.024481 *
## wealthindex2Middle -0.873 0.382877
## wealthindex3Rich 0.641 0.521678
## rank_order2nd/3rd birth rank, <=2 years 0.806 0.420440
## rank_order2nd/3rd birth rank, >2 years 1.000 0.317109
## rank_order4th birth rank, <=2 years -0.402 0.688009
## rank_order4th birth rank, >2 years -0.278 0.781224
## rank_order5+ birth order 0.360 0.718551
## fac_CS_deliveryFac without CS delivery 2.271 0.023166 *
## fac_CS_deliveryHome 1.049 0.294083
## sibling_aliveDead 1.308 0.190828
## residenceUrban -0.050 0.959897
## regionAfar 0.832 0.405357
## regionAmhara 0.476 0.634353
## regionBen-Gumuz -0.199 0.842085
## regionDire Dawa 0.678 0.497576
## regionGambella 0.510 0.610076
## regionHarari 0.249 0.803278
## regionOromia 0.825 0.409261
## regionSNNP 0.767 0.442913
## regionSomali 0.417 0.676465
## regionTigray 0.014 0.989069
## v201 2.277 0.022765 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef)
## yr2016 1.1466 0.8722
## mother_age20-29 0.3870 2.5839
## mother_age30-39 0.2484 4.0265
## mother_age40-49 0.2428 4.1181
## educ1Primary 1.1181 0.8944
## educ2Secodary and Higher 0.7354 1.3598
## occup1Unemployed 0.8637 1.1578
## occup2Agricultural/Manual(skilled/unskilled) 1.1286 0.8861
## marst1nevermarried 0.4363 2.2921
## marst3widowed/divorced/separated 1.6293 0.6138
## child_sexMale 1.7556 0.5696
## child_wanted2Later 0.7803 1.2816
## child_wanted3Not at all 0.9504 1.0522
## TT_pregnancy1TT 1.1416 0.8759
## TT_pregnancy2+TT 0.6838 1.4623
## contraceptive_useYes 0.6687 1.4955
## antenatal_visits1-4 visits 1.0144 0.9858
## antenatal_visits5+ Visits 0.8010 1.2484
## post_natal1No 2.1865 0.4574
## wealthindex2Middle 0.8349 1.1977
## wealthindex3Rich 1.1124 0.8989
## rank_order2nd/3rd birth rank, <=2 years 1.3051 0.7662
## rank_order2nd/3rd birth rank, >2 years 1.3149 0.7605
## rank_order4th birth rank, <=2 years 0.8085 1.2368
## rank_order4th birth rank, >2 years 0.8987 1.1127
## rank_order5+ birth order 1.1702 0.8546
## fac_CS_deliveryFac without CS delivery 2.4825 0.4028
## fac_CS_deliveryHome 1.3646 0.7328
## sibling_aliveDead 1.2198 0.8198
## residenceUrban 0.9855 1.0147
## regionAfar 1.5215 0.6573
## regionAmhara 1.2810 0.7806
## regionBen-Gumuz 0.9026 1.1079
## regionDire Dawa 1.4353 0.6967
## regionGambella 1.2878 0.7765
## regionHarari 1.1347 0.8813
## regionOromia 1.5138 0.6606
## regionSNNP 1.4682 0.6811
## regionSomali 1.2290 0.8137
## regionTigray 1.0069 0.9932
## v201 1.1383 0.8785
## lower .95 upper .95
## yr2016 0.8339 1.5766
## mother_age20-29 0.2229 0.6719
## mother_age30-39 0.1238 0.4984
## mother_age40-49 0.1082 0.5451
## educ1Primary 0.8261 1.5133
## educ2Secodary and Higher 0.2427 2.2280
## occup1Unemployed 0.5838 1.2779
## occup2Agricultural/Manual(skilled/unskilled) 0.7210 1.7667
## marst1nevermarried 0.1131 1.6835
## marst3widowed/divorced/separated 1.0359 2.5626
## child_sexMale 1.3490 2.2847
## child_wanted2Later 0.5370 1.1338
## child_wanted3Not at all 0.5960 1.5156
## TT_pregnancy1TT 0.7106 1.8341
## TT_pregnancy2+TT 0.4731 0.9885
## contraceptive_useYes 0.4623 0.9673
## antenatal_visits1-4 visits 0.7144 1.4403
## antenatal_visits5+ Visits 0.4320 1.4850
## post_natal1No 1.1059 4.3228
## wealthindex2Middle 0.5568 1.2521
## wealthindex3Rich 0.8030 1.5410
## rank_order2nd/3rd birth rank, <=2 years 0.6828 2.4945
## rank_order2nd/3rd birth rank, >2 years 0.7691 2.2480
## rank_order4th birth rank, <=2 years 0.2865 2.2814
## rank_order4th birth rank, >2 years 0.4231 1.9092
## rank_order5+ birth order 0.4978 2.7506
## fac_CS_deliveryFac without CS delivery 1.1325 5.4418
## fac_CS_deliveryHome 0.7635 2.4387
## sibling_aliveDead 0.9057 1.6429
## residenceUrban 0.5571 1.7433
## regionAfar 0.5662 4.0887
## regionAmhara 0.4617 3.5541
## regionBen-Gumuz 0.3293 2.4740
## regionDire Dawa 0.5052 4.0781
## regionGambella 0.4872 3.4041
## regionHarari 0.4199 3.0664
## regionOromia 0.5655 4.0525
## regionSNNP 0.5505 3.9161
## regionSomali 0.4666 3.2369
## regionTigray 0.3789 2.6756
## v201 1.0182 1.2725
##
## Concordance= 0.68 (se = 0.013 )
## Rsquare= NA (max possible= NA )
## Likelihood ratio test= NA on 41 df, p=NA
## Wald test = 110 on 41 df, p=3.191e-08
## Score (logrank) test = NA on 41 df, p=NAsummary(cox.child)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (1239) clusters.
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = edhs_all[is.na(edhs_all$weight) == F, ], nest = T)## Call:
## svycoxph(formula = Surv(death.age, d.event) ~ yr + mother_age +
## educ + occup + marst + child_sex + child_wanted + TT_pregnancy +
## educ + contraceptive_use + antenatal_visits + post_natal +
## wealthindex + rank_order + child_wanted + fac_CS_delivery +
## sibling_alive + residence + region + v201, design = des3)
##
## n= 12991, number of events= 534
## (8179 observations deleted due to missingness)
##
## coef exp(coef) se(coef)
## yr2016 0.12282 1.13069 0.15100
## mother_age20-29 -0.99808 0.36859 0.28740
## mother_age30-39 -1.39305 0.24832 0.34492
## mother_age40-49 -1.40294 0.24587 0.38630
## educ1Primary 0.11135 1.11778 0.14805
## educ2Secodary and Higher -0.35045 0.70437 0.55979
## occup1Unemployed -0.03614 0.96450 0.18213
## occup2Agricultural/Manual(skilled/unskilled) 0.10109 1.10637 0.21038
## marst1nevermarried -0.87643 0.41627 0.68113
## marst3widowed/divorced/separated 0.58651 1.79770 0.20249
## child_sexMale 0.52646 1.69294 0.12191
## child_wanted2Later -0.29895 0.74160 0.18157
## child_wanted3Not at all -0.12492 0.88257 0.22705
## TT_pregnancy1TT 0.07388 1.07668 0.23638
## TT_pregnancy2+TT -0.33062 0.71848 0.17110
## contraceptive_useYes -0.41027 0.66347 0.18040
## antenatal_visits1-4 visits -0.05364 0.94778 0.16886
## antenatal_visits5+ Visits -0.24119 0.78569 0.28798
## post_natal1No 0.71399 2.04212 0.33095
## wealthindex2Middle -0.12749 0.88031 0.18119
## wealthindex3Rich 0.03158 1.03208 0.15610
## rank_order2nd/3rd birth rank, <=2 years 0.16564 1.18014 0.32118
## rank_order2nd/3rd birth rank, >2 years 0.23360 1.26314 0.25589
## rank_order4th birth rank, <=2 years 0.34827 1.41662 0.40034
## rank_order4th birth rank, >2 years 0.07783 1.08093 0.33817
## rank_order5+ birth order 0.18621 1.20467 0.40284
## fac_CS_deliveryFac without CS delivery 0.90434 2.47029 0.39947
## fac_CS_deliveryHome 0.39375 1.48253 0.29098
## sibling_aliveDead 0.18767 1.20643 0.14405
## residenceUrban 0.07811 1.08124 0.28802
## regionAfar 0.53268 1.70348 0.48947
## regionAmhara 0.33086 1.39216 0.50490
## regionBen-Gumuz 0.09530 1.09999 0.49609
## regionDire Dawa 0.50633 1.65919 0.53413
## regionGambella 0.42630 1.53158 0.47884
## regionHarari 0.28020 1.32339 0.49133
## regionOromia 0.56690 1.76280 0.48971
## regionSNNP 0.45750 1.58012 0.49111
## regionSomali 0.23433 1.26406 0.48071
## regionTigray 0.24694 1.28010 0.48307
## v201 0.12169 1.12940 0.05423
## z Pr(>|z|)
## yr2016 0.813 0.415978
## mother_age20-29 -3.473 0.000515 ***
## mother_age30-39 -4.039 5.37e-05 ***
## mother_age40-49 -3.632 0.000282 ***
## educ1Primary 0.752 0.451983
## educ2Secodary and Higher -0.626 0.531293
## occup1Unemployed -0.198 0.842701
## occup2Agricultural/Manual(skilled/unskilled) 0.481 0.630870
## marst1nevermarried -1.287 0.198186
## marst3widowed/divorced/separated 2.897 0.003773 **
## child_sexMale 4.318 1.57e-05 ***
## child_wanted2Later -1.646 0.099675 .
## child_wanted3Not at all -0.550 0.582199
## TT_pregnancy1TT 0.313 0.754617
## TT_pregnancy2+TT -1.932 0.053318 .
## contraceptive_useYes -2.274 0.022955 *
## antenatal_visits1-4 visits -0.318 0.750750
## antenatal_visits5+ Visits -0.838 0.402293
## post_natal1No 2.157 0.030975 *
## wealthindex2Middle -0.704 0.481687
## wealthindex3Rich 0.202 0.839682
## rank_order2nd/3rd birth rank, <=2 years 0.516 0.606049
## rank_order2nd/3rd birth rank, >2 years 0.913 0.361288
## rank_order4th birth rank, <=2 years 0.870 0.384332
## rank_order4th birth rank, >2 years 0.230 0.817985
## rank_order5+ birth order 0.462 0.643913
## fac_CS_deliveryFac without CS delivery 2.264 0.023583 *
## fac_CS_deliveryHome 1.353 0.175990
## sibling_aliveDead 1.303 0.192648
## residenceUrban 0.271 0.786250
## regionAfar 1.088 0.276471
## regionAmhara 0.655 0.512273
## regionBen-Gumuz 0.192 0.847666
## regionDire Dawa 0.948 0.343157
## regionGambella 0.890 0.373320
## regionHarari 0.570 0.568488
## regionOromia 1.158 0.247012
## regionSNNP 0.932 0.351566
## regionSomali 0.487 0.625926
## regionTigray 0.511 0.609224
## v201 2.244 0.024846 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef)
## yr2016 1.1307 0.8844
## mother_age20-29 0.3686 2.7131
## mother_age30-39 0.2483 4.0271
## mother_age40-49 0.2459 4.0671
## educ1Primary 1.1178 0.8946
## educ2Secodary and Higher 0.7044 1.4197
## occup1Unemployed 0.9645 1.0368
## occup2Agricultural/Manual(skilled/unskilled) 1.1064 0.9039
## marst1nevermarried 0.4163 2.4023
## marst3widowed/divorced/separated 1.7977 0.5563
## child_sexMale 1.6929 0.5907
## child_wanted2Later 0.7416 1.3484
## child_wanted3Not at all 0.8826 1.1331
## TT_pregnancy1TT 1.0767 0.9288
## TT_pregnancy2+TT 0.7185 1.3918
## contraceptive_useYes 0.6635 1.5072
## antenatal_visits1-4 visits 0.9478 1.0551
## antenatal_visits5+ Visits 0.7857 1.2728
## post_natal1No 2.0421 0.4897
## wealthindex2Middle 0.8803 1.1360
## wealthindex3Rich 1.0321 0.9689
## rank_order2nd/3rd birth rank, <=2 years 1.1801 0.8474
## rank_order2nd/3rd birth rank, >2 years 1.2631 0.7917
## rank_order4th birth rank, <=2 years 1.4166 0.7059
## rank_order4th birth rank, >2 years 1.0809 0.9251
## rank_order5+ birth order 1.2047 0.8301
## fac_CS_deliveryFac without CS delivery 2.4703 0.4048
## fac_CS_deliveryHome 1.4825 0.6745
## sibling_aliveDead 1.2064 0.8289
## residenceUrban 1.0812 0.9249
## regionAfar 1.7035 0.5870
## regionAmhara 1.3922 0.7183
## regionBen-Gumuz 1.1000 0.9091
## regionDire Dawa 1.6592 0.6027
## regionGambella 1.5316 0.6529
## regionHarari 1.3234 0.7556
## regionOromia 1.7628 0.5673
## regionSNNP 1.5801 0.6329
## regionSomali 1.2641 0.7911
## regionTigray 1.2801 0.7812
## v201 1.1294 0.8854
## lower .95 upper .95
## yr2016 0.8410 1.5201
## mother_age20-29 0.2098 0.6474
## mother_age30-39 0.1263 0.4882
## mother_age40-49 0.1153 0.5242
## educ1Primary 0.8363 1.4941
## educ2Secodary and Higher 0.2351 2.1101
## occup1Unemployed 0.6750 1.3783
## occup2Agricultural/Manual(skilled/unskilled) 0.7325 1.6710
## marst1nevermarried 0.1095 1.5818
## marst3widowed/divorced/separated 1.2088 2.6735
## child_sexMale 1.3331 2.1499
## child_wanted2Later 0.5195 1.0586
## child_wanted3Not at all 0.5656 1.3773
## TT_pregnancy1TT 0.6775 1.7112
## TT_pregnancy2+TT 0.5138 1.0047
## contraceptive_useYes 0.4659 0.9449
## antenatal_visits1-4 visits 0.6807 1.3196
## antenatal_visits5+ Visits 0.4468 1.3816
## post_natal1No 1.0675 3.9065
## wealthindex2Middle 0.6172 1.2556
## wealthindex3Rich 0.7600 1.4015
## rank_order2nd/3rd birth rank, <=2 years 0.6289 2.2147
## rank_order2nd/3rd birth rank, >2 years 0.7650 2.0858
## rank_order4th birth rank, <=2 years 0.6464 3.1047
## rank_order4th birth rank, >2 years 0.5571 2.0973
## rank_order5+ birth order 0.5470 2.6532
## fac_CS_deliveryFac without CS delivery 1.1291 5.4047
## fac_CS_deliveryHome 0.8382 2.6223
## sibling_aliveDead 0.9097 1.6000
## residenceUrban 0.6148 1.9014
## regionAfar 0.6527 4.4460
## regionAmhara 0.5175 3.7451
## regionBen-Gumuz 0.4160 2.9084
## regionDire Dawa 0.5824 4.7266
## regionGambella 0.5992 3.9150
## regionHarari 0.5052 3.4667
## regionOromia 0.6751 4.6030
## regionSNNP 0.6035 4.1374
## regionSomali 0.4927 3.2430
## regionTigray 0.4966 3.2994
## v201 1.0155 1.2561
##
## Concordance= 0.674 (se = 0.013 )
## Rsquare= NA (max possible= NA )
## Likelihood ratio test= NA on 41 df, p=NA
## Wald test = 120.3 on 41 df, p=1.012e-09
## Score (logrank) test = NA on 41 df, p=NAfit.test.child<-cox.zph(cox.child)
fit.test.child## rho chisq p
## yr2016 -0.01788 6.16e-01 4.33e-01
## mother_age20-29 0.00248 1.08e-02 9.17e-01
## mother_age30-39 -0.00249 1.08e-02 9.17e-01
## mother_age40-49 0.00913 1.31e-01 7.17e-01
## educ1Primary 0.02538 8.85e-01 3.47e-01
## educ2Secodary and Higher -0.05939 1.19e+01 5.58e-04
## occup1Unemployed -0.02444 8.91e-01 3.45e-01
## occup2Agricultural/Manual(skilled/unskilled) -0.02832 1.32e+00 2.50e-01
## marst1nevermarried -0.05439 2.56e+00 1.10e-01
## marst3widowed/divorced/separated -0.03559 1.78e+00 1.82e-01
## child_sexMale -0.08568 1.06e+01 1.15e-03
## child_wanted2Later -0.01323 2.82e-01 5.96e-01
## child_wanted3Not at all -0.00137 2.75e-03 9.58e-01
## TT_pregnancy1TT 0.05294 5.54e+00 1.86e-02
## TT_pregnancy2+TT 0.09716 1.58e+01 6.92e-05
## contraceptive_useYes -0.05076 4.19e+00 4.07e-02
## antenatal_visits1-4 visits -0.08386 1.20e+01 5.35e-04
## antenatal_visits5+ Visits -0.01084 2.14e-01 6.44e-01
## post_natal1No -0.05477 3.91e+00 4.79e-02
## wealthindex2Middle 0.02422 8.96e-01 3.44e-01
## wealthindex3Rich -0.01158 1.76e-01 6.75e-01
## rank_order2nd/3rd birth rank, <=2 years -0.06896 8.60e+00 3.35e-03
## rank_order2nd/3rd birth rank, >2 years -0.11056 2.28e+01 1.79e-06
## rank_order4th birth rank, <=2 years 0.02655 9.34e-01 3.34e-01
## rank_order4th birth rank, >2 years -0.05810 5.63e+00 1.77e-02
## rank_order5+ birth order -0.06451 6.76e+00 9.33e-03
## fac_CS_deliveryFac without CS delivery -0.02377 1.16e+00 2.81e-01
## fac_CS_deliveryHome 0.06226 6.41e+00 1.13e-02
## sibling_aliveDead -0.05474 4.47e+00 3.46e-02
## residenceUrban 0.08498 1.31e+01 3.03e-04
## regionAfar 0.03739 1.37e+00 2.42e-01
## regionAmhara -0.02057 4.47e-01 5.04e-01
## regionBen-Gumuz 0.02525 5.64e-01 4.53e-01
## regionDire Dawa 0.01143 1.26e-01 7.23e-01
## regionGambella 0.01532 1.85e-01 6.67e-01
## regionHarari 0.00728 4.56e-02 8.31e-01
## regionOromia 0.00803 6.74e-02 7.95e-01
## regionSNNP 0.02429 5.75e-01 4.48e-01
## regionSomali -0.02214 4.64e-01 4.96e-01
## regionTigray 0.01710 2.50e-01 6.17e-01
## v201 0.00806 9.03e-02 7.64e-01
## GLOBAL NA 1.36e+02 3.77e-12par(mfrow=c(3,3))
plot(fit.test.child, df=2)
par(mfrow=c(1,1))
From the above cox proportional hazard model we see that some of the predictors appear to be not associated with the residuals, which suggests t suggest that they are having some proportionality issue. In order to resolve the proportionality issue, we stratify the varaible which were not proportional using type of place of residence.
## stratifying some varaibles in order to resolve the problem
cox.child<-svycoxph(Surv(death.age,d.event)~yr+mother_age+educ+occup+marst+child_sex+child_wanted+TT_pregnancy+educ+contraceptive_use+antenatal_visits+post_natal+wealthindex+rank_order+child_wanted+fac_CS_delivery+sibling_alive+strata(residence)+region+v201, design=des3)b) Cox Fraility Model
Next, we apply the cox fraility model in order to examine wheather the level of under-5 mortality varied by region or not.
fit.cox<-coxph(Surv(death.age, d.event)~yr+mother_age+educ+occup+marst+child_sex+child_wanted+TT_pregnancy+educ+contraceptive_use+antenatal_visits+post_natal+wealthindex+rank_order+child_wanted+fac_CS_delivery+sibling_alive+residence+region+v201, data=edhs_all, weights=weight)
fit.frail<-coxme(Surv(death.age, d.event)~yr+mother_age+educ+occup+marst+child_sex+child_wanted+TT_pregnancy+educ+contraceptive_use+antenatal_visits+post_natal+wealthindex+rank_order+child_wanted+fac_CS_delivery+sibling_alive+residence+(1|region)+v201, data=edhs_all, weights=weight)
summary(fit.frail)## Cox mixed-effects model fit by maximum likelihood
## Data: edhs_all
## events, n = 534, 12991 (8179 observations deleted due to missingness)
## Iterations= 5 33
## NULL Integrated Fitted
## Log-likelihood -5433.298 -5324.382 -5322.434
##
## Chisq df p AIC BIC
## Integrated loglik 217.83 32.00 0 153.83 16.86
## Penalized loglik 221.73 32.62 0 156.49 16.87
##
## Model: Surv(death.age, d.event) ~ yr + mother_age + educ + occup + marst + child_sex + child_wanted + TT_pregnancy + educ + contraceptive_use + antenatal_visits + post_natal + wealthindex + rank_order + child_wanted + fac_CS_delivery + sibling_alive + residence + (1 | region) + v201
## Fixed coefficients
## coef exp(coef)
## yr2016 0.12903394 1.1377287
## mother_age20-29 -0.99988256 0.3679226
## mother_age30-39 -1.40791922 0.2446518
## mother_age40-49 -1.42705446 0.2400149
## educ1Primary 0.12167367 1.1293855
## educ2Secodary and Higher -0.35576681 0.7006360
## occup1Unemployed -0.04069724 0.9601198
## occup2Agricultural/Manual(skilled/unskilled) 0.07703392 1.0800787
## marst1nevermarried -0.89165548 0.4099765
## marst3widowed/divorced/separated 0.57905143 1.7843450
## child_sexMale 0.52680836 1.6935186
## child_wanted2Later -0.28835839 0.7494929
## child_wanted3Not at all -0.11429667 0.8919933
## TT_pregnancy1TT 0.07590517 1.0788603
## TT_pregnancy2+TT -0.31815727 0.7274884
## contraceptive_useYes -0.41060841 0.6632466
## antenatal_visits1-4 visits -0.06770338 0.9345376
## antenatal_visits5+ Visits -0.27750069 0.7576750
## post_natal1No 0.72928934 2.0736064
## wealthindex2Middle -0.11364645 0.8925735
## wealthindex3Rich 0.04963759 1.0508902
## rank_order2nd/3rd birth rank, <=2 years 0.17121408 1.1867448
## rank_order2nd/3rd birth rank, >2 years 0.23610307 1.2663048
## rank_order4th birth rank, <=2 years 0.35910166 1.4320424
## rank_order4th birth rank, >2 years 0.07973685 1.0830020
## rank_order5+ birth order 0.19148703 1.2110491
## fac_CS_deliveryFac without CS delivery 0.87996442 2.4108139
## fac_CS_deliveryHome 0.40304754 1.4963780
## sibling_aliveDead 0.18453159 1.2026550
## residenceUrban 0.01931691 1.0195047
## v201 0.12370325 1.1316800
## se(coef) z p
## yr2016 0.09157610 1.41 1.6e-01
## mother_age20-29 0.18587925 -5.38 7.5e-08
## mother_age30-39 0.22088257 -6.37 1.8e-10
## mother_age40-49 0.26118873 -5.46 4.7e-08
## educ1Primary 0.10105528 1.20 2.3e-01
## educ2Secodary and Higher 0.32294219 -1.10 2.7e-01
## occup1Unemployed 0.12760199 -0.32 7.5e-01
## occup2Agricultural/Manual(skilled/unskilled) 0.13201017 0.58 5.6e-01
## marst1nevermarried 0.87566238 -1.02 3.1e-01
## marst3widowed/divorced/separated 0.13292742 4.36 1.3e-05
## child_sexMale 0.08688380 6.06 1.3e-09
## child_wanted2Later 0.11827539 -2.44 1.5e-02
## child_wanted3Not at all 0.14105723 -0.81 4.2e-01
## TT_pregnancy1TT 0.13812879 0.55 5.8e-01
## TT_pregnancy2+TT 0.10857907 -2.93 3.4e-03
## contraceptive_useYes 0.10687765 -3.84 1.2e-04
## antenatal_visits1-4 visits 0.10617289 -0.64 5.2e-01
## antenatal_visits5+ Visits 0.17953350 -1.55 1.2e-01
## post_natal1No 0.26414409 2.76 5.8e-03
## wealthindex2Middle 0.11139958 -1.02 3.1e-01
## wealthindex3Rich 0.10683047 0.46 6.4e-01
## rank_order2nd/3rd birth rank, <=2 years 0.22242395 0.77 4.4e-01
## rank_order2nd/3rd birth rank, >2 years 0.16755125 1.41 1.6e-01
## rank_order4th birth rank, <=2 years 0.30040761 1.20 2.3e-01
## rank_order4th birth rank, >2 years 0.22676083 0.35 7.3e-01
## rank_order5+ birth order 0.25589008 0.75 4.5e-01
## fac_CS_deliveryFac without CS delivery 0.24989206 3.52 4.3e-04
## fac_CS_deliveryHome 0.19248612 2.09 3.6e-02
## sibling_aliveDead 0.09261530 1.99 4.6e-02
## residenceUrban 0.17981862 0.11 9.1e-01
## v201 0.03562683 3.47 5.2e-04
##
## Random effects
## Group Variable Std Dev Variance
## region Intercept 0.086929609 0.007556757anova(fit.cox, fit.frail)## Analysis of Deviance Table
## Cox model: response is Surv(death.age, d.event)
## Model 1: ~yr + mother_age + educ + occup + marst + child_sex + child_wanted + TT_pregnancy + educ + contraceptive_use + antenatal_visits + post_natal + wealthindex + rank_order + child_wanted + fac_CS_delivery + sibling_alive + residence + region + v201
## Model 2: ~yr + mother_age + educ + occup + marst + child_sex + child_wanted + TT_pregnancy + educ + contraceptive_use + antenatal_visits + post_natal + wealthindex + rank_order + child_wanted + fac_CS_delivery + sibling_alive + residence + (1 | region) + v201
## loglik Chisq Df P(>|Chi|)
## 1 -5320.8
## 2 -5324.4 7.1505 9 0.6215AIC(fit.cox)## [1] 10723.61AIC(fit.frail)## [1] 10710.11So, we see that the frailty model has amost equal variance, and that the AICs does not show a significance diference to the ordinary Cox model, suggesting that both models fit almost equally. The model likelihood ratio test also confirms that the frailty model fits almost equally as the ordinary cox model.
E) Discrete-time Models
- Perform a person-period data set.
edhs_all$death.age<-ifelse(edhs_all$b5==1, ((((edhs_all$v008))+1900)-(((edhs_all$b3))+1900)),edhs_all$b7)+1
#censoring indicator for death by age 5, in months (<=60 months)
edhs_all$d.event<-ifelse(is.na(edhs_all$b7)==T|edhs_all$b7>60,0,1)
edhs_all$d.eventfac<-factor(edhs_all$d.event); levels(edhs_all$d.eventfac)<-c("Alive at Age 5", "Dead by Age 5")
table(edhs_all$d.eventfac)##
## Alive at Age 5 Dead by Age 5
## 19974 1196#make person period file
pp<-survSplit(Surv(death.age, d.event)~., data=edhs_all, cut=c(12, 24, 36, 48, 60), episode="year_age", subset = death.age!=0)
head(pp[, c("death.age", "d.event", "year_age", "agefb", "residence")], n=20)## death.age d.event year_age agefb residence
## 1 12 0 1 <20 Rural
## 2 24 0 2 <20 Rural
## 3 36 0 3 <20 Rural
## 4 48 0 4 <20 Rural
## 5 60 0 5 <20 Rural
## 6 76 0 6 <20 Rural
## 7 12 0 1 <20 Rural
## 8 24 0 2 <20 Rural
## 9 36 0 3 <20 Rural
## 10 48 0 4 <20 Rural
## 11 60 0 5 <20 Rural
## 12 68 0 6 <20 Rural
## 13 12 0 1 <20 Rural
## 14 24 0 2 <20 Rural
## 15 36 0 3 <20 Rural
## 16 44 0 4 <20 Rural
## 17 12 0 1 <20 Rural
## 18 24 0 2 <20 Rural
## 19 36 0 3 <20 Rural
## 20 43 0 4 <20 Rural#generate survey design
pp<-pp[pp$v201!=0,]
options(survey.lonely.psu = "adjust")
pp$weight<-pp$v005/1000000
options(survey.lonely.psu = "adjust")
pp$yrstratum<-paste(pp$yr, pp$v022, sep="-")
pp$yrpsu<-paste(pp$yr, pp$v021, sep="-")
des<-svydesign(ids=~yrpsu, strata=~yrstratum, weights=~weight, data = pp[is.na(pp$weight)==F,], nest=T)
#Fit the basic logistic model with ONLY time in the model
#I do -1 so that no intercept is fit in the model, and we get a hazard estimate for each time period
fit.0<-svyglm(d.event ~as.factor(year_age)-1,design=des , family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!summary(fit.0)##
## Call:
## svyglm(formula = d.event ~ as.factor(year_age) - 1, design = des,
## family = "binomial")
##
## Survey design:
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = pp[is.na(pp$weight) == F, ], nest = T)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## as.factor(year_age)1 -3.08708 0.05027 -61.41 <2e-16 ***
## as.factor(year_age)2 -4.92958 0.14525 -33.94 <2e-16 ***
## as.factor(year_age)3 -5.25432 0.17654 -29.76 <2e-16 ***
## as.factor(year_age)4 -5.34224 0.22444 -23.80 <2e-16 ***
## as.factor(year_age)5 -5.72396 0.33266 -17.21 <2e-16 ***
## as.factor(year_age)6 -5.95162 0.38200 -15.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.000014)
##
## Number of Fisher Scoring iterations: 8#Plot the hazard function on the probability scale
haz<-1/(1+exp(-coef(fit.0)))
time<-seq(1,6,1)
plot(haz~time, type="l", ylab="h(t)")
title(main="Discrete Time Hazard Function for Child Mortality")
#now we include a single predictor and examine the proportional hazards
fit.1<-svyglm(d.event~as.factor(year_age)+rural-1,design=des , family=binomial(link="cloglog"))## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!###fit.1<-svyglm(d.event~as.factor(year_age)+rural-1+educ+parity+age+agefb+antenatal, design=des, family=binomial(link="cloglog"))
summary(fit.1)##
## Call:
## svyglm(formula = d.event ~ as.factor(year_age) + rural - 1, design = des,
## family = binomial(link = "cloglog"))
##
## Survey design:
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = pp[is.na(pp$weight) == F, ], nest = T)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## as.factor(year_age)1 -3.3303 0.1127 -29.538 <2e-16 ***
## as.factor(year_age)2 -5.1502 0.1627 -31.648 <2e-16 ***
## as.factor(year_age)3 -5.4689 0.1917 -28.527 <2e-16 ***
## as.factor(year_age)4 -5.5494 0.2374 -23.372 <2e-16 ***
## as.factor(year_age)5 -5.9243 0.3375 -17.551 <2e-16 ***
## as.factor(year_age)6 -6.1464 0.3844 -15.988 <2e-16 ***
## rural 0.2633 0.1184 2.225 0.0263 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.981687)
##
## Number of Fisher Scoring iterations: 8t<-data.frame(cbind(y=1/(1+exp(-fit.1$linear.predictors)), year=pp$year_age,rural=pp$rural))
t0<-t[t$rural==0,]
t0<-unique(t0[order(t0$year),])
t1<-t[t$rural==1,]
t1<-unique(t1[order(t1$year),])
#Here we plot the hazards, I edhs_alltract .5 from each time so it looks like the hazard is at the midpoint of the year
plot(t0$year-.5, t0$y, type="b", ylab="Hazard", xlim=c(1,6), xlab="Year", col="red", main =c("Discrete time hazard model", "Effect of Rural Residence"))
points(t1$year-.5, t1$y, type="b", col="blue")
legend("topright", legend=c("Urban", "Rural"), col=c("red", "blue"),lty=1)
t<-data.frame(cbind(y=1/(1+exp(-fit.1$linear.predictors)), year=pp$year_age, post_natal_yes=pp$post_natal_yes))
t0<-t[t$post_natal_yes==0,]
t0<-unique(t0[order(t0$year),])
t1<-t[t$post_natal_yes==1,]
t1<-unique(t1[order(t1$year),])
#Here we plot the hazards, I edhs_alltract .5 from each time so it looks like the hazard is at the midpoint of the year
plot(t0$year-.5, t0$y, type="b", ylab="Hazard", ylim=c(0,0.02), xlim=c(1,6), xlab="Year", col="red", main =c("Discrete time hazard model", "Effect of Postnatal Care Visits"))
points(t1$year-.5, t1$y, type="b", col="blue")
legend("topright", legend=c("Post Natal Care No", "Post Natal Care-Yes"), col=c("red", "blue"),lty=1)
fit.2<-svyglm(d.event~as.factor(year_age)*Sec_educ-1,design=des , family=binomial(link="cloglog"))## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!summary(fit.2)##
## Call:
## svyglm(formula = d.event ~ as.factor(year_age) * Sec_educ - 1,
## design = des, family = binomial(link = "cloglog"))
##
## Survey design:
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = pp[is.na(pp$weight) == F, ], nest = T)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## as.factor(year_age)1 -3.07854 0.04966 -61.997 < 2e-16 ***
## as.factor(year_age)2 -4.86214 0.14590 -33.324 < 2e-16 ***
## as.factor(year_age)3 -5.20341 0.18060 -28.812 < 2e-16 ***
## as.factor(year_age)4 -5.25247 0.22332 -23.520 < 2e-16 ***
## as.factor(year_age)5 -5.62511 0.33215 -16.935 < 2e-16 ***
## as.factor(year_age)6 -5.85120 0.38249 -15.297 < 2e-16 ***
## Sec_educ -0.52988 0.30289 -1.749 0.08048 .
## as.factor(year_age)2:Sec_educ -1.64693 0.74761 -2.203 0.02779 *
## as.factor(year_age)3:Sec_educ -0.51162 0.57036 -0.897 0.36989
## as.factor(year_age)4:Sec_educ -12.75801 0.36624 -34.836 < 2e-16 ***
## as.factor(year_age)5:Sec_educ -12.37756 0.44267 -27.961 < 2e-16 ***
## as.factor(year_age)6:Sec_educ -3.06517 1.10711 -2.769 0.00572 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.9793705)
##
## Number of Fisher Scoring iterations: 17regTermTest(fit.2, ~as.factor(year_age):Sec_educ)## Wald test for as.factor(year_age):Sec_educ
## in svyglm(formula = d.event ~ as.factor(year_age) * Sec_educ - 1,
## design = des, family = binomial(link = "cloglog"))
## F = 356.6841 on 5 and 1180 df: p= < 2.22e-16#No effect of the interaction
dat<-expand.grid(year_age=1:6, Sec_educ=c(0,1))
dat$fitted<-predict(fit.2, type = "response", newdata=dat)
#plot the interaction with time
plot(dat$year[dat$Sec_educ==0]-.5, dat$fitted[dat$Sec_educ==0], type="b", ylab="Hazard", xlim=c(0,6), xlab="Year", col="red", main =c("Discrete time hazard model", "effect of High School or Above Education "))
points(dat$year[dat$Sec_educ==1]-.5, dat$fitted[dat$Sec_educ==1], type="b", col="blue")
legend("topright", legend=c("Primary or less", " Secondary and above"), col=c("red", "blue"),lty=1)
Which shows a lower risk of death for children from primary or less education.
fit.2<-svyglm(d.event~rural*Sec_educ-1,design=des , family=binomial(link="cloglog"))## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!summary(fit.2)##
## Call:
## svyglm(formula = d.event ~ rural * Sec_educ - 1, design = des,
## family = binomial(link = "cloglog"))
##
## Survey design:
## svydesign(ids = ~yrpsu, strata = ~yrstratum, weights = ~weight,
## data = pp[is.na(pp$weight) == F, ], nest = T)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## rural -3.99677 0.05272 -75.816 < 2e-16 ***
## Sec_educ -5.05019 0.34041 -14.835 < 2e-16 ***
## rural:Sec_educ 4.78971 0.61388 7.802 1.32e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.068043)
##
## Number of Fisher Scoring iterations: 7regTermTest(fit.2, ~rural:Sec_educ)## Wald test for rural:Sec_educ
## in svyglm(formula = d.event ~ rural * Sec_educ - 1, design = des,
## family = binomial(link = "cloglog"))
## F = 60.87652 on 1 and 1189 df: p= 1.3235e-14#There is interaction effect between education and residence.