#Cortical Scaling of the Neonatal Brain in Typical and Altered Development
# Alexandra F. Bonthrone 2024

library(ggplot2)
library(visreg)
library(dplyr)
library(knitr)
library(png)
library(gvlma)
library(png)
library(jtools)
library(interactions)
library(robustbase)
library(gvlma)
library(broom)
library(ggplot2)
library(data.table)
library(scales)
library(dplyr)
library(car)
library(insight)
library(lmerPerm)

######## import and organise data

data <- read.csv("Bonthrone_allometric_data.csv")
row.names(data) <- data$Baby_ID

# Creating our summary variables

processed.data <- setDT(as.data.frame(data))


processed.data[, log_Outer_Hull := log10(Total_Outer_Hull)]                                            
processed.data[, log_Mean_Corr_Thick := log10(Mean_Corr_Thick)]                                                           
processed.data[, log_Total_SA := log10(Total_SA)]    

processed.data[, supratentorial_volume := Cortical_grey_matter + white_matter + deep_grey_matter + Hippocampus_and_amygdala]

processed.data[, log_supratentorial_volume := log10(supratentorial_volume)]

processed.data[, wholebrain_GI := Total_SA/Total_Outer_Hull]

processed.data[, log_SA_by_thick0.5 := log10(Total_SA*(Mean_Corr_Thick^0.5))]    

processed.data[, wholebrain_offset_term:= log_Total_SA - (1.25*log_Outer_Hull)+ (0.25*log(Mean_Corr_Thick^2))]

processed.data[, wholebrain_isometric_term := log_Total_SA + log_Outer_Hull + log(Mean_Corr_Thick^2)]

processed.data[, wholebrain_shape_term := (1.5*log_Total_SA) + (0.75*log_Outer_Hull) - (2.25*log(Mean_Corr_Thick^2))]


processed.data[, wholebrain_offset_term_Z:= (wholebrain_offset_term - mean(subset(processed.data, Group==0)$wholebrain_offset_term, na.rm=TRUE))/sd(subset(processed.data, Group==0)$wholebrain_offset_term, na.rm=TRUE)]

processed.data[, wholebrain_isometric_term_Z := (wholebrain_isometric_term - mean(subset(processed.data, Group==0)$wholebrain_isometric_term, na.rm=TRUE))/sd(subset(processed.data, Group==0)$wholebrain_isometric_term, na.rm=TRUE)]

processed.data[, wholebrain_shape_term_Z := (wholebrain_shape_term - mean(subset(processed.data, Group==0)$wholebrain_shape_term, na.rm=TRUE))/sd(subset(processed.data, Group==0)$wholebrain_shape_term, na.rm=TRUE)]
 
 
processed.data <- as.data.frame(processed.data)

#### organising data

CHD_preterms <- subset(processed.data, Group==1 | Group ==2)


processed.data$Group <- as.factor(processed.data$Group)
CHD_preterms$Group <- as.factor(CHD_preterms$Group)

#### Getting coefficients for each group

term.SA <-(lmrob(log_Total_SA ~ PMA + I(PMA^2)+male+Multiple+log_supratentorial_volume, data=subset(processed.data, Group==0), settings="KS2014"))
preterm.SA <-(lm(log_Total_SA ~ PMA + I(PMA^2)+male+Multiple+log_supratentorial_volume, data=subset(processed.data, Group==1)))
CHD.SA <-(lm(log_Total_SA ~ PMA + I(PMA^2)+male+Multiple+log_supratentorial_volume, data=subset(processed.data, Group==2)))

term.Thick <-(lm(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+male+Multiple+log_supratentorial_volume, data=subset(processed.data, Group==0)))
preterm.Thick <-(lmrob(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+male+Multiple+log_supratentorial_volume, data=subset(processed.data, Group==1), settings="KS2014"))
CHD.Thick <-(lmrob(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+male+Multiple+log_supratentorial_volume, data=subset(processed.data, Group==2),settings="KS2014"))

term.fold <-(lmrob(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+male+Multiple+log_Outer_Hull, data=subset(processed.data, Group==0),setting="KS2014"))
preterm.fold <-(lmrob(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+male+Multiple+log_Outer_Hull, data=subset(processed.data, Group==1),setting="KS2014"))
CHD.fold <-(lm(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+male+Multiple+log_Outer_Hull, data=subset(processed.data, Group==2)))


