#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","tibble","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, cnt_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_er <- rowMeans(dat[c("ext01","ext03","ext04","ext06","ext08","opn01","opn02","opn05")], na.rm=TRUE) #dat$com_er <- rowMeans(dat[c("agr02","agr04","agr05","agr07","agr09","cns01","cns03","cns06","cns08")], 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) # ##N = 9,328,610 # ##select data #dat[is.na(dat)] <- NA #dat <- subset(dat, cnt_now != -999 & 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, cnt_frq = ave(seq(nrow(dat)), cnt_now, FUN=length)) #dat <- subset(dat, cnt_frq >= 300) # ##N = 2,672,820 # ##NOTE: there are no missings at the agency-communion level (checked as follows) #sum(is.na(dat$agy_er)) #sum(is.na(dat$com_er)) #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)) # ##psychometrics #dat2 <- aggregate(record_id ~ cnt_now, dat, mean) #dat2$cnt <- 1:nrow(dat2) #dat2$record_id <- NULL #dat <- merge(dat, dat2, by = "cnt_now", all = TRUE) #rm(dat2) # #cnt <- c(1:max(dat$cnt)) #agy_er_alpha <- c(1:max(dat$cnt)) #com_er_alpha <- c(1:max(dat$cnt)) #agy_ts_alpha <- c(1:max(dat$cnt)) #com_ts_alpha <- c(1:max(dat$cnt)) #agy_ac_alpha <- c(1:max(dat$cnt)) #com_ac_alpha <- c(1:max(dat$cnt)) #agy_bf_alpha <- c(1:max(dat$cnt)) #com_bf_alpha <- c(1:max(dat$cnt)) #agr_alpha <- c(1:max(dat$cnt)) #cns_alpha <- c(1:max(dat$cnt)) #ext_alpha <- c(1:max(dat$cnt)) #opn_alpha <- c(1:max(dat$cnt)) #neu_alpha <- c(1:max(dat$cnt)) # #dat2 <- data.frame(cnt,agy_er_alpha,com_er_alpha,agy_ts_alpha,com_ts_alpha,agy_ac_alpha,com_ac_alpha, # agy_bf_alpha,com_bf_alpha,agr_alpha,cns_alpha,ext_alpha,opn_alpha,neu_alpha) # #rm(cnt,agy_er_alpha,com_er_alpha,agy_ts_alpha,com_ts_alpha,agy_ac_alpha,com_ac_alpha, # agy_bf_alpha,com_bf_alpha,agr_alpha,cns_alpha,ext_alpha,opn_alpha,neu_alpha) # #i <- 1 #while (i <= max(dat$cnt)) { # dat3 <- subset(dat, cnt == i) # dat3$id <- seq.int(nrow(dat3)) # dat2$cnt[i] <- i # dat.agy_er <- subset(dat3, select = c(ext01,ext03,ext04,ext06,ext08,opn01,opn02,opn05)) # dat.com_er <- subset(dat3, select = c(agr02,agr04,agr05,agr07,agr09,cns01,cns03,cns06,cns08)) # dat.agy_ts <- subset(dat3, select = c(ext01,ext03,ext04,ext05r,ext06,ext07r,ext08,opn01)) # dat.com_ts <- subset(dat3, select = c(agr01r,agr02,agr03r,agr04,agr07,agr08r,agr09,cns07)) # dat.agy_ac <- subset(dat3, select = c(ext03,ext06,ext07r,opn04,opn07r,neu07)) # dat.com_ac <- subset(dat3, select = c(agr01r,agr03r,agr05,agr07,agr08r,agr09,cns01)) # dat.agy_bf <- subset(dat3, select = c(cns06,ext02r,ext05r,opn01,opn03,opn05,opn10,neu02)) # dat.com_bf <- subset(dat3, select = c(agr02,agr03r,agr04,agr06r,agr07,agr08r,agr09,cns01)) # dat.agr <- subset(dat3, select = c(agr01r,agr02,agr03r,agr04,agr05,agr06r,agr07,agr08r,agr09)) # dat.cns <- subset(dat3, select = c(cns01,cns02r,cns03,cns04r,cns05r,cns06,cns07,cns08,cns09r)) # dat.ext <- subset(dat3, select = c(ext01,ext02r,ext03,ext04,ext05r,ext06,ext07r,ext08)) # dat.opn <- subset(dat3, select = c(opn01,opn02,opn03,opn04,opn05,opn06,opn07r,opn08,opn09r,opn10)) # dat.neu <- subset(dat3, select = c(neu01,neu02r,neu03,neu04,neu05r,neu06,neu07r,neu08)) # dat2$agy_er_alpha[i] <- round(psych::alpha(dat.agy_er, na.rm=T)$total[1], digits = 2) # dat2$com_er_alpha[i] <- round(psych::alpha(dat.com_er, na.rm=T)$total[1], digits = 2) # dat2$agy_ts_alpha[i] <- round(psych::alpha(dat.agy_ts, na.rm=T)$total[1], digits = 2) # dat2$com_ts_alpha[i] <- round(psych::alpha(dat.com_ts, na.rm=T)$total[1], digits = 2) # dat2$agy_ac_alpha[i] <- round(psych::alpha(dat.agy_ac, na.rm=T)$total[1], digits = 2) # dat2$com_ac_alpha[i] <- round(psych::alpha(dat.com_ac, na.rm=T)$total[1], digits = 2) # dat2$agy_bf_alpha[i] <- round(psych::alpha(dat.agy_bf, na.rm=T)$total[1], digits = 2) # dat2$com_bf_alpha[i] <- round(psych::alpha(dat.com_bf, na.rm=T)$total[1], digits = 2) # dat2$agr_alpha[i] <- round(psych::alpha(dat.agr, na.rm=T)$total[1], digits = 2) # dat2$cns_alpha[i] <- round(psych::alpha(dat.cns, na.rm=T)$total[1], digits = 2) # dat2$ext_alpha[i] <- round(psych::alpha(dat.ext, na.rm=T)$total[1], digits = 2) # dat2$opn_alpha[i] <- round(psych::alpha(dat.opn, na.rm=T)$total[1], digits = 2) # dat2$neu_alpha[i] <- round(psych::alpha(dat.neu, na.rm=T)$total[1], digits = 2) # i <- i + 1 #} # #rm(dat3,i,dat.agy_er,dat.com_er,dat.agy_ts,dat.com_ts,dat.agy_ac,dat.com_ac,dat.agy_bf,dat.com_bf, # dat.agr,dat.cns,dat.ext,dat.opn,dat.neu) # ##save data #fwrite(dat2, file = "PCmatch_psychmet_MainText.csv") #rm(dat2) # ##measurement invariance #dat2 <- setNames(data.frame(matrix(ncol = 16, nrow = 6)), c("model","fit_index","agy_er","com_er","agy_ts","com_ts","agy_ac","com_ac", # "agy_bf","com_bf","agr","cns","ext","opn","neu","opn_4pcl")) #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_er <- 'agy =~ ext01 + ext03 + ext04 + ext06 + ext08 + opn01 + opn02 + opn05' #config_agy_er <- cfa(mod_agy_er, data = dat, missing = "ML", group = "cnt_now") #dat2$agy_er[1] <- fitMeasures(config_agy_er, "cfi") #dat2$agy_er[2] <- fitMeasures(config_agy_er, "rmsea") #dat2$agy_er[3] <- fitMeasures(config_agy_er, "srmr") #metric_agy_er <- cfa(mod_agy_er, data = dat, missing = "ML", group = "cnt_now", group.equal = c("loadings")) #dat2$agy_er[4] <- fitMeasures(metric_agy_er, "cfi") #dat2$agy_er[5] <- fitMeasures(metric_agy_er, "rmsea") #dat2$agy_er[6] <- fitMeasures(metric_agy_er, "srmr") #rm(mod_agy_er,config_agy_er,metric_agy_er) # #mod_com_er <- 'com =~ agr02 + agr04 + agr05 + agr07 + agr09 + cns01 + cns03 + cns06 + cns08' #config_com_er <- cfa(mod_com_er, data = dat, missing = "ML", group = "cnt_now") #dat2$com_er[1] <- fitMeasures(config_com_er, "cfi") #dat2$com_er[2] <- fitMeasures(config_com_er, "rmsea") #dat2$com_er[3] <- fitMeasures(config_com_er, "srmr") #metric_com_er <- cfa(mod_com_er, data = dat, missing = "ML", group = "cnt_now", group.equal = c("loadings")) #dat2$com_er[4] <- fitMeasures(metric_com_er, "cfi") #dat2$com_er[5] <- fitMeasures(metric_com_er, "rmsea") #dat2$com_er[6] <- fitMeasures(metric_com_er, "srmr") #rm(mod_com_er,config_com_er,metric_com_er) # #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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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 = "cnt_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) # ##parceling openness (following Little et al.) #dat.opn <- subset(dat, select = c(opn01,opn02,opn03,opn04,opn05,opn06,opn07r,opn08,opn09r,opn10)) #itot.opn <- item.total(dat.opn) #itot.opn <- itot.opn[order(-itot.opn$Item.Total),] #dat$opn_4pcl1 <- rowMeans(dat[c(as.character(itot.opn$Variable[1]),as.character(itot.opn$Variable[8]),as.character(itot.opn$Variable[9]))], na.rm=TRUE) #dat$opn_4pcl2 <- rowMeans(dat[c(as.character(itot.opn$Variable[2]),as.character(itot.opn$Variable[7]),as.character(itot.opn$Variable[10]))], na.rm=TRUE) #dat$opn_4pcl3 <- rowMeans(dat[c(as.character(itot.opn$Variable[3]),as.character(itot.opn$Variable[6]))], na.rm=TRUE) #dat$opn_4pcl4 <- rowMeans(dat[c(as.character(itot.opn$Variable[4]),as.character(itot.opn$Variable[5]))], na.rm=TRUE) #rm(itot.opn,dat.opn) # #mod_opn_4pcl <- 'opn.pcl =~ opn_4pcl1 + opn_4pcl2 + opn_4pcl3 + opn_4pcl4' #config_opn_4pcl <- cfa(mod_opn_4pcl, data = dat, missing = "ML", group = "cnt_now") #dat2$opn_4pcl[1] <- fitMeasures(config_opn_4pcl, "cfi") #dat2$opn_4pcl[2] <- fitMeasures(config_opn_4pcl, "rmsea") #dat2$opn_4pcl[3] <- fitMeasures(config_opn_4pcl, "srmr") #metric_opn_4pcl <- cfa(mod_opn_4pcl, data = dat, missing = "ML", group = "cnt_now", group.equal = c("loadings")) #dat2$opn_4pcl[4] <- fitMeasures(metric_opn_4pcl, "cfi") #dat2$opn_4pcl[5] <- fitMeasures(metric_opn_4pcl, "rmsea") #dat2$opn_4pcl[6] <- fitMeasures(metric_opn_4pcl, "srmr") #rm(mod_opn_4pcl,config_opn_4pcl,metric_opn_4pcl) # ##select variables #dat <- subset(dat, select = c(record_id, cnt_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_MainText.csv") #fwrite(dat2, file = "PCmatch_mi_MainText.csv") # #rm(dat2) #####################################Analyses with publicly available dataset##################################### #read data dat <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/330/4/PCmatch_csv-data_MainText.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)) #Table SOM-R2 dat$cntrel <- ave(dat$relig, dat$cnt_now, FUN = function(x) mean(x, na.rm=T)) dat2 <- subset(dat, select = c(cnt_now, age2, female, cntrel)) dat2 <- transform(dat2, cnt_frq = ave(seq(nrow(dat2)), cnt_now, FUN=length)) dat_s2_1 <- aggregate(cnt_frq ~ cnt_now, dat2, FUN = function(x) mean(x, na.rm=T)) dat_s2_2 <- aggregate(age2 ~ cnt_now, dat2, FUN = function(x) round(mean(x, na.rm=T),2)) dat_s2_3 <- aggregate(age2 ~ cnt_now, dat2, FUN = function(x) round(sd(x, na.rm=T),2)) dat_s2_4 <- aggregate(female ~ cnt_now, dat2, FUN = function(x) round(100*mean(x, na.rm=T),2)) dat_s2_5 <- aggregate(cntrel ~ cnt_now, dat2, FUN = function(x) round(mean(x, na.rm=T),2)) colnames(dat_s2_2)[2] <- "age_mean" colnames(dat_s2_3)[2] <- "age_sd" dat_s2_2$cnt_now <- NULL dat_s2_3$cnt_now <- NULL dat_s2_4$cnt_now <- NULL dat_s2_5$cnt_now <- NULL dat_s2 <- cbind(dat_s2_1,dat_s2_2,dat_s2_3,dat_s2_4,dat_s2_5) dat_s2_6 <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/330/1/cnt_labels.csv", header = T, sep = ";")) dat_s2 <- merge(dat_s2, dat_s2_6, by = "cnt_now") dat_s2 <- dat_s2[order(-dat_s2$cntrel),] dat_s2$cnt_now <- NULL dat_s2 <- data.frame(dat_s2[c("cnt_name","cnt_frq","age_mean","age_sd","female","cntrel")]) rm(dat2,dat_s2_1,dat_s2_2,dat_s2_3,dat_s2_4,dat_s2_5,dat_s2_6) col_keys <- c("cnt_name","cnt_frq","age_mean","age_sd","female","cntrel") head1 <- c("country","N","age","age","% women","country-level") head2 <- c("","","M","SD","","religiosity") head <- data.frame(col_keys,head1,head2, stringsAsFactors = FALSE) rm(col_keys,head1,head2) tbl <- flextable(dat_s2) tbl <- set_header_df(tbl, mapping=head, key="col_keys") tbl <- colformat_num(tbl, j=2, digits=0) tbl <- colformat_num(tbl, J=3:6, digits=2) tbl <- hline_top(tbl, j=1:6, border=fp_border(width=2), part="header") tbl <- merge_at(tbl, i=1, j=3:4, part="header") tbl <- hline(tbl, i=1, j=3:4, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=1:6, 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 <- italic(tbl, i=1, j=2, part="header") tbl <- italic(tbl, i=2, j=c(3,4), part="header") tbl <- align(tbl, align="center", part="all") tbl <- align(tbl, j=1, align="left", part="all") tbl <- width(tbl, j=1, width=1.4) tbl <- width(tbl, j=c(2,5,6), width=1) tbl <- width(tbl, j=c(3,4), width=0.5) tbl <- height_all(tbl, height=.1, part="all") tbl doc <- read_docx() doc <- body_add_flextable(doc, value = tbl) print(doc, target = "Table_SOM-R2.docx") rm(dat_s2,head,tbl,doc) #Table SOM-R3 dat1 <- aggregate(record_id ~ cnt_now, dat, mean) dat1$cnt <- 1:nrow(dat1) dat1$record_id <- NULL dat <- merge(dat, dat1, by = "cnt_now", all = TRUE) cnt <- c(1:max(dat$cnt)) m_agy_ts <- c(1:max(dat$cnt)) m_com_ts <- c(1:max(dat$cnt)) m_agy_ac <- c(1:max(dat$cnt)) m_com_ac <- c(1:max(dat$cnt)) m_agy_bf <- c(1:max(dat$cnt)) m_com_bf <- c(1:max(dat$cnt)) m_agr <- c(1:max(dat$cnt)) m_cns <- c(1:max(dat$cnt)) m_ext <- c(1:max(dat$cnt)) m_opn <- c(1:max(dat$cnt)) m_neu <- c(1:max(dat$cnt)) m_se <- c(1:max(dat$cnt)) m_relig <- c(1:max(dat$cnt)) sd_agy_ts <- c(1:max(dat$cnt)) sd_com_ts <- c(1:max(dat$cnt)) sd_agy_ac <- c(1:max(dat$cnt)) sd_com_ac <- c(1:max(dat$cnt)) sd_agy_bf <- c(1:max(dat$cnt)) sd_com_bf <- c(1:max(dat$cnt)) sd_agr <- c(1:max(dat$cnt)) sd_cns <- c(1:max(dat$cnt)) sd_ext <- c(1:max(dat$cnt)) sd_opn <- c(1:max(dat$cnt)) sd_neu <- c(1:max(dat$cnt)) sd_se <- c(1:max(dat$cnt)) sd_relig <- c(1:max(dat$cnt)) dat2 <- data.frame(cnt,m_se,sd_se,m_relig,sd_relig,m_agy_ts,sd_agy_ts,m_com_ts,sd_com_ts,m_agy_ac,sd_agy_ac,m_com_ac,sd_com_ac,m_agy_bf,sd_agy_bf,m_com_bf,sd_com_bf, m_agr,sd_agr,m_cns,sd_cns,m_ext,sd_ext,m_opn,sd_opn,m_neu,sd_neu) rm(cnt,m_se,sd_se,m_relig,sd_relig,m_agy_ts,sd_agy_ts,m_com_ts,sd_com_ts,m_agy_ac,sd_agy_ac,m_com_ac,sd_com_ac,m_agy_bf,sd_agy_bf,m_com_bf,sd_com_bf, m_agr,sd_agr,m_cns,sd_cns,m_ext,sd_ext,m_opn,sd_opn,m_neu,sd_neu) i <- 1 while (i <= max(dat$cnt)) { dat3 <- subset(dat, cnt == i) dat3$id <- seq.int(nrow(dat3)) dat2$cnt[i] <- i j <- 6 k <- 2 l <- 3 while (j <= 18) { dat2[i,k] <- round(mean(dat3[,j], na.rm=T),digits=2) dat2[i,l] <- round(sd(dat3[,j], na.rm=T),digits=2) j <- j + 1 k <- k + 2 l <- l + 2 } i <- i + 1 } dat2 <- merge(dat1, dat2, by = "cnt") dat3 <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/330/1/cnt_labels.csv", header = T, sep = ";")) dat3 <- merge(dat3, dat2, by = "cnt_now") dat3 <- dat3[order(dat3$cnt_name),] dat3 <- data.frame(dat3[c("cnt_name","m_se","sd_se","m_relig","sd_relig","m_agy_ts","sd_agy_ts","m_com_ts","sd_com_ts","m_agy_ac","sd_agy_ac","m_com_ac","sd_com_ac", "m_agy_bf","sd_agy_bf","m_com_bf","sd_com_bf","m_agr","sd_agr","m_cns","sd_cns","m_ext","sd_ext","m_opn","sd_opn","m_neu","sd_neu")]) dat3$se <- paste0(dat3$m_se, " (", dat3$sd_se, ")") dat3$relig <- paste0(dat3$m_relig, " (", dat3$sd_relig, ")") dat3$agy_ts <- paste0(dat3$m_agy_ts, " (", dat3$sd_agy_ts, ")") dat3$com_ts <- paste0(dat3$m_com_ts, " (", dat3$sd_com_ts, ")") dat3$agy_ac <- paste0(dat3$m_agy_ac, " (", dat3$sd_agy_ac, ")") dat3$com_ac <- paste0(dat3$m_com_ac, " (", dat3$sd_com_ac, ")") dat3$agy_bf <- paste0(dat3$m_agy_bf, " (", dat3$sd_agy_bf, ")") dat3$com_bf <- paste0(dat3$m_com_bf, " (", dat3$sd_com_bf, ")") dat3$agr <- paste0(dat3$m_agr, " (", dat3$sd_agr, ")") dat3$cns <- paste0(dat3$m_cns, " (", dat3$sd_cns, ")") dat3$ext <- paste0(dat3$m_ext, " (", dat3$sd_ext, ")") dat3$opn <- paste0(dat3$m_opn, " (", dat3$sd_opn, ")") dat3$neu <- paste0(dat3$m_neu, " (", dat3$sd_neu, ")") dat3$blank <- NA dat3$blank.1 <- NA dat3$blank.2 <- NA dat3$blank.3 <- NA dat3$blank.4 <- NA dat3$blank.5 <- NA dat3 <- dat3[c("cnt_name","blank","se","blank.1","relig","blank.2","agy_ts","com_ts","blank.3","agy_ac","com_ac","blank.4", "agy_bf","com_bf","blank.5","agr","cns","ext","opn","neu")] col_keys <- c("cnt_name","blank","se","blank.1","relig","blank.2","agy_ts","com_ts","blank.3","agy_ac","com_ac","blank.4", "agy_bf","com_bf","blank.5","agr","cns","ext","opn","neu") head1 <- c("Country","","Self-esteem","","Religiosity","","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(dat3) tbl <- set_header_df(tbl, mapping=head, key="col_keys") tbl <- colformat_num(tbl, j=c(3,5,7,8,10,11,13,14,16:20), digits=2) tbl <- colformat_lgl(tbl, j =~ blank + blank.1 + blank.2 + blank.3 + blank.4 + blank.5, na_str="") tbl <- merge_at(tbl, i=1, j=7:14, part="header") tbl <- merge_at(tbl, i=1, j=16:20, part="header") tbl <- merge_at(tbl, i=2, j=7:8, part="header") tbl <- merge_at(tbl, i=2, j=10:11, part="header") tbl <- merge_at(tbl, i=2, j=13:14, part="header") tbl <- hline_top(tbl, j=1:20, border=fp_border(width=2), part="header") tbl <- hline(tbl, i=1, j=7:14, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=1, j=16:20, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=7:8, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=10:11, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=2, j=13:14, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=3, j=1:20, 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=1, align="left", part="all") tbl <- width(tbl, j=1, width=0.8) tbl <- width(tbl, j=c(2,4,6,9,12,15), width=0.01) tbl <- height_all(tbl, height=.1, part="all") tbl doc <- read_docx() doc <- body_add_flextable(doc, value = tbl) print(doc, target = "Table_SOM-R3.docx") rm(dat1,dat2,i,j,k,l,dat3,head,tbl,doc) #Table SOM-R4 dat1 <- fread("https://madata.bib.uni-mannheim.de/330/15/scalelength.csv", header = T, sep = ";") dat2 <- fread("https://madata.bib.uni-mannheim.de/330/14/PCmatch_psychmet_MainText.csv", header = T, sep = ",") dat3 <- fread("https://madata.bib.uni-mannheim.de/330/9/PCmatch_mi_MainText.csv", header = T, sep = ",") dat3 <- dat3[, c(1:12,14,13,15)] dat3[,3:15] <- round(dat3[,3:15], digits=2) dat4 <- dat1 dat4$label1 <- "number of items" dat4$label2 <- NA dat4 <- dat4[, c(14,15,1:13)] dat4[] <- lapply(dat4, as.character) dat4 <- add_row(dat4) dat4 <- add_row(dat4) dat4[3,1] <- paste("Cronbach's", intToUtf8(945)) dat4[3,2] <- "(M/SD)" dat4[3,3] <- paste0(f_num(mean(dat2$agy_er_alpha), digits=2),"/",f_num(sd(dat2$agy_er_alpha), digits=2)) dat4[3,4] <- paste0(f_num(mean(dat2$com_er_alpha), digits=2),"/",f_num(sd(dat2$com_er_alpha), digits=2)) dat4[3,5] <- paste0(f_num(mean(dat2$agy_ts_alpha), digits=2),"/",f_num(sd(dat2$agy_ts_alpha), digits=2)) dat4[3,6] <- paste0(f_num(mean(dat2$com_ts_alpha), digits=2),"/",f_num(sd(dat2$com_ts_alpha), digits=2)) dat4[3,7] <- paste0(f_num(mean(dat2$agy_ac_alpha), digits=2),"/",f_num(sd(dat2$agy_ac_alpha), digits=2)) dat4[3,8] <- paste0(f_num(mean(dat2$com_ac_alpha), digits=2),"/",f_num(sd(dat2$com_ac_alpha), digits=2)) dat4[3,9] <- paste0(f_num(mean(dat2$agy_bf_alpha), digits=2),"/",f_num(sd(dat2$agy_bf_alpha), digits=2)) dat4[3,10] <- paste0(f_num(mean(dat2$com_bf_alpha), digits=2),"/",f_num(sd(dat2$com_bf_alpha), digits=2)) dat4[3,11] <- paste0(f_num(mean(dat2$agr_alpha), digits=2),"/",f_num(sd(dat2$agr_alpha), digits=2)) dat4[3,12] <- paste0(f_num(mean(dat2$cns_alpha), digits=2),"/",f_num(sd(dat2$cns_alpha), digits=2)) dat4[3,13] <- paste0(f_num(mean(dat2$opn_alpha), digits=2),"/",f_num(sd(dat2$opn_alpha), digits=2)) dat4[3,14] <- paste0(f_num(mean(dat2$ext_alpha), digits=2),"/",f_num(sd(dat2$ext_alpha), digits=2)) dat4[3,15] <- paste0(f_num(mean(dat2$neu_alpha), digits=2),"/",f_num(sd(dat2$neu_alpha), digits=2)) dat4 <- add_row(dat4) dat4 <- add_row(dat4) dat4[5,1] <- "MI" dat4[5,2] <- "model" dat4 <- add_row(dat4) dat4 <- add_row(dat4) dat4[7,1] <- "CFI" dat4[7,2] <- "config" dat4[7,3:15] <- mutate_all(data.frame(dat3[1,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4[8,2] <- "metric" dat4[8,3:15] <- mutate_all(data.frame(dat3[4,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4[9,2] <- intToUtf8(916) dat4[9,3:15] <- mutate_all(data.frame(dat3[1,3:15] - dat3[4,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4 <- add_row(dat4) dat4[11,1] <- "RMSEA" dat4[11,2] <- "config" dat4[11,3:15] <- mutate_all(data.frame(dat3[2,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4[12,2] <- "metric" dat4[12,3:15] <- mutate_all(data.frame(dat3[5,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4[13,2] <- intToUtf8(916) dat4[13,3:15] <- mutate_all(data.frame(dat3[2,3:15] - dat3[5,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4 <- add_row(dat4) dat4[15,1] <- "SRMR" dat4[15,2] <- "config" dat4[15,3:15] <- mutate_all(data.frame(dat3[3,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4[16,2] <- "metric" dat4[16,3:15] <- mutate_all(data.frame(dat3[6,3:15]), list(~f_num(unlist(.), digits=2))) dat4 <- add_row(dat4) dat4[17,2] <- intToUtf8(916) dat4[17,3:15] <- mutate_all(data.frame(dat3[3,3:15] - dat3[6,3:15]), list(~f_num(unlist(.), digits=2))) dat4$blank1 <- NA dat4$blank2 <- NA dat4$blank3 <- NA dat4$blank4 <- NA dat4 <- dat4[, c(1:4,16,5:6,17,7:8,18,9:10,19,11:15)] col_keys <- c("label1","label2","agy_er","com_er","blank1","agy_ts","com_ts","blank2","agy_ac","com_ac","blank3","agy_bf","com_bf","blank4","agr","cns","opn","ext","neu") head1 <- c("","","Big Two","Big Two","Big Two","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("","","ER","ER","","TS","TS","","AC","AC","","BF","BF","","","","","","") head3 <- c("","","agy","com","","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(dat4) tbl <- set_header_df(tbl, mapping=head, key="col_keys") tbl <- colformat_lgl(tbl, j =~ blank1 + blank2 + blank3 + blank4, na_str="") tbl <- hline_top(tbl, j=1:19, border=fp_border(width=2), part="header") tbl <- merge_at(tbl, i=1, j=3:13, part="header") tbl <- merge_at(tbl, i=1, j=15:19, 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 <- merge_at(tbl, i=2, j=12:13, part="header") tbl <- hline(tbl, i=1, j=3:13, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=1, j=15:19, 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=2, j=12:13, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=3, j=3:13, border=fp_border(width=1), part="header") tbl <- hline(tbl, i=3, j=15:19, 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 <- italic(tbl, i=2, j=2, part="body") tbl <- align(tbl, align="center", part="all") tbl <- align(tbl, j = c("label1","label2"), align="left", part="body") tbl <- width(tbl, j =~ label1, width=1) tbl <- width(tbl, j =~ label2, width=.7) tbl <- width(tbl, j =~ agy_er + com_er + 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 + blank3, width=.1) tbl <- width(tbl, j =~ blank4, width=.3) tbl <- height(tbl, i=c(2,4,6,10,14), height=.1, part="body") tbl doc <- read_docx() doc <- body_add_flextable(doc, value = tbl) print(doc, target = "Table_SOM-R4.docx") rm(dat1,dat2,dat3,dat4,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$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zagytsgrp <- scale(dat$agy_ts - ave(dat$agy_ts, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zcomtsgrp <- scale(dat$com_ts - ave(dat$com_ts, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zagyacgrp <- scale(dat$agy_ac - ave(dat$agy_ac, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zcomacgrp <- scale(dat$com_ac - ave(dat$com_ac, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zagybfgrp <- scale(dat$agy_bf - ave(dat$agy_bf, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zcombfgrp <- scale(dat$com_bf - ave(dat$com_bf, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zagrgrp <- scale(dat$agr - ave(dat$agr, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zcnsgrp <- scale(dat$cns - ave(dat$cns, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zextgrp <- scale(dat$ext - ave(dat$ext, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zopngrp <- scale(dat$opn - ave(dat$opn, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) dat$zneugrp <- scale(dat$neu - ave(dat$neu, dat$cnt_now, FUN = function(x) mean(x, na.rm=T))) ##compute country-level religiosiy (z-standardized) dat1 <- aggregate(relig ~ cnt_now, dat, mean) dat1$zcrel <- scale(dat1$relig) dat1$relig <- NULL dat <- merge(dat, dat1, by = "cnt_now", all = TRUE) ###correlation between computed country-level religiosiy and ###country-level religiosity index based on Gallup World Poll data (Joshanloo & Gebauer, 2020) dat2 <- as.data.frame(fread("https://madata.bib.uni-mannheim.de/330/17/crel_gallup.csv", header = T, sep = ",")) dat3 <- merge(dat1, dat2, by = "cnt_now") dat3 <- subset(dat3, select = c(zcrel,crel_gallup)) print(corr.test(dat3), short=F) rm(dat1,dat2,dat3) ##re-center personality traits and country-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$zcrelmax <- dat$zcrel - abs(max(dat$zcrel)) dat$zcrelmin <- dat$zcrel + abs(min(dat$zcrel)) ##select variables dat <- subset(dat, select = c(record_id, age, female, langu, cnt_now, zse:zcrelmin)) #####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$cnt_now <- as.factor(dat$cnt_now) julia_assign("dat", dat) ##BIG TWO TS julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zcrel * zcomtsgrp + zrelgrp * zcrel * zagytsgrp + (1 + zrelgrp + zcomtsgrp + zagytsgrp | cnt_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 * zcrelmax * zcomtsgrphigh + zrelgrp * zcrelmax * zagytsgrplow + (1 + zrelgrp + zcomtsgrphigh + zagytsgrplow | cnt_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 * zcrelmax * zcomtsgrpveryhigh + zrelgrp * zcrelmax * zagytsgrpverylow + (1 + zrelgrp + zcomtsgrpveryhigh + zagytsgrpverylow | cnt_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 * zcrelmax * zcomtsgrp + zrelgrp * zcrelmax * zagytsgrp + (1 + zrelgrp + zcomtsgrp + zagytsgrp | cnt_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 * zcrelmax * zcomtsgrplow + zrelgrp * zcrelmax * zagytsgrphigh + (1 + zrelgrp + zcomtsgrplow + zagytsgrphigh | cnt_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 * zcrelmax * zcomtsgrpverylow + zrelgrp * zcrelmax * zagytsgrpveryhigh + (1 + zrelgrp + zcomtsgrpverylow + zagytsgrpveryhigh | cnt_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 * zcrelmin * zcomtsgrphigh + zrelgrp * zcrelmin * zagytsgrplow + (1 + zrelgrp + zcomtsgrphigh + zagytsgrplow | cnt_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 * zcrelmin * zcomtsgrpveryhigh + zrelgrp * zcrelmin * zagytsgrpverylow + (1 + zrelgrp + zcomtsgrpveryhigh + zagytsgrpverylow | cnt_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 * zcrelmin * zcomtsgrp + zrelgrp * zcrelmin * zagytsgrp + (1 + zrelgrp + zcomtsgrp + zagytsgrp | cnt_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 * zcrelmin * zcomtsgrplow + zrelgrp * zcrelmin * zagytsgrphigh + (1 + zrelgrp + zcomtsgrplow + zagytsgrphigh | cnt_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 * zcrelmin * zcomtsgrpverylow + zrelgrp * zcrelmin * zagytsgrpveryhigh + (1 + zrelgrp + zcomtsgrpverylow + zagytsgrpveryhigh | cnt_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 * zcrel * zcomacgrp + zrelgrp * zcrel * zagyacgrp + (1 + zrelgrp + zcomacgrp + zagyacgrp | cnt_now)), dat)") BigTwoAC <- julia_eval("BigTwoAC = 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 * zcrelmax * zcomacgrphigh + zrelgrp * zcrelmax * zagyacgrplow + (1 + zrelgrp + zcomacgrphigh + zagyacgrplow | cnt_now)), dat)") BigTwoAC_cmax_asshigh <- julia_eval("BigTwoAC_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 * zcrelmax * zcomacgrpveryhigh + zrelgrp * zcrelmax * zagyacgrpverylow + (1 + zrelgrp + zcomacgrpveryhigh + zagyacgrpverylow | cnt_now)), dat)") BigTwoAC_cmax_assveryhigh <- julia_eval("BigTwoAC_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 * zcrelmax * zcomacgrp + zrelgrp * zcrelmax * zagyacgrp + (1 + zrelgrp + zcomacgrp + zagyacgrp | cnt_now)), dat)") BigTwoAC_cmax_average <- julia_eval("BigTwoAC_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 * zcrelmax * zcomacgrplow + zrelgrp * zcrelmax * zagyacgrphigh + (1 + zrelgrp + zcomacgrplow + zagyacgrphigh | cnt_now)), dat)") BigTwoAC_cmax_conthigh <- julia_eval("BigTwoAC_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 * zcrelmax * zcomacgrpverylow + zrelgrp * zcrelmax * zagyacgrpveryhigh + (1 + zrelgrp + zcomacgrpverylow + zagyacgrpveryhigh | cnt_now)), dat)") BigTwoAC_cmax_contveryhigh <- julia_eval("BigTwoAC_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 * zcrelmin * zcomacgrphigh + zrelgrp * zcrelmin * zagyacgrplow + (1 + zrelgrp + zcomacgrphigh + zagyacgrplow | cnt_now)), dat)") BigTwoAC_cmin_asshigh <- julia_eval("BigTwoAC_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 * zcrelmin * zcomacgrpveryhigh + zrelgrp * zcrelmin * zagyacgrpverylow + (1 + zrelgrp + zcomacgrpveryhigh + zagyacgrpverylow | cnt_now)), dat)") BigTwoAC_cmin_assveryhigh <- julia_eval("BigTwoAC_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 * zcrelmin * zcomacgrp + zrelgrp * zcrelmin * zagyacgrp + (1 + zrelgrp + zcomacgrp + zagyacgrp | cnt_now)), dat)") BigTwoAC_cmin_average <- julia_eval("BigTwoAC_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 * zcrelmin * zcomacgrplow + zrelgrp * zcrelmin * zagyacgrphigh + (1 + zrelgrp + zcomacgrplow + zagyacgrphigh | cnt_now)), dat)") BigTwoAC_cmin_conthigh <- julia_eval("BigTwoAC_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 * zcrelmin * zcomacgrpverylow + zrelgrp * zcrelmin * zagyacgrpveryhigh + (1 + zrelgrp + zcomacgrpverylow + zagyacgrpveryhigh | cnt_now)), dat)") BigTwoAC_cmin_contveryhigh <- julia_eval("BigTwoAC_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 BF julia_eval("fm = fit(LinearMixedModel, @formula(zse ~ 1 + zrelgrp * zcrel * zcombfgrp + zrelgrp * zcrel * zagybfgrp + (1 + zrelgrp + zcombfgrp + zagybfgrp | cnt_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 * zcrelmax * zcombfgrphigh + zrelgrp * zcrelmax * zagybfgrplow + (1 + zrelgrp + zcombfgrphigh + zagybfgrplow | cnt_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 * zcrelmax * zcombfgrpveryhigh + zrelgrp * zcrelmax * zagybfgrpverylow + (1 + zrelgrp + zcombfgrpveryhigh + zagybfgrpverylow | cnt_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 * zcrelmax * zcombfgrp + zrelgrp * zcrelmax * zagybfgrp + (1 + zrelgrp + zcombfgrp + zagybfgrp | cnt_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 * zcrelmax * zcombfgrplow + zrelgrp * zcrelmax * zagybfgrphigh + (1 + zrelgrp + zcombfgrplow + zagybfgrphigh | cnt_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 * zcrelmax * zcombfgrpverylow + zrelgrp * zcrelmax * zagybfgrpveryhigh + (1 + zrelgrp + zcombfgrpverylow + zagybfgrpveryhigh | cnt_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 * zcrelmin * zcombfgrphigh + zrelgrp * zcrelmin * zagybfgrplow + (1 + zrelgrp + zcombfgrphigh + zagybfgrplow | cnt_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 * zcrelmin * zcombfgrpveryhigh + zrelgrp * zcrelmin * zagybfgrpverylow + (1 + zrelgrp + zcombfgrpveryhigh + zagybfgrpverylow | cnt_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 * zcrelmin * zcombfgrp + zrelgrp * zcrelmin * zagybfgrp + (1 + zrelgrp + zcombfgrp + zagybfgrp | cnt_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 * zcrelmin * zcombfgrplow + zrelgrp * zcrelmin * zagybfgrphigh + (1 + zrelgrp + zcombfgrplow + zagybfgrphigh | cnt_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 * zcrelmin * zcombfgrpverylow + zrelgrp * zcrelmin * zagybfgrpveryhigh + (1 + zrelgrp + zcombfgrpverylow + zagybfgrpveryhigh | cnt_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 * zcrel * zagrgrp + zrelgrp * zcrel * zcnsgrp + zrelgrp * zcrel * zopngrp + zrelgrp * zcrel * zextgrp + zrelgrp * zcrel * zneugrp + (1 + zrelgrp + zagrgrp + zcnsgrp + zopngrp + zextgrp + zneugrp | cnt_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 * zcrelmax * zagrgrphigh + zrelgrp * zcrelmax * zcnsgrplow + zrelgrp * zcrelmax * zopngrplow + zrelgrp * zcrelmax * zextgrplow + zrelgrp * zcrelmax * zneugrphigh + (1 + zrelgrp + zagrgrphigh + zcnsgrplow + zopngrplow + zextgrplow + zneugrphigh | cnt_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 * zcrelmax * zagrgrpveryhigh + zrelgrp * zcrelmax * zcnsgrpverylow + zrelgrp * zcrelmax * zopngrpverylow + zrelgrp * zcrelmax * zextgrpverylow + zrelgrp * zcrelmax * zneugrpveryhigh + (1 + zrelgrp + zagrgrpveryhigh + zcnsgrpverylow + zopngrpverylow + zextgrpverylow + zneugrpveryhigh | cnt_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 * zcrelmax * zagrgrp + zrelgrp * zcrelmax * zcnsgrp + zrelgrp * zcrelmax * zopngrp + zrelgrp * zcrelmax * zextgrp + zrelgrp * zcrelmax * zneugrp + (1 + zrelgrp + zagrgrp + zcnsgrp + zopngrp + zextgrp + zneugrp | cnt_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 * zcrelmax * zagrgrplow + zrelgrp * zcrelmax * zcnsgrphigh + zrelgrp * zcrelmax * zopngrphigh + zrelgrp * zcrelmax * zextgrphigh + zrelgrp * zcrelmax * zneugrplow + (1 + zrelgrp + zagrgrplow + zcnsgrphigh + zopngrphigh + zextgrphigh + zneugrplow | cnt_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 * zcrelmax * zagrgrpverylow + zrelgrp * zcrelmax * zcnsgrpveryhigh + zrelgrp * zcrelmax * zopngrpveryhigh + zrelgrp * zcrelmax * zextgrpveryhigh + zrelgrp * zcrelmax * zneugrpverylow + (1 + zrelgrp + zagrgrpverylow + zcnsgrpveryhigh + zopngrpveryhigh + zextgrpveryhigh + zneugrpverylow | cnt_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 * zcrelmin * zagrgrphigh + zrelgrp * zcrelmin * zcnsgrplow + zrelgrp * zcrelmin * zopngrplow + zrelgrp * zcrelmin * zextgrplow + zrelgrp * zcrelmin * zneugrphigh + (1 + zrelgrp + zagrgrphigh + zcnsgrplow + zopngrplow + zextgrplow + zneugrphigh | cnt_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 * zcrelmin * zagrgrpveryhigh + zrelgrp * zcrelmin * zcnsgrpverylow + zrelgrp * zcrelmin * zopngrpverylow + zrelgrp * zcrelmin * zextgrpverylow + zrelgrp * zcrelmin * zneugrpveryhigh + (1 + zrelgrp + zagrgrpveryhigh + zcnsgrpverylow + zopngrpverylow + zextgrpverylow + zneugrpveryhigh | cnt_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 * zcrelmin * zagrgrp + zrelgrp * zcrelmin * zcnsgrp + zrelgrp * zcrelmin * zopngrp + zrelgrp * zcrelmin * zextgrp + zrelgrp * zcrelmin * zneugrp + (1 + zrelgrp + zagrgrp + zcnsgrp + zopngrp + zextgrp + zneugrp | cnt_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 * zcrelmin * zagrgrplow + zrelgrp * zcrelmin * zcnsgrphigh + zrelgrp * zcrelmin * zopngrphigh + zrelgrp * zcrelmin * zextgrphigh + zrelgrp * zcrelmin * zneugrplow + (1 + zrelgrp + zagrgrplow + zcnsgrphigh + zopngrphigh + zextgrphigh + zneugrplow | cnt_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 * zcrelmin * zagrgrpverylow + zrelgrp * zcrelmin * zcnsgrpveryhigh + zrelgrp * zcrelmin * zopngrpveryhigh + zrelgrp * zcrelmin * zextgrpveryhigh + zrelgrp * zcrelmin * zneugrpverylow + (1 + zrelgrp + zagrgrpverylow + zcnsgrpveryhigh + zopngrpveryhigh + zextgrpveryhigh + zneugrpverylow | cnt_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 1 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", "cnt rel", "com", "agy", "rel x cnt rel", "rel x com", "cnt rel x com", "rel x agy", "cnt rel x agy", "rel x cnt rel x com", "rel x cnt rel x agy", "", "", "", "", "", "", "", "", "", "", "", "") dat1[,11] <- c("(intercept)","rel", "cnt rel", "agr", "cns", "opn", "ext", "neu", "rel x cnt rel", "rel x agr", "cnt rel x agr", "rel x cns", "cnt rel x cns", "rel x opn", "cnt rel x opn", "rel x ext", "cnt rel x ext", "rel x neu", "cnt rel x neu", "rel x cnt rel x agr", "rel x cnt rel x cns", "rel x cnt rel x opn", "rel x cnt rel x ext", "rel x cnt 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 = "Table1.docx") rm(dat1,head,tbl,doc) #simple slopes for the religiosity x country-level religiosity interactions ##mean association between religiosity and self-esteem in least religious countries 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 countries: zPE = ", round(cmin_zpe, digits=2), " 95% CI [", round(cmin_lowerCI, digits=3), ", ", round(cmin_upperCI, digits=2), "]") ##mean association between religiosity and self-esteem in most religious countries 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 countries: 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 1 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 countries","AC/least religious countries","BF/least religious countries","TS/most religious countries","AC/most religious countries","BF/most religious countries", "TS/least religious countries","AC/least religious countries","BF/least religious countries","TS/most religious countries","AC/most religious countries","BF/most religious countries", "TS/least religious countries","AC/least religious countries","BF/least religious countries","TS/most religious countries","AC/most religious countries","BF/most religious countries", "TS/least religious countries","AC/least religious countries","BF/least religious countries","TS/most religious countries","AC/most religious countries","BF/most religious countries", "TS/least religious countries","AC/least religious countries","BF/least religious countries","TS/most religious countries","AC/most religious countries","BF/most religious countries") 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 1a. 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) cnt <- c("least religious countries","most religious countries","least religious countries","most religious countries","least religious countries","most religious countries", "least religious countries","most religious countries","least religious countries","most religious countries") 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(cnt,ass_cont,zPE,lowerCI,upperCI) rm(cnt,ass_cont,zPE,lowerCI,upperCI) dat1$cnt <- factor(dat1$cnt, levels=unique(as.character(dat1$cnt))) dat1$ass_cont <- factor(dat1$ass_cont, levels=unique(as.character(dat1$ass_cont))) BigFivePlot <- ggplot(data=dat1, aes(fill=cnt, x=ass_cont, y=zPE, color=cnt, group=cnt)) + 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 1b. Big Five") + ylab("religiosity-esteem association (zPE)") figure <- grid.arrange(BigTwoPlot,BigFivePlot) ggsave("Figure1.jpg", plot=figure, device="jpg", width=210, height=210, units="mm", dpi=720) ggsave("Figure1.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=3), " 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, "]")