Exposure to placental insufficiency alters postnatal growth trajectory in extremely low birth weight infants

0. Packages and functions

packages <- c("tidyverse",
              "lme4",
              "lmerTest",
              "psych",
              "vcd",
              "ggpubr",
              "optimx",
              "effsize",
              "PropCIs")

for (i in packages) {
  if (!(i %in% installed.packages())) {
    install.packages(i)
  }
}

for (i in packages) {
  lapply(i, library, character.only = TRUE)
}

## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

## Registered S3 method overwritten by 'rvest':
##   method            from
##   read_xml.response xml2

## ?? Attaching packages ??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????? tidyverse 1.2.1 ??

## ? ggplot2 3.1.1       ? purrr   0.3.2  
## ? tibble  2.1.1       ? dplyr   https://protect-eu.mimecast.com/s/uw71CWP0PfzpqRJuzlzZs?domain=0.8.0.1
## ? tidyr   0.8.3       ? stringr 1.4.0  
## ? readr   1.3.1       ? forcats 0.4.0

## ?? Conflicts ?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????? tidyverse_conflicts() ??
## ? dplyr::filter() masks stats::filter()
## ? dplyr::lag()    masks stats::lag()

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:tidyr':
## 
##     expand

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

## Loading required package: grid

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

## 
## Attaching package: 'effsize'

## The following object is masked from 'package:psych':
## 
##     cohen.d

require(ggplot2)

# check and create folder
folder.check <- function(folder.name = NULL, current.dir) {
  setwd(current.dir)
  if (file.exists(folder.name)) {
    return(setwd(file.path(current.dir, folder.name)))
  } else {
    dir.create(file.path(current.dir, folder.name))
    return(setwd(file.path(current.dir, folder.name)))
  }
}

cat.cont.graph <- function(cat.table, cont.table, cat.name, cont.name) {
  temp.df <- data.frame(cat.name = cat.table[cat.name], cont.name = cont.table[cont.name])
  temp.df <- na.omit(temp.df)
  formula <- as.formula(paste0(cont.name, " ~ ", cat.name))
  m <- aggregate(formula, temp.df, mean)
  colnames(m)[colnames(m) %in% cont.name] = "mean"
  sd <- aggregate(temp.df[cont.name], temp.df[cat.name], sd)
  colnames(sd)[colnames(sd) %in% cont.name] = "sd"
  table <- merge(m, sd)

  if (nrow(table) == 2) {
  ttest <- t.test(cont.table[,cont.name]~cat.table[,cat.name])
  pval <- round(ttest$p.value, 4)
  stat <- round(ttest$statistic, 2)
  df <- round(ttest$parameter, 2)
  } else {
    pval <- ""
    stat <- ""
    df <- ""
  }
  k <- ggplot2::ggplot(table, aes_(as.name(cat.name), as.name("mean"), fill = as.name(cat.name))) + 
    geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), width = 0.2) +
    labs(x = cat.name, y = cont.name, 
         title = paste0(cont.name, " ~ ", cat.name), 
         subtitle = paste0("t[df=",df,"] = ", stat, ", p-value = ", pval)) +
    theme_bw(base_size = 16) +
    theme(axis.text=element_text(size=14),
          axis.title=element_text(size=16, face="bold"),
          plot.margin = margin(5,5,5,5, "mm"))
  
  return(k)
}

1. Data upload to R

require(readr)
continuous.df <- data.frame(readr::read_csv("continuous_variables.csv"))

## Parsed with column specification:
## cols(
##   id = col_double(),
##   calories_7d = col_double(),
##   aminoAcid_7d = col_double(),
##   protein_fullEF = col_double(),
##   calories_fullEF = col_double(),
##   Abx_days = col_double(),
##   birth_wt = col_double()
## )

categorical.df <- data.frame(readr::read_csv("categorical_variables.csv", col_types = cols(
  id = col_double(),
  pla_insuff = col_factor(levels = c("No", "Yes")),
  Gender = col_factor(levels = c("male", "female")),
  BPD = col_factor(levels = c("No", "Yes")),
  IVH = col_factor(levels = c("No", "Yes")),
  Sepsis = col_factor(levels = c("No", "Yes")),
  ROP = col_factor(levels = c("No", "Yes")),
  sga = col_factor(levels = c("No", "Yes"))
)))

## Warning: The following named parsers don't match the column names: sga

weight <- data.frame(readr::read_csv("weight.csv"))

## Parsed with column specification:
## cols(
##   id = col_double(),
##   time_week = col_double(),
##   wt_z_score = col_double()
## )