linearHypothesis(term.SA, c("log_supratentorial_volume = 0.66667"))
linearHypothesis(preterm.SA, c("log_supratentorial_volume = 0.66667"))
linearHypothesis(CHD.SA, c("log_supratentorial_volume = 0.66667"))


linearHypothesis(term.Thick, c("log_supratentorial_volume = 0.33333"))
linearHypothesis(preterm.Thick, c("log_supratentorial_volume = 0.33333"))
linearHypothesis(CHD.Thick, c("log_supratentorial_volume = 0.33333"))

linearHypothesis(term.fold, c("log_Outer_Hull = 1.25"))
linearHypothesis(preterm.fold, c("log_Outer_Hull = 1.25"))
linearHypothesis(CHD.fold, c("log_Outer_Hull = 1.25"))




### Controls effect of demographics

thick.sex<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2) +Multiple+ GA+Bwzscore_UK90+male*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

thick.PMA<-lmerp(log_Mean_Corr_Thick ~ I(PMA^2) + male +Multiple+ GA+Bwzscore_UK90+PMA*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

thick.PMA2<-lmerp(log_Mean_Corr_Thick ~ PMA + male + GA+Multiple+Bwzscore_UK90+I(PMA^2)*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

thick.GA<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2) + male+Bwzscore_UK90 + Multiple+ GA*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

thick.multiple<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Bwzscore_UK90 + male + GA+Multiple*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

term.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Multiple+male+GA+log_supratentorial_volume*Bwzscore_UK90,subset(processed.data, Group==0), R=5000, minimum=5000)



SA.sex<-lmerp(log_Total_SA ~ PMA + I(PMA^2) + GA +Multiple+Bwzscore_UK90+ male*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

SA.PMA<-lmerp(log_Total_SA ~ I(PMA^2) + male + GA +Multiple+Bwzscore_UK90+ PMA*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

SA.PMA2<-lmerp(log_Total_SA ~ PMA + male + GA+Multiple+Bwzscore_UK90+ I(PMA^2)*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

SA.GA<-lmerp(log_Total_SA ~ PMA + I(PMA^2) + male+Bwzscore_UK90 +Multiple+ GA*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

SA.multiple<-lmerp(log_Total_SA ~ PMA + I(PMA^2) +Bwzscore_UK90+ male + GA+Multiple*log_supratentorial_volume,subset(processed.data, Group==0), R=5000, minimum=5000)

term.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Multiple+ GA+male+log_supratentorial_volume*Bwzscore_UK90,subset(processed.data, Group==0), R=5000, minimum=5000)


log_SA_by_thick0.5.sex<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ GA+male*log_Outer_Hull,subset(processed.data, Group==0), R=5000, minimum=5000)

log_SA_by_thick0.5.PMA<-lmerp(log_SA_by_thick0.5~ I(PMA^2)+Bwzscore_UK90 + male +Multiple+ GA+PMA*log_Outer_Hull,subset(processed.data, Group==0), R=5000, minimum=5000)

log_SA_by_thick0.5.PMA2<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + male +Multiple+ GA+I(PMA^2)*log_Outer_Hull,subset(processed.data, Group==0), R=5000, minimum=5000)

log_SA_by_thick0.5.GA<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ male + GA*log_Outer_Hull,subset(processed.data, Group==0), R=5000, minimum=5000)

log_SA_by_thick0.5.Multiple<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90 + male + GA+Multiple*log_Outer_Hull,subset(processed.data, Group==0), R=5000, minimum=5000)

term.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Multiple+ GA+male+log_Outer_Hull*Bwzscore_UK90,subset(processed.data, Group==0), R=5000, minimum=5000)


term.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90,subset(processed.data, Group==0), R=5000, minimum=5000)

term.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90,subset(processed.data, Group==0), R=5000, minimum=5000)

term.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2) + male +Multiple+ GA+Bwzscore_UK90,subset(processed.data, Group==0), R=5000, minimum=5000)



### preterm effect of demographics

