Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Master-Thesis-Analyses/AntepartumDepression_Cortisol.R
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
2817 lines (2675 sloc)
130 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## Thesis Results: The diurnal cortisol rythm during pregnancy | |
| #setwd("C:/Users/User/Desktop/Internship/RScripts/Cortisol/Master_Thesis") | |
| setwd("C:/Users/alici/Desktop/Git_Folder/ITU_cortisol_analyses/Master_Thesis") | |
| ############ Data preparation ########################## | |
| ## libraries | |
| library(dplyr) | |
| library(ggplot2) | |
| library(nlme) | |
| library(sjPlot) | |
| library(data.table) | |
| library(ggeffects) | |
| library(tibble) | |
| library(DescTools) | |
| library(reghelper) | |
| library(APAstyler) | |
| library(car) | |
| ### wide format ################################## | |
| load("Rdata/ITU_combined_cortisol_dates_times_wide_format.Rdata") | |
| load("Rdata/processed_register_data.Rdata") | |
| wide_cort_reg_data <- left_join(wide_cort, | |
| register_data) | |
| #exclude corticosteroid treatment | |
| no_cortico_treat <- wide_cort_reg_data[wide_cort_reg_data$Maternal_Corticosteroid_Treatment_during_Pregnancy == "no",] | |
| #exclude violattions against cortisol collection protocol | |
| comments_to_be_excluded <- c(88,6,12,3,15,9999,7,16,45,5,34,17,57,67,99) | |
| final_wide_cort <- no_cortico_treat[!(no_cortico_treat$notes %in% comments_to_be_excluded),] | |
| ##pregstage factoring | |
| final_wide_cort$pregstage <- factor(final_wide_cort$pregstage, ordered = FALSE) | |
| ##factoring of predictors | |
| final_wide_cort$pregstage <- factor(final_wide_cort$pregstage) | |
| final_wide_cort$participantID <- factor(final_wide_cort$participantID) | |
| final_wide_cort$ELISA_analysis_plate <- factor(final_wide_cort$ELISA_analysis_plate) | |
| final_wide_cort$Maternal_Diabetes_Disorders_anyVSnone <- factor(final_wide_cort$Maternal_Diabetes_Disorders_anyVSnone, | |
| levels = c(-999,0,1), | |
| labels = c("none","none", "yes")) | |
| final_wide_cort$Maternal_Hypertensive_Disorders_anyVSnone <- factor(final_wide_cort$Maternal_Hypertensive_Disorders_anyVSnone, | |
| levels = c(-999,0,1), | |
| labels = c("none","none", "yes")) | |
| #data preparation | |
| #boxplot(final_wide_cort$cort_d_AUC ~ final_wide_cort$pregstage) | |
| final_wide_cort$time_of_day_in_hours_since_midnight_S1 <- as.numeric(final_wide_cort$time_of_day_in_hours_since_midnight_S1) | |
| final_wide_cort$season <- factor(final_wide_cort$season) | |
| ########### add well_being data to it | |
| load("Rdata/processed_wellbeingduringpreg_completevars.Rdata") | |
| #for now only sleep, nausea, anxiety and CESD are included | |
| cols <- c("1", "2",grep("gestage", names(q_data_complete)), | |
| grep("GWdiff_cort_qclosest", names(q_data_complete)), | |
| grep("Cesd", names(q_data_complete)), | |
| # grep("puqe", names(q_data_complete)), | |
| grep("PSQI", names(q_data_complete)), | |
| # grep("ESS", names(q_data_complete)), | |
| grep("BAI",names(q_data_complete))) | |
| q_data_sub <- q_data_complete[,c(as.numeric(cols))] | |
| #combine with cortisol | |
| final_wide_cort <- left_join(final_wide_cort, | |
| q_data_sub, | |
| by = c("participantID", "pregstage", "gestage_weeks")) | |
| ## Sample stratification | |
| ## those with multiple assessments | |
| complete_assessments <- c() | |
| two_assessments <- c() | |
| at_least_two_assessments <- c() | |
| once <- c() | |
| final_wide_cort <- as.data.frame(final_wide_cort) | |
| IDs <- unique(final_wide_cort$participantID) | |
| IDs_I <- unique(final_wide_cort$participantID[final_wide_cort$pregstage == "I"]) | |
| IDs_II <- unique(final_wide_cort$participantID[final_wide_cort$pregstage == "II"]) | |
| IDs_III <- unique(final_wide_cort$participantID[final_wide_cort$pregstage == "III"]) | |
| for(i in 1:length(IDs)){ | |
| ID <- IDs[[i]] | |
| #print(ID) | |
| trims <- unique(final_wide_cort[final_wide_cort$participantID == ID, c("pregstage")]) | |
| #print(nrow(trims)) | |
| if(length(trims) == 3){ | |
| complete_assessments <- append(complete_assessments,ID) | |
| } | |
| if(length(trims) == 2){ | |
| two_assessments <- append(two_assessments, ID) | |
| } | |
| if(length(trims) >= 2){ | |
| at_least_two_assessments <- append(at_least_two_assessments, ID) | |
| } | |
| if(length(trims) == 1){ | |
| once <- append(once, ID) | |
| } | |
| } | |
| # IDs_T1_T2 <- which(two_assessments %in% IDs_I & two_assessments %in% IDs_II) | |
| # length(IDs_T1_T2) #68 | |
| # IDs_T1_T3 <- which(two_assessments %in% IDs_I & two_assessments %in% IDs_III) | |
| # length(IDs_T1_T3) #32 | |
| # IDs_T2_T3 <- which(two_assessments %in% IDs_II & two_assessments %in% IDs_III) | |
| # length(IDs_T2_T3) #145 | |
| # Q_df_wide <- subset(final_wide_cort, abs(GWdiff_cort_qclosest_to_cortGW) < 2) | |
| # Q_df_wide %>% | |
| # group_by(pregstage) %>% | |
| # summarise(n=n(), | |
| # GW = mean(gestage_weeks,na.rm=T), | |
| # min_GW = min(gestage_weeks,na.rm=T), | |
| # max_GW = max(gestage_weeks,na.rm=T), | |
| # | |
| # ) | |
| ### long format #### | |
| load("Rdata/ITU_combined_cortisol_dates_times_long_format.Rdata") | |
| # long_cort <- long_cort %>% | |
| # group_by(participantID) %>% | |
| # mutate(nr_assessments = max(row_number())) | |
| # Nrs <- long_cort[,c("participantID", "nr_assessments")] | |
| # Nrs <- subset(Nrs, !duplicated(Nrs)) | |
| # mean(Nrs$nr_assessments, na.rm = T) | |
| # | |
| # long_cort <- long_cort %>% | |
| # group_by(participantID, pregstage) %>% | |
| # mutate(nr_assessments = max(row_number())) | |
| # | |
| # #T1 | |
| # Nrs <- long_cort[long_cort$pregstage == "I",c("participantID", "nr_assessments")] | |
| # Nrs <- subset(Nrs, !duplicated(Nrs)) | |
| # mean(Nrs$nr_assessments, na.rm = T) | |
| # | |
| # #T2 | |
| # Nrs <- long_cort[long_cort$pregstage == "II",c("participantID", "nr_assessments")] | |
| # Nrs <- subset(Nrs, !duplicated(Nrs)) | |
| # mean(Nrs$nr_assessments, na.rm = T) | |
| # | |
| # #T3 | |
| # Nrs <- long_cort[long_cort$pregstage == "III",c("participantID", "nr_assessments")] | |
| # Nrs <- subset(Nrs, !duplicated(Nrs)) | |
| # mean(Nrs$nr_assessments, na.rm = T) | |
| # | |
| # length(unique(long_cort$participantID[long_cort$nr_assessments == 7])) | |
| #register data as prepared | |
| load("Rdata/processed_register_data.Rdata") | |
| final_cort <- left_join(long_cort, | |
| register_data) | |
| ## Exclusion based on register data: betamethasone treatment and protocol violations #### | |
| no_cortico_treat <- final_cort[final_cort$Maternal_Corticosteroid_Treatment_during_Pregnancy == "no",] | |
| #length(unique(final_cort$participantID[final_cort$Maternal_Corticosteroid_Treatment_during_Pregnancy == 1])) #N=15 | |
| #exclude violations against protocol | |
| comments_to_be_excluded <- c(88,6,12,3,15,9999,7,16,45,5,34,17,57,67,99) | |
| final_long_cort <- no_cortico_treat[!(no_cortico_treat$notes %in% comments_to_be_excluded),] | |
| ## merge with well-being questionnaires #### | |
| final_long_cort <- left_join(final_long_cort, | |
| as.data.frame(final_wide_cort[,c("participantID", | |
| "gestage_weeks_qclosest_to_cortGW", | |
| "pregstage", | |
| "Cesd_qclosest_to_cortGW", | |
| "Cesd_qclosest_to_cortGW_pregMean", | |
| "Cesd_qclosest_to_cortGW_pregMeanDeviat", | |
| "PSQI_qclosest_to_cortGW", | |
| "PSQI_qclosest_to_cortGW_pregMean", | |
| "PSQI_qclosest_to_cortGW_pregMeanDeviat", | |
| "BAI_qclosest_to_cortGW", | |
| "BAI_qclosest_to_cortGW_pregMean", | |
| "BAI_qclosest_to_cortGW_pregMeanDeviat", | |
| "qclosest_based_clin_Cesd", # "puqe_sumscore_qclosest_to_cortGW", | |
| # "puqe_sumscore_qclosest_to_cortGW_pregMean", | |
| # "puqe_sumscore_qclosest_to_cortGW_pregMeanDeviat", | |
| "GWdiff_cort_qclosest_to_cortGW" | |
| )]), | |
| copy=TRUE, | |
| by=c("participantID", "pregstage")) | |
| ##factoring of predictors | |
| final_long_cort$pregstage <- factor(final_long_cort$pregstage) | |
| final_long_cort$participantID <- factor(final_long_cort$participantID) | |
| final_long_cort$ELISA_analysis_plate <- factor(final_long_cort$ELISA_analysis_plate) | |
| #Gestational week centered at mean GW at first visit | |
| final_long_cort$gestage_weeks_cent <- | |
| final_long_cort$gestage_weeks - mean(final_long_cort$gestage_weeks[final_long_cort$pregstage == "I"], na.rm=T) | |
| # ##Questionnaires centered at T1 scores | |
| # final_long_cort$Cesd_T1 <- NA | |
| # final_long_cort$PSQI_T1 <- NA | |
| # final_long_cort$BAI_T1 <- NA | |
| # | |
| # | |
| # for(i in 1:length(final_long_cort$participantID)){ #for each participant ID | |
| # ID <- final_long_cort$participantID[[i]] | |
| # | |
| # cesd <- unique(final_long_cort$Cesd_qclosest_to_cortGW[final_long_cort$pregstage == "I" & | |
| # final_long_cort$participantID == unique(ID)]) | |
| # psqi <- unique(final_long_cort$PSQI_qclosest_to_cortGW[final_long_cort$pregstage == "I" & | |
| # final_long_cort$participantID == unique(ID)]) | |
| # bai <- unique(final_long_cort$BAI_qclosest_to_cortGW[final_long_cort$pregstage == "I" & | |
| # final_long_cort$participantID == unique(ID)]) | |
| # | |
| # if(length(cesd) > 0){ | |
| # final_long_cort$Cesd_T1[final_long_cort$participantID == ID] <- cesd | |
| # } | |
| # if(length(psqi) > 0){ | |
| # final_long_cort$PSQI_T1[final_long_cort$participantID == ID] <- psqi | |
| # } | |
| # if(length(bai) > 0){ | |
| # final_long_cort$BAI_T1[final_long_cort$participantID == ID] <- bai | |
| # } | |
| # } | |
| # #view(final_long_cort[,c("participantID", "Cesd_qclosest_to_cortGW","Cesd_T1")]) | |
| # | |
| # ## Change in psych symptoms from T1 | |
| # final_long_cort$Cesd_T1.WS <- final_long_cort$Cesd_T1 - final_long_cort$Cesd_qclosest_to_cortGW | |
| # final_long_cort$PSQI_T1.WS <- final_long_cort$PSQI_T1 - final_long_cort$PSQI_qclosest_to_cortGW | |
| # final_long_cort$BAI_T1.WS <- final_long_cort$BAI_T1 - final_long_cort$BAI_qclosest_to_cortGW | |
| final_long_cort$Maternal_Diabetes_Disorders_anyVSnone <- factor(final_long_cort$Maternal_Diabetes_Disorders_anyVSnone, | |
| levels = c(-999,0,1), | |
| labels = c("none","none", "yes")) | |
| final_long_cort$Maternal_Hypertensive_Disorders_anyVSnone <- factor(final_long_cort$Maternal_Hypertensive_Disorders_anyVSnone, | |
| levels = c(-999,0,1), | |
| labels = c("none","none", "yes")) | |
| final_long_cort$season <- factor(final_long_cort$season) | |
| ## Add Maternal Education #### | |
| maternal_edu <- read.delim("Rdata/ITU maternal education.dat") | |
| names(maternal_edu)[1] <- "participantID" | |
| names(maternal_edu)[3] <- "Maternal_Education" | |
| maternal_edu <- maternal_edu[,c("participantID", | |
| "Maternal_Education")] | |
| library(naniar) | |
| maternal_edu <- maternal_edu %>% replace_with_na(replace = list(Maternal_Education = -9)) | |
| maternal_edu$Maternal_Education <- factor(maternal_edu$Maternal_Education, | |
| levels = c(1,2,3), | |
| labels = c("primary", "applied university", "university")) | |
| final_long_cort <- left_join(final_long_cort, | |
| maternal_edu) | |
| ### add medication data #### | |
| medication <- read.delim("Rdata/ITU_psychotrophicmedication05July21_Maternal_CurrentPregnancy.dat") | |
| names(medication)[1] <- "participantID" | |
| final_long_cort <- left_join(final_long_cort, | |
| medication[,c("participantID", | |
| "ITUbroadpsychiatricmedication_18_KELA", | |
| "antidepressants_18_KELA")], | |
| copy=TRUE) | |
| final_long_cort$ITUbroadpsychiatricmedication_18_KELA <- factor(final_long_cort$ITUbroadpsychiatricmedication_18_KELA, | |
| levels = c(0,1), | |
| labels = c("no", "yes")) | |
| #tidy up | |
| rm(list= ls()[!(ls() %in% c("final_long_cort", "at_least_two_assessments"))]) | |
| ## Exclusion Criteria based on cortisol data and protocol violations #### | |
| ### exclude observations CAR protocol violation (i.e., stayed in bed | |
| #at S2 and S3, or took sample 1 hour after leaving the bed) | |
| #mark the observations | |
| final_long_cort$CAR_exclude <- NA | |
| for(i in 1:nrow(final_long_cort)){ | |
| tp <- final_long_cort$timeslot[[i]] | |
| CAR_time <- final_long_cort$hours_since_waking[[i]] | |
| note <- final_long_cort$notes[[i]] | |
| if(!is.na(tp) & !is.na(note)){ | |
| if((tp == "S2" | tp == "S3") & note == 2){ | |
| final_long_cort$CAR_exclude[[i]] <- "exc" | |
| } | |
| } | |
| if(!is.na(tp) & !is.na(CAR_time)){ | |
| if((tp == "S2" | tp == "S3") & (CAR_time > 0.75)){ | |
| final_long_cort$CAR_exclude[[i]] <- "exc" | |
| } | |
| } | |
| } | |
| # only include Q-data within 2 week radius of cortisol assessment | |
| final_long_cort <- subset(final_long_cort, is.na(CAR_exclude)) | |
| #length(unique(final_long_cort$participantID[!is.na(final_long_cort$Cesd_qclosest_to_cortGW_pregMean)])) | |
| #length(unique(final_long_cort$participantID[final_long_cort$GWdiff_cort_qclosest_to_cortGW <= 2])) | |
| #length(unique(final_long_cort$participantID[final_long_cort$GWdiff_cort_qclosest_to_cortGW <= 1])) | |
| #apply 2 week interval for Q and cortisol data assessment | |
| final_long_cort <- subset(final_long_cort, abs(GWdiff_cort_qclosest_to_cortGW) < 2) | |
| #summary(final_long_cort$GWdiff_cort_qclosest_to_cortGW) | |
| ### Centering #### | |
| final_long_cort$Gestational_Age_Weeks_cent <- c(scale(final_long_cort$Gestational_Age_Weeks, scale = F)) | |
| final_long_cort$waking_time_cent <- c(scale(final_long_cort$waking_time, scale = F)) | |
| final_long_cort$Maternal_Age_Years_cent <- c(scale(final_long_cort$Maternal_Age_Years, scale = F)) | |
| final_long_cort$Maternal_Body_Mass_Index_in_Early_Pregnancy_cent <- c(scale(final_long_cort$Maternal_Body_Mass_Index_in_Early_Pregnancy, scale = F)) | |
| final_long_cort$Weight_Gain_cent <- c(scale(final_long_cort$Weight_Gain, scale = F)) | |
| ### Depressive Symptom Severity #### | |
| final_long_cort$Depressive_Symptom_Severity <- NA | |
| for(i in 1:nrow(final_long_cort)){ | |
| m <- final_long_cort$Cesd_qclosest_to_cortGW_pregMean[[i]] | |
| if(!is.na(m)){ | |
| if(m<16){ | |
| final_long_cort$Depressive_Symptom_Severity[[i]] <- "Below_Clinical (CES-D < 16)" | |
| } | |
| if(m>=16){ | |
| final_long_cort$Depressive_Symptom_Severity[[i]] <- "Clinical (CES-D >= 16)" | |
| } | |
| } | |
| } | |
| final_long_cort$Depressive_Symptom_Severity <- factor(final_long_cort$Depressive_Symptom_Severity) | |
| ## Functions to save results #### | |
| round_df <- function(x, cols1, digits1, cols2, digits2){ | |
| x <- as.data.frame(x) | |
| # round all variables | |
| # x: data frame | |
| # cols: select which columns to round | |
| # digits: number of digits to round corresponding to cols | |
| for(i in 1:ncol(x[,c(cols1)])){ | |
| colname <- names(x[,c(cols1)])[[i]] | |
| col <- x[,c(colname)] | |
| for(j in 1:length(col)){ | |
| value <- as.numeric(col[[j]]) | |
| if(!is.na(value)){ | |
| # print(value) | |
| if(abs(value) > 0.005){ | |
| x[,c(cols1)][[j,i]] <- format(round(value, digits1)) | |
| } | |
| } | |
| } | |
| } | |
| x[,c(cols2)] <- round(x[,c(cols2)], digits2) | |
| return(x) | |
| } | |
| #calculate percent changes | |
| res_Bchange <- function(res_df){ ##tTable of a summary model object as dataframe | |
| res_df[["pchange"]] <- NA | |
| for(i in 1:nrow(res_df)){ | |
| p <- as.numeric(res_df[["p-value"]][[i]]) | |
| val <- res_df[["Value"]][[i]] | |
| if(p <= 0.05){ | |
| res_df[["pchange"]][[i]] <- ((exp(val))-1)*100 | |
| } | |
| } | |
| df <- res_df[,-c(3)] | |
| df[1,5] <- NA | |
| names(df) <- c("Coefficient", "SE", "t-statistic", "p", "Interpretation") | |
| df <- round_df(df, c(1:3,5), 2, 4,3) | |
| #to follow APA style, delete leading zero | |
| for(j in 1:nrow(df)){ | |
| P <- df[[j, "p"]] | |
| if(P < 1 & P > 0){ | |
| P <- snip(P, lead = 1) | |
| df[[j, "p"]] <- snip(P, lead = 1) | |
| } | |
| } | |
| return(df) | |
| } | |
| #extract the random part | |
| res_random <- function(model_summary){ #model summary object | |
| random <- VarCorr(model_summary)[(c(2:3,5:6)), c(1,3)] | |
| random <- as.data.frame(random) | |
| names(random) <- c("Variance_Component", "Covariance") | |
| random$Variance_Component <- as.numeric(random$Variance_Component) | |
| random$Covariance <- as.numeric(random$Covariance) | |
| #random <- round_df(random, 1:2,2, 1:2,2) | |
| n <- nobs(model_summary) | |
| random <- rbind(random, c(n, NA)) | |
| varests <- as.numeric(VarCorr(model_summary)[(c(2:3,5:6)), c(1)]) | |
| total_R_explained <- format(round(sum(varests[1:3])/sum(varests)*100,1)) | |
| random <- rbind(random, c(total_R_explained, NA)) | |
| row.names(random) <- c("Intercept (person-level)", | |
| "Time", | |
| "Intercept (pregstage-level)", | |
| "Residual", | |
| "Number of Observations", | |
| "Total R explained") | |
| return(random) | |
| } | |
| model_res <- function(model_summary){ | |
| tTable <- res_Bchange(as.data.frame(model_summary[["tTable"]])) | |
| #tTable <- round_df(tTable,c(2:4,6),2,5,3) | |
| random_part <- res_random(model_summary) | |
| random_part$empty1 <- NA | |
| random_part$empty2 <- NA | |
| random_part$empty3 <- NA | |
| tTable <- rbind(tTable,c("Variance Component", "Covariance", rep(NA,3))) | |
| names(random_part) <- names(tTable) | |
| Res <- rbind(tTable, random_part) | |
| Res <- setDT(Res, keep.rownames = "predictor") | |
| return(Res) | |
| } | |
| ############################################################################# | |
| ####### CAR Model #### | |
| CAR_times <- c("S1", "S2", "S3") | |
| CAR_df <- subset(final_long_cort, timeslot %in% CAR_times) | |
| ### Baseline Model: Random effects and variance-covariance structure #### | |
| #Including diurnal characteristics | |
| CAR_Form <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent | |
| ## Random Intercept only model | |
| M.CAR.RI <- lme(CAR_Form, | |
| random = list(~ 1|participantID, | |
| ~ 1|pregstage), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| # plot(M.CAR.RI) | |
| # qqnorm(M.CAR.RI, abline = c(0,1)) | |
| ## Does adding ELISA improve model fit? | |
| M.CAR.RI_wELISA <- update(M.CAR.RI, random = list(~1|ELISA_analysis_plate, | |
| ~1|participantID, | |
| ~1|pregstage)) | |
| tab1 <- anova(M.CAR.RI,M.CAR.RI_wELISA) #no improvement, in fact non-sig less model fit | |
| ## Does adding a random slope for CAR improve model fit? | |
| #i.e., allowing CAR to vary across participants | |
| M.CAR.RS <- update(M.CAR.RI, random = list(~hours_since_waking|participantID, | |
| ~1|pregstage)) | |
| res.RS <- anova(M.CAR.RI, M.CAR.RS) #Yes, it does improve fit | |
| tab2 <- rbind.data.frame(tab1, res.RS) | |
| #Does adding heterogenous WS variance for the effect of time improve fit? | |
| M.CAR.hetero <- lme(CAR_Form, | |
| random = list(~hours_since_waking|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| res.hetero <- anova(M.CAR.RS, M.CAR.hetero) #Does not improve fit! | |
| tab3 <- rbind.data.frame(tab2, res.hetero) | |
| #plot(M.CAR.hetero) #no fit | |
| #autocorrelation | |
| #Does adding auto-correlation improve fit? | |
| M.CAR.auto <- lme(CAR_Form, | |
| random = list(~hours_since_waking|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| res.auto <- anova(M.CAR.RS, M.CAR.auto) #Does improve fit! | |
| RES_CAR_Mselect <- rbind.data.frame(tab3, res.auto) | |
| #plot(M.CAR.auto) | |
| #qqnorm(M.CAR.auto, abline = c(0,1)) | |
| ### Baseline Model: Fixed Component #### | |
| M.CAR_ML1 <- lme(CAR_Form, | |
| random = list(~hours_since_waking|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "ML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| ## check whether quadratic term of time since waking indeed improves model fit | |
| M.CAR_ML2 <- update(M.CAR_ML1, .~. -I(hours_since_waking^2)) | |
| anova(M.CAR_ML1, M.CAR_ML2) #sig worse fit | |
| ## adding a quadratic term for GW | |
| M.CAR_ML3 <- update(M.CAR_ML1, .~. +I(gestage_weeks_cent^2)) | |
| anova(M.CAR_ML1, M.CAR_ML3) #sig better fit | |
| summary(M.CAR_ML3) | |
| ## adding a waking time and season | |
| M.CAR_ML4 <- update(M.CAR_ML3, .~. + | |
| waking_time_cent + | |
| waking_time_cent:hours_since_waking + | |
| season + | |
| season:hours_since_waking) | |
| baseline_CARres <- car::Anova(M.CAR_ML4) | |
| summary(M.CAR_ML3) # 4361.47 | |
| summary(M.CAR_ML4) #4341.392 | |
| #document waking time findings | |
| RES_CAR_B.ns_waking_main <- list(format(round(baseline_CARres$Chisq[5],2)), baseline_CARres$Df[5], format(round(baseline_CARres$`Pr(>Chisq)`[5],2))) | |
| names(RES_CAR_B.ns_waking_main) <- c("chi", "df", "p") | |
| RES_CAR_B.ns_waking_slp <- list(format(round(baseline_CARres$Chisq[8],2)), baseline_CARres$Df[8], format(round(baseline_CARres$`Pr(>Chisq)`[8],2))) | |
| names(RES_CAR_B.ns_waking_slp) <- c("chi", "df", "p") | |
| #season | |
| RES_CAR_B.ns_season_slp <- list(format(round(baseline_CARres$Chisq[9],2)), baseline_CARres$Df[9], format(round(baseline_CARres$`Pr(>Chisq)`[9],2))) | |
| names(RES_CAR_B.ns_season_slp) <- c("chi", "df", "p") | |
| ## re-evaluate when only including p<0.1 | |
| M.CAR_ML5 <- update(M.CAR_ML3, .~. + | |
| season) | |
| ##### Baseline Model: Re-evaluate random effects structure | |
| ### Final Baseline Model: Model Validation #### | |
| Morning_Cort_Form <- | |
| lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season | |
| M.Baseline_Morning_cort <- lme(Morning_Cort_Form, | |
| random = list(~hours_since_waking|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| #plot(M.Baseline_Morning_cort) | |
| #qqnorm(M.Baseline_Morning_cort, abline = c(0,1)) #overall normal distirbution of residuals can be assumed! | |
| #qqnorm(M.Baseline_Morning_cort, ~ranef(., level=1), abline=c(0,1)) #random effects not normal! | |
| #qqnorm(M.Baseline_Morning_cort, ~ranef(., level=2), abline=c(0,1)) #random effects not normal! | |
| ### (1) Baseline Model: Effects of Gestation #### | |
| car::Anova(M.Baseline_Morning_cort) | |
| #tab_model(M.Baseline_Morning_cort, show.r2 = TRUE, show.intercept = TRUE) | |
| RES_CAR_Baseline <- model_res(summary(M.Baseline_Morning_cort)) | |
| RES_CAR_Baseline <- mutate_all(RES_CAR_Baseline,funs(replace(., is.na(.), ""))) | |
| ## Total Variance explained | |
| varests <- as.numeric(VarCorr(M.Baseline_Morning_cort)[(c(2:3,5:6)), c(1)]) | |
| RES_CAR_total_R_explained <- format(round(sum(varests[1:3])/sum(varests)*100,1)) #R^2=54% | |
| CAR_residual <- sum(varests[4])/sum(varests) #45.6% | |
| #CAR_GW_variance <- varests[2]/sum(varests)*100 #0.0099% | |
| # ## vaiance explained in Null-Model | |
| # CAR_Null <- lg_cort_concentration_ug_l ~ 1 | |
| # | |
| # M.Baseline_Morning_cort <- lme(CAR_Null, | |
| # random = list(~1|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # varests <- as.numeric(VarCorr(M.Baseline_Morning_cort)[(c(2,4:5)), 1]) | |
| # RES_CAR_subject_R_explained <- format(round(sum(varests[1])/sum(varests)*100,2)) #R^2=26% | |
| # RES_CAR_pregstage_R_explained <- format(round(sum(varests[2])/sum(varests)*100,5)) #R^2=26% | |
| #plot_model(M.Baseline_Morning_cort, type = "re", level = 1) | |
| ##Extract Parameters for associations with child outcome | |
| CAR_parameters <- as.data.frame(coef(M.Baseline_Morning_cort)) | |
| save(CAR_parameters, file = "Rdata/CAR_parameters.Rdata") | |
| # RI of given participant indicates the difference of individual cortisol concentration | |
| # to overall sample mean concentration across pregnancy, and RS the individual differences | |
| # in CAR across prengnancy | |
| #head(ranef(M.Baseline_Morning_cort, level = 1)) | |
| #random effect of pregnancy stage, indicates the difference in cortisol concentration | |
| #within individuals from one measurement day to another | |
| #head(ranef(M.Baseline_Morning_cort, level = 2)) | |
| #head(coef(M.Baseline_Morning_cort)) | |
| #Intercept = cortisol concentration at awakening, at on average first GW of cortisol assessment period | |
| #CAR = linear CAR slope | |
| #CAR^2 = degree of curvilinearity | |
| #gestage_weeks_cent = change of average cortisol concentration from baseline | |
| #gestage_weeks_cent^2 = degree of curvilinearity | |
| #CAR:gestage_weeks_cent = effect of pregressing gestation on the diurnal slope | |
| # mean_II <- mean(final_long_cort$gestage_weeks_cent[final_long_cort$pregstage == "II"], na.rm=T) | |
| # mean_III <- mean(final_long_cort$gestage_weeks_cent[final_long_cort$pregstage == "III"], na.rm=T) | |
| # GW_values <- c(0,mean_II, mean_III) | |
| # CAR_across_preg <- ggpredict(M.Baseline_Morning_cort, terms = c("hours_since_waking", "gestage_weeks_cent[GW_values]")) | |
| # plot(CAR_across_preg) | |
| # Average_cort_across_preg <- ggpredict(M.Baseline_Morning_cort, terms = c("gestage_weeks_cent[GW_values]")) | |
| # #plot(Average_cort_across_preg) | |
| # season_across_preg <- ggpredict(M.Baseline_Morning_cort, terms = c("season")) | |
| # plot(season_across_preg) | |
| # Average_cort_across_preg <- ggpredict(M.Baseline_Morning_cort, terms = c("pregstage")) | |
| # plot(Average_cort_across_preg) | |
| # Average_cort_across_preg <- ggpredict(M.Baseline_Morning_cort, terms = c("CAR", "pregstage")) | |
| # plot(Average_cort_across_preg) | |
| ### (2) Covariates: Effects of maternal characteristics #### | |
| Baseline_CAR_Form.COV <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Maternal_Education + | |
| Parity + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| Weight_Gain_cent + | |
| Maternal_Hypertensive_Disorders_anyVSnone + | |
| Maternal_Diabetes_Disorders_anyVSnone + | |
| Maternal_Smoking_During_Pregnancy + | |
| Gestational_Age_Weeks_cent + | |
| Mother_Cohabiting + | |
| #interaction with time slope | |
| caseVScontrol:hours_since_waking + | |
| Maternal_Age_Years_cent:hours_since_waking + | |
| Maternal_Education:hours_since_waking + | |
| Parity:hours_since_waking + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| Weight_Gain_cent:hours_since_waking + | |
| Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| Gestational_Age_Weeks_cent:hours_since_waking + | |
| Mother_Cohabiting:hours_since_waking + | |
| #interaction with GW slope | |
| caseVScontrol:gestage_weeks_cent + | |
| Maternal_Age_Years_cent:gestage_weeks_cent + | |
| Maternal_Education:gestage_weeks_cent + | |
| Parity:gestage_weeks_cent + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent + | |
| Weight_Gain_cent:gestage_weeks_cent + | |
| Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent + | |
| Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent + | |
| Maternal_Smoking_During_Pregnancy:gestage_weeks_cent + | |
| Gestational_Age_Weeks_cent:gestage_weeks_cent + | |
| Mother_Cohabiting:gestage_weeks_cent | |
| M.Baseline_CAR_cort.COV <- lme(Baseline_CAR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_CAR_cort.COV) | |
| # ### instead investigating each covariates seperately #### | |
| # ## Case VS Control #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # caseVScontrol + | |
| # caseVScontrol:hours_since_waking + | |
| # caseVScontrol:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.CoCa <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.CoCa) | |
| # | |
| # # Repeat in subsample #### | |
| # M.Baseline_CAR_cort.COV.CoCa <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.CoCa) | |
| # | |
| # | |
| # ## Maternal Age #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Maternal_Age_Years_cent + | |
| # Maternal_Age_Years_cent:hours_since_waking + | |
| # Maternal_Age_Years_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.age <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.age) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.age <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.age) | |
| # | |
| # | |
| # ## maternal education #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Maternal_Education + | |
| # Maternal_Education:hours_since_waking + | |
| # Maternal_Education:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.edu <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.edu) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.edu <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.edu) | |
| # | |
| # ## Parity #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Parity + | |
| # Parity:hours_since_waking + | |
| # Parity:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.par <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.par) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.par <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.par) | |
| # | |
| # ## BMI #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.bmi <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.bmi) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.bmi <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.bmi) | |
| # | |
| # ## Weight Gain #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Weight_Gain_cent + | |
| # Weight_Gain_cent:hours_since_waking + | |
| # Weight_Gain_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.wg <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.wg) | |
| # | |
| # # repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.wg <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.wg) | |
| # | |
| # ## Hypertensive disorders #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Maternal_Hypertensive_Disorders_anyVSnone + | |
| # Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| # Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.hyp <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.hyp) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.hyp <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.hyp) | |
| # | |
| # ## Diabetes Disorders #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Maternal_Diabetes_Disorders_anyVSnone + | |
| # Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.diab <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.diab) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.diab <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.diab) | |
| # | |
| # ## Cohabitation #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Mother_Cohabiting + | |
| # Mother_Cohabiting:hours_since_waking + | |
| # Mother_Cohabiting:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.Co <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.Co) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.Co <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.Co) | |
| # | |
| # ## Smoking #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Maternal_Smoking_During_Pregnancy + | |
| # Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| # Maternal_Smoking_During_Pregnancy:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.Smoke <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.Smoke) | |
| # | |
| # # repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.Smoke <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.Smoke) | |
| # | |
| # ## Gestational Length #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Gestational_Age_Weeks_cent + | |
| # Gestational_Age_Weeks_cent:hours_since_waking + | |
| # Gestational_Age_Weeks_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # | |
| # | |
| # ## Child Sex #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Child_Sex + | |
| # Child_Sex:hours_since_waking + | |
| # Child_Sex:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # | |
| # ## PSQI #### | |
| Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season + | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent + | |
| PSQI_qclosest_to_cortGW_pregMean + | |
| PSQI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| PSQI_qclosest_to_cortGW_pregMean:gestage_weeks_cent | |
| M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # | |
| # ## BAI #### | |
| # Baseline_CAR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent + | |
| # BAI_qclosest_to_cortGW_pregMean + | |
| # BAI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| # BAI_qclosest_to_cortGW_pregMean:gestage_weeks_cent | |
| # | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_CAR_cort.COV.gl <- lme(Baseline_CAR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_CAR_cort.COV.gl) | |
| ### (3) Covariate Selection Result #### | |
| Baseline_CAR_Form.COV <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Maternal_Education + | |
| Parity + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| Weight_Gain_cent + | |
| #Maternal_Hypertensive_Disorders_anyVSnone + | |
| Maternal_Diabetes_Disorders_anyVSnone + | |
| #Maternal_Smoking_During_Pregnancy + | |
| #Gestational_Age_Weeks_cent + | |
| Mother_Cohabiting + | |
| PSQI_qclosest_to_cortGW_pregMean + | |
| #BAI_qclosest_to_cortGW_pregMean + | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| #BAI_qclosest_to_cortGW_pregMeanDeviat + | |
| ##interaction with time slope | |
| #caseVScontrol:hours_since_waking + | |
| #Maternal_Age_Years_cent:hours_since_waking + | |
| Maternal_Education:hours_since_waking + | |
| #Parity:hours_since_waking + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| #Weight_Gain_cent:hours_since_waking + | |
| #Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| #Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| #Gestational_Age_Weeks_cent:hours_since_waking + | |
| #Mother_Cohabiting:hours_since_waking + | |
| PSQI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| #BAI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| ##interaction with GW slope | |
| caseVScontrol:gestage_weeks_cent + | |
| #Maternal_Age_Years_cent:gestage_weeks_cent + | |
| #Maternal_Education:gestage_weeks_cent + | |
| Parity:gestage_weeks_cent + | |
| #Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent + | |
| #Weight_Gain_cent:gestage_weeks_cent + | |
| #Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent + | |
| #Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent + | |
| #Maternal_Smoking_During_Pregnancy:gestage_weeks_cent + | |
| #Gestational_Age_Weeks_cent:gestage_weeks_cent + | |
| Mother_Cohabiting:gestage_weeks_cent | |
| M.Baseline_CAR_cort.COV_wPSQI <- lme(Baseline_CAR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_CAR_cort.COV_wPSQI) | |
| RES_CAR_cov.select.PSQI <- model_res(summary(M.Baseline_CAR_cort.COV_wPSQI)) | |
| ##BAI | |
| Baseline_CAR_Form.COV <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Maternal_Education + | |
| Parity + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| Weight_Gain_cent + | |
| #Maternal_Hypertensive_Disorders_anyVSnone + | |
| Maternal_Diabetes_Disorders_anyVSnone + | |
| #Maternal_Smoking_During_Pregnancy + | |
| #Gestational_Age_Weeks_cent + | |
| Mother_Cohabiting + | |
| #PSQI_qclosest_to_cortGW_pregMean + | |
| BAI_qclosest_to_cortGW_pregMean + | |
| #PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| BAI_qclosest_to_cortGW_pregMeanDeviat + | |
| ##interaction with time slope | |
| #caseVScontrol:hours_since_waking + | |
| #Maternal_Age_Years_cent:hours_since_waking + | |
| Maternal_Education:hours_since_waking + | |
| #Parity:hours_since_waking + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| #Weight_Gain_cent:hours_since_waking + | |
| #Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| #Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| #Gestational_Age_Weeks_cent:hours_since_waking + | |
| #Mother_Cohabiting:hours_since_waking + | |
| #PSQI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| BAI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| ##interaction with GW slope | |
| caseVScontrol:gestage_weeks_cent + | |
| #Maternal_Age_Years_cent:gestage_weeks_cent + | |
| #Maternal_Education:gestage_weeks_cent + | |
| Parity:gestage_weeks_cent + | |
| #Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent + | |
| #Weight_Gain_cent:gestage_weeks_cent + | |
| #Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent + | |
| #Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent + | |
| #Maternal_Smoking_During_Pregnancy:gestage_weeks_cent + | |
| #Gestational_Age_Weeks_cent:gestage_weeks_cent + | |
| Mother_Cohabiting:gestage_weeks_cent | |
| M.Baseline_CAR_cort.COV_wBAI <- lme(Baseline_CAR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_CAR_cort.COV_wBAI) | |
| RES_CAR_cov.select.BAI <- model_res(summary(M.Baseline_CAR_cort.COV_wBAI)) | |
| ### (4.1) Hypothesis Testing: BS and WS effects #### | |
| CAR_CES_D <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season + | |
| Cesd_qclosest_to_cortGW_pregMean + | |
| Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| M1_CAR_CES_D <- lme(CAR_CES_D, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M1_CAR_CES_D) | |
| RES_CAR_CESD_M1 <- model_res(summary(M1_CAR_CES_D)) | |
| RES_CAR_CESD_M1 <- mutate_all(RES_CAR_CESD_M1,funs(replace(., is.na(.), ""))) | |
| #explained BS variance by depressive symptoms | |
| summary(M1_CAR_CES_D) #N=512, obs = 2812 --> 5.49 obs per participant on average | |
| summary(M.Baseline_Morning_cort) #N=514, obs = 2860 --> 5.56 obs per participant on average | |
| BS1 <- as.numeric(VarCorr(M1_CAR_CES_D)[2,1]) #with CES-D | |
| resid1 <- as.numeric(VarCorr(M1_CAR_CES_D)[6,1]) | |
| BS2 <- as.numeric(VarCorr(M.Baseline_Morning_cort)[2,1]) #without CES-D | |
| resid2 <- as.numeric(VarCorr(M.Baseline_Morning_cort)[6,1]) | |
| CAR_added_BS_var <- format(round(1 - ((BS1^2 + (resid1^2/5.49))/(BS2^2 + resid2^2/5.56)),2)) | |
| ## Model 2: including maternal characteristics as covariates #### | |
| CAR_CES_D2 <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| I(gestage_weeks_cent^2) + | |
| hours_since_waking:gestage_weeks_cent + | |
| season + | |
| caseVScontrol + | |
| Maternal_Education + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| Weight_Gain_cent + | |
| #Maternal_Diabetes_Disorders_anyVSnone + | |
| Mother_Cohabiting + | |
| ##interaction with time slope | |
| Maternal_Education:hours_since_waking + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| #Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| ##interaction with GW slope | |
| Mother_Cohabiting:gestage_weeks_cent + | |
| #caseVScontrol:gestage_weeks_cent + | |
| Cesd_qclosest_to_cortGW_pregMean + | |
| Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| M2_CAR_CES_D <- lme(CAR_CES_D2, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M2_CAR_CES_D) | |
| RES_CAR_CESD_M2 <- model_res(summary(M2_CAR_CES_D)) | |
| RES_CAR_CESD_M2 <- mutate_all(RES_CAR_CESD_M2,funs(replace(., is.na(.), ""))) | |
| # ## Model 3: including sleep and anxiety as coviates #### not included as | |
| # neither BAI nor PSQI showed significant effect on cortisol concentrations | |
| # CAR_CES_D3 <- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # caseVScontrol + | |
| # Maternal_Education + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| # Weight_Gain_cent + | |
| # Maternal_Diabetes_Disorders_anyVSnone + | |
| # Mother_Cohabiting + | |
| # ##interaction with time slope | |
| # Maternal_Education:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| # ##interaction with GW slope | |
| # Mother_Cohabiting:gestage_weeks_cent + | |
| # Cesd_qclosest_to_cortGW_pregMean + | |
| # Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent + | |
| # PSQI_qclosest_to_cortGW_pregMean + | |
| # PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| # BAI_qclosest_to_cortGW_pregMean + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat | |
| # M3_CAR_CES_D <- lme(CAR_CES_D3, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M3_CAR_CES_D) | |
| # RES_CAR_CESD_M3 <- model_res(summary(M3_CAR_CES_D)) | |
| # RES_CAR_CESD_M3 <- mutate_all(RES_CAR_CESD_M3,funs(replace(., is.na(.), ""))) | |
| RES_CAR_CESD <- rbindlist(list(RES_CAR_CESD_M1,RES_CAR_CESD_M2)) | |
| #RES_CAR_CESD[[c(14), "predictor"]] <- "Random Effects" | |
| #RES_CAR_CESD[[c(47), "predictor"]] <- "Random Effects" | |
| #RES_CAR_CESD[[c(83), "predictor"]] <- "Random Effects" | |
| RES_CAR_CESD$p <- as.numeric(RES_CAR_CESD$p) | |
| #to display in main part of thesis | |
| Cesd_rows <- c(as.character(grep("Cesd_", c(as.character(RES_CAR_CESD$predictor))))) | |
| random_part_rows <- c(rownames(RES_CAR_CESD)[c(14:20,47:53)]) #,83:89 | |
| RES_CAR_CESD_sub <- subset(RES_CAR_CESD, rownames(RES_CAR_CESD) %in% Cesd_rows | rownames(RES_CAR_CESD) %in% random_part_rows | p <= 0.05) | |
| RES_CAR_CESD_sub <- mutate_all(RES_CAR_CESD_sub,funs(replace(., is.na(.), ""))) | |
| diurnal_factors <- c(rownames(RES_CAR_CESD_sub)[c(1:5,8,18:23,30)]) #46:50 | |
| RES_CAR_CESD_sub <- subset(RES_CAR_CESD_sub, !(rownames(RES_CAR_CESD_sub) %in% diurnal_factors)) | |
| ## explained BS variance by adding CES-D, and when adding all covariates | |
| # model fit comparison with and without depressive symptoms | |
| # CAR_CES_D <- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Cesd_qclosest_to_cortGW_pregMean + | |
| # Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| # | |
| # M1_CAR_CES_D <- lme(CAR_CES_D, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "ML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # CAR_M1 <- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # I(gestage_weeks_cent^2) + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season | |
| # | |
| # M1_CAR <- lme(CAR_M1, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "ML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # anova(M1_CAR, M1_CAR_CES_D) | |
| # summary(M1_CAR) #3120.991 | |
| # summary(M1_CAR_CES_D) #3074.746s | |
| ### Sensitivity Analysis: Part I #### | |
| ##M1 | |
| M1_CAR_CES_D_sub1 <- lme(CAR_CES_D, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M1_CAR_CES_D_sub1) | |
| AppenB1_RES_CAR_CESD_M1 <- model_res(summary(M1_CAR_CES_D_sub1)) | |
| AppenB1_RES_CAR_CESD_M1 <- mutate_all(AppenB1_RES_CAR_CESD_M1,funs(replace(., is.na(.), ""))) | |
| ##M2 | |
| M2_CAR_CES_D_sub1 <- lme(CAR_CES_D2, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M2_CAR_CES_D_sub1) | |
| AppenB1_RES_CAR_CESD_M2 <- model_res(summary(M2_CAR_CES_D_sub1)) | |
| AppenB1_RES_CAR_CESD_M2 <- mutate_all(AppenB1_RES_CAR_CESD_M2,funs(replace(., is.na(.), ""))) | |
| # ##M3 | |
| # M3_CAR_CES_D_sub1 <- lme(CAR_CES_D3, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M3_CAR_CES_D_sub1) | |
| # AppenB1_RES_CAR_CESD_M3 <- model_res(summary(M3_CAR_CES_D_sub1)) | |
| # AppenB1_RES_CAR_CESD_M3 <- mutate_all(AppenB1_RES_CAR_CESD_M3,funs(replace(., is.na(.), ""))) | |
| #save all in one df | |
| AppenB1_RES_CAR <- rbindlist(list(AppenB1_RES_CAR_CESD_M1, | |
| AppenB1_RES_CAR_CESD_M2)) | |
| ## Sensitivity Analysis: Part II #### | |
| CAR_df_no_Meads <- subset(CAR_df, ITUbroadpsychiatricmedication_18_KELA == "no") | |
| ##M1 | |
| M1_CAR_CES_D_sub2 <- lme(CAR_CES_D, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df_no_Meads, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M1_CAR_CES_D_sub2) | |
| AppenB2_RES_CAR_CESD_M1 <- model_res(summary(M1_CAR_CES_D_sub2)) | |
| AppenB2_RES_CAR_CESD_M1 <- mutate_all(AppenB2_RES_CAR_CESD_M1,funs(replace(., is.na(.), ""))) | |
| ##M2 | |
| M2_CAR_CES_D_sub2 <- lme(CAR_CES_D2, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = CAR_df_no_Meads, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M2_CAR_CES_D_sub2) | |
| AppenB2_RES_CAR_CESD_M2 <- model_res(summary(M2_CAR_CES_D_sub2)) | |
| AppenB2_RES_CAR_CESD_M2 <- mutate_all(AppenB2_RES_CAR_CESD_M2,funs(replace(., is.na(.), ""))) | |
| ##M3 | |
| # M3_CAR_CES_D_sub2 <- lme(CAR_CES_D3, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df_no_Meads, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M3_CAR_CES_D_sub2) | |
| # AppenB2_RES_CAR_CESD_M3 <- model_res(summary(M3_CAR_CES_D_sub2)) | |
| #AppenB2_RES_CAR_CESD_M3 <- mutate_all(AppenB2_RES_CAR_CESD_M3,funs(replace(., is.na(.), ""))) | |
| #save all in one df | |
| AppenB2_RES_CAR <- rbindlist(list(AppenB2_RES_CAR_CESD_M1, | |
| AppenB2_RES_CAR_CESD_M2)) | |
| ### (4.2) Hypothesis Testing: Clinically relevant symptoms #### | |
| # CAR_CES_D <- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # Depressive_Symptom_Severity + | |
| # Depressive_Symptom_Severity:hours_since_waking + | |
| # Depressive_Symptom_Severity:gestage_weeks_cent + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| # | |
| # M1_CAR_CES_D_cat <- lme(CAR_CES_D, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M1_CAR_CES_D_cat) | |
| # RES_CAR_CESD_M1_cat <- summary(M1_CAR_CES_D) | |
| # RES_CAR_CESD_M1_cat <- RES_CAR_CESD_M1$tTable | |
| # | |
| # ## repeat in subscample #### | |
| # M1_CAR_CES_D_sub1_cat <- lme(CAR_CES_D, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M1_CAR_CES_D_sub1_cat) | |
| # # RES_CAR_CESD_M1 <- summary(M1_CAR_CES_D_sub1) | |
| # # RES_CAR_CESD_M1 <- M1_CAR_CES_D_sub1$tTable | |
| # | |
| # ##including maternal characteristics as covariates | |
| # CAR_CES_D2 <- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # caseVScontrol + | |
| # Maternal_Education + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| # Weight_Gain_cent + | |
| # Maternal_Diabetes_Disorders_anyVSnone + | |
| # Mother_Cohabiting + | |
| # ##interaction with time slope | |
| # Maternal_Education:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| # ##interaction with GW slope | |
| # Mother_Cohabiting:gestage_weeks_cent + | |
| # Depressive_Symptom_Severity + | |
| # Depressive_Symptom_Severity:hours_since_waking + | |
| # Depressive_Symptom_Severity:gestage_weeks_cent + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| # | |
| # M2_CAR_CES_D_cat <- lme(CAR_CES_D2, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M2_CAR_CES_D_cat) | |
| # # RES_CAR_CESD_M2 <- summary(M2_CAR_CES_D_cat) | |
| # # RES_CAR_CESD_M2 <- RES_CAR_CESD_M1$tTable | |
| # | |
| # ## repeat in subsample 1 #### | |
| # M2_CAR_CES_D_sub1 <- lme(CAR_CES_D2, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M2_CAR_CES_D_sub1) | |
| # | |
| # ## including sleep and anxiety as coviates | |
| # ##including maternal characteristics as covariates | |
| # CAR_CES_D3 <- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # season + | |
| # caseVScontrol + | |
| # Maternal_Education + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| # Weight_Gain_cent + | |
| # Maternal_Diabetes_Disorders_anyVSnone + | |
| # Mother_Cohabiting + | |
| # ##interaction with time slope | |
| # Maternal_Education:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| # ##interaction with GW slope | |
| # Mother_Cohabiting:gestage_weeks_cent + | |
| # Depressive_Symptom_Severity + | |
| # Depressive_Symptom_Severity:hours_since_waking + | |
| # Depressive_Symptom_Severity:gestage_weeks_cent + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent + | |
| # PSQI_qclosest_to_cortGW_pregMean + | |
| # PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| # BAI_qclosest_to_cortGW_pregMean + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat | |
| # | |
| # M3_CAR_CES_D_cat <- lme(CAR_CES_D3, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M3_CAR_CES_D_cat) | |
| # # RES_CAR_CESD_M3 <- summary(M3_CAR_CES_D) | |
| # # RES_CAR_CESD_M3 <- RES_CAR_CESD_M3$tTable | |
| # | |
| # ## repeat in subscample | |
| # M3_CAR_CES_D_cat <- lme(CAR_CES_D3, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = CAR_df[CAR_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M3_CAR_CES_D_cat) | |
| ################################################################################# | |
| ####### Diurnal Model #### | |
| DIURNAL_times <- c("S1", "S4", "S5", "S6", "S7") | |
| DIURNAL_df <- subset(final_long_cort, timeslot %in% DIURNAL_times) | |
| ### Baseline Model: Random effects and variance-covariance structure #### | |
| #Including diurnal characteristics | |
| DIURNAL_Form <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent | |
| DIURNAL_Null <- lg_cort_concentration_ug_l ~ 1 | |
| ## Random Intercept only model | |
| M.DIUR.RI <- lme(DIURNAL_Form, | |
| random = list(~ 1|participantID, | |
| ~ 1|pregstage), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| #plot(M.DIUR.RI) | |
| #qqnorm(M.DIUR.RI, abline = c(0,1)) | |
| ## Does adding ELISA improve model fit? | |
| M.DIUR.RI_wELISA <- update(M.DIUR.RI, random = list(~1|ELISA_analysis_plate, | |
| ~1|participantID, | |
| ~1|pregstage)) | |
| tab1 <- anova(M.DIUR.RI,M.DIUR.RI_wELISA) #no improvement, in fact non-sig less model fit | |
| ## Does adding a random slope for time since waking improve model fit? | |
| #i.e., allowing the diurnal decline to vary across participants | |
| # M.DIUR.RS1 <- update(M.DIUR.RI, random = list(~hours_since_waking|participantID, | |
| # ~1|pregstage)) #singular convergence | |
| #Does adding heterogenous WS variance for the effect of time improve fit? | |
| M.DIUR.hetero <- lme(DIURNAL_Form, | |
| random = list(~1|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| anova(M.DIUR.RI, M.DIUR.hetero) #Yes, it does improve fit! | |
| #plot(M.DIUR.hetero) | |
| ### Can the slope for cortisol decline now be added without convergence issues | |
| #Does adding heterogenous WS variance for the effect of time improve fit? | |
| M.DIUR.hetero2 <- lme(DIURNAL_Form, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| anova(M.DIUR.RI, M.DIUR.hetero, M.DIUR.hetero2) #Yes, it does improve fit! | |
| #plot(M.DIUR.hetero) | |
| #Does adding autocorrelation for WS variance for the effect of time improve fit? | |
| M.DIUR.hetero_auto <- lme(DIURNAL_Form, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| tab2 <- anova(M.DIUR.RI, M.DIUR.hetero, M.DIUR.hetero2,M.DIUR.hetero_auto) #Yes, it does improve fit! | |
| RES_DIUR_Mselect <- rbind.data.frame(tab1, tab2) | |
| # plot(M.DIUR.hetero) | |
| # qqnorm(M.DIUR.hetero, abline = c(0,1)) | |
| ### Baseline Model: Fixed Component #### | |
| M.DIUR_ML1 <- lme(DIURNAL_Form, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "ML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| ## check whether quadratic term of time since waking indeed improves model fit | |
| # M.DIUR_ML2 <- update(M.DIUR_ML1, .~. -I(hours_since_waking^2)) #leads to singular convergence | |
| # anova(M.DIUR_ML1, M.DIUR_ML2) #worse fit to exclude quadratic term | |
| # ## adding a quadratic term for GW | |
| # M.DIUR_ML3 <- update(M.DIUR_ML1, .~. +I(gestage_weeks_cent^2)) | |
| # anova(M.DIUR_ML1, M.DIUR_ML3) #no sig. improvement | |
| ## adding a waking time | |
| M.DIUR_ML4 <- update(M.DIUR_ML1, .~. + | |
| waking_time + | |
| waking_time:hours_since_waking + | |
| season + | |
| season:hours_since_waking) | |
| car::Anova(M.DIUR_ML4) | |
| ### Baseline Model: Model Validation #### | |
| Baseline_DIUR_Form <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent + | |
| waking_time + | |
| waking_time:hours_since_waking + | |
| season + | |
| season:hours_since_waking | |
| M.Baseline_DIUR_cort <- lme(Baseline_DIUR_Form, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| # plot(M.Baseline_DIUR_cort) | |
| # qqnorm(M.Baseline_DIUR_cort, abline = c(0,1)) #overall normal distirbution of residuals can be assumed! | |
| # qqnorm(M.Baseline_DIUR_cort, ~ranef(., level=1), abline=c(0,1)) #random effects not normal! | |
| # qqnorm(M.Baseline_DIUR_cort, ~ranef(., level=2), abline=c(0,1)) #random effects not normal! | |
| DIUR_parameters <- as.data.frame(coef(M.Baseline_DIUR_cort)) | |
| save(DIUR_parameters, file = "Rdata/Diurnal_Parameters.Rdata") | |
| ### (1) Cortisol Across Gestation #### | |
| car::Anova(M.Baseline_DIUR_cort) | |
| RES_DIUR_Baseline <- model_res(summary(M.Baseline_DIUR_cort)) | |
| RES_DIUR_Baseline <- mutate_all(RES_DIUR_Baseline,funs(replace(., is.na(.), ""))) | |
| ## Total Variance explained | |
| varests <- as.numeric(VarCorr(M.Baseline_DIUR_cort)[(c(2:3,5:6)), 1]) | |
| RES_DIUR_total_R_explained <- format(round(sum(varests[1:3])/sum(varests)*100,1)) #R^2=23.9% | |
| DIUR_residual <- sum(varests[4])/sum(varests) | |
| ### (2) Covariates: Effects of maternal characteristics #### | |
| Baseline_DIUR_Form.COV <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent + | |
| waking_time + | |
| waking_time:hours_since_waking + | |
| season + | |
| season:hours_since_waking + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Maternal_Education + | |
| Parity + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| Weight_Gain_cent + | |
| Maternal_Hypertensive_Disorders_anyVSnone + | |
| Maternal_Diabetes_Disorders_anyVSnone + | |
| Maternal_Smoking_During_Pregnancy + | |
| Gestational_Age_Weeks_cent + | |
| Mother_Cohabiting + | |
| #interaction with time slope | |
| caseVScontrol:hours_since_waking + | |
| Maternal_Age_Years_cent:hours_since_waking + | |
| Maternal_Education:hours_since_waking + | |
| Parity:hours_since_waking + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| Weight_Gain_cent:hours_since_waking + | |
| Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| Gestational_Age_Weeks_cent:hours_since_waking + | |
| Mother_Cohabiting:hours_since_waking + | |
| #interaction with GW slope | |
| caseVScontrol:gestage_weeks_cent + | |
| Maternal_Age_Years_cent:gestage_weeks_cent + | |
| Maternal_Education:gestage_weeks_cent + | |
| Parity:gestage_weeks_cent + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent + | |
| Weight_Gain_cent:gestage_weeks_cent + | |
| Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent + | |
| Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent + | |
| Maternal_Smoking_During_Pregnancy:gestage_weeks_cent + | |
| Gestational_Age_Weeks_cent:gestage_weeks_cent + | |
| Mother_Cohabiting:gestage_weeks_cent | |
| M.Baseline_DIUR_cort.COV <- lme(Baseline_DIUR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_DIUR_cort.COV) | |
| # ### instead investigating each covariates seperately #### | |
| # ## Case VS Control #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # caseVScontrol + | |
| # caseVScontrol:hours_since_waking + | |
| # caseVScontrol:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.CoCa <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.CoCa) | |
| # | |
| # # Repeat in subsample #### | |
| # M.Baseline_DIUR_cort.COV.CoCa <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.CoCa) | |
| # | |
| # | |
| # ## Maternal Age #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Maternal_Age_Years_cent + | |
| # Maternal_Age_Years_cent:hours_since_waking + | |
| # Maternal_Age_Years_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.age <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.age) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.age <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.age) | |
| # | |
| # ## maternal education #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Maternal_Education + | |
| # Maternal_Education:hours_since_waking + | |
| # Maternal_Education:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.edu <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.edu) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.edu <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.edu) | |
| # | |
| # ## Parity #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Parity + | |
| # Parity:hours_since_waking + | |
| # Parity:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.par <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.par) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.par <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.par) | |
| # | |
| # ## BMI #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.bmi <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.bmi) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.bmi <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.bmi) | |
| # | |
| # ## Weight Gain #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Weight_Gain_cent + | |
| # Weight_Gain_cent:hours_since_waking + | |
| # Weight_Gain_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.wg <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.wg) | |
| # | |
| # # repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.wg <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.wg) | |
| # ## Hypertensive disorders #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Maternal_Hypertensive_Disorders_anyVSnone + | |
| # Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| # Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.hyp <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.hyp) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.hyp <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.hyp) | |
| # | |
| # ## Diabetes Disorders #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone + | |
| # Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.diab <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.diab) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.diab <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.diab) | |
| # | |
| # ## Cohabitation #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Mother_Cohabiting + | |
| # Mother_Cohabiting:hours_since_waking + | |
| # Mother_Cohabiting:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.Co <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.Co) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.Co <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.Co) | |
| # | |
| # ## Smoking #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Maternal_Smoking_During_Pregnancy + | |
| # Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| # Maternal_Smoking_During_Pregnancy:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.Smoke <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.Smoke) | |
| # | |
| # # repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.Smoke <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.Smoke) | |
| # ## Gestational Length #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Gestational_Age_Weeks_cent + | |
| # Gestational_Age_Weeks_cent:hours_since_waking + | |
| # Gestational_Age_Weeks_cent:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # | |
| # ## Child Sex #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # Child_Sex + | |
| # Child_Sex:hours_since_waking + | |
| # Child_Sex:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # | |
| # ## BAI #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # BAI_qclosest_to_cortGW_pregMean + | |
| # BAI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| # BAI_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # BAI_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # | |
| # ## PSQI #### | |
| # Baseline_DIUR_Form.COV<- lg_cort_concentration_ug_l ~ | |
| # hours_since_waking + | |
| # I(hours_since_waking^2) + | |
| # gestage_weeks_cent + | |
| # hours_since_waking:gestage_weeks_cent + | |
| # waking_time + | |
| # waking_time:hours_since_waking + | |
| # season + | |
| # season:hours_since_waking + | |
| # PSQI_qclosest_to_cortGW_pregMean + | |
| # PSQI_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| # PSQI_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| # PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| # PSQI_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| # PSQI_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| # | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df, | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| # #repeat in sub #### | |
| # M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| # random = list(~(hours_since_waking)|participantID, | |
| # ~1|pregstage), | |
| # weights = varExp(form = ~hours_since_waking|pregstage), | |
| # correlation = corCAR1(form = ~ (hours_since_waking)), | |
| # method = "REML", | |
| # data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| # na.action = "na.omit", | |
| # control= lmeControl(msMaxIter = 200)) | |
| # car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| ### (3) Final Covariate Selection #### | |
| Baseline_DIUR_Form.COV <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent + | |
| waking_time + | |
| waking_time:hours_since_waking + | |
| season + | |
| season:hours_since_waking + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Maternal_Education + | |
| Parity + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent + | |
| Child_Sex + | |
| #Weight_Gain_cent + | |
| #Maternal_Hypertensive_Disorders_anyVSnone + | |
| Maternal_Diabetes_Disorders_anyVSnone + | |
| #Maternal_Smoking_During_Pregnancy + | |
| Gestational_Age_Weeks_cent + | |
| #Mother_Cohabiting + | |
| PSQI_qclosest_to_cortGW_pregMean + | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| ##interaction with time slope | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| caseVScontrol:hours_since_waking + | |
| #Maternal_Age_Years_cent:hours_since_waking + | |
| Maternal_Education:hours_since_waking + | |
| Parity:hours_since_waking + | |
| Child_Sex:hours_since_waking + | |
| # Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:hours_since_waking + | |
| # Weight_Gain_cent:hours_since_waking + | |
| # Maternal_Hypertensive_Disorders_anyVSnone:hours_since_waking + | |
| # Maternal_Diabetes_Disorders_anyVSnone:hours_since_waking + | |
| # Maternal_Smoking_During_Pregnancy:hours_since_waking + | |
| # Gestational_Age_Weeks_cent:hours_since_waking + | |
| # Mother_Cohabiting:hours_since_waking + | |
| ##interaction with GW slope | |
| # caseVScontrol:gestage_weeks_cent + | |
| # Maternal_Age_Years_cent:gestage_weeks_cent + | |
| # Maternal_Education:gestage_weeks_cent + | |
| # Parity:gestage_weeks_cent + | |
| Maternal_Body_Mass_Index_in_Early_Pregnancy_cent:gestage_weeks_cent + | |
| # Weight_Gain_cent:gestage_weeks_cent + | |
| # Maternal_Hypertensive_Disorders_anyVSnone:gestage_weeks_cent + | |
| Maternal_Diabetes_Disorders_anyVSnone:gestage_weeks_cent | |
| # Maternal_Smoking_During_Pregnancy:gestage_weeks_cent + | |
| # Gestational_Age_Weeks_cent:gestage_weeks_cent + | |
| # Mother_Cohabiting:gestage_weeks_cent | |
| M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| RES_DIUR_cov.select <- model_res(summary(M.Baseline_DIUR_cort.COV.gl)) | |
| RES_DIUR_cov.select <- mutate_all(RES_DIUR_cov.select,funs(replace(., is.na(.), ""))) | |
| #repeat in sub #### | |
| M.Baseline_DIUR_cort.COV.gl <- lme(Baseline_DIUR_Form.COV, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M.Baseline_DIUR_cort.COV.gl) | |
| ### (4) Hypothesis Testing #### | |
| DIUR_CES_D <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent + | |
| waking_time + | |
| waking_time:hours_since_waking + | |
| season + | |
| season:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean + | |
| Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| M1_DIUR_CES_D <- lme(DIUR_CES_D, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M1_DIUR_CES_D) | |
| RES_DIUR_CESD_M1 <- model_res(summary(M1_DIUR_CES_D)) | |
| RES_DIUR_CESD_M1 <- mutate_all(RES_DIUR_CESD_M1,funs(replace(., is.na(.), ""))) | |
| #explained BS variance by depressive symptoms | |
| summary(M1_DIUR_CES_D) #N=512, obs = 4574 --> 8.93 obs per participant on average | |
| summary(M.Baseline_DIUR_cort) #N=514, obs = 4656 --> 9.06 obs per participant on average | |
| BS1 <- as.numeric(VarCorr(M1_DIUR_CES_D)[2,1]) #with CES-D | |
| resid1 <- as.numeric(VarCorr(M1_DIUR_CES_D)[6,1]) | |
| BS2 <- as.numeric(VarCorr(M.Baseline_DIUR_cort)[2,1]) #without CES-D | |
| resid2 <- as.numeric(VarCorr(M.Baseline_DIUR_cort)[6,1]) | |
| DIUR_added_BS_var <- format(round(1 - ((BS1^2 + (resid1^2/8.93))/(BS2^2 + resid2^2/9.06)),3)) | |
| ## M2 adding Covariates | |
| DIUR_CES_D2 <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent + | |
| waking_time + | |
| season + | |
| season:hours_since_waking + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Parity + | |
| Gestational_Age_Weeks_cent + | |
| Child_Sex + | |
| ##interaction with time slope | |
| caseVScontrol:hours_since_waking + | |
| Parity:hours_since_waking + | |
| Child_Sex:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean + | |
| Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| M2_DIUR_CES_D <- lme(DIUR_CES_D2, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M2_DIUR_CES_D) | |
| RES_DIUR_CESD_M2 <- model_res(summary(M2_DIUR_CES_D)) | |
| RES_DIUR_CESD_M2 <- mutate_all(RES_DIUR_CESD_M2,funs(replace(., is.na(.), ""))) | |
| ## M3 Adding other questionnaires | |
| DIUR_CES_D3 <- lg_cort_concentration_ug_l ~ | |
| hours_since_waking + | |
| I(hours_since_waking^2) + | |
| gestage_weeks_cent + | |
| hours_since_waking:gestage_weeks_cent + | |
| waking_time + | |
| season + | |
| season:hours_since_waking + | |
| caseVScontrol + | |
| Maternal_Age_Years_cent + | |
| Parity + | |
| Gestational_Age_Weeks_cent + | |
| Child_Sex + | |
| PSQI_qclosest_to_cortGW_pregMean + | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat + | |
| ##interaction with time slope | |
| PSQI_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| caseVScontrol:hours_since_waking + | |
| Parity:hours_since_waking + | |
| Child_Sex:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean + | |
| Cesd_qclosest_to_cortGW_pregMean:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMean:gestage_weeks_cent + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:hours_since_waking + | |
| Cesd_qclosest_to_cortGW_pregMeanDeviat:gestage_weeks_cent | |
| M3_DIUR_CES_D <- lme(DIUR_CES_D3, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| car::Anova(M3_DIUR_CES_D) | |
| RES_DIUR_CESD_M3 <- model_res(summary(M3_DIUR_CES_D)) | |
| RES_DIUR_CESD_M3 <- mutate_all(RES_DIUR_CESD_M3,funs(replace(., is.na(.), ""))) | |
| RES_DIUR_CESD <- rbindlist(list(RES_DIUR_CESD_M1,RES_DIUR_CESD_M2,RES_DIUR_CESD_M3)) | |
| RES_DIUR_CESD[[c(16), "predictor"]] <- "Random Effects" | |
| RES_DIUR_CESD[[c(45), "predictor"]] <- "Random Effects" | |
| RES_DIUR_CESD[[c(77), "predictor"]] <- "Random Effects" | |
| RES_DIUR_CESD$p <- as.numeric(RES_DIUR_CESD$p) | |
| #to display in main part of thesis | |
| Cesd_rows <- c(as.character(grep("Cesd_", c(as.character(RES_DIUR_CESD$predictor))))) | |
| random_part_rows <- c(rownames(RES_DIUR_CESD)[c(16:22, 45:51, 77:83)]) | |
| RES_DIUR_CESD_sub <- subset(RES_DIUR_CESD, rownames(RES_DIUR_CESD) %in% Cesd_rows | rownames(RES_DIUR_CESD) %in% random_part_rows | p <= 0.05) | |
| RES_DIUR_CESD_sub <- mutate_all(RES_DIUR_CESD_sub,funs(replace(., is.na(.), ""))) | |
| diurnal_factors <- c(rownames(RES_DIUR_CESD_sub)[c(1:6, 9:10, 22:27, 33:34, 47:52, 58:59)]) | |
| RES_DIUR_CESD_sub <- subset(RES_DIUR_CESD_sub, !(rownames(RES_DIUR_CESD_sub) %in% diurnal_factors)) | |
| ### Sensitivity analysis part I #### | |
| #M1 | |
| M1_DIUR_CES_D_sub1 <- lme(DIUR_CES_D, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| AppenB1_RES_DIUR_CESD_M1 <- model_res(summary(M1_DIUR_CES_D_sub1)) | |
| AppenB1_RES_DIUR_CESD_M1 <- mutate_all(AppenB1_RES_DIUR_CESD_M1,funs(replace(., is.na(.), ""))) | |
| #M2 | |
| M2_DIUR_CES_D_sub1 <- lme(DIUR_CES_D2, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| AppenB1_RES_DIUR_CESD_M2 <- model_res(summary(M2_DIUR_CES_D_sub1)) | |
| AppenB1_RES_DIUR_CESD_M2 <- mutate_all(AppenB1_RES_DIUR_CESD_M2,funs(replace(., is.na(.), ""))) | |
| #M3 | |
| M3_DIUR_CES_D_sub1 <- lme(DIUR_CES_D3, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIURNAL_df[DIURNAL_df$participantID %in% at_least_two_assessments,], | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| AppenB1_RES_DIUR_CESD_M3 <- model_res(summary(M3_DIUR_CES_D_sub1)) | |
| AppenB1_RES_DIUR_CESD_M3 <- mutate_all(AppenB1_RES_DIUR_CESD_M3,funs(replace(., is.na(.), ""))) | |
| AppenB1_RES_DIUR_CESD <- rbindlist(list(AppenB1_RES_DIUR_CESD_M1, | |
| AppenB1_RES_DIUR_CESD_M2, | |
| AppenB1_RES_DIUR_CESD_M3)) | |
| ### Sensitivity analysis part II #### | |
| DIUR_no_meds <- subset(DIURNAL_df, !(ITUbroadpsychiatricmedication_18_KELA == "yes")) | |
| #M1 | |
| M1_DIUR_CES_D_sub2 <- lme(DIUR_CES_D, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIUR_no_meds, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| AppenB2_RES_DIUR_CESD_M1 <- model_res(summary(M1_DIUR_CES_D_sub2)) | |
| AppenB2_RES_DIUR_CESD_M1 <- mutate_all(AppenB2_RES_DIUR_CESD_M1,funs(replace(., is.na(.), ""))) | |
| #M2 | |
| M2_DIUR_CES_D_sub2 <- lme(DIUR_CES_D2, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIUR_no_meds, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| AppenB2_RES_DIUR_CESD_M2 <- model_res(summary(M2_DIUR_CES_D_sub2)) | |
| AppenB2_RES_DIUR_CESD_M2 <- mutate_all(AppenB2_RES_DIUR_CESD_M2,funs(replace(., is.na(.), ""))) | |
| #M3 | |
| M3_DIUR_CES_D_sub2 <- lme(DIUR_CES_D3, | |
| random = list(~(hours_since_waking)|participantID, | |
| ~1|pregstage), | |
| weights = varExp(form = ~hours_since_waking|pregstage), | |
| correlation = corCAR1(form = ~ (hours_since_waking)), | |
| method = "REML", | |
| data = DIUR_no_meds, | |
| na.action = "na.omit", | |
| control= lmeControl(msMaxIter = 200)) | |
| AppenB2_RES_DIUR_CESD_M3 <- model_res(summary(M3_DIUR_CES_D_sub2)) | |
| AppenB2_RES_DIUR_CESD_M3 <- mutate_all(AppenB2_RES_DIUR_CESD_M3,funs(replace(., is.na(.), ""))) | |
| #save results all in one df | |
| AppenB2_RES_DIUR_CESD <- rbindlist(list(AppenB2_RES_DIUR_CESD_M1, | |
| AppenB2_RES_DIUR_CESD_M2, | |
| AppenB2_RES_DIUR_CESD_M3)) | |
| ### Figures #### | |
| ## Plot diurnal cortisol curve #### | |
| final_long_cort$pregstage <- factor(final_long_cort$pregstage, | |
| levels = c("I", "II", "III"), | |
| labels = c("T1", "T2", "T3")) | |
| final_long_cort$Antenatal_Care <- factor(final_long_cort$caseVScontrol, | |
| levels = c("case", "control"), | |
| labels = c("Referral", "Routine")) | |
| v <- final_long_cort %>% | |
| group_by(pregstage, timeslot, Antenatal_Care) %>% | |
| summarise( | |
| n = n(), | |
| c = mean(cort_concentration_ug_l, na.rm = T), | |
| SE = sd(cort_concentration_ug_l, na.rm = T)/sqrt(n)) #, | |
| #t = mean(as.numeric(time_of_day_in_hours_since_midnight), na.rm = T)) | |
| library(jtools) #APA theme for ggplot | |
| #plot raw data per pregstage | |
| cortisol_by_pregstage <- ggplot(na.omit(v), | |
| aes(x = timeslot, | |
| y = c, | |
| color = Antenatal_Care)) + | |
| geom_line(aes(group = Antenatal_Care)) + | |
| geom_line(data = na.omit(v)) + | |
| geom_point() + | |
| facet_grid(. ~ pregstage) + | |
| geom_errorbar(aes(ymin=c-SE, ymax=c+SE), width=.1) + | |
| #jtools::theme_apa() + | |
| #theme_classic() + | |
| theme_bw() + | |
| theme(axis.line = element_line(colour = "black"), | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| #panel.border = element_blank(), | |
| panel.background = element_blank()) + | |
| xlab("\n Time of Sampling") + | |
| ylab("Cortisol Concentration (ug/l)\n") + | |
| theme(legend.position = "none") + | |
| scale_x_discrete(labels = c("7:15", | |
| "+15", | |
| "+30", | |
| "10:30", | |
| "12:00", | |
| "17:00", | |
| "23:05")) | |
| #mean time at awakening | |
| mean(final_long_cort$time_of_day_in_hours_since_midnight[final_long_cort$timeslot == "S1"], na.rm=T) | |
| mean(final_long_cort$time_of_day_in_hours_since_midnight[final_long_cort$timeslot == "S7"], na.rm=T) | |
| # library(export) | |
| # graph2ppt (cortisol_by_pregstage, file = "myGraph.pptx") | |
| # | |
| # # Create a new powerpoint document | |
| # xl <- read_xlsx() | |
| # xl <- add_sheet(xl,label = "LAB") | |
| # doc = addPlot( doc = doc , fun = print, x = cortisol_by_pregstage) | |
| # ## Plot morning concentrations by CES-D #### | |
| # CAR_df <- CAR_df%>% | |
| # group_by(participantID) %>% | |
| # mutate(mean_cort = mean(lg_cort_concentration_ug_l, na.rm=T)) | |
| # CAR_plot <- CAR_df[,c("participantID", "mean_cort","Depressive_Symptom_Severity")] | |
| # CAR_plot <- subset(CAR_plot, !duplicated(CAR_plot)) | |
| # | |
| # p <- ggplot(na.omit(CAR_plot), aes(x=Depressive_Symptom_Severity, y=mean_cort)) + | |
| # geom_boxplot() + | |
| # #facet_grid(.~pregstage) + | |
| # jtools::theme_apa() + | |
| # # theme_bw() + | |
| # # theme(axis.line = element_line(colour = "black"), | |
| # # panel.grid.major = element_blank(), | |
| # # panel.grid.minor = element_blank(), | |
| # # #panel.border = element_blank(), | |
| # # panel.background = element_blank()) + | |
| # ylab("Mean Morgning Cortisol Concentration [ln(ug/l)]") | |
| # p | |
| # ## depression severity by SD | |
| # CAR_df$CESD_Severity <- NA | |
| # SD <- sd(CAR_df$Cesd_qclosest_to_cortGW_pregMean, na.rm = T) | |
| # m <- mean(CAR_df$Cesd_qclosest_to_cortGW_pregMean, na.rm = T) | |
| # low_SD <- m-SD | |
| # high_SD <- m+SD | |
| # | |
| # for(i in 1:nrow(CAR_df)){ | |
| # cesd <- CAR_df$Cesd_qclosest_to_cortGW_pregMean[[i]] | |
| # if(!is.na(cesd)){ | |
| # if(cesd <= low_SD){ | |
| # CAR_df$CESD_Severity[[i]] <- -1 | |
| # } | |
| # if(cesd >= high_SD){ | |
| # CAR_df$CESD_Severity[[i]] <- 1 | |
| # } | |
| # } | |
| # } | |
| # | |
| # CAR_df$CESD_Severity <- factor(CAR_df$CESD_Severity) | |
| # | |
| # ##plot | |
| # CAR_df <- CAR_df%>% | |
| # group_by(participantID) %>% | |
| # mutate(mean_cort = mean(lg_cort_concentration_ug_l, na.rm=T)) | |
| # CAR_plot <- CAR_df[,c("participantID", "mean_cort","CESD_Severity")] | |
| # CAR_plot <- subset(CAR_plot, !duplicated(CAR_plot)) | |
| # | |
| # p <- ggplot(na.omit(CAR_plot), aes(x=CESD_Severity, y=mean_cort)) + | |
| # geom_boxplot() + | |
| # #facet_grid(.~pregstage) + | |
| # jtools::theme_apa() + | |
| # # theme_bw() + | |
| # # theme(axis.line = element_line(colour = "black"), | |
| # # panel.grid.major = element_blank(), | |
| # # panel.grid.minor = element_blank(), | |
| # # #panel.border = element_blank(), | |
| # # panel.background = element_blank()) + | |
| # ylab("Mean Morgning Cortisol Concentration [ln(ug/l)]") | |
| # | |
| # ### by pregstage values #### | |
| # CAR_df <- CAR_df%>% | |
| # group_by(participantID, pregstage) %>% | |
| # mutate(mean_cort = mean(lg_cort_concentration_ug_l, na.rm=T)) | |
| # CAR_plot <- CAR_df[,c("participantID", "pregstage","mean_cort","qclosest_based_clin_Cesd")] | |
| # CAR_plot$qclosest_based_clin_Cesd <- factor(CAR_plot$qclosest_based_clin_Cesd, | |
| # levels = c("non_clinical", "clinical"), | |
| # labels = c("below_clinical", "clinical")) | |
| # p <- ggplot(na.omit(CAR_plot), aes(x=qclosest_based_clin_Cesd, y=mean_cort)) + | |
| # geom_boxplot() + | |
| # facet_grid(.~pregstage) + | |
| # #jtools::theme_apa() + | |
| # theme_bw() + | |
| # theme(axis.line = element_line(colour = "black"), | |
| # panel.grid.major = element_blank(), | |
| # panel.grid.minor = element_blank(), | |
| # #panel.border = element_blank(), | |
| # panel.background = element_blank()) + | |
| # ylab("Mean Morgning Cortisol Concentration [ln(ug/l)]\n") + | |
| # xlab("\nCES-D score below versus above clinical cut-off") | |
| # p | |