maternal_demographics.df <- data.frame(readr::read_csv("maternal_characteristics.csv"))

## Parsed with column specification:
## cols(
##   id = col_double(),
##   maternal_age = col_double(),
##   maternal_race = col_character(),
##   antenatal_steroid = col_double(),
##   delivery_method = col_double(),
##   maternal_obesity = col_double(),
##   maternal_smoking = col_double()
## )

neonatal_demographics.df <- data.frame(readr::read_csv("neonatal_characteristics.csv"))

## Parsed with column specification:
## cols(
##   id = col_double(),
##   gender = col_character(),
##   apgar_1min = col_double(),
##   apgar_5min = col_double(),
##   ga_birth = col_double(),
##   birth_weight = col_double(),
##   small_for_ga = col_double()
## )

#data summary
summary(continuous.df)

##        id          calories_7d     aminoAcid_7d   protein_fullEF 
##  Min.   :  1.00   Min.   :25.43   Min.   :1.540   Min.   :1.360  
##  1st Qu.: 27.50   1st Qu.:46.43   1st Qu.:2.950   1st Qu.:3.385  
##  Median : 50.00   Median :53.29   Median :3.100   Median :3.580  
##  Mean   : 50.65   Mean   :53.47   Mean   :3.056   Mean   :3.479  
##  3rd Qu.: 75.50   3rd Qu.:61.07   3rd Qu.:3.245   3rd Qu.:3.735  
##  Max.   :100.00   Max.   :77.00   Max.   :3.650   Max.   :4.330  
##  calories_fullEF    Abx_days        birth_wt      
##  Min.   : 72.4   Min.   : 0.00   Min.   :-2.3500  
##  1st Qu.:100.8   1st Qu.: 4.00   1st Qu.:-0.9950  
##  Median :104.3   Median : 9.00   Median :-0.3300  
##  Mean   :104.6   Mean   :11.63   Mean   :-0.4263  
##  3rd Qu.:109.8   3rd Qu.:17.00   3rd Qu.: 0.1350  
##  Max.   :120.6   Max.   :57.00   Max.   : 1.7800

summary(categorical.df)

##        id         pla_insuff    Gender    BPD      IVH     Sepsis  
##  Min.   :  1.00   No :57     male  :40   No :30   No :83   No :73  
##  1st Qu.: 27.50   Yes:34     female:51   Yes:61   Yes: 8   Yes:18  
##  Median : 50.00                                                    
##  Mean   : 50.65                                                    
##  3rd Qu.: 75.50                                                    
##  Max.   :100.00                                                    
##   ROP         SGA           
##  No :71   Length:91         
##  Yes:20   Class :character  
##           Mode  :character  
##                             
##                             
## 

summary(maternal_demographics.df)

##        id          maternal_age   maternal_race      antenatal_steroid
##  Min.   :  1.00   Min.   :15.00   Length:91          Min.   :0.000    
##  1st Qu.: 27.50   1st Qu.:23.00   Class :character   1st Qu.:1.000    
##  Median : 50.00   Median :27.00   Mode  :character   Median :1.000    
##  Mean   : 50.65   Mean   :27.93                      Mean   :0.978    
##  3rd Qu.: 75.50   3rd Qu.:32.00                      3rd Qu.:1.000    
##  Max.   :100.00   Max.   :43.00                      Max.   :1.000    
##  delivery_method  maternal_obesity maternal_smoking
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :1.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.7802   Mean   :0.0989   Mean   :0.3297  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000

summary(neonatal_demographics.df)

##        id            gender            apgar_1min      apgar_5min   
##  Min.   :  1.00   Length:91          Min.   :1.000   Min.   :2.000  
##  1st Qu.: 27.50   Class :character   1st Qu.:3.000   1st Qu.:6.000  
##  Median : 50.00   Mode  :character   Median :5.000   Median :7.000  
##  Mean   : 50.65                      Mean   :4.725   Mean   :6.857  
##  3rd Qu.: 75.50                      3rd Qu.:6.000   3rd Qu.:8.000  
##  Max.   :100.00                      Max.   :9.000   Max.   :9.000  
##     ga_birth      birth_weight      small_for_ga   
##  Min.   :23.14   Min.   :-2.3500   Min.   :0.0000  
##  1st Qu.:25.14   1st Qu.:-0.9950   1st Qu.:0.0000  
##  Median :26.29   Median :-0.3250   Median :0.0000  
##  Mean   :26.37   Mean   :-0.4261   Mean   :0.1868  
##  3rd Qu.:27.50   3rd Qu.: 0.1345   3rd Qu.:0.0000  
##  Max.   :31.57   Max.   : 1.7780   Max.   :1.0000