prem.thick.sex<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2) +Multiple+ Bwzscore_UK90+GA+male*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.thick.PMA<-lmerp(log_Mean_Corr_Thick ~ I(PMA^2) + male +Multiple+Bwzscore_UK90+ GA+PMA*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.thick.PMA2<-lmerp(log_Mean_Corr_Thick ~ PMA + male +Bwzscore_UK90+ GA+Multiple+I(PMA^2)*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.thick.GA<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2) +Bwzscore_UK90+ male + Multiple+ GA*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.thick.multiple<-lmerp(log_Mean_Corr_Thick ~ PMA +Bwzscore_UK90+ I(PMA^2) + male + GA+Multiple*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Multiple+male+GA+log_supratentorial_volume*Bwzscore_UK90,subset(processed.data, Group==1), R=5000, minimum=5000)


prem.SA.sex<-lmerp(log_Total_SA ~ PMA + I(PMA^2) + GA+Bwzscore_UK90 +Multiple+ male*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.SA.PMA<-lmerp(log_Total_SA ~ I(PMA^2) + male + GA+Bwzscore_UK90 +Multiple+ PMA*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.SA.PMA2<-lmerp(log_Total_SA ~ PMA + male+Bwzscore_UK90 + GA+Multiple+ I(PMA^2)*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.SA.GA<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Bwzscore_UK90 + male +Multiple+ GA*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.SA.multiple<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Bwzscore_UK90 + male + GA+Multiple*log_supratentorial_volume,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Multiple+ GA+male+log_supratentorial_volume*Bwzscore_UK90,subset(processed.data, Group==1), R=5000, minimum=5000)


prem.log_SA_by_thick0.5.sex<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ GA+male*log_Outer_Hull,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.log_SA_by_thick0.5.PMA<-lmerp(log_SA_by_thick0.5~ I(PMA^2)+Bwzscore_UK90 + male +Multiple+ GA+PMA*log_Outer_Hull,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.log_SA_by_thick0.5.PMA2<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + male +Multiple+ GA+I(PMA^2)*log_Outer_Hull,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.log_SA_by_thick0.5.GA<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + I(PMA^2) +Multiple+ male + GA*log_Outer_Hull,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.log_SA_by_thick0.5.Multiple<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + I(PMA^2) + male + GA+Multiple*log_Outer_Hull,subset(processed.data, Group==1), R=5000, minimum=5000)

prem.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Multiple+ GA+male+log_Outer_Hull*Bwzscore_UK90,subset(processed.data, Group==1), R=5000, minimum=5000)


preterm.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90,subset(processed.data, Group==1), R=5000, minimum=5000)

preterm.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90,subset(processed.data, Group==1), R=5000, minimum=5000)

preterm.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2) + male +Multiple+ GA+Bwzscore_UK90,subset(processed.data, Group==1), R=5000, minimum=5000)


### CHD effect of demographics

CHD.thick.sex<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2) +Multiple+Bwzscore_UK90+ GA+male*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.thick.PMA<-lmerp(log_Mean_Corr_Thick ~ I(PMA^2) + male +Multiple+Bwzscore_UK90+ GA+PMA*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.thick.PMA2<-lmerp(log_Mean_Corr_Thick ~ PMA + male + GA+Multiple+Bwzscore_UK90+I(PMA^2)*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.thick.GA<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2) + male +Bwzscore_UK90+ Multiple+ GA*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.thick.multiple<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Bwzscore_UK90 + male + GA+Multiple*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Multiple+male+GA+log_supratentorial_volume*Bwzscore_UK90,subset(processed.data, Group==2), R=5000, minimum=5000)


CHD.SA.sex<-lmerp(log_Total_SA ~ PMA + I(PMA^2) + GA +Multiple+Bwzscore_UK90+ male*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.SA.PMA<-lmerp(log_Total_SA ~ I(PMA^2) + male + GA +Multiple+Bwzscore_UK90+ PMA*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.SA.PMA2<-lmerp(log_Total_SA ~ PMA + male + GA+Multiple+Bwzscore_UK90+ I(PMA^2)*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.SA.GA<-lmerp(log_Total_SA ~ PMA + I(PMA^2) + male +Multiple+Bwzscore_UK90+ GA*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.SA.multiple<-lmerp(log_Total_SA ~ PMA + I(PMA^2) + male+Bwzscore_UK90 + GA+Multiple*log_supratentorial_volume,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Multiple+ GA+male+log_supratentorial_volume*Bwzscore_UK90,subset(processed.data, Group==2), R=5000, minimum=5000)


CHD.log_SA_by_thick0.5.sex<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ GA+male*log_Outer_Hull,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.log_SA_by_thick0.5.PMA<-lmerp(log_SA_by_thick0.5~ I(PMA^2)+Bwzscore_UK90 + male +Multiple+ GA+PMA*log_Outer_Hull,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.log_SA_by_thick0.5.PMA2<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + male +Multiple+ GA+I(PMA^2)*log_Outer_Hull,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.log_SA_by_thick0.5.GA<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + I(PMA^2) +Multiple+ male + GA*log_Outer_Hull,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.log_SA_by_thick0.5.Multiple<-lmerp(log_SA_by_thick0.5 ~ PMA+Bwzscore_UK90 + I(PMA^2) + male + GA+Multiple*log_Outer_Hull,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Multiple+ GA+male+log_Outer_Hull*Bwzscore_UK90,subset(processed.data, Group==2), R=5000, minimum=5000)


CHD.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2) + male +Multiple+ GA+Bwzscore_UK90,subset(processed.data, Group==2), R=5000, minimum=5000)




### preterm risk factors

respsupport.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90+respsupportdays,subset(processed.data, Group==1), R=5000, minimum=5000)

respsupport.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90+respsupportdays,subset(processed.data, Group==1), R=5000, minimum=5000)

respsupport.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2) + male +Multiple+ GA+Bwzscore_UK90+respsupportdays,subset(processed.data, Group==1), R=5000, minimum=5000)

