#Code for analyses included in: #Dodd et al. #Trajectories of anxiety when children start school: the role of Behavioural Inhibition and attention bias to angry and happy faces #Published in Journal of Abnormal Psychology library(lme4) library(papaja) library(psych) library(emmeans) library(corrplot) library(mice) library(mitml) library(Hmisc) library(MASS) library(dplyr) library(readr) library(tidyverse) # Read in data from multiple source files #### #setwd("~/Dropbox/Work master folder/Research projects/ESRC FRL") df_cov <- read_csv(file.path("R", "data", "Parent baseline questionnaire.csv")) df_cov %>% select(child_ID, aBIQtot) %>% str(give.attr = FALSE) df_labtab <- read_csv(file.path("R", "data", "LabTAB observation of BI.csv")) df_cov2 <- merge(df_cov, df_labtab, id=child_ID, all = TRUE) df_cov2$ZBIScore_BIQ <- scale(df_cov2$aBIQtot, center=TRUE, scale = TRUE) df_cov2$ZBIScore_LABTAB <- scale (df_cov2$BIScore_LABTAB, center=TRUE, scale = TRUE) df_cov2$CombinedBI <- rowMeans(df_cov2[c('ZBIScore_BIQ', 'ZBIScore_LABTAB')], na.rm=FALSE) df_cov2 %>% select(ZBIScore_BIQ, ZBIScore_LABTAB, CombinedBI) %>% summary() Agedata <- read_csv(file.path("R", "data", "Ages at each stage.csv")) Agedata2<-Agedata [, c("child_ID", "As_age", "Sch_age", "Time_as_sch")] df_cov3 <- merge(df_cov2, Agedata2, id=child_ID, all = TRUE) anx_av <- read_csv(file.path("R", "data", "Text message anxiety ratings data.csv")) anx_av2 <- anx_av [, c("child_ID", "av_anxiety")] df_cov4 <- merge(df_cov3, anx_av2, id=child_ID, all = TRUE) followup_df <- read_csv(file.path("R", "data", "Parent follow-up 1 questionnaires.csv")) followup_df2 <- followup_df [, c("child_ID", "bPAStot")] df_cov5 <- merge(df_cov4, followup_df2, id=child_ID, all = TRUE) ET_df <- read_csv(file.path("R", "data", "Eyetracking bias data.csv")) ET_df2 <- ET_df [, c("child_ID", "orientationbias_angry", "orientationbias_happy", "dwellbias_angry", "dwellbias_happy")] df_cov6 <- merge(df_cov5, ET_df2, id=child_ID, all = TRUE) df_cov <- df_cov6 rm (df_cov2, df_cov3, df_cov4, df_cov5, df_cov6, df_labtab, Agedata, Agedata2, anx_av, anx_av2, ET_df, ET_df2, followup_df, followup_df2) df_cov <- df_cov %>% rename(DwellBias_A = dwellbias_angry) df_cov <- df_cov %>% rename(DwellBias_H = dwellbias_happy ) df_cov <- df_cov %>% rename(OrientationBias_Angry = orientationbias_angry) df_cov <- df_cov %>% rename(OrientationBias_Happy = orientationbias_happy) df_cov$combined_BI_cent <- scale(df_cov$CombinedBI, scale = FALSE)[,1] df_cov$ori_bias_ang_cent <- scale(df_cov$OrientationBias_Angry, scale = FALSE)[,1] df_cov$ori_bias_hap_cent <- scale(df_cov$OrientationBias_Happy, scale = FALSE)[,1] df_cov$dwell_bias_ang_cent <- scale(df_cov$DwellBias_A, scale = FALSE)[,1] df_cov$dwell_bias_hap_cent <- scale(df_cov$DwellBias_H, scale = FALSE)[,1] ## Create a correlation matrix between variables #### #Create dataframe with only variables I want in correlation matrix. correlation_df <- df_cov %>% select(DwellBias_A, DwellBias_H, aPAStot, bPAStot, CombinedBI, av_anxiety) #rename variables as I want them for the correlation matrix plot names(correlation_df)[names(correlation_df) == "DwellBias_A"] <- "Dwell Bias Angry Faces" names(correlation_df)[names(correlation_df) == "DwellBias_H"] <- "Dwell Bias Happy Faces" names(correlation_df)[names(correlation_df) == "aPAStot"] <- "Baseline PAS" names(correlation_df)[names(correlation_df) == "bPAStot"] <- "Follow-up PAS" names(correlation_df)[names(correlation_df) == "CombinedBI"] <- "BI" names(correlation_df)[names(correlation_df) == "av_anxiety"] <- "Average Daily Anxiety Rating" library(apaTables) apa.cor.table(correlation_df, filename="Table1.doc", table.number=1) #Explore missingness df_cov %>% select(DwellBias_A, DwellBias_H, CombinedBI, bPAStot, aPAStot) %>% md.pattern() # run MI for regressions #### df_cov_analysis<-df_cov [, c("child_ID", "ori_bias_ang_cent", "ori_bias_hap_cent", "dwell_bias_ang_cent", "dwell_bias_hap_cent", "combined_BI_cent", "aPAStot", "bPAStot")] A<-is.na(df_cov_analysis) describe(A) imp <- mice(data = df_cov_analysis, maxit=0) pred <- imp$pred pred[c("ori_bias_ang_cent", "ori_bias_hap_cent", "dwell_bias_ang_cent", "dwell_bias_hap_cent", "combined_BI_cent", "aPAStot", "bPAStot"), "child_ID"] <- 0 pred library(miceadds) # run MI df_cov_imp <- mice(data = df_cov_analysis, pred = pred, m = 20) summary(df_cov_imp) df_cov_imp_stacked <- mice::complete(df_cov_imp, action='long', include=TRUE) df_cov_imp_stacked$.imp_factor <- factor(df_cov_imp_stacked$.imp) #Rename variables that are to be centred df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(ori_bias_ang = ori_bias_ang_cent) df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(ori_bias_hap = ori_bias_hap_cent) df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(dwell_bias_ang = dwell_bias_ang_cent) df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(dwell_bias_hap = dwell_bias_hap_cent) df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(bPAStotRaw = bPAStot) df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(combined_BI = combined_BI_cent) df_cov_imp_stacked <- df_cov_imp_stacked %>% rename(aPAStotRaw = aPAStot) #Centre variables in new stacked long file so that new variables have names used in analysis code below df_cov_imp_stacked_1 <- transform(df_cov_imp_stacked, ori_bias_ang_cent=ave(ori_bias_ang, .imp_factor, FUN=scale)) df_cov_imp_stacked_2 <- transform(df_cov_imp_stacked_1, ori_bias_hap_cent=ave(ori_bias_hap, .imp_factor, FUN=scale)) df_cov_imp_stacked_3 <- transform(df_cov_imp_stacked_2, dwell_bias_ang_cent=ave(dwell_bias_ang, .imp_factor, FUN=scale)) df_cov_imp_stacked_4 <- transform(df_cov_imp_stacked_3, dwell_bias_hap_cent=ave(dwell_bias_hap, .imp_factor, FUN=scale)) df_cov_imp_stacked_5 <- transform(df_cov_imp_stacked_4, bPAStot=ave(bPAStotRaw, .imp_factor, FUN=scale)) df_cov_imp_stacked_6 <- transform(df_cov_imp_stacked_5, combined_BI_cent=ave(combined_BI, .imp_factor, FUN=scale)) df_cov_imp_stacked_7 <- transform(df_cov_imp_stacked_6, aPAStot=ave(aPAStotRaw, .imp_factor, FUN=scale)) df_cov_imp <- as.mids(df_cov_imp_stacked_7) # Cross-sectional regressions #### #Note: Orientation bias analyses not included in paper due to poor reliability of bias score. #Editor asked for these analyses to be taken out of the final version. They are commented out. #reg_mod_1 <- # with(data=df_cov_imp, exp=lm(aPAStot ~ combined_BI_cent * ori_bias_ang_cent)) #pool(reg_mod_1) #summary(pool(reg_mod_1), conf.int=TRUE) #pool.r.squared(reg_mod_1, adjusted = FALSE) reg_mod_2 <- with(data=df_cov_imp, exp=lm(aPAStot ~ combined_BI_cent * dwell_bias_ang_cent)) pool(reg_mod_2) summary(pool(reg_mod_2), conf.int=TRUE) pool.r.squared(reg_mod_2, adjusted = FALSE) #reg_mod_3 <- # with(data=df_cov_imp, exp=lm(aPAStot ~ combined_BI_cent * ori_bias_hap_cent)) #pool(reg_mod_3) #summary(pool(reg_mod_3), conf.int=TRUE) #pool.r.squared(reg_mod_3, adjusted = FALSE) reg_mod_4 <- with(data=df_cov_imp, exp=lm(aPAStot ~ combined_BI_cent * dwell_bias_hap_cent)) pool(reg_mod_4) summary(pool(reg_mod_4), conf.int=TRUE) pool.r.squared(reg_mod_4, adjusted = FALSE) #Note: Analysis added following reviewer comments. reg_mod_5 <- with(data=df_cov_imp, exp=lm(aPAStot ~ combined_BI_cent * dwell_bias_hap_cent * dwell_bias_ang_cent)) pool(reg_mod_5) summary(pool(reg_mod_5), conf.int=TRUE) pool.r.squared(reg_mod_5, adjusted = FALSE) # Regressions for predicting follow-up PAS #### #reg_mod_1 <- with(data=df_cov_imp, exp=lm(bPAStot ~ combined_BI_cent * ori_bias_ang_cent)) #pool(reg_mod_1) #summary(pool(reg_mod_1), conf.int=TRUE) #pool.r.squared(reg_mod_1, adjusted = FALSE) reg_mod_2 <- with(data=df_cov_imp, exp=lm(bPAStot ~ combined_BI_cent * dwell_bias_ang_cent)) pool(reg_mod_2) summary(pool(reg_mod_2), conf.int=TRUE) pool.r.squared(reg_mod_2, adjusted = FALSE) #reg_mod_3 <- # with(data=df_cov_imp, exp=lm(bPAStot ~ combined_BI_cent * ori_bias_hap_cent)) #pool(reg_mod_3) #summary(pool(reg_mod_3), conf.int=TRUE) #pool.r.squared(reg_mod_3, adjusted = FALSE) reg_mod_4 <- with(data=df_cov_imp, exp=lm(bPAStot ~ combined_BI_cent * dwell_bias_hap_cent)) pool(reg_mod_4) summary(pool(reg_mod_4), conf.int=TRUE) pool.r.squared(reg_mod_4, adjusted = FALSE) #Analysis added after reviewer comments reg_mod_5 <- with(data=df_cov_imp, exp=lm(bPAStot ~ combined_BI_cent * dwell_bias_hap_cent * dwell_bias_ang_cent)) pool(reg_mod_5) summary(pool(reg_mod_5), conf.int=TRUE) pool.r.squared(reg_mod_5, adjusted = FALSE) # Preparing for GCA #### #Read in text message data df_anx <- read_csv(file.path("R", "data", "Text message anxiety ratings data.csv")) df_anx$av_anxiety<-NULL df_anx$X45<-NULL df_anx$X46<-NULL df_anx$X47<-NULL df_anx$X48<-NULL df_anx$X49<-NULL df_anx$X50<-NULL df_anx %>% select(`child_ID`, ends_with("Anxiety")) %>% str(give.attr = FALSE) df_anx <- df_anx %>% mutate_if(is.integer, as.numeric) df_anx_long <- df_anx %>% gather(key = occasion, value = rating, -`child_ID`, -`Usual anxiety/pre survey`) %>% extract(occasion, c("day", "type"), "[^\\d]*(\\d+)\\s(.*)") %>% spread(key = type, value = rating) head(df_anx_long, n = 20) df_anx_long$day <- as.integer(df_anx_long$day) summary(df_anx_long$Anxiety) length(which(df_anx_long$Anxiety<1)) df_cov$child_ID %>% unique %>% sort df_anx_long$`child_ID` %>% unique %>% sort df_cov %>% select(child_ID, OrientationBias_Angry, OrientationBias_Happy, DwellBias_A, DwellBias_H, CombinedBI, aPAStot, bPAStot) setdiff(unique(df_cov$child_ID), unique(df_anx_long$`child_ID`)) setdiff(unique(df_anx_long$`child_ID`), unique(df_cov$child_ID)) df_analysis <- left_join(select(df_anx_long, child_ID, day, Anxiety), select(df_cov, child_ID, OrientationBias_Angry, OrientationBias_Happy, DwellBias_A, DwellBias_H, CombinedBI, aPAStot, bPAStot), by = "child_ID") df_analysis <- df_analysis %>% arrange(child_ID, day) head(df_analysis, n = 20) summary(df_analysis$Anxiety) # set maximum number of polynomials (degrees) degrees <- 5 # create orthogonal polynomials days_orthogonal_polynomials <- poly(df_analysis$day, degrees) # shape output so that we can add it into our dataframe for analysis days_orthogonal_polynomials <- as_tibble(days_orthogonal_polynomials[,]) names(days_orthogonal_polynomials) <- paste0("day_", 1:degrees) # add it to our dataframe for analysis df_analysis <- bind_cols(df_analysis, days_orthogonal_polynomials) #Centre variables df_analysis$combined_BI_cent <- scale(df_analysis$CombinedBI, scale = FALSE)[,1] df_analysis$ori_bias_ang_cent <- scale(df_analysis$OrientationBias_Angry, scale = FALSE)[,1] df_analysis$ori_bias_hap_cent <- scale(df_analysis$OrientationBias_Happy, scale = FALSE)[,1] df_analysis$dwell_bias_ang_cent <- scale(df_analysis$DwellBias_A, scale = FALSE)[,1] df_analysis$dwell_bias_hap_cent <- scale(df_analysis$DwellBias_H, scale = FALSE)[,1] # Multiple imputation for GCA data #### df_analysis['const']='1' df_analysis$const <- as.numeric(df_analysis$const) df_analysis_MI<-df_analysis [, c("child_ID", "day", "day_1", "day_2", "day_3", "Anxiety", "ori_bias_ang_cent", "ori_bias_hap_cent", "dwell_bias_ang_cent", "dwell_bias_hap_cent", "combined_BI_cent", "aPAStot", "bPAStot", "const")] df_analysis_MI$child_ID <- as.integer(df_analysis_MI$child_ID) #Create a matrix for 'where' function in MI below. Set to FALSE variables I don't want imputed B<-is.na(df_analysis_MI) describe(B) B[, c(1, 2, 3, 4, 5, 12, 14) ]<-FALSE describe(B) #Set predictors in pred matrix. Some changed to 0 because errors showed models overspecified. #Can't get it to work for Anxiety as well imp <- mice(data = df_analysis_MI, maxit=0) imp pred <- imp$pred pred pred["Anxiety",] <- c(-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 2) pred["combined_BI_cent",] <- c(-2, 0, 0, 0, 0, 3, 3, 3, 3, 3, 0, 3, 0, 2) pred["ori_bias_ang_cent",] <- c(-2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3, 3, 3, 2) pred["ori_bias_hap_cent",] <- c(-2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3, 3, 3, 2) pred["dwell_bias_ang_cent",] <- c(-2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3, 3, 3, 2) pred["dwell_bias_hap_cent",] <- c(-2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3, 3, 3, 2) pred["bPAStot",] <- c(-2, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 0, 2) pred library(miceadds) #Run MI using mice df_analysis_MI_imp <- mice(data = df_analysis_MI, pred = pred, method = "2l.norm", where = B, m = 20, seed = 123) summary(df_analysis_MI_imp) attributes(df_analysis_MI_imp) plot(df_analysis_MI_imp) stripplot(df_analysis_MI_imp, combined_BI_cent~.imp, pch=20, cex=2) stripplot(df_analysis_MI_imp, Anxiety~.imp, pch=20, cex=2) stripplot(df_analysis_MI_imp, ori_bias_ang_cent~.imp, pch=20, cex=2) stripplot(df_analysis_MI_imp, ori_bias_hap_cent~.imp, pch=20, cex=2) stripplot(df_analysis_MI_imp, dwell_bias_ang_cent~.imp, pch=20, cex=2) stripplot(df_analysis_MI_imp, dwell_bias_hap_cent~.imp, pch=20, cex=2) stripplot(df_analysis_MI_imp, bPAStot~.imp, pch=20, cex=2) df_analysis_imp_stacked <- mice::complete(df_analysis_MI_imp, action='long', include=TRUE) df_analysis_imp_stacked$.imp_factor <- factor(df_analysis_imp_stacked$.imp) #Rename variables that are to be centred df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(ori_bias_ang = ori_bias_ang_cent) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(ori_bias_hap = ori_bias_hap_cent) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(dwell_bias_ang = dwell_bias_ang_cent) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(dwell_bias_hap = dwell_bias_hap_cent) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(bPAStotRaw = bPAStot) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(combined_BI = combined_BI_cent) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(aPAStotRaw = aPAStot) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(day_1_1 = day_1) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(day_2_1 = day_2) df_analysis_imp_stacked <- df_analysis_imp_stacked %>% rename(day_3_1 = day_3) #Centre variables in new stacked long file, new variables have names used in analysis df_analysis_imp_stacked_1 <- transform(df_analysis_imp_stacked, ori_bias_ang_cent=ave(ori_bias_ang, .imp_factor, FUN=scale)) df_analysis_imp_stacked_2 <- transform(df_analysis_imp_stacked_1, ori_bias_hap_cent=ave(ori_bias_hap, .imp_factor, FUN=scale)) df_analysis_imp_stacked_3 <- transform(df_analysis_imp_stacked_2, dwell_bias_ang_cent=ave(dwell_bias_ang, .imp_factor, FUN=scale)) df_analysis_imp_stacked_4 <- transform(df_analysis_imp_stacked_3, dwell_bias_hap_cent=ave(dwell_bias_hap, .imp_factor, FUN=scale)) df_analysis_imp_stacked_5 <- transform(df_analysis_imp_stacked_4, bPAStot=ave(bPAStotRaw, .imp_factor, FUN=scale)) df_analysis_imp_stacked_6 <- transform(df_analysis_imp_stacked_5, combined_BI_cent=ave(combined_BI, .imp_factor, FUN=scale)) df_analysis_imp_stacked_7 <- transform(df_analysis_imp_stacked_6, aPAStot=ave(aPAStotRaw, .imp_factor, FUN=scale)) df_analysis_imp_stacked_8 <- transform(df_analysis_imp_stacked_7, day_1=ave(day_1_1, .imp_factor, FUN=scale)) df_analysis_imp_stacked_9 <- transform(df_analysis_imp_stacked_8, day_2=ave(day_2_1, .imp_factor, FUN=scale)) df_analysis_imp_stacked_10 <- transform(df_analysis_imp_stacked_9, day_3=ave(day_3_1, .imp_factor, FUN=scale)) #Centre variables in df_analysis df df_analysis_1 <- transform(df_analysis, ori_bias_ang_cent=ave(ori_bias_ang_cent, FUN=scale)) df_analysis_2 <- transform(df_analysis_1, ori_bias_hap_cent=ave(ori_bias_hap_cent, FUN=scale)) df_analysis_3 <- transform(df_analysis_2, dwell_bias_ang_cent=ave(dwell_bias_ang_cent, FUN=scale)) df_analysis_4 <- transform(df_analysis_3, dwell_bias_hap_cent=ave(dwell_bias_hap_cent, FUN=scale)) df_analysis_5 <- transform(df_analysis_4, bPAStot=ave(bPAStot, FUN=scale)) df_analysis_6 <- transform(df_analysis_5, combined_BI_cent=ave(combined_BI_cent, FUN=scale)) df_analysis_7 <- transform(df_analysis_6, aPAStot=ave(aPAStot, FUN=scale)) df_analysis_8 <- transform(df_analysis_7, day_1=ave(day_1, FUN=scale)) df_analysis_9 <- transform(df_analysis_8, day_2=ave(day_2, FUN=scale)) df_analysis_10 <- transform(df_analysis_9, day_3=ave(day_3, FUN=scale)) df_analysis_MI_imp <- as.mids(df_analysis_imp_stacked_10) # create a list of completed data sets implist2 <- mids2mitml.list(df_analysis_MI_imp) # GCA Models for BI and Attention bias predicting anxiety #### contr <- lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=100000000)) fit0 <- with(implist2, lmer(Anxiety ~ (1 | child_ID))) summary(est <- pool(fit0)) # Angry dwell bias model #### #fit3_dwell_a <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_ang_cent * (day_1 + day_2 + day_3) + (1 + day_1 + day_2 + day_3 | child_ID), # ), control = contr) #fit3_dwell_a <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_ang_cent * (day_1 + day_2 + day_3) + (1 | child_ID), # ), control = contr) #summary(est <- pool(fit3_dwell_a)) #fit2_dwell_a <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_ang_cent * (day_1 + day_2 ) + (1 + day_1 + day_2 | child_ID), # ), control = contr) #Final model fit2_dwell_a <- with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_ang_cent * (day_1 + day_2 ) + (1 | child_ID), ), control = contr) round(summary(pool(fit2_dwell_a), conf.int=TRUE), 3) testEstimates(as.mitml.result(fit2_dwell_a), var.comp = T)$var.comp multilevelR2(fit2_dwell_a) #Non-imputed version of final model (for plot) fit2_dwell_a <- lmer(Anxiety ~ combined_BI_cent * dwell_bias_ang_cent * (day_1 + day_2 ) + (1 | child_ID), data = df_analysis_10, control = contr) summary(fit2_dwell_a) #Plot angry dwell bias interaction n_days <- max(df_analysis_10$day) day_1_vals <- df_analysis_10$day_1[1:n_days] day_2_vals <- df_analysis_10$day_2[1:n_days] low_dwell_bias <- - sd(df_analysis_10$dwell_bias_ang_cent, na.rm = TRUE) high_dwell_bias <- sd(df_analysis_10$dwell_bias_ang_cent, na.rm = TRUE) low_bi <- - sd(df_analysis_10$combined_BI_cent, na.rm = TRUE) high_bi <- sd(df_analysis_10$combined_BI_cent, na.rm = TRUE) df_yhat_ldb_lbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_ang_cent = low_dwell_bias, combined_BI_cent = low_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "Low dwell bias", bi_group = "Low behavioural inhibtion") df_yhat_ldb_hbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_ang_cent = low_dwell_bias, combined_BI_cent = high_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "Low dwell bias", bi_group = "High behavioural inhibtion") df_yhat_hdb_lbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_ang_cent = high_dwell_bias, combined_BI_cent = low_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "High dwell bias", bi_group = "Low behavioural inhibtion") df_yhat_hdb_hbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_ang_cent = high_dwell_bias, combined_BI_cent = high_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "High dwell bias", bi_group = "High behavioural inhibtion") df_yhat <- bind_rows(df_yhat_ldb_lbi, df_yhat_ldb_hbi, df_yhat_hdb_lbi, df_yhat_hdb_hbi) # predictions df_yhat$Anxiety <- predict(fit2_dwell_a, newdata = df_yhat, allow.new.levels = TRUE) # add CI's by getting the standard errors from the effects function. mes <- effects::effect('combined_BI_cent * dwell_bias_ang_cent * day_2', fit2_dwell_a, xlevels = list(combined_BI_cent = c(low_bi, high_bi), dwell_bias_ang_cent = c(low_dwell_bias, high_dwell_bias), day_2 = day_2_vals), se=TRUE, confidence.level=.95, typical=mean) mes <- as.data.frame(mes) mes <- mes %>% arrange(dwell_bias_ang_cent, combined_BI_cent) df_yhat$lower <- df_yhat$Anxiety - (1.96 * mes$se) df_yhat$upper <- df_yhat$Anxiety + (1.96 * mes$se) # plot ggplot(df_yhat, aes(x = day, y = Anxiety, group = ob_group, colour = ob_group, fill = ob_group)) + geom_ribbon(aes(ymin = df_yhat$lower, ymax = df_yhat$upper), alpha = 0.25, linetype = 0, fill = "grey70") + geom_line() + facet_wrap("~ bi_group") + scale_colour_grey(start = 0, end = 0.7, name = "Dwell Bias Angry Faces") + scale_x_continuous(name = "Day") + scale_fill_discrete(guide = FALSE) + theme_light() # Happy dwell bias model #### #fit3_dwell_h <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 + day_2 + day_3) + (1 + day_1 + day_2 + day_3 | child_ID), # ), control = contr) #fit3_dwell_h <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 + day_2 + day_3) + (1 | child_ID), # ), control = contr) #summary(est <- pool(fit3_dwell_h)) #fit2_dwell_h <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 + day_2 ) + (1 + day_1 + day_2 | child_ID), # ), control = contr) #fit2_dwell_h <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 + day_2 ) + (1 | child_ID), # ), control = contr) #summary(est <- pool(fit2_dwell_h)) #fit1_dwell_h <- # with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 ) + (1 + day_1 | child_ID), # ), control = contr) #summary(est <- pool(fit1_dwell_h)) #Final model fit1_dwell_h_noslopes <- with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 ) + (1 | child_ID), ), control = contr) round(summary(pool(fit1_dwell_h_noslopes), conf.int=TRUE), 3) testEstimates(as.mitml.result(fit1_dwell_h_noslopes), var.comp = T) multilevelR2(fit1_dwell_h_noslopes) #Non-imputed version of final model for plot fit1_dwell_h <- lmer(Anxiety ~ combined_BI_cent * dwell_bias_hap_cent * (day_1 ) + (1 | child_ID), data = df_analysis_10, control = contr) summary(fit1_dwell_h) #Plot happy dwell bias interaction n_days <- max(df_analysis_10$day) day_1_vals <- df_analysis_10$day_1[1:n_days] day_2_vals <- df_analysis_10$day_2[1:n_days] low_dwell_bias <- - sd(df_analysis_10$dwell_bias_hap_cent, na.rm = TRUE) high_dwell_bias <- sd(df_analysis_10$dwell_bias_hap_cent, na.rm = TRUE) low_bi <- - sd(df_analysis_10$combined_BI_cent, na.rm = TRUE) high_bi <- sd(df_analysis_10$combined_BI_cent, na.rm = TRUE) df_yhat_ldb_lbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_hap_cent = low_dwell_bias, combined_BI_cent = low_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "Low dwell bias", bi_group = "Low behavioural inhibtion") df_yhat_ldb_hbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_hap_cent = low_dwell_bias, combined_BI_cent = high_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "Low dwell bias", bi_group = "High behavioural inhibtion") df_yhat_hdb_lbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_hap_cent = high_dwell_bias, combined_BI_cent = low_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "High dwell bias", bi_group = "Low behavioural inhibtion") df_yhat_hdb_hbi <- tibble(child_ID = 999, day = 1:n_days, dwell_bias_hap_cent = high_dwell_bias, combined_BI_cent = high_bi, day_1 = day_1_vals, day_2 = day_2_vals, ob_group = "High dwell bias", bi_group = "High behavioural inhibtion") df_yhat <- bind_rows(df_yhat_ldb_lbi, df_yhat_ldb_hbi, df_yhat_hdb_lbi, df_yhat_hdb_hbi) # predictions df_yhat$Anxiety <- predict(fit1_dwell_h, newdata = df_yhat, allow.new.levels = TRUE) # add CI's by getting the standard errors from the effects function. mes <- effects::effect('combined_BI_cent * dwell_bias_hap_cent * day_1', fit1_dwell_h, xlevels = list(combined_BI_cent = c(low_bi, high_bi), dwell_bias_hap_cent = c(low_dwell_bias, high_dwell_bias), day_1 = day_1_vals), se=TRUE, confidence.level=.95, typical=mean) mes <- as.data.frame(mes) mes <- mes %>% arrange(dwell_bias_hap_cent, combined_BI_cent) df_yhat$lower <- df_yhat$Anxiety - (1.96 * mes$se) df_yhat$upper <- df_yhat$Anxiety + (1.96 * mes$se) # plot ggplot(df_yhat, aes(x = day, y = Anxiety, group = ob_group, colour = ob_group, fill = ob_group)) + geom_ribbon(aes(ymin = df_yhat$lower, ymax = df_yhat$upper), alpha = 0.25, linetype = 0, fill = "grey70") + geom_line() + facet_wrap("~ bi_group") + scale_colour_grey(start = 0, end = 0.7, name = "Dwell Bias Happy Faces") + scale_x_continuous(name = "Day") + scale_fill_discrete(guide = FALSE) + theme_light() # Angry by happy bias (check asked for by reviewers) #### fit2_dwell_a_h <- with(implist2, lmer(Anxiety ~ combined_BI_cent * dwell_bias_ang_cent * dwell_bias_hap_cent * (day_1 + day_2 ) + (1 | child_ID), ), control = contr) round(summary(pool(fit2_dwell_a_h), conf.int=TRUE), 3) testEstimates(as.mitml.result(fit2_dwell_a_h), var.comp = T) multilevelR2(fit2_dwell_a_h)