summary(weight)

##        id           time_week       wt_z_score    
##  Min.   :  1.00   Min.   :24.00   Min.   :-4.530  
##  1st Qu.: 28.00   1st Qu.:30.86   1st Qu.:-2.060  
##  Median : 51.00   Median :35.00   Median :-1.420  
##  Mean   : 50.73   Mean   :35.50   Mean   :-1.544  
##  3rd Qu.: 74.00   3rd Qu.:39.43   3rd Qu.:-0.900  
##  Max.   :100.00   Max.   :50.00   Max.   : 1.430

2. Correlation analysis

require(psych)
require(vcd)

# set up folder for graph storage
folder.check(folder.name = "Results", getwd())
folder.check(folder.name = "correlations", getwd())

# continuous variables
pdf(file = paste0(getwd(), "/continuous_correlation.pdf"), width = 24, height = 24)
psych::pairs.panels(continuous.df[-1], method = "spearman", bg = "yellow", pch = 21, stars = F)
dev.off()

## png 
##   2

#for html markdown
psych::pairs.panels(continuous.df[-1], method = "spearman", bg = "yellow", pch = 21, stars = F)

# categorical variables
pdf(file = paste0(getwd(), "/categorical_correlation.pdf"), width = 10, height = 7, onefile = TRUE)
categorical.df.temp <- categorical.df[-1]
for (i in 1:(ncol(categorical.df.temp)-1)) {
    for (j in (i+1):ncol(categorical.df.temp)) {
    d <- c(colnames(categorical.df.temp)[i], colnames(categorical.df.temp[j]))
    temp.df <- categorical.df.temp[, d]
    temp.df <- temp.df[complete.cases(temp.df),]
    temp.df.table <- table(temp.df)
    if(min(matrix(temp.df.table)) < 5) {
      fishert <- fisher.test(temp.df.table)
    }
    mosaic.plot <- vcd::mosaic(temp.df.table,
                          main = paste0(colnames(categorical.df.temp)[i], " vs ", colnames(categorical.df.temp)[j]),
                          sub = paste0("Fisher exact test p-value: ",ifelse(exists("fishert"), round(fishert$p.value, digits = 3),"")),
                          border = "black",
                          gp = shading_hsv, 
                          gp_args = list(interpolate = 1:4),
                          legend=TRUE
                          )
  }

}
dev.off()

## png 
##   2

#for html markdown
categorical.df.temp <- categorical.df[-1]
for (i in 1:(ncol(categorical.df.temp)-1)) {
    for (j in (i+1):ncol(categorical.df.temp)) {
    d <- c(colnames(categorical.df.temp)[i], colnames(categorical.df.temp[j]))
    temp.df <- categorical.df.temp[, d]
    temp.df <- temp.df[complete.cases(temp.df),]
    temp.df.table <- table(temp.df)
    if(min(matrix(temp.df.table)) < 5) {
      fishert <- fisher.test(temp.df.table)
    }
    mosaic.plot <- vcd::mosaic(temp.df.table,
                          main = paste0(colnames(categorical.df.temp)[i], " vs ", colnames(categorical.df.temp)[j]),
                          sub = paste0("Fisher exact test p-value: ",ifelse(exists("fishert"), round(fishert$p.value, digits = 3),"")),
                          border = "black",
                          gp = shading_hsv, 
                          gp_args = list(interpolate = 1:4),
                          legend=TRUE
                          )
  }

}

#categorical vs. continuous variables
categorical.df.temp <- categorical.df[-1]
continuous.df.temp <- continuous.df[-1]
cat.col <- colnames(categorical.df.temp)
cont.col <- colnames(continuous.df.temp)

pdf(file = paste0(getwd(), "/continuous_categorical_correlation.pdf"), width = 18, height = 10, onefile = TRUE)
for (i in cat.col) {
  p <- list()
  for (j in cont.col) {
    p[[j]] <- cat.cont.graph(categorical.df.temp, continuous.df.temp, i, j)
  }

plot <- ggpubr::ggarrange(plotlist = p, ncol = 3, nrow = 2)
plot <- ggpubr::annotate_figure(plot,
                        top = text_grob(paste("Categorical variable:", i), color = "maroon", size = 32)
                        )
print(plot)
}
dev.off()

## png 
##   2

#for html markdown
print(p)

## $calories_7d

## 
## $aminoAcid_7d