respsupport.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Multiple+male+GA+Bwzscore_UK90+log_supratentorial_volume*respsupportdays,subset(processed.data, Group==1), R=5000, minimum=5000)

respsupport.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Multiple+ GA+Bwzscore_UK90+male+log_supratentorial_volume*respsupportdays,subset(processed.data, Group==1), R=5000, minimum=5000)

respsupport.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Multiple+ GA+Bwzscore_UK90+male+log_Outer_Hull*respsupportdays,subset(processed.data, Group==1), R=5000, minimum=5000)



parenteralnutrition.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2) +Multiple+Bwzscore_UK90+ male + GA+parenteralnutritiondays,subset(processed.data, Group==1), R=5000, minimum=5000)

parenteralnutrition.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2) +Multiple+ male + GA+Bwzscore_UK90+parenteralnutritiondays,subset(processed.data, Group==1), R=5000, minimum=5000)

parenteralnutrition.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2) + male +Multiple+ GA+Bwzscore_UK90+parenteralnutritiondays,subset(processed.data, Group==1), R=5000, minimum=5000)

parenteralnutrition.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Multiple+male+GA+Bwzscore_UK90+log_supratentorial_volume*parenteralnutritiondays,subset(processed.data, Group==1), R=5000, minimum=5000)

parenteralnutrition.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Multiple+ GA+male+Bwzscore_UK90+log_supratentorial_volume*parenteralnutritiondays,subset(processed.data, Group==1), R=5000, minimum=5000)

parenteralnutrition.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Multiple+ GA+Bwzscore_UK90+male+log_Outer_Hull*parenteralnutritiondays,subset(processed.data, Group==1), R=5000, minimum=5000)


### CHD risk factors

CHD_group.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Bwzscore_UK90+Multiple+male+GA+log_supratentorial_volume*CHD_group,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD_group.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Bwzscore_UK90+Multiple+ GA+male+log_supratentorial_volume*CHD_group,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD_group.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90+Multiple+ GA+male+log_Outer_Hull*CHD_group,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD_group.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ male + GA+CHD_group,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD_group.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ male + GA+CHD_group,subset(processed.data, Group==2), R=5000, minimum=5000)

CHD_group.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2)+Bwzscore_UK90 + male +Multiple+ GA+CHD_group,subset(processed.data, Group==2), R=5000, minimum=5000)



CDO2.thick.bw<-lmerp(log_Mean_Corr_Thick ~ PMA + I(PMA^2)+Multiple+Bwzscore_UK90+male+GA+log_supratentorial_volume*CDO2,subset(processed.data, Group==2), R=5000, minimum=5000)

CDO2.SA.bw<-lmerp(log_Total_SA ~ PMA + I(PMA^2)+Multiple+Bwzscore_UK90+ GA+male+log_supratentorial_volume*CDO2,subset(processed.data, Group==2), R=5000, minimum=5000)

CDO2.log_SA_by_thick0.5.bw<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2)+Bwzscore_UK90+Multiple+ GA+male+log_Outer_Hull*CDO2,subset(processed.data, Group==2), R=5000, minimum=5000)

CDO2.wholebrain.offset <- lmerp(wholebrain_offset_term ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ male + GA+CDO2,subset(processed.data, Group==2), R=5000, minimum=5000)

CDO2.wholebrain.isometric <- lmerp(wholebrain_isometric_term  ~ PMA + I(PMA^2)+Bwzscore_UK90 +Multiple+ male + GA+CDO2,subset(processed.data, Group==2), R=5000, minimum=5000)

