##Martin et al. Birth mode morbidity Qom

##Accompanying dataset is pre-processed, including dataframe merges, data cleaning, and imputed categorical and missing data. 


###(1) Import file and load required packages, 
###(2) Mixed-models of cesarean on breastfeeding and WAZ (Supplementary Tables 1-3, Figure 1)
###(3) Mixed-effects models of GI and RI risk (Table 3)

####(1)#####

bmorb <- read.csv("./QomBabiesMorbidity_share.csv") #set wd and import .csv file

#load packages
require(lme4)
require(lmerTest)
require(ggplot2)
require(sjPlot)

##pre-processed variable names/descriptions
---------------
SubID= child identifier; date= date of interview; age.yrs/.wks/.days = child age in years, weeks, days; 
bage.grp = child age group (0-5 months, 6-11, 12-23, 24-35, 35+ months);
gi= gi symptoms observed (1=yes, 0=no), see Supplementary Text for description of symptoms);
ri= ri symptoms observed (1=yes, 0=no);
weight= child weight (kg); height = child length/height (cm);
#WHO igrowup computed variables: cbmi = BMI; zlen = length-for-age z-score (LAZ/HAZ); zwei = WAZ; zwfl = WHZ; zbmi= BMIZ;
dob= date of birth; sex= male/female; birthweight = birthweight (kg); gestage= gestational age (weeks);
cscet = cesarean (1=yes, 0=no); mom.age= birth mother age (years); mom.age.grp = mother's age group (12-19, 20-29, 30-45 years);
primipar2= primiparous (1=yes, 0=no); bf = breastfeeding; ebf = exclusively breastfeeding (yes/no);
bfstatus= bf, ebf, not bf; respondent = caregiver interviewed; ObsNo= sequential observation number/subject; 
comments = text comments noted by interviewer from original data files
------------------

#note summary statistics performed previously on a different, non-pre-processed dataset from all subjects with available birth data (n=91)
#summary statistics performed using R base commands (e.g. summary()) and summary functions in "pastecs" and "doBy"
#univariate analysis performed with chisq.test() and glm() functions in base package

####(2)cesarean-WAZ and breastfeeding models####

##variable pre-processing

#recategorize baby age group
bmorb$bage.grp2 <- NA
bmorb$bage.grp2[bmorb$bage.grp=="0-5 months"] <- "0-5 months"
bmorb$bage.grp2[bmorb$bage.grp=="6-11 months"] <- "6-11 months"
bmorb$bage.grp2[bmorb$bage.grp=="12-23 months"] <- "12-23 months"
bmorb$bage.grp2[bmorb$bage.grp=="24-36 months"] <- "24+ months"
bmorb$bage.grp2[bmorb$bage.grp=="36+ months"] <- "24+ months"
bmorb$bage.grp2 <- as.factor(bmorb$bage.grp2)

bmorb$wean <- as.factor(ifelse(bmorb$bf=="yes",0,1)) #Relabel as numeric category

bmorb$mom.age.grp2 <- relevel(bmorb$mom.age.grp, ref="20-29") #relevel reference category

bmorb$pre.early <- as.factor(ifelse(bmorb$gestage < 39,1,0)) #bin into < $ >= 39 weeks 

#sort by subject and date 
bmorb <- bmorb[order(bmorb$SubID, bmorb$date),]


####(2a)mixed-effects linear regression of cesarean and WAZ (Supplementary Table 1 Figure 1)####

#rename birthmodel for plot
bmorb$birthmode <- ifelse(bmorb$csect==0, "Vaginal", "Cesarean")
bmorb$birthmode <- factor(bmorb$birthmode, levels=c("Vaginal","Cesarean"))
summary(bmorb$birthmode)
summary(bmorb$csect)

#best fit model (compared to random intercept only, interaction significant and lowers AIC)
w1b <- lmer(zwei~age.yrs*birthmode + (age.yrs | SubID), data=bmorb, REML=FALSE)

summary(w1b) #summary data for Supp Table 1

cc <- confint(w1b, parm="beta_", method="Wald") 
ctab <- cbind(est=fixef(w1b),cc)
rtab <- as.data.frame(exp(ctab))

set_theme(
    base = theme_classic(),
    axis.title.size = 1.3,
    axis.textsize = 1.3,
    legend.size = 1.5,
    legend.title.size = 1.3,
    geom.label.size = 3,
    legend.inside=TRUE,
    legend.pos = "top"
)

plot_model(w1b, type="int", title="", #sjPlot function, takes a while
           axis.title = c("Age (years)", "WAZ"),
           axis.lable= "",
           colors=c("coral", "cadetblue"))

####(2b)mixed-effect logistic regression of cesarean and 'not breastfeeding'####

#best fit model (compared to random intercept only and nested age group/ID model)
b3 <- glmer(wean~age.yrs + csect  + (age.yrs | SubID),
            data=bmorb, family = binomial("logit"))

