--- title: "SHS Daily Resilience Stress and Affect Reactivity Paper" date: "08 Mar 2021" author(s): Natasha Tung (yan.tung@monash.edu), Joshua F. Wiley (joshua.wiley@monash.edu) --- ```{r knitroptions, include = FALSE} knitr::opts_chunk$set(error = TRUE) ``` # Background This document employs literate programming to include notes and documentation along with the actual source code used for analyses. Comments and notes can be typed throughout using plain text. Hash tags are used to create headings, with one hashtag creating a Level 1 heading, two hashtags a Level 2 heading, etc. `R` code is written in "chunks" special parts of the document separated by a triple back tick, curly braces and `r`. The end of a chunk is indicated by triple back ticks on their own line. # Load R Packages A number of R packages are used for data management and analyses. First if not already installed, uncomment and run the code to install the packages. That should only need to be done once (or when packages need to be updated). ```{r packages} ## ## uncomment and run these once only # install.packages("knitr", dependencies = TRUE) # install.packages("rmarkdown", dependencies = TRUE) # install.packages("pander", dependencies = TRUE) # install.packages("data.table", dependencies = TRUE) # install.packages("chron", dependencies = TRUE) # install.packages("bit64", dependencies = TRUE) # install.packages("devtools", dependencies = TRUE) # install.packages("extraoperators", dependencies = TRUE) # install.packages("ggplot2", dependencies = TRUE) # install.packages("cowplot", dependencies = TRUE) # install.packages("visreg", dependencies = TRUE) # install.packages("emmeans", dependencies = TRUE) # install.packages("iccMixed", dependencies = TRUE) # devtools::install_github("JWiley/JWileymisc") Sys.setenv(TZ = "Australia/Melbourne") library(rmarkdown) library(data.table) library(chron) library(extraoperators) library(ggplot2) library(cowplot) library(JWileymisc) library(lmerTest) library(multilevelTools) library(visreg) library(emmeans) library(tufte) library(haven) library(ggplot2) library(ggpubr) library(lme4) ``` # Load Data This code block reads the data into `R`. The baseline data are stored in an `R` object, `db`. ```{r, results = 'hide', cache=FALSE} ## set the base location for data if (Sys.info()[["user"]] %in% c("jwile")) { loc.data <- "g:/Shared drives/SHS Research Interns/Data" } else { loc.data <- "/Volumes/GoogleDrive/Shared drives/SHS Research Interns/Data" } ## read in data dd <- readRDS(file.path(loc.data, "shs_all.RDS")) ## create a binary after (1) / before (0) covid lock down start in Victoria ## dd[, COVID := as.numeric(BaselineStartDateTime > "2020-03-08")] ## if baseline prepost covid wanted dd[, COVID := as.numeric(SurveyDay > "2020-03-08")] ## binary variable if any stress source reported dd[, ANYSTRESSOR := as.integer(STRESSDUE_Count >= 1)] ## calculate which month each survey completed in as a covariate dd[, SurveyMonth := months(SurveyDay)] ``` # View Variable Names and Data ```{r, results = 'hide'} names(dd) View(dd) ``` # Descriptives for Demographic and Daily Variables ```{r} ### descriptives for demographic and daily ### egltable(c("Age", "Sex", "RELA", "SES_1","Race", "BornAUS", "TimeInMelb","LeaveHome", "SASH_Total", "LANG_1", "Student", "Employment", "Raised", "RFQ_Total", "BRS_Total", "ANX_TScore", "DEP_TScore", "PSS14_Total" ,"STRESS", "NegAff", "PosAff", "COVID"), data = dd [!duplicated(ID)], strict = FALSE) ### differences in demographics by group egltable(c("Age", "Sex", "RELA", "SES_1","Race", "BornAUS", "TimeInMelb","LeaveHome", "SASH_Total", "LANG_1", "Student", "Employment","Raised", "RFQ_Total", "BRS_Total", "ANX_TScore", "DEP_TScore", "PSS14_Total", "COVID"), g="Group", data = dd [!duplicated(ID)], strict = FALSE) ## continuous outcomes by group m <- lm(Age~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(SES_1~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(TimeInMelb~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(SASH_Total~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(RFQ_Total~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(BRS_Total~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(ANX_TScore~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(DEP_TScore~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") m <- lm(PSS14_Total~ Group, data = dd[!duplicated(ID)]) pairs(emmeans(m, ~ Group), adjust = "none") ## binary outcomes m <- glm(Sex~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(RELA~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(BornAUS~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(LeaveHome~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(LANG_1~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(Student~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(Employment~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(Raised~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") m <- glm(COVID~ Group, data = dd[!duplicated(ID)], family = binomial()) pairs(emmeans(m, ~ Group), adjust = "none") ``` # Comparisons of stress levels by group and survey ``` {r} ## comparisons of stress levels by group ## STRESS <- emmeans(lmer(STRESS ~ Survey * Group + (1 | ID), data = dd), ~ Group | Survey) STRESSd <- as.data.frame(STRESS) ## means & CIs sprintf("%s %s: %s [%s, %s]", STRESSd$Survey, STRESSd$Group, format(round(STRESSd$emmean, 1), digits = 1, nsmall = 1), format(round(STRESSd$asymp.LCL, 1), digits = 1, nsmall = 1), format(round(STRESSd$asymp.UCL, 1), digits = 1, nsmall = 1)) ## pairwise comparisons confint(pairs(STRESS, adj = "none")) ``` # Comparisons of stressor types by group ``` {r} ## comparisons between stressor types by group ## m <- glmer(STRESSDUE_Work ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Heal ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Fina ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Argu ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Home ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Disc ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Rela~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") m <- glmer(STRESSDUE_Else ~ Group + (1 | ID), data = dd, family = binomial(), nAGQ = 20) em <- emmeans(m, ~ Group, tran = "logit", type = "response") emd <- as.data.frame(em) ## probabilities & CIs sprintf("%s: %s [%s, %s]", emd$Group, format(round(emd$prob * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.LCL * 100, 1), digits = 1, nsmall = 1), format(round(emd$asymp.UCL * 100, 1), digits = 1, nsmall = 1)) ## pairwise comparisons pairs(em, adj = "none") ``` ## differences in PA between groups by surveys ## ```{r} posaff <- emmeans(lmer(PosAff ~ Survey * Group + (1 | ID), data = dd), ~ Group | Survey) posaffd <- as.data.frame(posaff) ## means & CIs sprintf("%s %s: %s [%s, %s]", posaffd$Survey, posaffd$Group, format(round(posaffd$emmean, 1), digits = 1, nsmall = 1), format(round(posaffd$asymp.LCL, 1), digits = 1, nsmall = 1), format(round(posaffd$asymp.UCL, 1), digits = 1, nsmall = 1)) ## pairwise comparisons confint(pairs(posaff, adj = "none")) ``` ## differences in NA between groups by surveys ## ```{r} negaff <- emmeans(lmer(NegAff ~ Survey * Group + (1 | ID), data = dd), ~ Group | Survey) negaffd <- as.data.frame(negaff) ## means & CIs sprintf("%s %s: %s [%s, %s]", negaffd$Survey, negaffd$Group, format(round(negaffd$emmean, 1), digits = 1, nsmall = 1), format(round(negaffd$asymp.LCL, 1), digits = 1, nsmall = 1), format(round(negaffd$asymp.UCL, 1), digits = 1, nsmall = 1)) ## pairwise comparisons confint(pairs(negaff, adj = "none")) ``` #Stress and Affect Levels by Surveys # ``` {r} #Stress and Affect Levels by Surveys # dd [Survey == "Evening", ESTRESS := STRESS, by = ID] dd [Survey == "Wake", WakeSTRESS := STRESS, by = ID] dd [Survey == "Morning", MSTRESS := STRESS, by = ID] dd [Survey == "Afternoon", ASTRESS := STRESS, by = ID] dd [Survey == "Evening", EPosAff := PosAff, by = ID] dd [Survey == "Wake", WakePosAff := PosAff, by = ID] dd [Survey == "Morning", MPosAff := PosAff, by = ID] dd [Survey == "Afternoon", APosAff := PosAff, by = ID] dd [Survey == "Evening", ENegAff := NegAff, by = ID] dd [Survey == "Wake", WakeNegAff := NegAff, by = ID] dd [Survey == "Morning", MNegAff := NegAff, by = ID] dd [Survey == "Afternoon", ANegAff := NegAff, by = ID] egltable(c( "WakeSTRESS", "MSTRESS", "ASTRESS", "ESTRESS", "WakePosAff", "MPosAff", "APosAff", "EPosAff", "WakeNegAff", "MNegAff", "ANegAff", "ENegAff"), data = dd [!duplicated(ID)], strict = FALSE) ``` # surveys completed and ICC ``` {r} ## surveys completed, ICC ## dd [!is.na(NegAff),.N, by=Survey] dd [!is.na(PosAff),.N, by=Survey] dd [!is.na(STRESS),.N, by=Survey] dd [!is.na(STRESS)|!is.na(PosAff)|!is.na(NegAff),.N, by=Survey] dd [!is.na(ID),.N, by=Survey] iccMixed("NegAff", "ID", data=dd [Survey=="Wake"]) iccMixed("NegAff", "ID", data=dd [Survey=="Morning"]) iccMixed("NegAff", "ID", data=dd [Survey=="Afternoon"]) iccMixed("NegAff", "ID", data=dd [Survey=="Evening"]) iccMixed("PosAff", "ID", data=dd [Survey=="Wake"]) iccMixed("PosAff", "ID", data=dd [Survey=="Morning"]) iccMixed("PosAff", "ID", data=dd [Survey=="Afternoon"]) iccMixed("PosAff", "ID", data=dd [Survey=="Evening"]) iccMixed("STRESS", "ID", data=dd [Survey=="Wake"]) iccMixed("STRESS", "ID", data=dd [Survey=="Morning"]) iccMixed("STRESS", "ID", data=dd [Survey=="Afternoon"]) iccMixed("STRESS", "ID", data=dd [Survey=="Evening"]) ``` # Bivariate correlations ``` {r} ## create between person stressor number for correlations dd[, BSTRESSDUE_Count := mean(STRESSDUE_Count, na.rm = TRUE), by = .(ID, SurveyDay)] dd[, BSTRESSDUE_Count := mean(BSTRESSDUE_Count, na.rm = TRUE), by = .(ID)] # Correlation matrix # cors <- SEMSummary(~ BSTRESS + BPosAff + BNegAff + Age + Sex + COVID + SES_1 + SASH_Total + TimeInMelb + BornAUS + BSTRESSDUE_Count, data = dd[!duplicated(ID)], use = "pairwise.complete.obs") plot(cors) APAStyler(cors, type = "cor", stars=TRUE) ``` # Examine distributions of variables These plots help assess the distribution of a variable. By default it tests whether a distribution is normally distributed. The top plot is a density plot. The blue line is the density curve for a normal distribution. Individual values are shown as a "rug" plot, lines, under the density plot. Below that main plot is a QQ plot, but to save space instead of perfect fit being the diagonal, it is rotated to be horizontal. ``` {r} # model diagnostics# plot (modelDiagnostics (negaffbygrp), ncol = 2, nrow = 2, ask = FALSE) plot (modelDiagnostics (posaffbygrp), ncol = 2, nrow = 2, ask = FALSE) plot (modelDiagnostics (STRESSbygrp), ncol = 2, nrow = 2, ask = FALSE) plot (modelDiagnostics (negint), ncol = 2, nrow = 2, ask = FALSE) plot (modelDiagnostics (posint), ncol = 2, nrow = 2, ask = FALSE) `geom_smooth()` using formula 'y ~ x' ``` # mutlilevel analyses using lme4 package: main effects of Group on Stress, PA, NA ```{r} #mutlilevel analyses using lme4 package: main effects of Group on Stress, PA, NA# summary(STRESSbygrp <- lmer( STRESS ~ 1 + Group + STRESSDUE_Count + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + TimeInMelb + COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data=dd)) APAStyler(modelTest (STRESSbygrp)) pairs(emmeans(STRESSbygrp, ~ Group), adjust = "none") summary(negaffbygrp <- lmer( NegAff ~ 1 + Group + STRESSDUE_Count + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + TimeInMelb + COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data=dd)) APAStyler(modelTest (negaffbygrp)) pairs(emmeans(negaffbygrp, ~ Group), adjust = "none") summary(posaffbygrp <- lmer( PosAff ~ 1 + Group + STRESSDUE_Count + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + TimeInMelb+ COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data=dd)) APAStyler (modelTest (posaffbygrp)) pairs(emmeans(posaffbygrp, ~ Group), adjust = "none") ## plot graph visreg(negaffbygrp, xvar = "Group", partial = FALSE, rug = FALSE) visreg(posaffbygrp, xvar = "Group", partial = FALSE, rug = FALSE) ## mean differences STRESSbygrpmeans<- emmeans(STRESSbygrp, "Group") pairs(STRESSbygrpmeans) negaffbygrpmeans <- emmeans(negaffbygrp, "Group") pairs(negaffbygrpmeans) posaffbygrpmeans <- emmeans(posaffbygrp, "Group") pairs(posaffbygrpmeans) ``` ## Negative Affect Reactivity Stress x Group interaction predicting negative affect, adjusted for covariates ``` {r} ## interaction of stress and group on affect ## summary(negint <- lmer( NegAff ~ 1 + STRESS * Group + STRESSDUE_Count + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + TimeInMelb + COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data = dd)) negNOint <- lmer( NegAff ~ 1 + STRESS + Group + STRESSDUE_Count + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + TimeInMelb + COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data = dd) ## get marginal f2 and overall p-value for the interaction modelCompare(negint, negNOint) ## calculate simple effects (group means) at each level of stress: 0 - 10 simple.negint <- emmeans(negint, ~ Group | STRESS, at = list(STRESS = 0:10)) ## now calculate pairwise differences at each level of stress pairs(simple.negint, adjust = "none") ## from these simple effects we can tell when each group differs significantly ## from other groups and use this to construct shading to indicate it shading.negint <- data.table( Contrast = c("C-R", "C-V", "R-V"), ## these are labels for the comparions xmin = c(NA, 1, 1), ## these are the lowest stress level where the contrast is significant xmax = c(NA, 10, 10), ## highest stress level where the contrast is sign ymin = c(1.0, 1.501, 2.001), ## the start of the y axis region to shade ymax = c(1.5, 2.0, 2.5)) ## end of the y axis region to shade ## calculate simple slopes of negative affect on stress by group slopes.negint <- emtrends(negint, ~ Group, var = "STRESS") ## show slopes with confidence intervals and significance tests summary(slopes.negint, infer = TRUE) ## whether stress slope differs between groups pairs(slopes.negint, adjust = "none") ## make the final plot plot.negint <- visreg(negint, xvar = "STRESS", by = "Group", overlay=TRUE, partial = FALSE, rug = FALSE, gg = TRUE) + geom_rect(aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = Contrast), data = shading.negint, alpha = .2, inherit.aes = FALSE) + annotate("text", x = 0.1, y = 2.25, label = "R-V", hjust = 0) + annotate("text", x = 0.1, y = 1.75, label = "C-V", hjust = 0) + annotate("text", x = 0.1, y = 1.25, label = "C-R", hjust = 0) + scale_x_continuous("Stress", breaks = c(0, 2, 4, 6, 8, 10)) + ylab ("Negative Affect") + theme_pubr() + scale_fill_manual(values = c("Control" = "white", "Resilient" = "white", "Vulnerable" = "white", "C-R" = "grey80", "C-V" = "grey50", "R-V" = "grey20"), guide = FALSE) + scale_colour_manual("", values = c("Control" = "#b2abd2", "Resilient" = "#5e3c99", "Vulnerable" = "#e66101")) + theme(axis.title.x = element_text(family = NULL, face = NULL, colour = NULL, size = NULL, hjust = NULL, vjust = NULL, angle = NULL, lineheight = NULL, color = NULL, margin = NULL, debug = NULL, inherit.blank = FALSE), legend.key.width = unit(1, "cm")) + coord_cartesian(xlim = c(0,10.5), ylim = c(1,3), expand = FALSE) ``` ## Positive Affect Reactivity Stress x Group interaction predicting negative affect, adjusted for covariates ``` {r} ## interaction of stress and group on affect ## summary(posint <- lmer( PosAff ~ 1 + STRESS * Group + STRESSDUE_Count + TimeInMelb + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data = dd)) posNOint <- lmer( PosAff ~ 1 + STRESS + Group + STRESSDUE_Count + TimeInMelb + Age + BornAUS + Sex + SES_1 + RACE3G + SASH_Total + COVID + DayofWeek + StudyDay + SurveyMonth + (1 | ID), data = dd) ## get marginal f2 and overall p-value for the interaction modelCompare(posint, posNOint) ## calculate simple effects (group means) at each level of stress: 0 - 10 simple.posint <- emmeans(posint, ~ Group | STRESS, at = list(STRESS = 0:10)) ## now calculate pairwise differences at each level of stress pairs(simple.posint, adjust = "none") ## from these simple effects we can tell when each group differs significantly ## from other groups and use this to construct shading to indicate it shading.posint <- data.table( Contrast = c("C-R", "C-V", "R-V"), ## these are labels for the comparions xmin = c(NA, NA, NA), ## these are the lowest stress level where the contrast is significant xmax = c(NA, NA, NA), ## highest stress level where the contrast is sign ymin = c(1.5, 2.001, 2.501), ## the start of the y axis region to shade ymax = c(2.0, 2.5, 3.0)) ## end of the y axis region to shade ## calculate simple slopes of negative affect on stress by group slopes.posint <- emtrends(posint, ~ Group, var = "STRESS") ## show slopes with confidence intervals and significance tests summary(slopes.posint, infer = TRUE) ## whether stress slope differs between groups pairs(slopes.posint, adjust = "none") ## make the final plot plot.posint <- visreg(posint, xvar = "STRESS", by = "Group", overlay=TRUE, partial = FALSE, rug = FALSE, gg = TRUE) + ## geom_rect(aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = Contrast), ## data = shading.posint, alpha = .2, inherit.aes = FALSE) + annotate("text", x = 0.1, y = 2.75, label = "R-V", hjust = 0) + annotate("text", x = 0.1, y = 2.25, label = "C-V", hjust = 0) + annotate("text", x = 0.1, y = 1.75, label = "C-R", hjust = 0) + scale_x_continuous("Stress", breaks = c(0, 2, 4, 6, 8, 10)) + ylab ("Positive Affect") + theme_pubr() + scale_fill_manual(values = c("Control" = "white", "Resilient" = "white", "Vulnerable" = "white", "C-R" = "grey80", "C-V" = "grey50", "R-V" = "grey20"), guide = FALSE) + scale_colour_manual("", values = c("Control" = "#b2abd2", "Resilient" = "#5e3c99", "Vulnerable" = "#e66101")) + theme(axis.title.x = element_text(family = NULL, face = NULL, colour = NULL, size = NULL, hjust = NULL, vjust = NULL, angle = NULL, lineheight = NULL, color = NULL, margin = NULL, debug = NULL, inherit.blank = FALSE), legend.key.width = unit(1, "cm")) + coord_cartesian(xlim = c(0,10.5), ylim = c(1.5,3.5), expand = FALSE) ``` ## Combined Plot of Negative and Positive Affect Interactions ```{r} ## combine into one and store plot.negposint <- ggarrange(plot.negint, plot.posint, common.legend = TRUE, labels = c("A", "B")) ## view the plot print(plot.negposint) ## save a PDF of the plot, good high quality record to keep ggsave("stress_group_interaction.pdf", plot.negposint, width = 20, height = 10, units = "cm") ```