CDO2.wholebrain.shape <- lmerp(wholebrain_shape_term ~ PMA + I(PMA^2)+Bwzscore_UK90 + male +Multiple+ GA+CDO2,subset(processed.data, Group==2), R=5000, minimum=5000)



##Comparing scaling coefficients between groups


processed.data_young <- subset(processed.data, PMA <=42.29)


preterm_thick<-lmerp(log_Mean_Corr_Thick ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group)*log_supratentorial_volume, data=subset(processed.data, Group==0 | Group==1), R=5000, minimum=5000)


CHD_thick<-lmerp(log_Mean_Corr_Thick ~ PMA+I(PMA^2)+as.factor(Multiple)+GA+as.factor(male)+GA+Bwzscore_UK90+as.factor(Group)*log_supratentorial_volume, data=subset(processed.data_young, Group==0 | Group==2), R=5000, minimum=5000)


preterm_SA<-lmerp(log_Total_SA ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group)*log_supratentorial_volume, data=subset(processed.data, Group==0 | Group==1), R=5000, minimum=5000)


CHD_SA<-lmerp(log_Total_SA ~ PMA+I(PMA^2)+as.factor(Multiple)+GA+as.factor(male)+Bwzscore_UK90+as.factor(Group)*log_supratentorial_volume, data=subset(processed.data_young, Group==0 | Group==2), R=5000, minimum=5000)


preterm_hull<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2) + as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+log_Outer_Hull*as.factor(Group),subset(processed.data, Group == 0 | Group==1), R=5000, minimum=5000)

CHD_hull<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2) + as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+GA+log_Outer_Hull*as.factor(Group),subset(processed.data_young, Group ==0 | Group==2), R=5000, minimum=5000)


preterm_CHD_SA<-lmerp(log_Total_SA ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group)*log_supratentorial_volume, data=subset(processed.data, Group == 1 | Group==2), R=5000, minimum=5000)


preterm_CHD_thick<-lmerp(log_Mean_Corr_Thick ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group)*log_supratentorial_volume, data=subset(processed.data, Group == 1 | Group==2), R=5000, minimum=5000)

preterm_CHD_hull<-lmerp(log_SA_by_thick0.5 ~ PMA + I(PMA^2) + as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+log_Outer_Hull*as.factor(Group),subset(processed.data, Group == 1 | Group==2), R=5000, minimum=5000)


preterm_offset<-lmerp(wholebrain_offset_term ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data, Group==0 | Group==1), R=5000, minimum=5000)

preterm_isometric<-lmerp(wholebrain_isometric_term ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data, Group==0 | Group==1), R=5000, minimum=5000)

preterm_shape<-lmerp(wholebrain_shape_term ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data, Group==0 | Group==1), R=5000, minimum=5000)


CHD_offset<-lmerp(wholebrain_offset_term ~ PMA+I(PMA^2)+as.factor(Multiple)+GA+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data_young, Group==0 | Group==2), R=5000, minimum=5000)

CHD_isometric<-lmerp(wholebrain_isometric_term ~ PMA+I(PMA^2)+as.factor(Multiple)+GA+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data_young, Group==0 | Group==2), R=5000, minimum=5000)

CHD_shape<-lmerp(wholebrain_shape_term ~ PMA+I(PMA^2)+as.factor(Multiple)+GA+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data_young, Group==0 | Group==2), R=5000, minimum=5000)




CHD.preterm_offset<-lmerp(wholebrain_offset_term ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data, Group==2 | Group==1), R=5000, minimum=5000)

CHD.preterm_isometric<-lmerp(wholebrain_isometric_term ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data, Group==2 | Group==1), R=5000, minimum=5000)

CHD.preterm_shape<-lmerp(wholebrain_shape_term ~ PMA+I(PMA^2)+as.factor(Multiple)+as.factor(male)+Bwzscore_UK90+as.factor(Group), data=subset(processed.data, Group==2 | Group==1), R=5000, minimum=5000)




##comparing raw values between groups

for (i in c("Total_SA","Mean_Corr_Thick","supratentorial_volume","wholebrain_GI")) {
formula2 <- as.formula(sprintf("%s ~PMA+I(PMA^2)+male+Multiple+Bwzscore_UK90+Group", i))
result2 <- lmerp(formula=formula2,subset(processed.data, Group==0 | Group==1), R=5000, minimum=5000)
new_name2 <- paste0("raw_control_prem_",i)
print(new_name2)
assign(new_name2, result2)
}

