Survival Analysis of Twin Under Five Mortality in Ethiopia: Using Gamma Frailty Modelling Approach

Twins are relatively rare events, but several studies confirm that they contribute substantially to mortality in both neonatal and post-neonatal periods. The excess infant and child mortality rates among twins calls for a need to identify the main causes behind it. In studies involving multiple individuals from the same family, the independence assumption is not plausible unless all important familial factors were measured and controlled for in the model. Hence, it is important to include random effect or frailty in the model to account for the existence of correlation or unmeasured covariates. This study, therefore, employed survival analysis using gamma frailty models. Information on twin under-five mortality was found from the birth history of women who were included in the 2016 Ethiopia Demographic and Health survey. The results of the study show that sex of child, between twins birth order, preceding birth interval, and succeeding birth interval are significantly associated with under-five twin mortality. In addition, the result of the study confirmed the significance of shared /family random effect.


Introduction
The rate of occurrence of infant and child mortality rates reflects the country's level of socio economic development and quality of life. According to the WHO's report, 5.9 million Children under age five died in 2015 and in average 16 000 children under age five dies every day. The hazard of death for child under-five years is still highest in African Region (81 per 1000 live births), this is approximately 7.37 times the risk in European Region (11 per 1000 live births). Fortunately, under-five mortality rate decreased by 53%, from an estimated rate of 91 deaths per 1000 live births in 1990 to 43 deaths per 1000 live births in 2015 (WHO, 2016).
There are huge differences in child mortality among low and middle income countries and the industrial world with Sub-Saharan Africa and South East Asia carrying the highest burden of under-five mortality (You, New, & Wardlaw, 2014;WHO, 2006). The highest rates of child mortality are still in Sub-Saharan Africa-where 1 in 8 children dies before age 5, more than 17 times the average for developed regions (1 in 143)-and Southern Asia (1 in 15) (UNICEF, 2011). Ethiopia being one these country under five mortality still remain high, 88/1000 live birth (CSA, 2012) .
Sub-Saharan Africa not only has the highest rate of twin births in the world (Smits & Monden, 2011) but also the world's highest rate of infant and child mortality (WHO, 2016). Twin births contribute substantially to infant and child mortality in both neonatal and post-neonatal periods (Alam, Van Ginneken, & Bosch, 2007).
In Sub-Saharan African countries, there are few studies that evaluate the mortality of twins after the neonatal period. All of these studies showed the excess mortality of twin compared to singleton (Guo & Grummer-Strawn, 1993;Pison, 1992;Justesen & Kunst, 2000;Olalekan, Mubashir, & Ismail, 2000). Close contact between twin babies increases the chance of cross-infections (Justesen & Kunst, 2000). In addition, worldwide there are several studies on infant and child mortality showed that infant and child mortality is mainly associated twin birth (Koissi & Hgns, 2001;Munyamahoro, 2017;Bereka, Habtewold, & Nebi, 2017) among others. Furthermore, in studies involving multiple individuals from the same family (e.g. twins), the presence of correlation is expected, and obviously independence assumption is not plausible unless all important familial factors were measured and controlled for in the model. This is because, children belong to the same family share certain unobserved characteristics (or heterogeneity), which may not be sufficiently described by the covariates included in the models. The ignorance of such family-level correlation may lead to biased parameters estimates (Guo & Rodriguez, 1992).
The persistence of high level of infant and child mortality rates among twins call for a need to identify the determinants of twins mortality separately. Identifying all significant determinants of twin birth child mortality is essential to form policies and strategies to accelerate the reduction of child mortality and also to meet the United Nations the sustainable development goal (SDGs) (i.e., aiming to reduce under-five mortality to at least as low as

Study design
The study design was a retrospective study design.

Variables of the study
The response or outcome variable for this study is the survival time of a pair of under-five twins measured in months. Among the various covariates that are thought to have effect on survival of a twin child, those time invariant factors such as mothers age at child birth (below18 years, between 18-35 years and above 35 years) , sex of the child, preceding birth interval (1 st born ,below 18 months, between 18-24 months and above 24 months) , succeeding birth interval (below 18 months, between 18-24 months, above 24 months and last born), between twins birth order (first born, second born) and less likely to be time variant factor such as residence (urban, rural) were included in this study.

Method of Data analysis
Descriptive statistics tools such as frequency, percentage and cross tabulation was used to present the data. The log rank test was performed to test if there are statistically significant differences among the survival experience of the different groups of the covariates. When proceed to model fitting, first, the usual Cox Proportional Hazard (PH) model was fit including all the potential risk factors that were significant in log rank test. Then, univariate and shared gamma frailty variable (random effect) was introduced on the full multivariable Cox PH model to assess the presences and significances of individual and shared frailty. Likelihood ratio test was employed to test the significance of the random effect. The level of significance for all the inferential statistics was set at p<0.05. Statistical R software was used.

Results
A total of 908 or (454 pairs of) twin child deliveries were recorded in the 2016 Ethiopia Demographic and Health survey. The data on 622 children or (311 pairs of twin) born before 2011 were analyzed in this study. The overall information on censoring and covariates that are included in this study are presented below. 2) As can be seen from Table 1 , regarding age of mother at child birth, in the first category below 18 years, there were 6.8% of the study population, (of whom 69% of the children have died and 31% are alive); 6.1% of the children were born from mothers whose age exceeds 35 years, (of whom 60.5% of the children have died and 39.5% are alive); and the remaining 87.1% children belong to mothers whose age at birth was between 18 -35 years (of whom 42.3% of have died and 57.7% are alive). Regarding place of residence, 16.7% were from urban area (of which 55.8% are alive and 44.2% are dead); the rest 83.3% were from rural area (of which 54.6% are alive and 45.4% are dead). Regarding the sex of children, 55.5% were male (of which 51.6% are alive and 48.4% are dead); the rest 44.5% were female (of which 58.8% are alive and 41.2% are dead). Among the first born twin children, 61.1% of them are alive and 38.9% are dead; while from the second born twin children, 51.4% are dead and 48.6% are alive.
Regarding preceding child's birth order, 18.7% of the children were first born (of which 51.7% are alive and 48.3% are dead); 11.9% were born before their older sibling reaches the age 18 months (of which 25.7% are alive and 74.3% are dead) ; 18.8% of the children were born when their older sibling's age was between 18 -24 months (of which 57.3% are alive and 42.7% are dead); and the remaining 50.6% of the study population have older sibling whose age exceed 24 months (of which 61.9% of the children are alive and 38.1% are dead).
Regarding the succeeding birth interval, 15.1% were last born or no child after (of which 67% are alive and 33% are dead); 16.9% were having a succeeding birth interval below 18 months (of which 29.5% are alive and 70.5% are dead); 17.7% were having a succeeding birth interval between 18 months and 24 months (of which 57.3% are alive and 42.7% are dead); and the rest, 50.3% were having a succeeding birth interval of more than 24 months (of which 58.8% are alive and 41.2% are dead). .001* To test if there are statistically significant differences among the survival experience of the different groups of the covariates, the log rank test is performed. As one can see from Table 2, the results of the log rank test indicate that there are significant differences between the survival experiences of different groups for all the covariates except residence.
The log rank test was also used as a univariate analysis to check the significance of incorporating each one of the independent variables in a model that contains several (multiple) independent variables. If a given independent variable is insignificant in log rank test, then it is unnecessary to incorporate the variable in model fitting that contains multiple independent variables such as Cox PH model. The Cox PH model that employed in several child mortality studies assume a homogeneous population. In other words, Cox PH assume that all individuals sampled into that study are subject in principle under the same risk. Although Cox PH model's independence assumption is not plausible, for making comparison we fit Cox PH model. All independent variables that were found to be significant in log rank test were included in the model. The results are given in Table 3 Table 3 -1738.679 The results of the analysis showed that all included covariates such as age of mother at child birth, sex of child, twin birth order, preceding birth interval and succeeding birth interval were found to be significantly associated with twin under-five mortality. In order to check the validity of proportionality assumption, we fit Cox PH with all covariates and time-dependent interaction terms for all covariates simultaneously. The result showed that the coefficients of all time-dependent interaction terms are insignificant. This indicates that the proportional hazards assumption is not violated for all covariates.  -1738.679 The Cox's regression model cannot take the presence of the unobserved heterogeneity in the population into account. As a result parameter estimates are often biased when some influential factors are not observed .To account for unobserved heterogeneity this model should be modified to include hidden frailty together with observed covariate (Iachine, 1995). Each child has a proper susceptibility to infection, independently of his family membership. In addition, inside the common global family behavioural factor, parents may adopt a slightly different prenatal and neonatal attitude from one child to the next (Childs, Moxon, & Winkelstein, 1992).
Thus, univariate semi-parametric gamma frailty model was fit to assess the presence of unobserved heterogeneity at individual (child) level. The results of the univariate gamma frailty model were found to be exactly similar to the results of the Cox PH model given in Table 3. This is because the estimate of the heterogeneity parameter became 5e-07 that is almost zero. And the results of the hypothesis test about the significance of the variance of the random effect i.e.
: : found to be insignificant. This indicates that, there is no significant unobserved heterogeneity in the study population (at individual level). As a result, the estimates of covariate effect are exactly the same to the estimates of the Cox PH model given Table 3 Then shared gamma frailty is introduced on the full multivariable Cox PH model to assess the random effect shared by the twins. The result is given in Table 4 Comparison of the two models was made to decide which model fit the data better. As can be seen from Table 5, the shared gamma frailty model has the highest log likelihood and minimum AIC and BIC values, indicating that this model fits the data better than the Cox PH model which did not take in to account the shared random effect. Furthermore, as can be seen from Table 4 an estimate of the heterogeneity parameter is 1.969. Since Cox PH model is nested within shared gamma frailty model we conduct the likelihood ratio test for the significance of the heterogeneity parameter i.e. : : . As can be seen from Table 4,  : is rejected with p value <0.001 indicates that one cannot ignore the correlation between the pair of twins. Hence, shared gamma frailty model fit the data better than the Cox PH model which did not take in to account the presence of correlation. The results of the analysis using shared gamma frailty model given in Table 4 showed that covariates such as sex of child, twin birth order, preceding birth interval and succeeding birth interval were found to be significantly associated with twin under-five mortality. However, the incorporation of a shared gamma frailty variable make the variable age of mother at child birth insignificant at 5% level of significance.
As can be seen from Table 4, the estimated hazard ratio of a male twin child (i.e. a child of twin birth) is 1.733 (95% CI: 1.220 -2.462) implying that the risk of dying for a male twin child is 73.3% higher than a female twin child controlling for the other covariates in the model. This figure can be as low as 22.0% and as high as 146% with 95 % confidence.
The estimated hazard ratio of a first born twin child is 0.644 (95% CI: 0.503-0.826) implying that the risk of dying for a first born twin child is 35.6% less likely than a second born twin child (reference group) controlling for the other covariates in the model.
The estimated hazard ratio a child of twin birth who was born before the older sibling reaches the age of 18 months is 2.575 (95% CI: 1.332-4.977). Thus, the hazard rate of a child of twin birth who was born before the older sibling reaches the age of 18 months is 2.575 times than a child of twin birth who born when the older sibling's age exceeds 24 months (reference group) controlling for other covariates in the model. The confidence interval indicates that the risk of dying can be as low as 1.332 times and as high as 4.977 times than that of the reference group.
The estimated hazard ratio a child of twin birth who has a younger sibling that born within 18 months is 1.856 (95% CI: 1.054-3.266). This implies that the risk of dying for a child of twin birth who has a younger sibling that born within 18 months is 85.6% higher than a child of twin birth who has a younger sibling that born after 24 months controlling for the other covariates in the model. This figure can be as low as 5.4% and as high as 226% with 95 % confidence.

