#set working directory to the folder to which files created in R should be saved (e.g., "C:/PCmatch") setwd("C:/PCmatch") #install packages install.packages(c("psych","data.table","car","ggplot2","gridExtra","flextable","officer","tidyverse","dplyr","numform")) #load packages library(psych) library(data.table) library(car) library(ggplot2) library(gridExtra) library(flextable) library(officer) library(tidyverse) library(tibble) library(dplyr) library(numform) #####################################Analyses with full OOS Dataset (not publicly available)##################################### # ##install additional packages #install.packages(c("lavaan","semTools","multilevel")) # ##load additional packages #library(lavaan) #library(semTools) #library(multilevel) # ##read data #dat <- as.data.frame(fread("04_OOS_MASTER DATASET_2016_03_03.csv", header = T, sep = ',')) # ##select variables #dat <- subset(dat, select = c(record_id, sta_now, age, female, langu, relig, se, ext01:opn10)) # ##recode items to compute agy_ac and agy_bf #dat$neu02 <- car::recode(dat$neu02r, "1=5; 2=4; 3=3; 4=2; 5=1;") #dat$neu07 <- car::recode(dat$neu07r, "1=5; 2=4; 3=3; 4=2; 5=1;") # ##compute personality traits #dat$agr <- rowMeans(dat[c("agr01r","agr02","agr03r","agr04","agr05","agr06r","agr07","agr08r","agr09")], na.rm=TRUE) #dat$cns <- rowMeans(dat[c("cns01","cns02r","cns03","cns04r","cns05r","cns06","cns07","cns08","cns09r")], na.rm=TRUE) #dat$ext <- rowMeans(dat[c("ext01","ext02r","ext03","ext04","ext05r","ext06","ext07r","ext08")], na.rm=TRUE) #dat$opn <- rowMeans(dat[c("opn01","opn02","opn03","opn04","opn05","opn06","opn07r","opn08","opn09r","opn10")], na.rm=TRUE) #dat$neu <- rowMeans(dat[c("neu01","neu02r","neu03","neu04","neu05r","neu06","neu07r","neu08")], na.rm=TRUE) #dat$agy_ts <- rowMeans(dat[c("ext01","ext03","ext04","ext05r","ext06","ext07r","ext08","opn01")], na.rm=TRUE) #dat$com_ts <- rowMeans(dat[c("agr01r","agr02","agr03r","agr04","agr07","agr08r","agr09","cns07")], na.rm=TRUE) #dat$agy_ac <- rowMeans(dat[c("ext03","ext06","ext07r","opn04","opn07r","neu07")], na.rm=TRUE) #dat$com_ac <- rowMeans(dat[c("agr01r","agr03r","agr05","agr07","agr08r","agr09","cns01")], na.rm=TRUE) #dat$agy_bf <- rowMeans(dat[c("cns06","ext02r","ext05r","opn01","opn03","opn05","opn10","neu02")], na.rm=TRUE) #dat$com_bf <- rowMeans(dat[c("agr02","agr03r","agr04","agr06r","agr07","agr08r","agr09","cns01")], na.rm=TRUE) # #dat$sta_now <- car::recode(dat$sta_now, "1=1; 2=2; 4=4; 5=5; 7=7; 8=8; 9=9; 10=10; 11=11; # 14=14; 17=17; 19=19; 20=20; 21=21; 22=22; 23=23; 24=24; 25=25; # 27=27; 28=28; 30=30; 31=31; 32=32; 33=33; 34=34; 35=35; 36=36; # 37=37; 38=38; 41=41; 42=42; 43=43; 44=44; 47=47; 48=48; 52=52; # 53=53; 55=55; 57=57; 61=61; 63=63; 64=64; 65=65; 66=66; 67=67; # 69=69; 70=70; 71=71; 72=72; 73=73; 74=74; else=NA") # ##N = 9,328,610 # ##select data #dat[is.na(dat)] <- NA #dat <- subset(dat, is.na(sta_now)==F & is.na(relig)==F & is.na(se)==F & is.na(agr)==F & is.na(cns)==F & is.na(ext)==F & is.na(opn)==F & is.na(neu)==F) #dat <- transform(dat, sta_frq = ave(seq(nrow(dat)), sta_now, FUN=length)) #dat <- subset(dat, sta_frq >= 300) # ##N = 1,142,224 # ##NOTE: there are no missings at the agency-communion level (checked as follows) #sum(is.na(dat$agy_ts)) #sum(is.na(dat$com_ts)) #sum(is.na(dat$agy_ac)) #sum(is.na(dat$com_ac)) #sum(is.na(dat$agy_bf)) #sum(is.na(dat$com_bf)) # ##measurement invariance #dat2 <- setNames(data.frame(matrix(ncol = 13, nrow = 6)), c("model","fit_index","agy_ts","com_ts","agy_ac","com_ac", # "agy_bf","com_bf","agr","cns","ext","opn","neu")) #dat2$model[1:3] <- "config" #dat2$model[4:6] <- "metric" #dat2$fit_index[c(1,4)] <- "cfi" #dat2$fit_index[c(2,5)] <- "rmsea" #dat2$fit_index[c(3,6)] <- "srmr" # ##BIG TWO #mod_agy_ts <- 'agy =~ ext01 + ext03 + ext04 + ext05r + ext06 + ext07r + ext08 + opn01 # meth =~ 1*ext05r + 1*ext07r # agy ~~ 0*meth' #config_agy_ts <- cfa(mod_agy_ts, data = dat, missing = "ML", group = "sta_now") #dat2$agy_ts[1] <- fitMeasures(config_agy_ts, "cfi") #dat2$agy_ts[2] <- fitMeasures(config_agy_ts, "rmsea") #dat2$agy_ts[3] <- fitMeasures(config_agy_ts, "srmr") #metric_agy_ts <- cfa(mod_agy_ts, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$agy_ts[4] <- fitMeasures(metric_agy_ts, "cfi") #dat2$agy_ts[5] <- fitMeasures(metric_agy_ts, "rmsea") #dat2$agy_ts[6] <- fitMeasures(metric_agy_ts, "srmr") #rm(mod_agy_ts,config_agy_ts,metric_agy_ts) # #mod_com_ts <- 'com =~ agr01r + agr02 + agr03r + agr04 + agr07 + agr08r + agr09 + cns07 # meth =~ 1*agr01r + 1*agr03r + 1*agr08r # com ~~ 0*meth' #config_com_ts <- cfa(mod_com_ts, data = dat, missing = "ML", group = "sta_now") #dat2$com_ts[1] <- fitMeasures(config_com_ts, "cfi") #dat2$com_ts[2] <- fitMeasures(config_com_ts, "rmsea") #dat2$com_ts[3] <- fitMeasures(config_com_ts, "srmr") #metric_com_ts <- cfa(mod_com_ts, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$com_ts[4] <- fitMeasures(metric_com_ts, "cfi") #dat2$com_ts[5] <- fitMeasures(metric_com_ts, "rmsea") #dat2$com_ts[6] <- fitMeasures(metric_com_ts, "srmr") #rm(mod_com_ts,config_com_ts,metric_com_ts) # #mod_agy_ac <- 'agy =~ ext03 + ext06 + ext07r + opn04 + opn07r + neu07 # meth =~ 1*ext07r + 1*opn07r # agy ~~ 0*meth' #config_agy_ac <- cfa(mod_agy_ac, data = dat, missing = "ML", group = "sta_now") #dat2$agy_ac[1] <- fitMeasures(config_agy_ac, "cfi") #dat2$agy_ac[2] <- fitMeasures(config_agy_ac, "rmsea") #dat2$agy_ac[3] <- fitMeasures(config_agy_ac, "srmr") #metric_agy_ac <- cfa(mod_agy_ac, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$agy_ac[4] <- fitMeasures(metric_agy_ac, "cfi") #dat2$agy_ac[5] <- fitMeasures(metric_agy_ac, "rmsea") #dat2$agy_ac[6] <- fitMeasures(metric_agy_ac, "srmr") #rm(mod_agy_ac,config_agy_ac,metric_agy_ac) # #mod_com_ac <- 'com =~ agr01r + agr03r + agr05 + agr07 + agr08r + agr09 + cns01 # meth =~ 1*agr01r + 1*agr03r + 1*agr08r # com ~~ 0*meth' #config_com_ac <- cfa(mod_com_ac, data = dat, missing = "ML", group = "sta_now") #dat2$com_ac[1] <- fitMeasures(config_com_ac, "cfi") #dat2$com_ac[2] <- fitMeasures(config_com_ac, "rmsea") #dat2$com_ac[3] <- fitMeasures(config_com_ac, "srmr") #metric_com_ac <- cfa(mod_com_ac, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$com_ac[4] <- fitMeasures(metric_com_ac, "cfi") #dat2$com_ac[5] <- fitMeasures(metric_com_ac, "rmsea") #dat2$com_ac[6] <- fitMeasures(metric_com_ac, "srmr") #rm(mod_com_ac,config_com_ac,metric_com_ac) # #mod_agy_bf <- 'agy =~ cns06 + ext02r + ext05r + opn01 + opn03 + opn05 + opn10 + neu02 # meth =~ 1*ext02r + 1*ext05r # agy ~~ 0*meth' #config_agy_bf <- cfa(mod_agy_bf, data = dat, missing = "ML", group = "sta_now") #dat2$agy_bf[1] <- fitMeasures(config_agy_bf, "cfi") #dat2$agy_bf[2] <- fitMeasures(config_agy_bf, "rmsea") #dat2$agy_bf[3] <- fitMeasures(config_agy_bf, "srmr") #metric_agy_bf <- cfa(mod_agy_bf, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$agy_bf[4] <- fitMeasures(metric_agy_bf, "cfi") #dat2$agy_bf[5] <- fitMeasures(metric_agy_bf, "rmsea") #dat2$agy_bf[6] <- fitMeasures(metric_agy_bf, "srmr") #rm(mod_agy_bf,config_agy_bf,metric_agy_bf) # #mod_com_bf <- 'com =~ agr02 + agr03r + agr04 + agr06r + agr07 + agr08r + agr09 + cns01 # meth =~ 1*agr03r + 1*agr06r + 1*agr08r # com ~~ 0*meth' #config_com_bf <- cfa(mod_com_bf, data = dat, missing = "ML", group = "sta_now") #dat2$com_bf[1] <- fitMeasures(config_com_bf, "cfi") #dat2$com_bf[2] <- fitMeasures(config_com_bf, "rmsea") #dat2$com_bf[3] <- fitMeasures(config_com_bf, "srmr") #metric_com_bf <- cfa(mod_com_bf, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$com_bf[4] <- fitMeasures(metric_com_bf, "cfi") #dat2$com_bf[5] <- fitMeasures(metric_com_bf, "rmsea") #dat2$com_bf[6] <- fitMeasures(metric_com_bf, "srmr") #rm(mod_com_bf,config_com_bf,metric_com_bf) # ##BIG FIVE #mod_agr <- 'agr.lat =~ agr01r + agr02 + agr03r + agr04 + agr05 + agr06r + agr07 + agr08r + agr09 # meth =~ 1*agr01r + 1*agr03r + 1*agr06r + 1*agr08r # agr.lat ~~ 0*meth' #config_agr <- cfa(mod_agr, data = dat, missing = "ML", group = "sta_now") #dat2$agr[1] <- fitMeasures(config_agr, "cfi") #dat2$agr[2] <- fitMeasures(config_agr, "rmsea") #dat2$agr[3] <- fitMeasures(config_agr, "srmr") #metric_agr <- cfa(mod_agr, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$agr[4] <- fitMeasures(metric_agr, "cfi") #dat2$agr[5] <- fitMeasures(metric_agr, "rmsea") #dat2$agr[6] <- fitMeasures(metric_agr, "srmr") #rm(mod_agr,config_agr,metric_agr) # #mod_cns <- 'cns.lat =~ cns01 + cns02r + cns03 + cns04r + cns05r + cns06 + cns07 + cns08 + cns09r # meth =~ 1*cns02r + 1*cns04r + 1*cns05r + 1*cns09r # cns.lat ~~ 0*meth' #config_cns <- cfa(mod_cns, data = dat, missing = "ML", group = "sta_now") #dat2$cns[1] <- fitMeasures(config_cns, "cfi") #dat2$cns[2] <- fitMeasures(config_cns, "rmsea") #dat2$cns[3] <- fitMeasures(config_cns, "srmr") #metric_cns <- cfa(mod_cns, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$cns[4] <- fitMeasures(metric_cns, "cfi") #dat2$cns[5] <- fitMeasures(metric_cns, "rmsea") #dat2$cns[6] <- fitMeasures(metric_cns, "srmr") #rm(mod_cns,config_cns,metric_cns) # #mod_ext <- 'ext.lat =~ ext01 + ext02r + ext03 + ext04 + ext05r + ext06 + ext07r + ext08 # meth =~ 1*ext02r + 1*ext05r + 1*ext07r # ext.lat ~~ 0*meth' #config_ext <- cfa(mod_ext, data = dat, missing = "ML", group = "sta_now") #dat2$ext[1] <- fitMeasures(config_ext, "cfi") #dat2$ext[2] <- fitMeasures(config_ext, "rmsea") #dat2$ext[3] <- fitMeasures(config_ext, "srmr") #metric_ext <- cfa(mod_ext, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$ext[4] <- fitMeasures(metric_ext, "cfi") #dat2$ext[5] <- fitMeasures(metric_ext, "rmsea") #dat2$ext[6] <- fitMeasures(metric_ext, "srmr") #rm(mod_ext,config_ext,metric_ext) # #mod_opn <- 'opn.lat =~ opn01 + opn02 + opn03 + opn04 + opn05 + opn06 + opn07r + opn08 + opn09r + opn10 # meth =~ 1*opn07r + 1*opn09r # opn.lat ~~ 0*meth' #config_opn <- cfa(mod_opn, data = dat, missing = "ML", group = "sta_now") #dat2$opn[1] <- fitMeasures(config_opn, "cfi") #dat2$opn[2] <- fitMeasures(config_opn, "rmsea") #dat2$opn[3] <- fitMeasures(config_opn, "srmr") #metric_opn <- cfa(mod_opn, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$opn[4] <- fitMeasures(metric_opn, "cfi") #dat2$opn[5] <- fitMeasures(metric_opn, "rmsea") #dat2$opn[6] <- fitMeasures(metric_opn, "srmr") #rm(mod_opn,config_opn,metric_opn) # #mod_neu <- 'neu.lat =~ neu01 + neu02r + neu03 + neu04 + neu05r + neu06 + neu07r + neu08 # meth =~ 1*neu02r + 1*neu05r + 1*neu07r # neu.lat ~~ 0*meth' #config_neu <- cfa(mod_neu, data = dat, missing = "ML", group = "sta_now") #dat2$neu[1] <- fitMeasures(config_neu, "cfi") #dat2$neu[2] <- fitMeasures(config_neu, "rmsea") #dat2$neu[3] <- fitMeasures(config_neu, "srmr") #metric_neu <- cfa(mod_neu, data = dat, missing = "ML", group = "sta_now", group.equal = c("loadings")) #dat2$neu[4] <- fitMeasures(metric_neu, "cfi") #dat2$neu[5] <- fitMeasures(metric_neu, "rmsea") #dat2$neu[6] <- fitMeasures(metric_neu, "srmr") #rm(mod_neu,config_neu,metric_neu) # ##select variables #dat <- subset(dat, select = c(record_id, sta_now, age, female, langu, se, relig, agy_ts, com_ts, agy_ac, com_ac, agy_bf, com_bf, agr, cns, ext, opn, neu)) # ##save data #fwrite(dat, file = "PCmatch_csv-data_SOM-U2.csv") #fwrite(dat2, file = "PCmatch_mi_SOM-U2.csv") # #rm(dat2) #####################################Analyses with publicly available dataset##################################### #read data dat <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/330/6/PCmatch_csv-data_SOM-U2.csv", header = T, sep = ",")) #descripitives dat$female <- car::recode(dat$female, "-999=NA; 0=0; 1=1") 100*prop.table(table(dat$female)) dat$age <- as.numeric(dat$age) dat$age2 <- ifelse(dat$age > 150, NA, ifelse(dat$age < 0, NA, dat$age)) round(mean(dat$age2, na.rm=T),digits=2) round(sd(dat$age2, na.rm=T),digits=2) 100*prop.table(table(dat$langu)) #Measurement invariance Table dat1 <- fread("https://madata.bib.uni-mannheim.de/330/12/PCmatch_mi_SOM-U2.csv") dat1 <- dat1[, c(2,1,3:10,12,11,13)] dat1 <- add_row(dat1) dat1 <- add_row(dat1) dat1 <- dat1[c(1,4,7,10,2,5,8,11,3,6,9), ] dat1[,3:13] <- round(dat1[,3:13], digits=2) dat1[1,1] <- "CFI" dat1[2,1] <- "" dat1[3,1] <- "" dat1[3,2] <- intToUtf8(916) dat1[3,3:13] <- dat1[1,3:13] - dat1[2,3:13] dat1[5,1] <- "RMSEA" dat1[6,1] <- "" dat1[7,1] <- "" dat1[7,2] <- intToUtf8(916) dat1[7,3:13] <- dat1[5,3:13] - dat1[6,3:13] dat1[9,1] <- "SRMR" dat1[10,1] <- "" dat1[11,1] <- "" dat1[11,2] <- intToUtf8(916) dat1[11,3:13] <- dat1[9,3:13] - dat1[10,3:13] dat1[,3:13] <- mutate_all(data.frame(dat1[,3:13]), list(~f_num(unlist(.), digits=2))) dat1$blank1 <- NA dat1$blank2 <- NA dat1$blank3 <- NA dat1 <- dat1[, c(1:4,14,5:6,15,7:8,16,9:13)] col_keys <- c("fit_index","model","agy_ts","com_ts","blank1","agy_ac","com_ac","blank2","agy_bf","com_bf","blank3","agr","cns","opn","ext","neu") head1 <- c("","","Big Two","Big Two","Big Two","Big Two","Big Two","Big Two","Big Two","Big Two","","Big Five","Big Five","Big Five","Big Five","Big Five") head2 <- c("","","TS","TS","","AC","AC","","BF","BF","","","","","","") head3 <- c("","","agy","com","","agy","com","","agy","com","","agr","cns","opn","ext","neu") head <- data.frame(col_keys,head1,head2,head3, stringsAsFactors = FALSE) rm(col_keys,head1,head2,head3) tbl <- flextable(dat1) tbl <- set_header_df(tbl, mapping=head, key="col_keys") tbl <- colformat_lgl(tbl, j =~ blank1 + blank2 + blank3, na_str="") tbl <- hline_top(tbl, j=1:16, border=fp_border(width=2), part="header") tbl <- merge_at(tbl, i=1, j=3:10, part="header") tbl <- merge_at(tbl, i=1, j=12:16, part="header") tbl <- merge_at(tbl, i=2, j=3:4, part="header") tbl <- merge_at(tbl, i=2, j=6:7, part="header") tbl <- merge_at(tbl, i=2, j=9:10, part="header") tbl <- hline(tbl, i=1, j=3:10, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=1, j=12:16, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=3:4, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=6:7, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=9:10, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=3, j=3:10, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=3, j=12:16, border=fp_border(width=1), part="header") tbl <- font(tbl, fontname="Times", part="all") tbl <- fontsize(tbl, size=9, part="all") tbl <- bold(tbl, part="header") tbl <- align(tbl, align="center", part="all") tbl <- align(tbl, j = c("fit_index","model"), align="left", part="body") tbl <- width(tbl, j =~ fit_index, width=1) tbl <- width(tbl, j =~ model, width=.7) tbl <- width(tbl, j =~ agy_ts + com_ts + agy_ac + com_ac + agy_bf + com_bf + agr + cns + opn + ext + neu, width=.5) tbl <- width(tbl, j =~ blank1 + blank2, width=.1) tbl <- width(tbl, j =~ blank3, width=.3) tbl <- height(tbl, i=c(4,8), height=.1, part="body") tbl doc <- read_docx() doc <- body_add_flextable(doc, value = tbl) print(doc, target = "s2_mi_table.docx") rm(dat1,head,tbl,doc) #prepare data for mixed-effects models ##z-standardize self-esteem dat$zse <- scale(dat$se) ##group-mean center and z-standardize religiosity and personality traits dat$zrelgrp <- scale(dat$relig - ave(dat$relig, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zagytsgrp <- scale(dat$agy_ts - ave(dat$agy_ts, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zcomtsgrp <- scale(dat$com_ts - ave(dat$com_ts, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zagyacgrp <- scale(dat$agy_ac - ave(dat$agy_ac, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zcomacgrp <- scale(dat$com_ac - ave(dat$com_ac, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zagybfgrp <- scale(dat$agy_bf - ave(dat$agy_bf, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zcombfgrp <- scale(dat$com_bf - ave(dat$com_bf, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zagrgrp <- scale(dat$agr - ave(dat$agr, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zcnsgrp <- scale(dat$cns - ave(dat$cns, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zextgrp <- scale(dat$ext - ave(dat$ext, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zopngrp <- scale(dat$opn - ave(dat$opn, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) dat$zneugrp <- scale(dat$neu - ave(dat$neu, dat$sta_now, FUN = function(x) mean(x, na.rm=T))) ##compute state-level religiosiy (z-standardized) dat1 <- aggregate(relig ~ sta_now, dat, mean) dat1$zsrel <- scale(dat1$relig) dat1$relig <- NULL dat <- merge(dat, dat1, by = "sta_now", all = TRUE) #correlation between computed state-level religiosiy and #state-level religiosity index based on Gallup U.S. Poll data (Diener, Tay, & Myers, 2011) dat2 <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/330/18/srel_gallup.csv", header = T, sep = ',')) dat3 <- merge(dat1, dat2, by = "sta_now") dat3 <- subset(dat3, select = c(zsrel,srel_gallup)) print(corr.test(dat3), short=F) rm(dat1,dat2,dat3) ##re-center personality traits and state-level religiosity for simple slope analyses dat$zagytsgrpverylow <- dat$zagytsgrp + 2 dat$zagytsgrplow <- dat$zagytsgrp + 1 dat$zagytsgrphigh <- dat$zagytsgrp - 1 dat$zagytsgrpveryhigh <- dat$zagytsgrp - 2 dat$zcomtsgrpverylow <- dat$zcomtsgrp + 2 dat$zcomtsgrplow <- dat$zcomtsgrp + 1 dat$zcomtsgrphigh <- dat$zcomtsgrp - 1 dat$zcomtsgrpveryhigh <- dat$zcomtsgrp - 2 dat$zagyacgrpverylow <- dat$zagyacgrp + 2 dat$zagyacgrplow <- dat$zagyacgrp + 1 dat$zagyacgrphigh <- dat$zagyacgrp - 1 dat$zagyacgrpveryhigh <- dat$zagyacgrp - 2 dat$zcomacgrpverylow <- dat$zcomacgrp + 2 dat$zcomacgrplow <- dat$zcomacgrp + 1 dat$zcomacgrphigh <- dat$zcomacgrp - 1 dat$zcomacgrpveryhigh <- dat$zcomacgrp - 2 dat$zagybfgrpverylow <- dat$zagybfgrp + 2 dat$zagybfgrplow <- dat$zagybfgrp + 1 dat$zagybfgrphigh <- dat$zagybfgrp - 1 dat$zagybfgrpveryhigh <- dat$zagybfgrp - 2 dat$zcombfgrpverylow <- dat$zcombfgrp + 2 dat$zcombfgrplow <- dat$zcombfgrp + 1 dat$zcombfgrphigh <- dat$zcombfgrp - 1 dat$zcombfgrpveryhigh <- dat$zcombfgrp - 2 dat$zagrgrpverylow <- dat$zagrgrp + 2 dat$zagrgrplow <- dat$zagrgrp + 1 dat$zagrgrphigh <- dat$zagrgrp - 1 dat$zagrgrpveryhigh <- dat$zagrgrp - 2 dat$zcnsgrpverylow <- dat$zcnsgrp + 2 dat$zcnsgrplow <- dat$zcnsgrp + 1 dat$zcnsgrphigh <- dat$zcnsgrp - 1 dat$zcnsgrpveryhigh <- dat$zcnsgrp - 2 dat$zextgrpverylow <- dat$zextgrp + 2 dat$zextgrplow <- dat$zextgrp + 1 dat$zextgrphigh <- dat$zextgrp - 1 dat$zextgrpveryhigh <- dat$zextgrp - 2 dat$zopngrpverylow <- dat$zopngrp + 2 dat$zopngrplow <- dat$zopngrp + 1 dat$zopngrphigh <- dat$zopngrp - 1 dat$zopngrpveryhigh <- dat$zopngrp - 2 dat$zneugrpverylow <- dat$zneugrp + 2 dat$zneugrplow <- dat$zneugrp + 1 dat$zneugrphigh <- dat$zneugrp - 1 dat$zneugrpveryhigh <- dat$zneugrp - 2 dat$zsrelmax <- dat$zsrel - abs(max(dat$zsrel)) dat$zsrelmin <- dat$zsrel + abs(min(dat$zsrel)) #select variables dat <- subset(dat, select = c(record_id, age, female, langu, sta_now, zse:zsrelmin)) #####Before you continue with this R-Syntax, download and install Julia (https://julialang.org/downloads/)##### #install JuliaCall install.packages("JuliaCall") #load JuliaCall library(JuliaCall) #start Julia julia <- julia_setup() #set working directory to the Julia folder indicated in julia_setup (e.g., "C:/AppData/Local/JULIA-~1.1/bin") setwd("C:/AppData/Local/JULIA-~1.1/bin") #install Julia packages julia_install_package("MixedModels") julia_install_package("DataFrames") julia_install_package("StatsBase") #load Julia packages julia_library("MixedModels") julia_library("DataFrames") julia_library("StatsBase") #conduct mixed-effects models in Julia dat$sta_now <- as.factor(dat$sta_now) julia_assign("dat", dat) ##BIG TWO TS julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrel * zcomtsgrp + zrelgrp * zsrel * zagytsgrp + (1 + zrelgrp + zcomtsgrp + zagytsgrp | sta_now)), dat)") BigTwoTS <- julia_eval("BigTwoTS = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomtsgrphigh + zrelgrp * zsrelmax * zagytsgrplow + (1 + zrelgrp + zcomtsgrphigh + zagytsgrplow | sta_now)), dat)") BigTwoTS_cmax_asshigh <- julia_eval("BigTwoTS_cmax_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomtsgrpveryhigh + zrelgrp * zsrelmax * zagytsgrpverylow + (1 + zrelgrp + zcomtsgrpveryhigh + zagytsgrpverylow | sta_now)), dat)") BigTwoTS_cmax_assveryhigh <- julia_eval("BigTwoTS_cmax_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomtsgrp + zrelgrp * zsrelmax * zagytsgrp + (1 + zrelgrp + zcomtsgrp + zagytsgrp | sta_now)), dat)") BigTwoTS_cmax_average <- julia_eval("BigTwoTS_cmax_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomtsgrplow + zrelgrp * zsrelmax * zagytsgrphigh + (1 + zrelgrp + zcomtsgrplow + zagytsgrphigh | sta_now)), dat)") BigTwoTS_cmax_conthigh <- julia_eval("BigTwoTS_cmax_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomtsgrpverylow + zrelgrp * zsrelmax * zagytsgrpveryhigh + (1 + zrelgrp + zcomtsgrpverylow + zagytsgrpveryhigh | sta_now)), dat)") BigTwoTS_cmax_contveryhigh <- julia_eval("BigTwoTS_cmax_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomtsgrphigh + zrelgrp * zsrelmin * zagytsgrplow + (1 + zrelgrp + zcomtsgrphigh + zagytsgrplow | sta_now)), dat)") BigTwoTS_cmin_asshigh <- julia_eval("BigTwoTS_cmin_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomtsgrpveryhigh + zrelgrp * zsrelmin * zagytsgrpverylow + (1 + zrelgrp + zcomtsgrpveryhigh + zagytsgrpverylow | sta_now)), dat)") BigTwoTS_cmin_assveryhigh <- julia_eval("BigTwoTS_cmin_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomtsgrp + zrelgrp * zsrelmin * zagytsgrp + (1 + zrelgrp + zcomtsgrp + zagytsgrp | sta_now)), dat)") BigTwoTS_cmin_average <- julia_eval("BigTwoTS_cmin_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomtsgrplow + zrelgrp * zsrelmin * zagytsgrphigh + (1 + zrelgrp + zcomtsgrplow + zagytsgrphigh | sta_now)), dat)") BigTwoTS_cmin_conthigh <- julia_eval("BigTwoTS_cmin_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomtsgrpverylow + zrelgrp * zsrelmin * zagytsgrpveryhigh + (1 + zrelgrp + zcomtsgrpverylow + zagytsgrpveryhigh | sta_now)), dat)") BigTwoTS_cmin_contveryhigh <- julia_eval("BigTwoTS_cmin_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") ##BIG TWO AC julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrel * zcomacgrp + zrelgrp * zsrel * zagyacgrp + (1 + zrelgrp + zcomacgrp + zagyacgrp | sta_now)), dat)") BigTwoAC <- julia_eval("BigTwoAC = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomacgrphigh + zrelgrp * zsrelmax * zagyacgrplow + (1 + zrelgrp + zcomacgrphigh + zagyacgrplow | sta_now)), dat)") BigTwoAC_cmax_asshigh <- julia_eval("BigTwoAC_cmax_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomacgrpveryhigh + zrelgrp * zsrelmax * zagyacgrpverylow + (1 + zrelgrp + zcomacgrpveryhigh + zagyacgrpverylow | sta_now)), dat)") BigTwoAC_cmax_assveryhigh <- julia_eval("BigTwoAC_cmax_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomacgrp + zrelgrp * zsrelmax * zagyacgrp + (1 + zrelgrp + zcomacgrp + zagyacgrp | sta_now)), dat)") BigTwoAC_cmax_average <- julia_eval("BigTwoAC_cmax_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomacgrplow + zrelgrp * zsrelmax * zagyacgrphigh + (1 + zrelgrp + zcomacgrplow + zagyacgrphigh | sta_now)), dat)") BigTwoAC_cmax_conthigh <- julia_eval("BigTwoAC_cmax_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcomacgrpverylow + zrelgrp * zsrelmax * zagyacgrpveryhigh + (1 + zrelgrp + zcomacgrpverylow + zagyacgrpveryhigh | sta_now)), dat)") BigTwoAC_cmax_contveryhigh <- julia_eval("BigTwoAC_cmax_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=7)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomacgrphigh + zrelgrp * zsrelmin * zagyacgrplow + (1 + zrelgrp + zcomacgrphigh + zagyacgrplow | sta_now)), dat)") BigTwoAC_cmin_asshigh <- julia_eval("BigTwoAC_cmin_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomacgrpveryhigh + zrelgrp * zsrelmin * zagyacgrpverylow + (1 + zrelgrp + zcomacgrpveryhigh + zagyacgrpverylow | sta_now)), dat)") BigTwoAC_cmin_assveryhigh <- julia_eval("BigTwoAC_cmin_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomacgrp + zrelgrp * zsrelmin * zagyacgrp + (1 + zrelgrp + zcomacgrp + zagyacgrp | sta_now)), dat)") BigTwoAC_cmin_average <- julia_eval("BigTwoAC_cmin_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomacgrplow + zrelgrp * zsrelmin * zagyacgrphigh + (1 + zrelgrp + zcomacgrplow + zagyacgrphigh | sta_now)), dat)") BigTwoAC_cmin_conthigh <- julia_eval("BigTwoAC_cmin_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcomacgrpverylow + zrelgrp * zsrelmin * zagyacgrpveryhigh + (1 + zrelgrp + zcomacgrpverylow + zagyacgrpveryhigh | sta_now)), dat)") BigTwoAC_cmin_contveryhigh <- julia_eval("BigTwoAC_cmin_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=6)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=6)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=6))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") ##BIG TWO BF julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrel * zcombfgrp + zrelgrp * zsrel * zagybfgrp + (1 + zrelgrp + zcombfgrp + zagybfgrp | sta_now)), dat)") BigTwoBF <- julia_eval("BigTwoBF = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcombfgrphigh + zrelgrp * zsrelmax * zagybfgrplow + (1 + zrelgrp + zcombfgrphigh + zagybfgrplow | sta_now)), dat)") BigTwoBF_cmax_asshigh <- julia_eval("BigTwoBF_cmax_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcombfgrpveryhigh + zrelgrp * zsrelmax * zagybfgrpverylow + (1 + zrelgrp + zcombfgrpveryhigh + zagybfgrpverylow | sta_now)), dat)") BigTwoBF_cmax_assveryhigh <- julia_eval("BigTwoBF_cmax_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcombfgrp + zrelgrp * zsrelmax * zagybfgrp + (1 + zrelgrp + zcombfgrp + zagybfgrp | sta_now)), dat)") BigTwoBF_cmax_average <- julia_eval("BigTwoBF_cmax_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcombfgrplow + zrelgrp * zsrelmax * zagybfgrphigh + (1 + zrelgrp + zcombfgrplow + zagybfgrphigh | sta_now)), dat)") BigTwoBF_cmax_conthigh <- julia_eval("BigTwoBF_cmax_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zcombfgrpverylow + zrelgrp * zsrelmax * zagybfgrpveryhigh + (1 + zrelgrp + zcombfgrpverylow + zagybfgrpveryhigh | sta_now)), dat)") BigTwoBF_cmax_contveryhigh <- julia_eval("BigTwoBF_cmax_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcombfgrphigh + zrelgrp * zsrelmin * zagybfgrplow + (1 + zrelgrp + zcombfgrphigh + zagybfgrplow | sta_now)), dat)") BigTwoBF_cmin_asshigh <- julia_eval("BigTwoBF_cmin_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcombfgrpveryhigh + zrelgrp * zsrelmin * zagybfgrpverylow + (1 + zrelgrp + zcombfgrpveryhigh + zagybfgrpverylow | sta_now)), dat)") BigTwoBF_cmin_assveryhigh <- julia_eval("BigTwoBF_cmin_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcombfgrp + zrelgrp * zsrelmin * zagybfgrp + (1 + zrelgrp + zcombfgrp + zagybfgrp | sta_now)), dat)") BigTwoBF_cmin_average <- julia_eval("BigTwoBF_cmin_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcombfgrplow + zrelgrp * zsrelmin * zagybfgrphigh + (1 + zrelgrp + zcombfgrplow + zagybfgrphigh | sta_now)), dat)") BigTwoBF_cmin_conthigh <- julia_eval("BigTwoBF_cmin_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zcombfgrpverylow + zrelgrp * zsrelmin * zagybfgrpveryhigh + (1 + zrelgrp + zcombfgrpverylow + zagybfgrpveryhigh | sta_now)), dat)") BigTwoBF_cmin_contveryhigh <- julia_eval("BigTwoBF_cmin_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") ##BIG FIVE julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrel * zagrgrp + zrelgrp * zsrel * zcnsgrp + zrelgrp * zsrel * zopngrp + zrelgrp * zsrel * zextgrp + zrelgrp * zsrel * zneugrp + (1 + zrelgrp + zagrgrp + zcnsgrp + zopngrp + zextgrp + zneugrp | sta_now)), dat)") BigFive <- julia_eval("BigFive = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zagrgrphigh + zrelgrp * zsrelmax * zcnsgrplow + zrelgrp * zsrelmax * zopngrplow + zrelgrp * zsrelmax * zextgrplow + zrelgrp * zsrelmax * zneugrphigh + (1 + zrelgrp + zagrgrphigh + zcnsgrplow + zopngrplow + zextgrplow + zneugrphigh | sta_now)), dat)") BigFive_cmax_asshigh <- julia_eval("BigFive_cmax_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zagrgrpveryhigh + zrelgrp * zsrelmax * zcnsgrpverylow + zrelgrp * zsrelmax * zopngrpverylow + zrelgrp * zsrelmax * zextgrpverylow + zrelgrp * zsrelmax * zneugrpveryhigh + (1 + zrelgrp + zagrgrpveryhigh + zcnsgrpverylow + zopngrpverylow + zextgrpverylow + zneugrpveryhigh | sta_now)), dat)") BigFive_cmax_assveryhigh <- julia_eval("BigFive_cmax_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zagrgrp + zrelgrp * zsrelmax * zcnsgrp + zrelgrp * zsrelmax * zopngrp + zrelgrp * zsrelmax * zextgrp + zrelgrp * zsrelmax * zneugrp + (1 + zrelgrp + zagrgrp + zcnsgrp + zopngrp + zextgrp + zneugrp | sta_now)), dat)") BigFive_cmax_average <- julia_eval("BigFive_cmax_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zagrgrplow + zrelgrp * zsrelmax * zcnsgrphigh + zrelgrp * zsrelmax * zopngrphigh + zrelgrp * zsrelmax * zextgrphigh + zrelgrp * zsrelmax * zneugrplow + (1 + zrelgrp + zagrgrplow + zcnsgrphigh + zopngrphigh + zextgrphigh + zneugrplow | sta_now)), dat)") BigFive_cmax_conthigh <- julia_eval("BigFive_cmax_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmax * zagrgrpverylow + zrelgrp * zsrelmax * zcnsgrpveryhigh + zrelgrp * zsrelmax * zopngrpveryhigh + zrelgrp * zsrelmax * zextgrpveryhigh + zrelgrp * zsrelmax * zneugrpverylow + (1 + zrelgrp + zagrgrpverylow + zcnsgrpveryhigh + zopngrpveryhigh + zextgrpveryhigh + zneugrpverylow | sta_now)), dat)") BigFive_cmax_contveryhigh <- julia_eval("BigFive_cmax_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zagrgrphigh + zrelgrp * zsrelmin * zcnsgrplow + zrelgrp * zsrelmin * zopngrplow + zrelgrp * zsrelmin * zextgrplow + zrelgrp * zsrelmin * zneugrphigh + (1 + zrelgrp + zagrgrphigh + zcnsgrplow + zopngrplow + zextgrplow + zneugrphigh | sta_now)), dat)") BigFive_cmin_asshigh <- julia_eval("BigFive_cmin_asshigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zagrgrpveryhigh + zrelgrp * zsrelmin * zcnsgrpverylow + zrelgrp * zsrelmin * zopngrpverylow + zrelgrp * zsrelmin * zextgrpverylow + zrelgrp * zsrelmin * zneugrpveryhigh + (1 + zrelgrp + zagrgrpveryhigh + zcnsgrpverylow + zopngrpverylow + zextgrpverylow + zneugrpveryhigh | sta_now)), dat)") BigFive_cmin_assveryhigh <- julia_eval("BigFive_cmin_assveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zagrgrp + zrelgrp * zsrelmin * zcnsgrp + zrelgrp * zsrelmin * zopngrp + zrelgrp * zsrelmin * zextgrp + zrelgrp * zsrelmin * zneugrp + (1 + zrelgrp + zagrgrp + zcnsgrp + zopngrp + zextgrp + zneugrp | sta_now)), dat)") BigFive_cmin_average <- julia_eval("BigFive_cmin_average = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zagrgrplow + zrelgrp * zsrelmin * zcnsgrphigh + zrelgrp * zsrelmin * zopngrphigh + zrelgrp * zsrelmin * zextgrphigh + zrelgrp * zsrelmin * zneugrplow + (1 + zrelgrp + zagrgrplow + zcnsgrphigh + zopngrphigh + zextgrphigh + zneugrplow | sta_now)), dat)") BigFive_cmin_conthigh <- julia_eval("BigFive_cmin_conthigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zsrelmin * zagrgrpverylow + zrelgrp * zsrelmin * zcnsgrpveryhigh + zrelgrp * zsrelmin * zopngrpveryhigh + zrelgrp * zsrelmin * zextgrpveryhigh + zrelgrp * zsrelmin * zneugrpverylow + (1 + zrelgrp + zagrgrpverylow + zcnsgrpveryhigh + zopngrpveryhigh + zextgrpveryhigh + zneugrpverylow | sta_now)), dat)") BigFive_cmin_contveryhigh <- julia_eval("BigFive_cmin_contveryhigh = rename!(DataFrame([(coeftable(fm).rownms) (round.(fixef(fm, false),digits=5)) (round.(((fixef(fm, false))-1.96*(stderror(fm))),digits=5)) (round.(((fixef(fm, false))+1.96*(stderror(fm))),digits=5))]), :x1 => :predictor, :x2 => :zPe, :x3 => :lowerCI, :x4 => :upperCI)") #set working directory to the folder to which files created in R should be saved (e.g., "C:/PCmatch") setwd("C:/PCmatch") #Table SOM-U2 BigTwoTS[, c(2:4)] <- sapply(BigTwoTS[, c(2:4)], as.numeric) BigTwoAC[, c(2:4)] <- sapply(BigTwoAC[, c(2:4)], as.numeric) BigTwoBF[, c(2:4)] <- sapply(BigTwoBF[, c(2:4)], as.numeric) BigFive[, c(2:4)] <- sapply(BigFive[, c(2:4)], as.numeric) dat1 <- data.frame(BigTwoTS,BigTwoAC,BigTwoBF,BigFive) dat1[,5] <- NULL dat1[,8] <- NULL dat1[,1] <- c("(intercept)","rel", "sta rel", "com", "agy", "rel x sta rel", "rel x com", "sta rel x com", "rel x agy", "sta rel x agy", "rel x sta rel x com", "rel x sta rel x agy", "", "", "", "", "", "", "", "", "", "", "", "") dat1[,11] <- c("(intercept)","rel", "sta rel", "agr", "cns", "opn", "ext", "neu", "rel x sta rel", "rel x agr", "sta rel x agr", "rel x cns", "sta rel x cns", "rel x opn", "sta rel x opn", "rel x ext", "sta rel x ext", "rel x neu", "sta rel x neu", "rel x sta rel x agr", "rel x sta rel x cns", "rel x sta rel x opn", "rel x sta rel x ext", "rel x sta rel x neu") trunc0e <- function(x){ ifelse(abs(x)< 0.0005, formatC(x, format = "e", digits = 0), f_num(unlist(x), digits=3)) } dat1 <- dat1 %>% mutate_if(is.numeric, trunc0e) dat1$ci <- paste0("[", dat1$lowerCI, ", ", dat1$upperCI, "]") dat1$ci.1 <- paste0("[", dat1$lowerCI.1, ", ", dat1$upperCI.1, "]") dat1$ci.2 <- paste0("[", dat1$lowerCI.2, ", ", dat1$upperCI.2, "]") dat1$ci.3 <- paste0("[", dat1$lowerCI.3, ", ", dat1$upperCI.3, "]") dat1$blank <- NA dat1$blank.1 <- NA dat1$blank.2 <- NA dat1 <- dat1[c("predictor","zPe","ci","blank","zPe.1","ci.1","blank.1","zPe.2","ci.2","blank.2","predictor.3","zPe.3","ci.3")] dat1[13:24,1:10] <- NA col_keys <- c("predictor","zPe","ci","blank","zPe.1","ci.1","blank.1","zPe.2","ci.2","blank.2","predictor.3","zPe.3","ci.3") head1 <- c("Big Two","Big Two","Big Two","","Big Two","Big Two","","Big Two","Big Two","","Big Five","Big Five","Big Five") head2 <- c("","TS","TS","","AC","AC","","BF","BF","","","","") head3 <- c(""," zPE"," 95% CI",""," zPE"," 95% CI",""," zPE"," 95% CI","",""," zPE"," 95% CI") head <- data.frame(col_keys,head1,head2,head3, stringsAsFactors = FALSE) rm(col_keys,head1,head2,head3) tbl <- flextable(dat1) tbl <- set_header_df(tbl, mapping=head, key="col_keys") tbl <- hline_top(tbl, j=1:13, border=fp_border(width=2), part="header") tbl <- merge_at(tbl, i=1, j=1:9, part="header") tbl <- merge_at(tbl, i=1, j=11:13, part="header") tbl <- merge_at(tbl, i=2, j=2:3, part="header") tbl <- merge_at(tbl, i=2, j=5:6, part="header") tbl <- merge_at(tbl, i=2, j=8:9, part="header") tbl <- hline(tbl, i=1, j=1:3, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=1, j=5:6, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=1, j=8:9, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=1, j=11:13, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=2, j=2:3, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=2, j=5:6, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=2, j=8:9, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=3, j=1:9, border=fp_border(width=1.2), part="header") tbl <- hline(tbl, i=3, j=11:13, border=fp_border(width=1.2), part="header") tbl <- font(tbl, fontname="Times", part="all") tbl <- fontsize(tbl, size=11, part="all") tbl <- bold(tbl, part="header") tbl <- italic(tbl, i=3, j=c("zPe","zPe.1","zPe.2","zPe.3"), part="header") tbl <- align(tbl, align="center", part="all") tbl <- align(tbl, j = c("ci","ci.1","ci.2","ci.3"), align="left", part="header") tbl <- align(tbl, j = c("predictor","predictor.3","ci","ci.1","ci.2","ci.3"), align="left", part="body") tbl <- align(tbl, j = c("zPe","zPe.1","zPe.2","zPe.3"), align="right", part="body") tbl <- width(tbl, j =~ predictor + predictor.3, width=1.3) tbl <- width(tbl, j =~ zPe + zPe.1 + zPe.2 + zPe.3, width=.67) tbl <- width(tbl, j =~ ci + ci.1 + ci.2 + ci.3, width=1.1) tbl <- width(tbl, j =~ blank + blank.1 + blank.2, width=.1) tbl <- colformat_lgl(tbl, j =~ blank + blank.1 + blank.2, na_str="") tbl <- height_all(tbl, height=.1, part="all") tbl doc <- read_docx() doc <- body_add_flextable(doc, value = tbl) print(doc, target = "Table_SOM-U2.docx") rm(dat1,head,tbl,doc) #simple slopes for the religiosity x state-level religiosity interactions ##mean association between religiosity and self-esteem in least religious states BigTwoTS_cmin_average[, c(2:4)] <- sapply(BigTwoTS_cmin_average[, c(2:4)], as.numeric) BigTwoAC_cmin_average[, c(2:4)] <- sapply(BigTwoAC_cmin_average[, c(2:4)], as.numeric) BigTwoBF_cmin_average[, c(2:4)] <- sapply(BigTwoBF_cmin_average[, c(2:4)], as.numeric) BigFive_cmin_average[, c(2:4)] <- sapply(BigFive_cmin_average[, c(2:4)], as.numeric) cmin_zpe <- (BigTwoTS_cmin_average[2,2] + BigTwoAC_cmin_average[2,2] + BigTwoBF_cmin_average[2,2] + BigFive_cmin_average[2,2]) / 4 cmin_lowerCI <- (BigTwoTS_cmin_average[2,3] + BigTwoAC_cmin_average[2,3] + BigTwoBF_cmin_average[2,3] + BigFive_cmin_average[2,3]) / 4 cmin_upperCI <- (BigTwoTS_cmin_average[2,4] + BigTwoAC_cmin_average[2,4] + BigTwoBF_cmin_average[2,4] + BigFive_cmin_average[2,4]) / 4 paste0("least religious states: zPE = ", round(cmin_zpe, digits=2), " 95% CI [", round(cmin_lowerCI, digits=4), ", ", round(cmin_upperCI, digits=2), "]") ##mean association between religiosity and self-esteem in most religious states BigTwoTS_cmax_average[, c(2:4)] <- sapply(BigTwoTS_cmax_average[, c(2:4)], as.numeric) BigTwoAC_cmax_average[, c(2:4)] <- sapply(BigTwoAC_cmax_average[, c(2:4)], as.numeric) BigTwoBF_cmax_average[, c(2:4)] <- sapply(BigTwoBF_cmax_average[, c(2:4)], as.numeric) BigFive_cmax_average[, c(2:4)] <- sapply(BigFive_cmax_average[, c(2:4)], as.numeric) cmax_zpe <- (BigTwoTS_cmax_average[2,2] + BigTwoAC_cmax_average[2,2] + BigTwoBF_cmax_average[2,2] + BigFive_cmax_average[2,2]) / 4 cmax_lowerCI <- (BigTwoTS_cmax_average[2,3] + BigTwoAC_cmax_average[2,3] + BigTwoBF_cmax_average[2,3] + BigFive_cmax_average[2,3]) / 4 cmax_upperCI <- (BigTwoTS_cmax_average[2,4] + BigTwoAC_cmax_average[2,4] + BigTwoBF_cmax_average[2,4] + BigFive_cmax_average[2,4]) / 4 paste0("most religious states: zPE = ", round(cmax_zpe, digits=2), " 95% CI [", round(cmax_lowerCI, digits=2), ", ", round(cmax_upperCI, digits=2), "]") #mean size of the person-culture match effect power_culture <- round(cmax_zpe, digits=2) - round(cmin_zpe, digits=2) power_lowerCI <- round(cmax_lowerCI, digits=2) - round(cmin_upperCI, digits=2) power_upperCI <- round(cmax_upperCI, digits=2) - round(cmin_lowerCI, digits=3) paste0("power of culture: delta_zPE = ", round(power_culture, digits=2), " 95% CI [", round(power_lowerCI, digits=2), ", ", round(power_upperCI, digits=2), "]") rm(cmin_zpe,cmin_lowerCI,cmin_upperCI,cmax_zpe,cmax_lowerCI,cmax_upperCI,power_culture,power_lowerCI,power_upperCI) #Figure SOM-U2 BigTwoTS_cmax_assveryhigh[, c(2:4)] <- sapply(BigTwoTS_cmax_assveryhigh[, c(2:4)], as.numeric) BigTwoTS_cmax_asshigh[, c(2:4)] <- sapply(BigTwoTS_cmax_asshigh[, c(2:4)], as.numeric) BigTwoTS_cmax_average[, c(2:4)] <- sapply(BigTwoTS_cmax_average[, c(2:4)], as.numeric) BigTwoTS_cmax_conthigh[, c(2:4)] <- sapply(BigTwoTS_cmax_conthigh[, c(2:4)], as.numeric) BigTwoTS_cmax_contveryhigh[, c(2:4)] <- sapply(BigTwoTS_cmax_contveryhigh[, c(2:4)], as.numeric) BigTwoTS_cmin_assveryhigh[, c(2:4)] <- sapply(BigTwoTS_cmin_assveryhigh[, c(2:4)], as.numeric) BigTwoTS_cmin_asshigh[, c(2:4)] <- sapply(BigTwoTS_cmin_asshigh[, c(2:4)], as.numeric) BigTwoTS_cmin_average[, c(2:4)] <- sapply(BigTwoTS_cmin_average[, c(2:4)], as.numeric) BigTwoTS_cmin_conthigh[, c(2:4)] <- sapply(BigTwoTS_cmin_conthigh[, c(2:4)], as.numeric) BigTwoTS_cmin_contveryhigh[, c(2:4)] <- sapply(BigTwoTS_cmin_contveryhigh[, c(2:4)], as.numeric) BigTwoAC_cmax_assveryhigh[, c(2:4)] <- sapply(BigTwoAC_cmax_assveryhigh[, c(2:4)], as.numeric) BigTwoAC_cmax_asshigh[, c(2:4)] <- sapply(BigTwoAC_cmax_asshigh[, c(2:4)], as.numeric) BigTwoAC_cmax_average[, c(2:4)] <- sapply(BigTwoAC_cmax_average[, c(2:4)], as.numeric) BigTwoAC_cmax_conthigh[, c(2:4)] <- sapply(BigTwoAC_cmax_conthigh[, c(2:4)], as.numeric) BigTwoAC_cmax_contveryhigh[, c(2:4)] <- sapply(BigTwoAC_cmax_contveryhigh[, c(2:4)], as.numeric) BigTwoAC_cmin_assveryhigh[, c(2:4)] <- sapply(BigTwoAC_cmin_assveryhigh[, c(2:4)], as.numeric) BigTwoAC_cmin_asshigh[, c(2:4)] <- sapply(BigTwoAC_cmin_asshigh[, c(2:4)], as.numeric) BigTwoAC_cmin_average[, c(2:4)] <- sapply(BigTwoAC_cmin_average[, c(2:4)], as.numeric) BigTwoAC_cmin_conthigh[, c(2:4)] <- sapply(BigTwoAC_cmin_conthigh[, c(2:4)], as.numeric) BigTwoAC_cmin_contveryhigh[, c(2:4)] <- sapply(BigTwoAC_cmin_contveryhigh[, c(2:4)], as.numeric) BigTwoBF_cmax_assveryhigh[, c(2:4)] <- sapply(BigTwoBF_cmax_assveryhigh[, c(2:4)], as.numeric) BigTwoBF_cmax_asshigh[, c(2:4)] <- sapply(BigTwoBF_cmax_asshigh[, c(2:4)], as.numeric) BigTwoBF_cmax_average[, c(2:4)] <- sapply(BigTwoBF_cmax_average[, c(2:4)], as.numeric) BigTwoBF_cmax_conthigh[, c(2:4)] <- sapply(BigTwoBF_cmax_conthigh[, c(2:4)], as.numeric) BigTwoBF_cmax_contveryhigh[, c(2:4)] <- sapply(BigTwoBF_cmax_contveryhigh[, c(2:4)], as.numeric) BigTwoBF_cmin_assveryhigh[, c(2:4)] <- sapply(BigTwoBF_cmin_assveryhigh[, c(2:4)], as.numeric) BigTwoBF_cmin_asshigh[, c(2:4)] <- sapply(BigTwoBF_cmin_asshigh[, c(2:4)], as.numeric) BigTwoBF_cmin_average[, c(2:4)] <- sapply(BigTwoBF_cmin_average[, c(2:4)], as.numeric) BigTwoBF_cmin_conthigh[, c(2:4)] <- sapply(BigTwoBF_cmin_conthigh[, c(2:4)], as.numeric) BigTwoBF_cmin_contveryhigh[, c(2:4)] <- sapply(BigTwoBF_cmin_contveryhigh[, c(2:4)], as.numeric) BFIBigTwo <- c("TS/least religious states","AC/least religious states","BF/least religious states","TS/most religious states","AC/most religious states","BF/most religious states", "TS/least religious states","AC/least religious states","BF/least religious states","TS/most religious states","AC/most religious states","BF/most religious states", "TS/least religious states","AC/least religious states","BF/least religious states","TS/most religious states","AC/most religious states","BF/most religious states", "TS/least religious states","AC/least religious states","BF/least religious states","TS/most religious states","AC/most religious states","BF/most religious states", "TS/least religious states","AC/least religious states","BF/least religious states","TS/most religious states","AC/most religious states","BF/most religious states") ass_cont <- c("strong\ncontrasters","strong\ncontrasters","strong\ncontrasters","strong\ncontrasters","strong\ncontrasters","strong\ncontrasters", "contrasters","contrasters","contrasters","contrasters","contrasters","contrasters", "average","average","average","average","average","average", "assimilators","assimilators","assimilators","assimilators","assimilators","assimilators", "strong\nassimilators","strong\nassimilators","strong\nassimilators","strong\nassimilators","strong\nassimilators","strong\nassimilators") zPE <- c(BigTwoTS_cmin_contveryhigh$zPe[2],BigTwoAC_cmin_contveryhigh$zPe[2],BigTwoBF_cmin_contveryhigh$zPe[2],BigTwoTS_cmax_contveryhigh$zPe[2],BigTwoAC_cmax_contveryhigh$zPe[2],BigTwoBF_cmax_contveryhigh$zPe[2], BigTwoTS_cmin_conthigh$zPe[2],BigTwoAC_cmin_conthigh$zPe[2],BigTwoBF_cmin_conthigh$zPe[2],BigTwoTS_cmax_conthigh$zPe[2],BigTwoAC_cmax_conthigh$zPe[2],BigTwoBF_cmax_conthigh$zPe[2], BigTwoTS_cmin_average$zPe[2],BigTwoAC_cmin_average$zPe[2],BigTwoBF_cmin_average$zPe[2],BigTwoTS_cmax_average$zPe[2],BigTwoAC_cmax_average$zPe[2],BigTwoBF_cmax_average$zPe[2], BigTwoTS_cmin_asshigh$zPe[2],BigTwoAC_cmin_asshigh$zPe[2],BigTwoBF_cmin_asshigh$zPe[2],BigTwoTS_cmax_asshigh$zPe[2],BigTwoAC_cmax_asshigh$zPe[2],BigTwoBF_cmax_asshigh$zPe[2], BigTwoTS_cmin_assveryhigh$zPe[2],BigTwoAC_cmin_assveryhigh$zPe[2],BigTwoBF_cmin_assveryhigh$zPe[2],BigTwoTS_cmax_assveryhigh$zPe[2],BigTwoAC_cmax_assveryhigh$zPe[2],BigTwoBF_cmax_assveryhigh$zPe[2]) lowerCI <- c(BigTwoTS_cmin_contveryhigh$lowerCI[2],BigTwoAC_cmin_contveryhigh$lowerCI[2],BigTwoBF_cmin_contveryhigh$lowerCI[2],BigTwoTS_cmax_contveryhigh$lowerCI[2],BigTwoAC_cmax_contveryhigh$lowerCI[2],BigTwoBF_cmax_contveryhigh$lowerCI[2], BigTwoTS_cmin_conthigh$lowerCI[2],BigTwoAC_cmin_conthigh$lowerCI[2],BigTwoBF_cmin_conthigh$lowerCI[2],BigTwoTS_cmax_conthigh$lowerCI[2],BigTwoAC_cmax_conthigh$lowerCI[2],BigTwoBF_cmax_conthigh$lowerCI[2], BigTwoTS_cmin_average$lowerCI[2],BigTwoAC_cmin_average$lowerCI[2],BigTwoBF_cmin_average$lowerCI[2],BigTwoTS_cmax_average$lowerCI[2],BigTwoAC_cmax_average$lowerCI[2],BigTwoBF_cmax_average$lowerCI[2], BigTwoTS_cmin_asshigh$lowerCI[2],BigTwoAC_cmin_asshigh$lowerCI[2],BigTwoBF_cmin_asshigh$lowerCI[2],BigTwoTS_cmax_asshigh$lowerCI[2],BigTwoAC_cmax_asshigh$lowerCI[2],BigTwoBF_cmax_asshigh$lowerCI[2], BigTwoTS_cmin_assveryhigh$lowerCI[2],BigTwoAC_cmin_assveryhigh$lowerCI[2],BigTwoBF_cmin_assveryhigh$lowerCI[2],BigTwoTS_cmax_assveryhigh$lowerCI[2],BigTwoAC_cmax_assveryhigh$lowerCI[2],BigTwoBF_cmax_assveryhigh$lowerCI[2]) upperCI <- c(BigTwoTS_cmin_contveryhigh$upperCI[2],BigTwoAC_cmin_contveryhigh$upperCI[2],BigTwoBF_cmin_contveryhigh$upperCI[2],BigTwoTS_cmax_contveryhigh$upperCI[2],BigTwoAC_cmax_contveryhigh$upperCI[2],BigTwoBF_cmax_contveryhigh$upperCI[2], BigTwoTS_cmin_conthigh$upperCI[2],BigTwoAC_cmin_conthigh$upperCI[2],BigTwoBF_cmin_conthigh$upperCI[2],BigTwoTS_cmax_conthigh$upperCI[2],BigTwoAC_cmax_conthigh$upperCI[2],BigTwoBF_cmax_conthigh$upperCI[2], BigTwoTS_cmin_average$upperCI[2],BigTwoAC_cmin_average$upperCI[2],BigTwoBF_cmin_average$upperCI[2],BigTwoTS_cmax_average$upperCI[2],BigTwoAC_cmax_average$upperCI[2],BigTwoBF_cmax_average$upperCI[2], BigTwoTS_cmin_asshigh$upperCI[2],BigTwoAC_cmin_asshigh$upperCI[2],BigTwoBF_cmin_asshigh$upperCI[2],BigTwoTS_cmax_asshigh$upperCI[2],BigTwoAC_cmax_asshigh$upperCI[2],BigTwoBF_cmax_asshigh$upperCI[2], BigTwoTS_cmin_assveryhigh$upperCI[2],BigTwoAC_cmin_assveryhigh$upperCI[2],BigTwoBF_cmin_assveryhigh$upperCI[2],BigTwoTS_cmax_assveryhigh$upperCI[2],BigTwoAC_cmax_assveryhigh$upperCI[2],BigTwoBF_cmax_assveryhigh$upperCI[2]) dat1 <- data.frame(BFIBigTwo,ass_cont,zPE,lowerCI,upperCI) rm(BFIBigTwo,ass_cont,zPE,lowerCI,upperCI) dat1$BFIBigTwo <- factor(dat1$BFIBigTwo, levels=unique(as.character(dat1$BFIBigTwo))) dat1$ass_cont <- factor(dat1$ass_cont, levels=unique(as.character(dat1$ass_cont))) BigTwoPlot <- ggplot(data=dat1, aes(fill=BFIBigTwo, x=ass_cont, y=zPE, group=BFIBigTwo, color=BFIBigTwo)) + geom_bar(stat="identity", position=position_dodge(), width=0.9, color="black") + geom_errorbar(aes(ymin=lowerCI, ymax=upperCI), position=position_dodge(width=0.9), width=0.3, size=0.5, color="black") + ylim(-0.075, 0.33) + scale_fill_manual(values=c("#f4a460", "#cd853f", "#8b4513", "#AAEEFF", "#2E9AFE", "#2E64FE")) + theme(axis.title.x = element_blank(), legend.position=c(0.17, 0.8), legend.title=element_blank(), axis.text = element_text(color="black", size=12), legend.text = element_text(size=12), axis.title.y = element_text(size=14), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(colour = "black")) + geom_hline(yintercept=0) + ggtitle("Figure SOM-U2a. Big Two") + ylab("religiosity-esteem association (zPE)") BigFive_cmax_assveryhigh[, c(2:4)] <- sapply(BigFive_cmax_assveryhigh[, c(2:4)], as.numeric) BigFive_cmax_asshigh[, c(2:4)] <- sapply(BigFive_cmax_asshigh[, c(2:4)], as.numeric) BigFive_cmax_average[, c(2:4)] <- sapply(BigFive_cmax_average[, c(2:4)], as.numeric) BigFive_cmax_conthigh[, c(2:4)] <- sapply(BigFive_cmax_conthigh[, c(2:4)], as.numeric) BigFive_cmax_contveryhigh[, c(2:4)] <- sapply(BigFive_cmax_contveryhigh[, c(2:4)], as.numeric) BigFive_cmin_assveryhigh[, c(2:4)] <- sapply(BigFive_cmin_assveryhigh[, c(2:4)], as.numeric) BigFive_cmin_asshigh[, c(2:4)] <- sapply(BigFive_cmin_asshigh[, c(2:4)], as.numeric) BigFive_cmin_average[, c(2:4)] <- sapply(BigFive_cmin_average[, c(2:4)], as.numeric) BigFive_cmin_conthigh[, c(2:4)] <- sapply(BigFive_cmin_conthigh[, c(2:4)], as.numeric) BigFive_cmin_contveryhigh[, c(2:4)] <- sapply(BigFive_cmin_contveryhigh[, c(2:4)], as.numeric) sta <- c("least religious states","most religious states","least religious states","most religious states","least religious states","most religious states", "least religious states","most religious states","least religious states","most religious states") ass_cont <- c("strong\ncontrasters","strong\ncontrasters","contrasters","contrasters","average","average","assimilators","assimilators","strong\nassimilators","strong\nassimilators") zPE <- c(BigFive_cmin_contveryhigh$zPe[2],BigFive_cmax_contveryhigh$zPe[2],BigFive_cmin_conthigh$zPe[2],BigFive_cmax_conthigh$zPe[2],BigFive_cmin_average$zPe[2],BigFive_cmax_average$zPe[2], BigFive_cmin_asshigh$zPe[2],BigFive_cmax_asshigh$zPe[2],BigFive_cmin_assveryhigh$zPe[2],BigFive_cmax_assveryhigh$zPe[2]) lowerCI <- c(BigFive_cmin_contveryhigh$lowerCI[2],BigFive_cmax_contveryhigh$lowerCI[2],BigFive_cmin_conthigh$lowerCI[2],BigFive_cmax_conthigh$lowerCI[2],BigFive_cmin_average$lowerCI[2],BigFive_cmax_average$lowerCI[2], BigFive_cmin_asshigh$lowerCI[2],BigFive_cmax_asshigh$lowerCI[2],BigFive_cmin_assveryhigh$lowerCI[2],BigFive_cmax_assveryhigh$lowerCI[2]) upperCI <- c(BigFive_cmin_contveryhigh$upperCI[2],BigFive_cmax_contveryhigh$upperCI[2],BigFive_cmin_conthigh$upperCI[2],BigFive_cmax_conthigh$upperCI[2],BigFive_cmin_average$upperCI[2],BigFive_cmax_average$upperCI[2], BigFive_cmin_asshigh$upperCI[2],BigFive_cmax_asshigh$upperCI[2],BigFive_cmin_assveryhigh$upperCI[2],BigFive_cmax_assveryhigh$upperCI[2]) dat1 <- data.frame(sta,ass_cont,zPE,lowerCI,upperCI) rm(sta,ass_cont,zPE,lowerCI,upperCI) dat1$sta <- factor(dat1$sta, levels=unique(as.character(dat1$sta))) dat1$ass_cont <- factor(dat1$ass_cont, levels=unique(as.character(dat1$ass_cont))) BigFivePlot <- ggplot(data=dat1, aes(fill=sta, x=ass_cont, y=zPE, color=sta, group=sta)) + geom_bar(stat="identity", position=position_dodge(), width=0.9, color="black") + geom_errorbar(aes(ymin=lowerCI, ymax=upperCI), position=position_dodge(width=0.9), width=0.2, size=0.5, color="black") + ylim(-0.075, 0.33) + scale_fill_manual(values=c("#cd853f", "#2E9AFE")) + theme(axis.title.x = element_blank(), legend.position=c(0.15, 0.8), legend.title=element_blank(), axis.text = element_text(color="black", size=12), legend.text = element_text(size=12), axis.title.y = element_text(size=14), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(colour = "black")) + geom_hline(yintercept=0) + ggtitle("Figure SOM-U2b. Big Five") + ylab("religiosity-esteem association (zPE)") figure <- grid.arrange(BigTwoPlot,BigFivePlot) ggsave("Figure_SOM-U2.jpg", plot=figure, device="jpg", width=210, height=210, units="mm", dpi=720) ggsave("Figure_SOM-U2.pdf", plot=figure, device="pdf", width=210, height=210, units="mm", dpi=720) rm(dat1,BigTwoPlot,BigFivePlot,figure) #size of the person-culture match effect in each model (power of culture) delta_zPE_b2ts <- BigTwoTS_cmax_average[2,2] - BigTwoTS_cmin_average[2,2] lowerCI_b2ts <- BigTwoTS_cmax_average[2,3] - BigTwoTS_cmin_average[2,4] upperCI_b2ts <- BigTwoTS_cmax_average[2,4] - BigTwoTS_cmin_average[2,3] delta_zPE_b2ac <- BigTwoAC_cmax_average[2,2] - BigTwoAC_cmin_average[2,2] lowerCI_b2ac <- BigTwoAC_cmax_average[2,3] - BigTwoAC_cmin_average[2,4] upperCI_b2ac <- BigTwoAC_cmax_average[2,4] - BigTwoAC_cmin_average[2,3] delta_zPE_b2bf <- BigTwoBF_cmax_average[2,2] - BigTwoBF_cmin_average[2,2] lowerCI_b2bf <- BigTwoBF_cmax_average[2,3] - BigTwoBF_cmin_average[2,4] upperCI_b2bf <- BigTwoBF_cmax_average[2,4] - BigTwoBF_cmin_average[2,3] delta_zPE_b5 <- BigFive_cmax_average[2,2] - BigFive_cmin_average[2,2] lowerCI_b5 <- BigFive_cmax_average[2,3] - BigFive_cmin_average[2,4] upperCI_b5 <- BigFive_cmax_average[2,4] - BigFive_cmin_average[2,3] paste0("Big Two TS Model: delta_zPE = ", round(delta_zPE_b2ts, digits=2), " 95% CI [", round(lowerCI_b2ts, digits=2), ", ", round(upperCI_b2ts, digits=2), "]") paste0("Big Two AC Model: delta_zPE = ", round(delta_zPE_b2ac, digits=2), " 95% CI [", round(lowerCI_b2ac, digits=2), ", ", round(upperCI_b2ac, digits=2), "]") paste0("Big Two BF Model: delta_zPE = ", round(delta_zPE_b2bf, digits=2), " 95% CI [", round(lowerCI_b2bf, digits=2), ", ", round(upperCI_b2bf, digits=2), "]") paste0("Big Five Model: delta_zPE = ", round(delta_zPE_b5, digits=2), " 95% CI [", round(lowerCI_b5, digits=2), ", ", round(upperCI_b5, digits=2), "]") rm(delta_zPE_b2ts,lowerCI_b2ts,upperCI_b2ts,delta_zPE_b2ac,lowerCI_b2ac,upperCI_b2ac,delta_zPE_b2bf,lowerCI_b2bf,upperCI_b2bf,delta_zPE_b5,lowerCI_b5,upperCI_b5) #extent to which personality alters the person-culture match effect in each model (power of personality) ##BIG TWO TS delta_zPE_ass <- BigTwoTS_cmax_assveryhigh$zPe[2] - BigTwoTS_cmin_assveryhigh$zPe[2] lowerCI_ass <- BigTwoTS_cmax_assveryhigh$lowerCI[2] - BigTwoTS_cmin_assveryhigh$upperCI[2] upperCI_ass <- BigTwoTS_cmax_assveryhigh$upperCI[2] - BigTwoTS_cmin_assveryhigh$lowerCI[2] delta_zPE_cont <- BigTwoTS_cmax_contveryhigh$zPe[2] - BigTwoTS_cmin_contveryhigh$zPe[2] lowerCI_cont <- BigTwoTS_cmax_contveryhigh$lowerCI[2] - BigTwoTS_cmin_contveryhigh$upperCI[2] upperCI_cont <- BigTwoTS_cmax_contveryhigh$upperCI[2] - BigTwoTS_cmin_contveryhigh$lowerCI[2] paste0("Big Two TS Model: strong assimilators: delta_zPE = ", round(delta_zPE_ass, digits=2), " 95% CI [", round(lowerCI_ass, digits=2), ", ", round(upperCI_ass, digits=2), "]") paste0("Big Two TS Model: strong contrasters: delta_zPE = ", round(delta_zPE_cont, digits=2), " 95% CI [", round(lowerCI_cont, digits=2), ", ", round(upperCI_cont, digits=2), "]") power_personality <- round(delta_zPE_ass, digits=2) - round(delta_zPE_cont, digits=2) power_lowerCI <- round(lowerCI_ass, digits=2) - round(upperCI_cont, digits=2) power_upperCI <- round(upperCI_ass, digits=2) - round(lowerCI_cont, digits=2) paste0("Power of Big Two TS: delta_zPE = ", power_personality, " 95% CI [", power_lowerCI, ", ", power_upperCI, "]") ##BIG TWO AC delta_zPE_ass <- BigTwoAC_cmax_assveryhigh$zPe[2] - BigTwoAC_cmin_assveryhigh$zPe[2] lowerCI_ass <- BigTwoAC_cmax_assveryhigh$lowerCI[2] - BigTwoAC_cmin_assveryhigh$upperCI[2] upperCI_ass <- BigTwoAC_cmax_assveryhigh$upperCI[2] - BigTwoAC_cmin_assveryhigh$lowerCI[2] delta_zPE_cont <- BigTwoAC_cmax_contveryhigh$zPe[2] - BigTwoAC_cmin_contveryhigh$zPe[2] lowerCI_cont <- BigTwoAC_cmax_contveryhigh$lowerCI[2] - BigTwoAC_cmin_contveryhigh$upperCI[2] upperCI_cont <- BigTwoAC_cmax_contveryhigh$upperCI[2] - BigTwoAC_cmin_contveryhigh$lowerCI[2] paste0("Big Two AC Model: strong assimilators: delta_zPE = ", round(delta_zPE_ass, digits=2), " 95% CI [", round(lowerCI_ass, digits=2), ", ", round(upperCI_ass, digits=2), "]") paste0("Big Two AC Model: strong contrasters: delta_zPE = ", round(delta_zPE_cont, digits=6), " 95% CI [", round(lowerCI_cont, digits=2), ", ", round(upperCI_cont, digits=2), "]") power_personality <- round(delta_zPE_ass, digits=2) - round(delta_zPE_cont, digits=2) power_lowerCI <- round(lowerCI_ass, digits=2) - round(upperCI_cont, digits=2) power_upperCI <- round(upperCI_ass, digits=2) - round(lowerCI_cont, digits=2) paste0("Power of Big Two AC: delta_zPE = ", power_personality, " 95% CI [", power_lowerCI, ", ", power_upperCI, "]") ##BIG TWO BF delta_zPE_ass <- BigTwoBF_cmax_assveryhigh$zPe[2] - BigTwoBF_cmin_assveryhigh$zPe[2] lowerCI_ass <- BigTwoBF_cmax_assveryhigh$lowerCI[2] - BigTwoBF_cmin_assveryhigh$upperCI[2] upperCI_ass <- BigTwoBF_cmax_assveryhigh$upperCI[2] - BigTwoBF_cmin_assveryhigh$lowerCI[2] delta_zPE_cont <- BigTwoBF_cmax_contveryhigh$zPe[2] - BigTwoBF_cmin_contveryhigh$zPe[2] lowerCI_cont <- BigTwoBF_cmax_contveryhigh$lowerCI[2] - BigTwoBF_cmin_contveryhigh$upperCI[2] upperCI_cont <- BigTwoBF_cmax_contveryhigh$upperCI[2] - BigTwoBF_cmin_contveryhigh$lowerCI[2] paste0("Big Two BF Model: strong assimilators: delta_zPE = ", round(delta_zPE_ass, digits=2), " 95% CI [", round(lowerCI_ass, digits=2), ", ", round(upperCI_ass, digits=2), "]") paste0("Big Two BF Model: strong contrasters: delta_zPE = ", round(delta_zPE_cont, digits=2), " 95% CI [", round(lowerCI_cont, digits=2), ", ", round(upperCI_cont, digits=2), "]") power_personality <- round(delta_zPE_ass, digits=2) - round(delta_zPE_cont, digits=2) power_lowerCI <- round(lowerCI_ass, digits=2) - round(upperCI_cont, digits=2) power_upperCI <- round(upperCI_ass, digits=2) - round(lowerCI_cont, digits=2) paste0("Power of Big Two BF: delta_zPE = ", power_personality, " 95% CI [", power_lowerCI, ", ", power_upperCI, "]") ##BIG FIVE delta_zPE_ass <- BigFive_cmax_assveryhigh$zPe[2] - BigFive_cmin_assveryhigh$zPe[2] lowerCI_ass <- BigFive_cmax_assveryhigh$lowerCI[2] - BigFive_cmin_assveryhigh$upperCI[2] upperCI_ass <- BigFive_cmax_assveryhigh$upperCI[2] - BigFive_cmin_assveryhigh$lowerCI[2] delta_zPE_cont <- BigFive_cmax_contveryhigh$zPe[2] - BigFive_cmin_contveryhigh$zPe[2] lowerCI_cont <- BigFive_cmax_contveryhigh$lowerCI[2] - BigFive_cmin_contveryhigh$upperCI[2] upperCI_cont <- BigFive_cmax_contveryhigh$upperCI[2] - BigFive_cmin_contveryhigh$lowerCI[2] paste0("Big Five Model: strong assimilators: delta_zPE = ", round(delta_zPE_ass, digits=2), " 95% CI [", round(lowerCI_ass, digits=2), ", ", round(upperCI_ass, digits=2), "]") paste0("Big Five Model: strong contrasters: delta_zPE = ", round(delta_zPE_cont, digits=2), " 95% CI [", round(lowerCI_cont, digits=2), ", ", round(upperCI_cont, digits=2), "]") power_personality <- round(delta_zPE_ass, digits=2) - round(delta_zPE_cont, digits=2) power_lowerCI <- round(lowerCI_ass, digits=2) - round(upperCI_cont, digits=2) power_upperCI <- round(upperCI_ass, digits=2) - round(lowerCI_cont, digits=2) paste0("Power of Big Five: delta_zPE = ", power_personality, " 95% CI [", power_lowerCI, ", ", power_upperCI, "]")