for (i in c("Total_SA","Mean_Corr_Thick","supratentorial_volume","wholebrain_GI")) {
formula2 <- as.formula(sprintf("%s ~PMA+I(PMA^2)+male+GA+Multiple+Bwzscore_UK90+Group", i))
result2 <- lmerp(formula=formula2,subset(processed.data, Group==0 | Group==2), R=5000, minimum=5000)
new_name2 <- paste0("raw_control_CHD_",i)
print(new_name2)
assign(new_name2, result2)
}

for (i in c("Total_SA","Mean_Corr_Thick","supratentorial_volume","wholebrain_GI")) {
formula2 <- as.formula(sprintf("%s ~PMA+I(PMA^2)+male+Multiple+Bwzscore_UK90+Group", i))
result2 <- lmerp(formula=formula2,subset(processed.data, Group==1 | Group==2), R=5000, minimum=5000)
new_name2 <- paste0("raw_prem_CHD_",i)
print(new_name2)
assign(new_name2, result2)
}



#getting volumetric z-scores adjusting for covs- zscoring to control mean and SD

vol.variables <- c("log_Total_SA","log_Mean_Corr_Thick")

Vol.zscores_CHD_preterms<-matrix(, nrow = 180, ncol = 2)
dimnames(Vol.zscores_CHD_preterms) <- list(CHD_preterms$Baby_ID, vol.variables)

Vol.zscores_terms<-matrix(, nrow = 345, ncol = 2)
dimnames(Vol.zscores_terms) <- list(subset(processed.data, Group==0)$Baby_ID, vol.variables)


for (i in vol.variables) {
keeps <- c(i,"PMA","male","log_supratentorial_volume","Multiple","Bwzscore_UK90")
processed.data <- subset(processed.data, Group==0)[keeps]
processed.data2 <- CHD_preterms[keeps]
print(colnames(processed.data)) 
fit <- lmrob(processed.data[,1] ~ log_supratentorial_volume+PMA+I(PMA^2)+as.factor(male)+Bwzscore_UK90+as.factor(Multiple), settings="KS2014", data=processed.data,k.max=500)
preds <- predict(fit, newdata=processed.data2, type="response",se=TRUE)
Vol.zscores_CHD_preterms[,i] <- ((preds$fit-processed.data2[,1])-mean(fit$resid))/sd(fit$resid)
Vol.zscores_terms[,i] <- (fit$resid-mean(fit$resid))/sd(fit$resid)
}

hull.variables <- c("log_SA_by_thick0.5")

Hull.zscores_CHD_preterms<-matrix(, nrow = 180, ncol = 1)
dimnames(Hull.zscores_CHD_preterms) <- list(CHD_preterms$Baby_ID, hull.variables)

Hull.zscores_terms<-matrix(, nrow = 345, ncol = 1)
dimnames(Hull.zscores_terms) <- list(subset(processed.data, Group==0)$Baby_ID, hull.variables)

for (i in hull.variables) {
keeps <- c(i,"PMA","male","log_Outer_Hull","Multiple","Bwzscore_UK90")
processed.data <- subset(processed.data, Group==0)[keeps]
processed.data2 <- CHD_preterms[keeps]
print(colnames(processed.data)) 
fit <- lmrob(processed.data[,1]~log_Outer_Hull+PMA+I(PMA^2)+as.factor(male)+as.factor(Multiple)+Bwzscore_UK90, settings="KS2014", data=processed.data, k.max=1000)
preds <- predict(fit, newdata=processed.data2, type="response",se=TRUE)
Hull.zscores_CHD_preterms[,i] <- ((preds$fit-processed.data2[,1])-mean(fit$resid))/sd(fit$resid)
Hull.zscores_terms[,i] <- (fit$resid-mean(fit$resid))/sd(fit$resid)
}

offset.variables <- c("wholebrain_offset_term","wholebrain_isometric_term","wholebrain_shape_term")


offset.zscores_CHD_preterms1<-matrix(, nrow = 180, ncol = 3)
dimnames(offset.zscores_CHD_preterms1) <- list(CHD_preterms$Baby_ID, offset.variables)

offset.zscores_terms1<-matrix(, nrow = 345, ncol =3)
dimnames(offset.zscores_terms1) <- list(subset(processed.data, Group==0)$Baby_ID, offset.variables)