## 
## $protein_fullEF

## 
## $calories_fullEF

## 
## $Abx_days

## 
## $birth_wt

3. Spaghetti plot

require(tidyverse)
require(ggpubr)

# retrieve fixed-effects variable names
categorical.fixed.variables <- colnames(categorical.df)[-1]
long.df <- data.frame()
long.df <- inner_join(weight, categorical.df, by = "id")
long.df <- inner_join(long.df, continuous.df, by = "id")
long.df <- inner_join(long.df, maternal_demographics.df, by = "id")
long.df <- inner_join(long.df, neonatal_demographics.df, by = "id")
long.df[["id"]] <- as.factor(long.df[["id"]])

# set up folder for plot storage
folder.check(folder.name = "Results", getwd())
folder.check(folder.name = "spaghetti_plot", getwd())

pdf(file = paste0(getwd(), "/spaghetti plots.pdf"), width = 18, height = 18, onefile = TRUE)

p <- list()

for (i in categorical.fixed.variables) {

  p[[i]] <- ggplot2::ggplot(long.df, aes(x = time_week, y = wt_z_score, group = id)) +
    geom_line(aes_(color = as.name(i))) +
    theme_bw() +
    theme(legend.position = "bottom", 
          legend.title = element_text(size=12),
          legend.text = element_text(size=12),
          axis.text=element_text(size=12),
          axis.title=element_text(size=14, face="bold", margin = margin(10,10,10,10,"mm")),
          plot.margin = margin(5,5,5,5, "mm")) +
    labs(x = "week", y = "weight z-score", 
         title = "Trend of weight z-score",
         subtitle = paste0("Color-coded by ", i))
}

plot <- ggpubr::ggarrange(plotlist = p, ncol = 3, nrow =3)
plot <- ggpubr::annotate_figure(plot,
                                top = text_grob("Spaghetti Plots", color = "maroon", size = 32))
print(plot)
dev.off()

## png 
##   2

#for html markdown
print(p)

## $pla_insuff

## 
## $Gender

## 
## $BPD

## 
## $IVH

## 
## $Sepsis

## 
## $ROP

## 
## $SGA

4. Grid search for time_week shift constant

GA <- unique(weight$time_week)
c.grid <- seq(min(GA), max(GA), 0.1)
cor.terms <- rep(NA, length(c.grid))

for (i in 1:length(c.grid)) {
  cor.terms[i] <- cor((GA-c.grid[i]), (GA-c.grid[i])^2)
}
shift.constant <- c.grid[which.min(abs(cor.terms))]
cat("Shift time_week by", shift.constant, "will minimize correlation between (time_week) and (time_week)^2")

## Shift time_week by 37.3 will minimize correlation between (time_week) and (time_week)^2

cat("\nPearson r before shift:", cor(GA, GA^2))

## 
## Pearson r before shift: 0.9960433

cat("\nPearson r after shift:", cor((GA-shift.constant), (GA-shift.constant)^2))

## 
## Pearson r after shift: 0.009998993

5. Mixed-effects modeling with lme4 without considering covariates

require(lme4)
require(lmerTest)
require(optimx)
require(kableExtra)

## Loading required package: kableExtra

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

long.df.scale <- long.df
long.df.scale[,"time_week"] <- as.numeric(long.df.scale[,"time_week"])
long.df.scale[,"time_week"] <- long.df.scale[,"time_week"] - shift.constant

formula_list <- list()

formula_list[[1]] <- "wt_z_score ~ 1 + time_week + (time_week | id)"
formula_list[[2]] <- "wt_z_score ~ 1 + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[3]] <- "wt_z_score ~ pla_insuff + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[4]] <- "wt_z_score ~ pla_insuff + pla_insuff:time_week + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[5]] <- "wt_z_score ~ pla_insuff + pla_insuff:I(time_week^2) + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[6]] <- "wt_z_score ~ pla_insuff + pla_insuff:time_week + pla_insuff:I(time_week^2) + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[7]] <- "wt_z_score ~ birth_weight + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[8]] <- "wt_z_score ~ birth_weight + birth_weight:time_week + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[9]] <- "wt_z_score ~ birth_weight + birth_weight:I(time_week^2) + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[10]] <- "wt_z_score ~ birth_weight + birth_weight:time_week + birth_weight:I(time_week^2) + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[11]] <- "wt_z_score ~ SGA + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[12]] <- "wt_z_score ~ SGA + SGA:time_week + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[13]] <- "wt_z_score ~ SGA + SGA:I(time_week^2) + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"
formula_list[[14]] <- "wt_z_score ~ SGA + SGA:time_week + SGA:I(time_week^2) + time_week + I(time_week^2) + (time_week + I(time_week^2) | id)"