#summary data for Supplementary Table 2 and Table 3
summary(b3) #AIC = 734.3, no significant relationship with not breastfeeding (age*csect not significant, AIC = 735.8)

cc <- confint(b3, parm="beta_", method="Wald") 
ctab <- cbind(est=fixef(b3),cc)
rtab <- as.data.frame(exp(ctab))

tab <- xtabs(~wean + bage.grp2 + csect, data=bmorb)
ftable(tab)

####(3) GI/RI risk models#####


####(3a)GI####

##Model 1: baseline "current risk model" (first determine best treatment of child age)
g1 <- glmer(gi ~ wean + zwei + age.yrs + (1 | SubID),
              data=bmorb, family = binomial("logit"))
summary(g1)#AIC = 2335, weaning and WAZ reduce risk, sex and age NS (adding sex raises AIC to 2337)

g1b <- glmer(gi ~ wean + zwei + age.yrs + (age.yrs | SubID),
            data=bmorb, family = binomial("logit"))
summary(g1b) #AIC = 2318 (age*waz interaction not significant, AIC= 2320)

##best fit baseline model
g1c <- glmer(gi ~ wean + zwei  + bage.grp2 + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g1c) #AIC = 2287.8  (model does not converge with zwei*bage.grp2 interaction)

g1d <- glmer(gi ~ wean + zwei + (1 | bage.grp2 ) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g1d) #AIC 2298


## Model 2: "current + gestational risk model

g2a <- glmer(gi ~ wean + zwei + age.yrs + pre.early + mom.age + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g2a) #AIC = 2333.9, random slope model does not converge

g2b <- glmer(gi ~ wean + zwei + bage.grp2 + pre.early + mom.age + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g2b) #model did not converge

#best fitting baseline model w/additive gestational risk factors (compare to AIC 2298)
g2c <- glmer(gi ~ wean + zwei  + pre.early + mom.age.grp + (1 | bage.grp2) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g2c) #AIC = 2295.8 with mom.age #AIC = 2294.4 with mom.age.grp
#AIC 2294.6 with mom.grp and removing pre.early; AIC 2297.7 with pre.early and no mom 

cc <- confint(g1e, parm="beta_", method="Wald") 
ctab <- cbind(est=fixef(g1e),cc)
rtab <- as.data.frame(exp(ctab))

## Model 3 "current + gestational risk + birth model model"
g2e <- glmer(gi ~ wean + zwei + pre.early + mom.age.grp2 + csect + (1 | bage.grp2) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g2e) #AIC = 2296.3

##try all one-to-one interactions between cesarean and other model terms
g2f <- glmer(gi ~ wean + zwei + pre.early + mom.age.grp2*csect + (1 | bage.grp2) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g2f) # AIC = 2298 for wean*csect NS, waz*csect NS
#AIC 2295.8 for pre.early*csect, but NS
#AIC 2299 for mom.age.grp*csect

#extra
g2x <- glmer(gi ~ wean + zwei + csect + (1 | bage.grp2) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(g2x) #2298.7


####(3b) RI risk####

#baseline "current risk"
r1 <- glmer(ri ~ wean + zwei + age.yrs + (1 | SubID),
            data=bmorb, family = binomial("logit"))
summary(r1)#AIC = 2520.5, nothing signficiant

r1b <- glmer(ri ~ wean + zwei + age.yrs + (age.yrs | SubID),
             data=bmorb, family = binomial("logit"))
summary(r1b) #AIC = 2520.7, reduced odds with age

r1c <- glmer(ri ~ wean + zwei  + bage.grp2 + (1 | SubID), #best fit but does not converge in additive models)
             data=bmorb, family = binomial("logit"))
summary(r1c) #AIC = 2507.8, higher odds at 12-23 and 6-11 months

r1d <- glmer(ri ~ wean + zwei + (1 | bage.grp2 ) + (1 | SubID), #model used as comparative baseline
             data=bmorb, family = binomial("logit"))
summary(r1d) #AIC 2513.6

cc <- confint(r1d, parm="beta_", method="Wald") 
ctab <- cbind(est=fixef(r1d),cc)
rtab <- as.data.frame(exp(ctab))

## "current risk + gestational"

r2a <- glmer(ri ~ wean + zwei  + bage.grp2 + pre.early + mom.age.grp + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(r2a) #failed to converge

#best fit
r2b <- glmer(ri ~ wean + zwei  + pre.early + mom.age.grp2 + (1 | bage.grp2) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(r2b) #AIC 2513.7 with pre.early & mom.age.grp 2; pre.early reduces odds
#AIC for gestage + mom.age.grp = 2514.8
#AIC for pre.early only = 2511; AIC for mom.age.grp only = 2516.2

## "current risk + gestational + birth mode"

r2c <- glmer(ri ~ wean + zwei  + pre.early + mom.age.grp2 + csect + (1 | bage.grp2) + (1 | SubID),
             data=bmorb, family = binomial("logit"))
summary(r2c) #AIC 2513.9, all NS