for (i in offset.variables) {
keeps <- c(i,"PMA","male","log_Outer_Hull","Multiple","Bwzscore_UK90")
processed.data <- subset(processed.data, Group==0)[keeps]
processed.data2 <- CHD_preterms[keeps]
print(colnames(processed.data)) 
fit <- lmrob(processed.data[,1]~PMA+I(PMA^2)+as.factor(male)+as.factor(Multiple)+Bwzscore_UK90, settings="KS2014", data=processed.data,k.max=500)
preds <- predict(fit, newdata=processed.data2, type="response",se=TRUE)
offset.zscores_CHD_preterms1[,i] <- ((preds$fit-processed.data2[,1])-mean(fit$resid))/sd(fit$resid)
offset.zscores_terms1[,i] <- (fit$resid-mean(fit$resid))/sd(fit$resid)
}

zscore.data2 <- cbind(Vol.zscores_CHD_preterms, Hull.zscores_CHD_preterms,offset.zscores_CHD_preterms1)
zscore.data3 <- cbind(Vol.zscores_terms, Hull.zscores_terms,offset.zscores_terms1)
zscore.data <- rbind(zscore.data2, zscore.data3)
zscore.data <- as.data.frame(zscore.data)

rownames(processed.data) <- processed.data$Baby_ID

results <- merge(zscore.data, zscore.data4, by = 'row.names', all = TRUE,incomparables = NA) 
rownames(results) <- results$Row.names

processed.data_zadjustingforcovs <- merge(processed.data, results,
                          by = 'row.names', all = TRUE,incomparables = NA) 

processed.data_zadjustingforcovs.bw <- merge(processed.data, zscore.data,
                          by = 'Row.names', all = TRUE,incomparables = NA) 

rownames(processed.data_zadjustingforcovs)<- processed.data_zadjustingforcovs$Row.names


write.csv(processed.data_zadjustingforcovs.bw, "processed.data_zscores_adjustingforcovs_bw.csv")



#Reloading Data
processed.data_zadjustingforcovs_final <- read.csv("processed.data_zscores_adjustingforcovs_bw.csv")
row.names(processed.data_zadjustingforcovs) <- processed.data_zadjustingforcovs$Baby_ID

processed.data_zadjustingforcovs_final$Group <- as.factor(processed.data_zadjustingforcovs_final$Group)

processed.data_zadjustingforcovs$CHD_group <- relevel(processed.data_zadjustingforcovs$CHD_group, ref = "Right")

deviance_scores<-c("log_Total_SA.y", "log_Mean_Corr_Thick.y", "log_SA_by_thick0.5.y", "wholebrain_offset_term.y","wholebrain_isometric_term.y","wholebrain_shape_term.y")


#Outcome Interactions

for (i in deviance_scores){
formula2 <- as.formula(sprintf("cog_comp ~IMD_rank+%s*Group", i))
result2 <- lmerp(formula=formula2,processed.data_zadjustingforcovs.bw, R=5000, minimum=5000)
new_name2 <- paste0("cog_comp_interact_",i)
print(new_name2)
assign(new_name2, result2)
formula3 <- as.formula(sprintf("language_comp ~IMD_rank+parent_firstlanguage_notenglish+%s*Group", i))
result3 <- lmerp(formula=formula3,processed.data_zadjustingforcovs.bw, R=5000, minimum=5000)
new_name3 <- paste0("language_comp_interact2_",i)
print(new_name3)
assign(new_name3, result3)
formula5 <- as.formula(sprintf("motor_comp ~IMD_rank+%s*Group", i))
result5 <- lmerp(formula=formula5,processed.data_zadjustingforcovs.bw, R=5000, minimum=5000)
new_name5 <- paste0("motor_comp_interact2_",i)
print(new_name5)
assign(new_name5, result5)
}

formula <- cog_comp~IMD_rank+wholebrain_shape_term.y
rsq_function <- function(formula, data, indices) {
  d <- data[indices,] #allows boot to select sample
  fit <- lm(formula, data=d)
  return(summary(fit)$r.square)
}
reps <- boot(data=subset(main_analysis, Group==1), statistic=rsq_function, R=5000, formula=formula)


### outcomes and raw metrics