summary <- list()
basic_analysis <- list()

for (i in seq_along(formula_list)) {
  formula <- as.formula(formula_list[[i]])
  lmer.analysis <- lmerTest::lmer(formula, 
                                  data = long.df.scale, 
                                  REML = FALSE,
                                  control = lme4::lmerControl(optimizer = "optimx", 
                                                              calc.derivs = FALSE,
                                                              optCtrl = list(method = "nlminb"))
 )
  summary[[i]] <- list(as.character(formula(lmer.analysis)),
                       as.numeric(summary(lmer.analysis)$AICtab["AIC"]),
                       as.numeric(summary(lmer.analysis)$AICtab["BIC"]),
                       as.numeric(summary(lmer.analysis)$AICtab["logLik"]),
                       as.numeric(summary(lmer.analysis)$AICtab["deviance"]),
                       as.numeric(summary(lmer.analysis)$AICtab["df.resid"]))
  basic_analysis[[i]] <- lmer.analysis
}

basic_AIC <- data.frame()
for (i in seq_along(1:length(formula_list))) {
  summary[[i]][[1]] <- as.character(summary[[i]][[1]])
  summary[[i]][[1]] <- paste(summary[[i]][[1]][2], summary[[i]][[1]][1], summary[[i]][[1]][3])
  for (j in seq_along(1:length(summary[[i]]))) {
    basic_AIC[i,j] <- summary[[i]][[j]]
  }
}

colnames(basic_AIC) <- c("formula", "AIC", "BIC", "logLik", "deviance", "df_resid")

folder.check(folder.name = "Results", getwd())
folder.check(folder.name = "mixed_effects_modeling", getwd())

write_csv(basic_AIC, "basic_analysis.csv")

#for html markdown
kable(basic_AIC[c(1,2,3,4,5,6),c(1:3)], row.names = FALSE) %>%
  add_header_above(c("Placental insufficiency" = 3)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T)

kable(basic_AIC[c(1,2,7:10),c(1:3)], row.names = FALSE) %>%
  add_header_above(c("Birth weight" = 3)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T)

kable(basic_AIC[c(1,2,11:14),c(1:3)], row.names = FALSE) %>%
  add_header_above(c("Small-for-gestataional age" = 3)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T)

6. Mixed-effects modeling with lme4 considering covariates

# packages
require(foreach)

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

require(doParallel)

## Loading required package: doParallel

## Loading required package: iterators

## Loading required package: parallel

require(optimx)
require(parallel)
require(kableExtra)

basic_variables_pla <- c("id", "time_week", "wt_z_score", "pla_insuff")
basic_variables_bw <- c("time_week", "wt_z_score", "birth_wt")
nutrition_variables <- c("calories_7d", "aminoAcid_7d", "protein_fullEF", "calories_fullEF")
comorbidity_variables <- c("BPD", "IVH", "Sepsis", "ROP", "Abx_days")
demographics_variables <- c("Gender", "maternal_obesity", "maternal_smoking", "maternal_race")

long.df.scale_nutrition <- long.df.scale[,c(basic_variables_pla, nutrition_variables)]
long.df.scale_comorbidity <- long.df.scale[,c(basic_variables_pla, comorbidity_variables)]
long.df.scale_demographics <- long.df.scale[,c(basic_variables_pla, demographics_variables)]

fm <- formula_list[[4]]

## construct all possible combinations of covariates
var_list_demographics <- data.frame()
var_list_comorbidity <- data.frame()
var_list_nutrition <- data.frame()