Discussion
The main aim of this study was obtaining the estimates of covariate effect on survival of under-5 twins. Since this study was about survival of under-5 twin children, it is expected to employ a model that can take the presence of correlation and unobserved heterogeneity into account. Thus, this study employs univariate and shared gamma frailty models.
The results of the univariate gamma frailty model showed that, the variance of the random effect is estimated to be almost zero, and the results of the Hypothesis test about the significance of the variance of the random effect found to be insignificant. This indicates that, there is no significant unobserved heterogeneity in the study population (at individual level). As a result, the estimates of covariate effect are almost similar to the estimates of the usual Cox PH model. Absence of heterogeneity at individual level does not indicate absence of correlation within groups. The estimate of the variance of the frailty from univariate data may has nothing to do with association. Univariate frailty variance is interpreted as a measure of unobserved heterogeneity in the study population (Wienke, 2011). The estimates of the regression parameters are adjusted for the unobserved heterogeneity in the population and only conditional given the same frailty.
As showed in Table 4, the variance of the random effect estimated from the shared gamma frailty model is significant at 5% level of significance. This indicates that, one cannot ignore the correlation between twins. Comparison of the results of Cox PH model and shared gamma frailty model showed there are differences in the estimates of covariate effects on magnitude and significances. Almost all coefficient's magnitudes are higher in shared gamma frailty model. Age of mother at birth was significant in Cox PH model but not significant in shared gamma frailty model. These kind of differences was also observed in a study by Flinn and Heckman (1982).
This study showed that the risk of dying for a male twin child is higher than a female twin child. In agreement with this result, several studies on child mortality showed that the risk of dying for a male child is higher than a female child (Bereka, Habtewold, & Nebi, 2017;Samuel, 2011;Tibebeu, 2011) .
This study result revealed that the effect of being first born children of twin birth is insignificant. Similar to this result, a study in Burkina Faso showed that the variable is statistically insignificant (Becher, Muller, Jahn, Gbangou, Kynast-Wolf, & Kouyate, 2004). However, a child mortality study showed that first born children experience low survival compared to those children who have higher birth order (Balk, Tom, Adam, Fern, & Melissa, 2003). Although to the best of the researchers knowledge not seen before, this study showed that that the risk of dying for a first born twin child is less likely than a second born twin child.
Unlike a study that showed children residing in urban areas have a better chance of survival than those residing in rural areas (Balk, Tom, Adam, Fern, & Melissa, 2003), this study found that the variable residence is not significant. In addition, disagreement with studies that showed children born from women at youngest and oldest age are subject to high risk of death (Guo & Rodriguez, 1992;Sastry, 1997), this study found the variable is not significantly associated with survival of under-5 twins.
The current study also showed that children of twin birth who were born before the older sibling reaches the age of 18 months experience low survival compared to those children of twin birth who have more than 24 month spacing. Several studies that are conducted on determinants of child mortality found similar results (Becher, Muller, Jahn, Gbangou, Kynast-Wolf, & Kouyate, 2004;Koissi & Hgns, 2001;Guo & Rodriguez, 1992).
The study also showed that the risk of dying for a child of twin birth who has a younger sibling that born within 18 months is higher than a child of twin birth who has a younger sibling that born after 24 months. Several studies that are conducted on determinants of child mortality found similar results (Becher, Muller, Jahn, Gbangou, Kynast-Wolf, & Kouyate, 2004;Koissi & Hgns, 2001;Guo & Rodriguez, 1992).

Conclusions and recommendations 5.1. Conclusion
The results of this study showed that the main factors associated with twin under-five mortality are sex of child, between twin birth order, preceding birth interval, and succeeding birth interval. Male twins have lower chance of survival compared to female twin children. The risk of dying for a first born twin child is less likely than a second born twin child. Those twins who were born before their sibling reaches 18 month experience low survival compare to those twins who have more than 24 month spacing. And those twins who had a younger sibling with age less than 18 months experience low survival compare to those twins who have a younger sibling with age more than 24 months. Since the parameter of the heterogeneity variable is significant, one can conclude that there is a random effect shared by twins in a pair and one cannot ignore the correlation between the pair of twins..

Recommendations
In order to reduce the rate of child mortality among twins, this study recommends the following: The significant effect of the birth spacing of the previous and succeeding sibling on the survival chance of children of twin birth indicate that efforts have to be exerted to educate the public about family planning and birth spacing.
The significant effect of sex of child and between twin birth order may indicate that it is important to conduct further studies perhaps biological studies; in order to explain why and how these variables are significant.
The modeling practice in this study revealed the significance of including a random effect in the model. Hence, researches in this area should incorporate and assess such random effects. Funding: The authors received no direct funding for this research.