for (i in c("Total_SA","Mean_Corr_Thick","supratentorial_volume","wholebrain_GI")){
formula <- as.formula(sprintf("cog_comp ~IMD_rank+PMA+I(PMA^2)+male+Multiple+Bwzscore_UK90+GA+%s", i))
formula2 <- as.formula(sprintf("language_comp ~ IMD_rank+PMA+I(PMA^2)+male+Multiple+Bwzscore_UK90+GA+parent_firstlanguage_notenglish+%s", i))
formula3 <- as.formula(sprintf("motor_comp ~IMD_rank+PMA+I(PMA^2)+male+Multiple+Bwzscore_UK90+GA+%s", i))
result <- lmerp(formula=formula, subset(processed.data, Group==0), R=5000, minimum=5000)
new_name <- paste0("cog_comp_controls_",i)
print(new_name)
assign(new_name, result)
result2 <- lmerp(formula=formula, subset(processed.data, Group==1), R=5000, minimum=5000)
new_name2 <- paste0("cog_comp_preterm_",i)
print(new_name2)
assign(new_name2, result2)
result3 <- lmerp(formula=formula, subset(processed.data, Group==2), R=5000, minimum=5000)
new_name3 <- paste0("cog_comp_CHD_",i)
print(new_name3)
assign(new_name3, result3)

result4 <- lmerp(formula=formula2, subset(processed.data, Group==0), R=5000, minimum=5000)
new_name4 <- paste0("lang_comp_controls_",i)
print(new_name4)
assign(new_name4, result4)
result5 <- lmerp(formula=formula2, subset(processed.data, Group==1), R=5000, minimum=5000)
new_name5 <- paste0("lang_comp_preterm_",i)
print(new_name5)
assign(new_name5, result5)
result6 <- lmerp(formula=formula2, subset(processed.data, Group==2), R=5000, minimum=5000)
new_name6 <- paste0("lang_comp_CHD_",i)
print(new_name6)
assign(new_name6, result6)

result7 <- lmerp(formula=formula3, subset(processed.data, Group==0), R=5000, minimum=5000)
new_name7 <- paste0("motor_comp_controls_",i)
print(new_name7)
assign(new_name7, result7)
result8 <- lmerp(formula=formula3, subset(processed.data, Group==1), R=5000, minimum=5000)
new_name8 <- paste0("motor_comp_preterm_",i)
print(new_name8)
assign(new_name8, result8)
result9 <- lmerp(formula=formula3, subset(processed.data, Group==2), R=5000, minimum=5000)
new_name9 <- paste0("motor_comp_CHD_",i)
print(new_name9)
assign(new_name9, result9)
}


#Histogram Figures

SA_hist <- ggplot(processed.data, aes(Total_SA, colour = Group, fill=Group)) + geom_histogram(alpha=0.5,position="identity") + xlab("Total Surface Area") + ylab("Number of Infants") + theme_bw() + scale_colour_manual(values = c("#CC508C","#56B4E9", "#009E8C"),labels = c("Control", "Preterm", "CHD")) + scale_fill_manual(values = c("#CC508C","#56B4E9", "#009E8C"),labels = c("Control", "Preterm", "CHD")) + theme(panel.border = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"),legend.position="right")


CT_hist <- ggplot(processed.data, aes(Mean_Corr_Thick, colour = Group, fill=Group)) + geom_histogram(alpha=0.5,position="identity") + xlab("Cortical Thickness") + ylab("Number of Infants") + theme_bw() + scale_colour_manual(values = c("#CC508C","#56B4E9", "#009E8C"),labels = c("Control", "Preterm", "CHD")) + scale_fill_manual(values = c("#CC508C","#56B4E9", "#009E8C"),labels = c("Control", "Preterm", "CHD")) + theme(panel.border = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"),legend.position="right")


Hull_hist <- ggplot(processed.data, aes(Total_Outer_Hull, colour = Group, fill=Group)) + geom_histogram(alpha=0.5,position="identity") + xlab("Outer Surface Area") + ylab("Number of infants") + theme_bw() + scale_colour_manual(values = c("#CC508C","#56B4E9", "#009E8C"),labels = c("Control", "Preterm", "CHD")) + scale_fill_manual(values = c("#CC508C","#56B4E9", "#009E8C"),labels = c("Control", "Preterm", "CHD")) + theme(panel.border = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"),legend.position="right")

png("SA_hist.png", units="cm", width=16, height=12, pointsize=12, res=1200)
print(SA_hist)
dev.off()

png("CT_hist.png", units="cm", width=16, height=12, pointsize=12, res=1200)
print(CT_hist)
dev.off()

png("Hull_hist.png", units="cm", width=16, height=12, pointsize=12, res=1200)
print(Hull_hist)
dev.off()