fit <- function(variables, data) {
  var_list <- data.frame()
  for (i in seq_along(1:length(variables))) {
    v <- list()
    v <- combn(variables, i, simplify = TRUE)
    v <- t(v)
    v <- as.data.frame(v)
    v <- unite(v, variables, sep = " + ")
    var_list <- rbind(var_list, v)
  }

  for (i in seq_along(1:length(var_list$variables))) {
    formula <- str_replace(fm, "pla_insuff \\+ pla_insuff\\:time_week \\+ ", paste0(var_list$variables[i], " + pla_insuff \\+ pla_insuff\\:time_week \\+ "))
   var_list$formula[i] <- formula
  }

  var_list <- rbind(c("", fm), var_list)

  cl <- makePSOCKcluster(detectCores())
  registerDoParallel(cl)

  fit.all <- list()
  fit.all <- foreach (i = seq_along(1:nrow(var_list)),
                    .combine = rbind,
                    .packages = c("lmerTest", "lme4", "optimx")) %dopar% {
                      formula <- var_list$formula[i]
                      formula.formula <- as.formula(formula)
                      fit <- lmerTest::lmer(formula.formula, 
                                            data=data, 
                                            REML=F,
                                            control = lmerControl(optimizer = "optimx", calc.derivs = FALSE, 
                                                                  optCtrl = list(method = "nlminb")))
                      return(list(formula(fit),
                                  as.numeric(summary(fit)$AICtab["AIC"]),
                                  as.numeric(summary(fit)$AICtab["BIC"]),
                                  as.numeric(summary(fit)$AICtab["logLik"]),
                                  as.numeric(summary(fit)$AICtab["deviance"]),
                                  as.numeric(summary(fit)$AICtab["df.resid"]),
                                  summary(fit)$coefficient,
                                  fit))
                    }

for (i in seq_along(1:nrow(var_list))) {
  fit.all[[i]] <- as.character(fit.all[[i]])
  fit.all[[i]] <- paste(fit.all[[i]][2], fit.all[[i]][1],fit.all[[i]][3])
}
fit.all <- as.data.frame(fit.all)
colnames(fit.all) <- c("formula", "AIC", "BIC", "logLik", "deviance", "df_resid", "coefficient", "results")
fit.all.results <- fit.all$results
fit.all$results <- NULL
fit.all$formula <- as.character(fit.all$formula)

var_list <- left_join(var_list, fit.all, by = c("formula" = "formula"))
var_list$AIC <- unlist(var_list$AIC)
var_list$BIC <- unlist(var_list$BIC)
var_list$logLik <- unlist(var_list$logLik)
var_list$deviance <- unlist(var_list$deviance)
var_list$df_resid <- unlist(var_list$df_resid)

stopCluster(cl)

var_list <- rownames_to_column(var_list, var = "id")
var_list$id <- as.numeric(var_list$id)

return(list(var_list, fit.all.results))
}

fit_nutrition <- fit(nutrition_variables, long.df.scale_nutrition)
fit_comorbidity <- fit(comorbidity_variables, long.df.scale_comorbidity)
fit_demographics <- fit(demographics_variables, long.df.scale_demographics)

nutrition_AIC <- fit_nutrition[[1]]
comorbidity_AIC <- fit_comorbidity[[1]]
demographics_AIC <- fit_demographics[[1]]

nutrition_results <- fit_nutrition[[2]]
comorbidity_results <- fit_comorbidity[[2]]
demographics_restuls <- fit_demographics[[2]]

rm(fit_nutrition)
rm(fit_comorbidity)
rm(fit_demographics)

folder.check(folder.name = "Results", getwd())
folder.check(folder.name = "mixed_effects_modeling", getwd())

write_csv(nutrition_AIC[,c(1:8)], "nutrition_covariates.csv")
write_csv(comorbidity_AIC[,c(1:8)], "comorbidity_covariates.csv")
write_csv(demographics_AIC[,c(1:8)], "demographics_covariates.csv")

#for html Markdown
kableExtra::kable(nutrition_AIC[order(nutrition_AIC$AIC),c(2,4,5)]) %>% 
  add_header_above(c("Variable combination - order by increasing AIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T) %>%
  scroll_box(width = "100%", height = "300px")

kableExtra::kable(nutrition_AIC[order(nutrition_AIC$BIC),c(2,4,5)]) %>% 
  add_header_above(c("Variable combination - order by increasing BIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T) %>%
  scroll_box(width = "100%", height = "300px")

kableExtra::kable(comorbidity_AIC[order(comorbidity_AIC$AIC),c(2,4,5)]) %>% 
  add_header_above(c("Variable combination - order by increasing AIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T) %>%
  scroll_box(width = "100%", height = "300px")

kableExtra::kable(comorbidity_AIC[order(comorbidity_AIC$BIC),c(2,4,5)]) %>% 
  add_header_above(c("Variable combination - order by increasing BIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T) %>%
  scroll_box(width = "100%", height = "300px")

kableExtra::kable(demographics_AIC[order(demographics_AIC$AIC),c(2,4,5)]) %>%
  add_header_above(c("Variable combination - order by increasing AIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T) %>%
  scroll_box(width = "100%", height = "300px")

kableExtra::kable(demographics_AIC[order(demographics_AIC$BIC),c(2,4,5)]) %>%
  add_header_above(c("Variable combination - order by increasing BIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T) %>%
  scroll_box(width = "100%", height = "300px")

7. Assessing interaction term between covariates and time_week

formula_cov <- list()
formula_cov[[1]] <- "wt_z_score ~ calories_7d + Abx_days + pla_insuff + time_week + I(time_week^2) + pla_insuff:time_week + (time_week + I(time_week^2) | id)"
formula_cov[[2]] <- "wt_z_score ~ calories_7d + calories_7d:time_week + Abx_days + pla_insuff + time_week + I(time_week^2) + pla_insuff:time_week + (time_week + I(time_week^2) | id)"
formula_cov[[3]] <- "wt_z_score ~ calories_7d + Abx_days + Abx_days:time_week + pla_insuff + time_week + I(time_week^2) + pla_insuff:time_week + (time_week + I(time_week^2) | id)"

summary_cov <- list()
cov_interaction_analysis <- list()

for (i in seq_along(formula_cov)) {
  formula <- as.formula(formula_cov[[i]])
  lmer.analysis <- lmerTest::lmer(formula, 
                                  data = long.df.scale, 
                                  REML = FALSE,
                                  control = lme4::lmerControl(optimizer = "optimx", 
                                                              calc.derivs = FALSE,
                                                              optCtrl = list(method = "nlminb"))
 )
  summary_cov[[i]] <- list(as.character(formula(lmer.analysis)),
                       as.numeric(summary(lmer.analysis)$AICtab["AIC"]),
                       as.numeric(summary(lmer.analysis)$AICtab["BIC"]),
                       as.numeric(summary(lmer.analysis)$AICtab["logLik"]),
                       as.numeric(summary(lmer.analysis)$AICtab["deviance"]),
                       as.numeric(summary(lmer.analysis)$AICtab["df.resid"]))
  cov_interaction_analysis[[i]] <- lmer.analysis
}

cov_AIC <- data.frame()
for (i in seq_along(1:length(formula_cov))) {
  summary_cov[[i]][[1]] <- as.character(summary_cov[[i]][[1]])
  summary_cov[[i]][[1]] <- paste(summary_cov[[i]][[1]][2], summary_cov[[i]][[1]][1], summary_cov[[i]][[1]][3])
  for (j in seq_along(1:length(summary_cov[[i]]))) {
    cov_AIC[i,j] <- summary_cov[[i]][[j]]
  }
}

colnames(cov_AIC) <- c("formula", "AIC", "BIC", "logLik", "deviance", "df_resid")

folder.check(folder.name = "Results", getwd())
folder.check(folder.name = "mixed_effects_modeling", getwd())

write_csv(cov_AIC, "covariate_interaction_analysis.csv")

kableExtra::kable(cov_AIC[order(cov_AIC$AIC),c(1:3)]) %>% 
  add_header_above(c("Various combination of covariate-by-time interaction - order by increasing AIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T)

kableExtra::kable(cov_AIC[order(cov_AIC$BIC),c(1:3)]) %>% 
    add_header_above(c("Various combination of covariate-by-time interaction - order by increasing BIC" = 4)) %>%
  kable_styling(bootstrap_options = "striped", full_width = T)

8. Table and plot

require(sjPlot)

## Loading required package: sjPlot

## Learn more about sjPlot with 'browseVignettes("sjPlot")'.

folder.check(folder.name = "Results", getwd())
folder.check(folder.name = "Tables and plots", getwd())
fit.plot.pla <- basic_analysis[[4]]
fit.plot.bw <- basic_analysis[[8]]

# Plotting placental insufficiency-by-time interaction
pdf(file = "plot_pla.pdf", onefile = TRUE, width = 8, height = 6)
k <- plot_model(fit.plot.pla, show.values = TRUE,
                title = "Mixed-effects modeling - coefficient plot", 
                auto.label = FALSE,
                value.offset = 0.3,
                axis.title = c("estimates", "variable"),
                vline.color = "darkorchid",
                order.terms = c(1,4,2,3),
                digits = 3)
l <- plot_model(fit.plot.pla, type = "pred",
                terms = c("time_week", "pla_insuff"),
                title = "Predicted values of wt z score by time")

## Model contains polynomial or cubic / quadratic terms. Consider using `terms="time_week [all]"` to get smooth plots. See also package-vignette 'Marginal Effects at Specific Values'.

print(k + theme_bw())
print(l + theme_bw())
dev.off()

## png 
##   2

# Plotting birth weight-by-time interaction
pdf(file = "plot_bw.pdf", onefile = TRUE, width = 8, height = 6)
m <- plot_model(fit.plot.bw, show.values = TRUE,
                title = "Mixed-effects modeling - coefficient plot", 
                auto.label = FALSE,
                value.offset = 0.3,
                axis.title = c("estimates", "variable"),
                vline.color = "darkorchid",
                order.terms = c(1,4,2,3),
                digits = 3)
n <- plot_model(fit.plot.bw, type = "pred",
                terms = c("time_week", "birth_weight"),
                title = "Predicted values of wt z score by time")

## Model contains polynomial or cubic / quadratic terms. Consider using `terms="time_week [all]"` to get smooth plots. See also package-vignette 'Marginal Effects at Specific Values'.

print(m + theme_bw())
print(n + theme_bw())
dev.off()

## png 
##   2

require(sjPlot)
# for html Markdown
tab_model(fit.plot.pla, show.icc = FALSE, auto.label = FALSE,
          title = "Mixed-effects modeling - coefficient table",
          file = "table_pla.html",
          digits = 3)

print(k + theme_bw())

print(l + theme_bw())

require(sjPlot)
#for html Markdown
tab_model(fit.plot.bw, show.icc = FALSE, auto.label = FALSE,
          title = "Mixed-effects modeling - coefficient table",
          file = "table_bw.html",
          digits = 3)

print(m + theme_bw())

print(n + theme_bw())

9. Session Info

folder.check(folder.name = "Results", getwd())
sink(file = paste0(getwd(), "/seesionInfo.txt"))
sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] sjPlot_2.6.3      doParallel_1.0.14 iterators_1.0.10 
##  [4] foreach_1.4.4     kableExtra_1.1.0  PropCIs_0.3-0    
##  [7] effsize_0.7.4     optimx_2018-7.10  ggpubr_0.2       
## [10] magrittr_1.5      vcd_1.4-4         psych_1.8.12     
## [13] lmerTest_3.1-0    lme4_1.1-21       Matrix_1.2-17    
## [16] forcats_0.4.0     stringr_1.4.0     dplyr_0.8.0.1    
## [19] purrr_0.3.2       readr_1.3.1       tidyr_0.8.3      
## [22] tibble_2.1.1      ggplot2_3.1.1     tidyverse_1.2.1  
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-139       lubridate_1.7.4    RColorBrewer_1.1-2
##  [4] insight_0.3.0      webshot_0.5.1      httr_1.4.0        
##  [7] numDeriv_2016.8-1  tools_3.6.0        TMB_1.7.15        
## [10] backports_1.1.4    R6_2.4.0           sjlabelled_1.0.17 
## [13] lazyeval_0.2.2     colorspace_1.4-1   withr_2.1.2       
## [16] tidyselect_0.2.5   gridExtra_2.3      mnormt_1.5-5      
## [19] emmeans_1.3.4      compiler_3.6.0     performance_0.1.0 
## [22] cli_1.1.0          rvest_0.3.3        xml2_1.2.0        
## [25] labeling_0.3       bayestestR_0.1.0   scales_1.0.0      
## [28] mvtnorm_1.0-10     lmtest_0.9-37      digest_0.6.18     
## [31] foreign_0.8-70     minqa_1.2.4        rmarkdown_1.12    
## [34] pkgconfig_2.0.2    htmltools_0.3.6    highr_0.8         
## [37] rlang_0.3.4        readxl_1.3.1       rstudioapi_0.10   
## [40] generics_0.0.2     zoo_1.8-5          jsonlite_1.6      
## [43] Rcpp_1.0.1         munsell_0.5.0      stringi_1.4.3     
## [46] yaml_2.2.0         snakecase_0.9.2    MASS_7.3-51.1     
## [49] plyr_1.8.4         sjmisc_2.7.9       crayon_1.3.4      
## [52] lattice_0.20-38    ggeffects_0.10.0   haven_2.1.0       
## [55] cowplot_0.9.4      splines_3.6.0      sjstats_0.17.4    
## [58] hms_0.4.2          knitr_1.22         pillar_1.3.1      
## [61] boot_1.3-20        estimability_1.3   codetools_0.2-16  
## [64] glue_1.3.1         evaluate_0.13      modelr_0.1.4      
## [67] nloptr_1.2.1       cellranger_1.1.0   gtable_0.3.0      
## [70] assertthat_0.2.1   xfun_0.6           xtable_1.8-3      
## [73] broom_0.5.2        coda_0.19-2        viridisLite_0.3.0 
## [76] glmmTMB_0.2.3      ellipsis_0.1.0

sink()