---
title: "Data analysis script - Zmigrod et al 2021 Philosophical Transactions"
output:
  html_document:
    df_print: paged
  word_document: default
  pdf_document: default
---

## Primary data analyses for Ideology/Ontology project

v20201201

Initial R code written by Leor Zmigrod

Converted to RMarkdown by Russ Poldrack


#### Load packages

```{r loadPackages,echo=FALSE,warning=FALSE,message=FALSE}
library(readxl)
library(ggpubr)
library(ggplot2)
library(dplyr)
library(knitr)
library(visreg)
library(mgcv)
library(cowplot)
library(jtools)
library(ggthemes)
library(corrplot)
library(Hmisc)
library(readxl)
library(psych)
library(lm.beta)
library(broom)
library(fmsb)
library(tibble)
library(GPArotation)
library(tidyr)
library(BayesFactor)
library(foreach)
library(RColorBrewer)
library(viridis)
library(BMA)
library(gridExtra)
library(tidyselect)
```

### Load Data

```{r loadData}
combined_data <- read_excel("Data - Zmigrod et al 2021 Philosophical Transactions.xlsx")
```

### Exploratory factor analysis

Run EFA on ideology data.

```{r EFA}
# relevant data for EFA on ideological orientations
efa_data <- data.frame(
        secs_social = combined_data$SocialConservatism, 
        secs_econ = combined_data$EconConservatism,
        nation = combined_data$Nationalism, 
        auth = combined_data$Authoritarianism,
        SJ = combined_data$SJ, 
        sdo = combined_data$SDO,
        patriot = combined_data$Patriotism,
        progrSum = combined_data$IngroupAttitude, #Note that this variable is at times called Ingroup Attitude / Progroup Attitude / or Extreme Pro-Group Actions.
        dogma = combined_data$Dogmatism,
        IH1 = -1*combined_data$IH1, #IH factors reversed because conceptually inverse of dogmatism
        IH2 = -1*combined_data$IH2,
        IH3 = -1*combined_data$IH3,
        IH4 = -1*combined_data$IH4,
        rel_atten = combined_data$ReligiousAttendance,
        rel_pray = combined_data$ReligiousPrayer,
        rel_imp = combined_data$ReligionImportance
      )

# EFA analysis
scree(efa_data) # scree plot
fa.parallel(efa_data, fm = "ml") # parallel analysis
Ideol_EFA <- fa(efa_data, 
                nfactors=3, 
                rotate="oblimin")     # EFA oblimin for 3 factors

# obtain factor loadings and add them to data
factor_scores <- as.data.frame(Ideol_EFA$scores)

combined_data <- add_column(
  combined_data, 
  Ideol_Factor1_Politics = factor_scores$MR1, 
  Ideol_Factor2_Dogmatism = factor_scores$MR2, 
  Ideol_Factor3_Religion = factor_scores$MR3  
  )

```

Print correlogram of EFA factor loadings for supplementary figure

```{r correlogramFigure, fig.width=8, fig.height=10}
## Correlogram of EFA factor loadings - Supplementary Figure

corrplot(Ideol_EFA$loadings, method="color", 
         addCoef.col = "black", number.cex = .75, 
         tl.col = "black", cl.pos = "b", 
         cl.length=4)

# save figure to file

dev.copy(pdf,"correlogram.pdf", width=8, height=10)
dev.off()

```


Print correlogram of ideological orientations for supplementary figure

```{r correlogramIdeolFigure, fig.width=8, fig.height=10}
## Correlogram of ideological orientations - Supplementary Figure

M <- cor(efa_data, use= "pairwise.complete.obs")
colnames(M) <- c("Social Cons.", "Econ. Cons.", "Nationalism", "Author.", "SJ", "SDO", "Patriotism", "Ingroup", "Dogmatism", "IH1", "IH2", "IH3", "IH4", "Rel. Attend.", "Rel. Prayer", "Rel. Import."  )
rownames(M) <- colnames(M)
P <- cor.mtest(efa_data)

corrplot(M, method="color", 
         addCoef.col = "black", number.cex = .75, 
         tl.col = "black", cl.pos = "b", 
         cl.length=4, p.mat = P$p, insig="blank")

# save figure to file

dev.copy(pdf,"correlogram_ideol.pdf", width=10, height=10)
dev.off()



```





#### Regression Data

First create data frame for regression analyses.

```{r regressionData}
data_regressions <- data.frame(
          SubjectID=combined_data$SubjectID,
          
          # demographic variables
          Age = combined_data$Age,
          Edu = combined_data$Edu,
          Gender = combined_data$Gender,
          Income = combined_data$Income,
          
          ## relevant ideology variables:
          secs_social = combined_data$SocialConservatism, 
          secs_econ = combined_data$EconConservatism,
          nation = combined_data$Nationalism, 
          auth = combined_data$Authoritarianism,
          SJ = combined_data$SJ, 
          sdo = combined_data$SDO,
          patriot_scale = combined_data$Patriotism,
          progrSum = combined_data$IngroupAttitude,
          
          dogma = combined_data$Dogmatism,
          IH1 = combined_data$IH1,
          IH2 = combined_data$IH2,
          IH3 = combined_data$IH3,
          IH4 = combined_data$IH4,
          
          rel_atten = combined_data$ReligiousAttendance,
          rel_pray = combined_data$ReligiousPrayer,
          rel_imp = combined_data$ReligionImportance,
          
          EFA_MR1_Pol = combined_data$Ideol_Factor1_Politics,
          EFA_MR2_Dog = combined_data$Ideol_Factor2_Dogmatism,
          EFA_MR3_Rel = combined_data$Ideol_Factor3_Religion,
          
          #policy attitudes
          abortion = combined_data$SECS_Abortion_reversed,
          business =combined_data$SECS_Business_Warmth,
          familyunit = combined_data$SECS_FamilyUnit_Warmth,
          fiscalres = combined_data$SECS_FiscalResponsibility_Warmth,
          gunown = combined_data$SECS_GunOwnership_Warmth,
          limgovt = combined_data$SECS_LimitedGovt_Warmth,
          security = combined_data$SECS_MilitaryNationalSecurity_Warmth,
          patriot = combined_data$SECS_Patriotism_Warmth,
          religion = combined_data$SECS_Religion_Warmth,
          marriage = combined_data$SECS_TraditionalMarriage_Warmth,
          values = combined_data$SECS_TraditionalValues_Warmth,
          welfare = combined_data$SECS_WelfareBen_reversed,
          
          samesex = combined_data$SECS_SameSexMarriage_reversed,
          immigr = combined_data$SECS_RestrictingImmigration_Warmth,
          parisagre = combined_data$SECS_ParisAgreement_reversed,
          maga = combined_data$SECS_MAGA_Warmth,
          
          ## survey ontology factors
          agree = combined_data$agree,
          fin_risk = combined_data$fin_risk,
          soc_risk = combined_data$soc_risk,
          ethic_risk = combined_data$ethic_risk,
          eat_cont= combined_data$eat_cont,
          risk_per = combined_data$risk_per,
          goal_dir = combined_data$goal_dir,
          rew_sens = combined_data$rew_sens,
          emot_cont = combined_data$emot_cont,
          impuls = combined_data$impuls,
          mindful = combined_data$mindful,
          senseseek = combined_data$senseseek,
          
          ## cognitive ontology factors
          SpeededIP = combined_data$SpeededIP,
          StrategicIP = combined_data$StrategicIP,
          Discount = combined_data$Discount,
          Caution = combined_data$Caution,
          Percep = combined_data$Percep
        )


```

#### Transform data to address skewness

```{r transformData}

# Zeynep Enkavi's skewness function (with permission)
require(psych)
require(tidyverse)
require(RCurl)
helper_func_path = 'https://raw.githubusercontent.com/zenkavi/SRO_Retest_Analyses/master/code/helper_functions/'
eval(parse(text = getURL(paste0(helper_func_path,'get_numeric_cols.R'), ssl.verifypeer = FALSE)))

"%w/o%" <- function(x, y) x[!x %in% y]

pos_log <- function(column){
  col_min = min(column, na.rm=T)
  a = 1-col_min
  column = column+a
  return(log(column))
}

neg_log <- function(column){
  col_max = max(column, na.rm=T)
  column = col_max+1-column
  return(log(column))
}

sqrt_const <- function(column, const){
  if(missing(const)|is.null(const)){
    col_min = min(column, na.rm=T)
    const = 1-col_min
  }
  column = column+const
  return(sqrt(column))
}

asin_count <- function(column, tpc){
  #tpc = total possible counts
  if(missing(tpc)|is.null(tpc)){
    tpc = max(column, na.rm=T)
  }
  column = sqrt(column/tpc)
  return(asin(column))
}

#transform: "log", "sqrt", "asin", "custom"
#if transform is custom: 
#mandatory argument "func" containing function to be applied must be provided
#optional argument "suffix" containing string to be added to transformed cols can be provided; else default string is .customTr

transform_remove_skew = function(data, columns=get_numeric_cols(df1=data, df2=data), threshold = 1, drop=FALSE, transform='log', ...){
  
  tmp = as.data.frame(apply(data[,columns],2,skew))
  names(tmp) = c("skew")
  tmp$dv = row.names(tmp)
  tmp = tmp %>% 
    filter(abs(skew)>threshold)
  
  skewed_variables = tmp$dv
  skew_subset = data[, skewed_variables]
  
  if(length(tmp$dv[tmp$skew>0])>0){
    positive_subset = data.frame(data[,tmp$dv[tmp$skew>0]])
    names(positive_subset) = tmp$dv[tmp$skew>0]
  }
  if(length(tmp$dv[tmp$skew<0]>0)){
    negative_subset = data.frame(data[,tmp$dv[tmp$skew<0]])
    names(negative_subset) = tmp$dv[tmp$skew<0]
  }
  
  cat(rep('*', 40))
  cat('\n')
  cat(paste0('Using transformation:', transform))
  cat('\n')
  
  if(exists('positive_subset')){ 
    if(!is.null(ncol(positive_subset)) & ncol(positive_subset)>0){
      
      if(transform == "log"){
        suffix = '.logTr'
        positive_subset = as.data.frame(apply(positive_subset, 2, pos_log))
      }
      
      if(transform == "sqrt"){
        suffix = '.sqrtTr'
        positive_subset = as.data.frame(apply(positive_subset, 2, sqrt_const, const=list(...)$const))
      }
      
      if(transform == "asin"){
        suffix = '.asinTr'
        positive_subset = as.data.frame(apply(positive_subset, 2, asin_count, tpc=list(...)$tpc))
      }
      
      if(transform == "custom"){
        if(!is.null(list(...)$suffix)){
          suffix = list(...)$suffix
        } else {
          suffix = '.customTr'
        }
        positive_subset = as.data.frame(apply(positive_subset, 2, list(...)$func))
      }
      
      successful_transforms = as.data.frame(apply(positive_subset, 2, skew))
      names(successful_transforms) = c('skew')
      successful_transforms$dv = row.names(successful_transforms)
      successful_transforms = successful_transforms %>% filter(abs(skew)<threshold)
      successful_transforms = successful_transforms$dv
      successful_transforms = positive_subset[,successful_transforms]
      dropped_vars = names(positive_subset) %w/o% names(successful_transforms)
      
      cat(rep('*', 40))
      cat('\n')
      cat(paste0(length(names(positive_subset)) ,' data positively skewed data were transformed:'))
      cat('\n')
      cat(names(positive_subset), sep = '\n')
      cat(rep('*', 40))
      cat('\n')
      
      # replace transformed variables
      data = data[,-c(which(names(data) %in% names(positive_subset)))]
      
      if(drop == TRUE){
        names(successful_transforms) = paste0(names(successful_transforms))
        cat(rep('*', 40))
        cat('\n')
        cat(paste0('Dropping ', length(dropped_vars) ,' positively skewed data that could not be transformed successfully:'))
        cat('\n')
        cat(dropped_vars, sep = '\n')
        cat(rep('*', 40))
        cat('\n')
        data = cbind(data, successful_transforms)
      } else{
        names(positive_subset) = paste0(names(positive_subset))
        cat(rep('*', 40))
        cat('\n')
        cat(paste0(length(dropped_vars) ,' positively skewed data could not be transformed successfully:'))
        cat('\n')
        cat(dropped_vars, sep = '\n')
        cat(rep('*', 40))
        cat('\n')
        data = cbind(data, positive_subset)
      }
    } 
  }else{
    cat(rep('*', 40))
    cat('\n')
    cat('No positively skewed variables found.')
    cat(rep('*', 40))
    cat('\n')
  }
  
  
  if(exists('negative_subset')){
    if(!is.null(ncol(negative_subset)) & ncol(negative_subset)>0){
      
      if(transform == "log"){
        suffix = '.ReflogTr'
        negative_subset = as.data.frame(apply(negative_subset, 2, neg_log))
      }
      
      if(transform == "sqrt"){
        suffix = '.sqrtTr'
        negative_subset = as.data.frame(apply(negative_subset, 2, sqrt_const, const=list(...)$const))
      }
      
      if(transform == "asin"){
        suffix = '.asinTr'
        negative_subset = as.data.frame(apply(negative_subset, 2, asin_count, tpc=list(...)$tpc))
      }
      if(transform == "custom"){
        if(!is.null(list(...)$suffix)){
          suffix = list(...)$suffix
        } else {
          suffix = '.customTr'
        }
        negative_subset = as.data.frame(apply(negative_subset, 2, list(...)$func))
      }
      
      successful_transforms = as.data.frame(apply(negative_subset, 2, skew))
      names(successful_transforms) = c('skew')
      successful_transforms$dv = row.names(successful_transforms)
      successful_transforms = successful_transforms %>% filter(abs(skew)<1)
      successful_transforms = successful_transforms$dv
      successful_transforms = negative_subset[,successful_transforms]
      dropped_vars = names(negative_subset) %w/o% names(successful_transforms)
      
      cat(rep('*', 40))
      cat('\n')
      cat(paste0(length(names(negative_subset)) ,' data negatively skewed data were transformed:'))
      cat('\n')
      cat(names(negative_subset), sep = '\n')
      cat(rep('*', 40))
      cat('\n')
      
      # replace transformed variables
      data = data[,-c(which(names(data) %in% names(negative_subset)))]
      if(drop == TRUE){
        names(successful_transforms) = paste0(names(successful_transforms))
        cat(rep('*', 40))
        cat('\n')
        cat(paste0('Dropping ', length(dropped_vars) ,' negatively skewed data that could not be transformed successfully:'))
        cat('\n')
        cat(dropped_vars, sep = '\n')
        cat(rep('*', 40))
        cat('\n')
        data = cbind(data, successful_transforms)
      } else{
        names(negative_subset) = paste0(names(negative_subset))
        cat(rep('*', 40))
        cat('\n')
        cat(paste0(length(dropped_vars) ,' negatively skewed data could not be transformed successfully:'))
        cat('\n')
        cat(dropped_vars, sep = '\n')
        cat(rep('*', 40))
        cat('\n')
        data = cbind(data, negative_subset)
      } 
    } 
  } else{
    cat(rep('*', 40))
    cat('\n')
    cat('No negatively skewed variables found.')
    cat(rep('*', 40))
    cat('\n')
  }
  
  return(data)
}


# address skewness in data
# removing IH2 variable from skewness transformation because it doesn't transform effectively 
# (i.e. transform does not work with this function and changes values in way
# that is incongruent with other IH variables)

data <- as_tibble(data_regressions)
data_to_be_transformed <- dplyr::select(data, -matches("IH2"))
transformed_data <- transform_remove_skew(data=data_regressions, 5:ncol(data_regressions), drop=FALSE)
data_regressions <- mutate(transformed_data, IH2 = data_regressions$IH2)


```

#### Set up functions to run and summarize regressions

```{r regressionFunctions,echo=FALSE}
# function to run regressions and return results

run_regression <- function(modelname,data,
                           formula,print.summary=FALSE){
  # modelname: label for model
  # varnames: human-readable variable names
  # data: data frame with data
  # formula: a formula object for the model
  # print.summary: should we print a summary?
  
  # set up list to store results
  results <- list(modelname=modelname)
  
  # run regression and extract statistics
  results$lm <- lm(formula,data=data)
  results$lm.beta <- lm.beta(results$lm)
  results$std <- data.frame(t(results$lm.beta$standardized.coefficients))
  
  p <- tidy(results$lm) %>% 
      select(p.value) %>% 
      t() %>% 
      as.data.frame()
  
  names(p) <- foreach(i=names(results$std)) %do% sprintf('pval:%s',i)
  results$pval <- p
  # now fix std names
  names(results$std) <- foreach(i=names(results$std)) %do% sprintf('stdbeta:%s',i)
  
  results$rs <- glance(results$lm)
  
  # optionally print summary
  if (print.summary){
    summary(results$reg)
  }
  return(results)
}  

#extract necessary elements from results object and return as data frame

summarize_results <- function(results){
  r <- cbind(results$std, results$pval,  results$rs)
  row.names(r) <- results$modelname
  return(r)
}

```

#### Set up regressor sets for analyses

```{r regressionVars}

varsets <- c('C'='SpeededIP + StrategicIP + Discount + Caution + Percep',
            'S'='agree + senseseek + risk_per + goal_dir + rew_sens +  soc_risk + fin_risk + ethic_risk + eat_cont + mindful + impuls  + emot_cont',
            'D'='Age + Edu + Gender + Income')
varsets['SD'] <- paste(varsets['S'],varsets['D'],sep=' + ')
varsets['CD'] <- paste(varsets['C'],varsets['D'],sep=' + ')
varsets['CS'] <- paste(varsets['S'], varsets['C'], sep=' + ')
varsets['CSD'] <- paste(varsets['C'],varsets['S'],varsets['D'],sep=' + ')


# set up target variables

targetvars <- c('SECS_Social'='secs_social',
                'SECS_Econ'='secs_econ',
                'Nationalism'='nation',
                'Authoritarianism'='auth',
                'System Justification'='SJ',
                'SDO'='sdo',
                'Dogmatism'='dogma',
                'IH1'='IH1',
                'IH2'='IH2',
                'IH3'='IH3',
                'IH4'='IH4',
                'Progroup_Actions'='progrSum',
                'Abortion'='abortion',
                'Business'='business',
                'FamilyUnit'='familyunit',
                'FiscalRespons'='fiscalres',
                'GunOwnership'='gunown',
                'LimitedGovt'='limgovt',
                'Security'='security',
                'Patriotism'='patriot',
                'Religion'='religion',
                'TraditionalMarriage'='marriage',
                'TraditionalValues'='values',
                'SameSexMarriage'='samesex',
                'RestrictingImmigration'='immigr',
                'ParisAgreement'='parisagre',
                'MAGA'='maga',
                'Religious_Attendance'='rel_atten',
                'Religious_Prayer'='rel_pray',
                'Religion_Importance'='rel_imp',
                'Patriotism scale'='patriot_scale',
                'welfare'='welfare',
                'EFA_MR1_Pol'='EFA_MR1_Pol',
                'EFA_MR2_Dog'='EFA_MR2_Dog',
                'EFA_MR3_Rel'='EFA_MR3_Rel'
)

```

#### Loop through all models for all target variables, estimate and save results to csv.  

In addition, all model results are saved to `results` list and all models saved to `models`.

```{r runRegressionModels}
all_results <- list()
results <- list()
models <- list()

for (v in names(varsets)){
  all_results[[v]]=c()
  for (t in names(targetvars)){
    modelname <- paste(t,v,sep='_')
    models[[modelname]] <- sprintf('%s ~ %s',targetvars[t],varsets[v])
    # cat(sprintf('running model: %s\n',modelname))
    results[[modelname]] <- run_regression(modelname,data=data_regressions,
                              formula=formula(models[[modelname]]))
    if (is.null(all_results[[v]])){
      all_results[[v]] <- summarize_results(results[[modelname]])
      
    } else {
      all_results[[v]] <- rbind(all_results[[v]],summarize_results(results[[modelname]]))
    }
  }
  write.csv(all_results[[v]], sprintf("linear_reg_all_outcomes_%s.csv",v))
}

```


### Figure 1 - Linear regression results

```{r mkFigure1, fig.width=10, fig.height=8}

# 1) Extracting data for R-squared values of all ideological orientations
R2Vals <- data.frame(
          Targets = targetvars,
          D = all_results$D$r.squared,
          C = all_results$C$r.squared,
          S = all_results$S$r.squared,
          CSD = all_results$CSD$r.squared)

# 2) Ideological orientations according to domain        

R2Values_Politics <- R2Vals %>%
          filter(Targets %in% c("secs_social",  
                                "secs_econ", 
                                "nation", 
                                "patriot_scale", 
                                "auth", 
                                "sdo", 
                                "SJ", 
                                "progrSum", 
                                "EFA_MR1_Pol")) %>% 
          gather(key = "Category", 
                 value = "R2", c("D", "C", "S", "CSD"), 
                 factor_key=TRUE)

R2Values_Dogmatism_Religiosity <- R2Vals %>%
          filter(Targets %in% c("dogma",  
                                "IH1", 
                                "IH2", 
                                "IH3", 
                                "IH4", 
                                "EFA_MR2_Dog", 
                                "rel_atten",  
                                "rel_imp", 
                                "rel_pray", 
                                "EFA_MR3_Rel")) %>% 
          gather(key = "Category", 
                 value = "R2", c("D", "C", "S", "CSD"), 
                 factor_key=TRUE)
        

R2Values_Dogmatism <- R2Vals %>%
          filter(Targets %in% c("dogma",  
                                "IH1", 
                                "IH2", 
                                "IH3", 
                                "IH4", 
                                "EFA_MR2_Dog")) %>% 
          gather(key = "Category", 
                 value = "R2", c("D", "C", "S", "CSD"), 
                 factor_key=TRUE)
        
R2Values_Religiosity <- R2Vals %>%
          filter(Targets %in% c("rel_atten",  
                                "rel_imp", 
                                "rel_pray", 
                                "EFA_MR3_Rel")) %>% 
          gather(key = "Category", 
                 value = "R2", c("D", "C", "S", "CSD"), 
                 factor_key=TRUE)

R2Values_Policies <- R2Vals %>%
          filter(Targets %in% c("abortion", 
                                "business", 
                                "familyunit", 
                                "fiscalres", 
                                "gunown",
                                "immigr", 
                                "limgovt", 
                                "maga", 
                                "marriage", 
                                "parisagre", 
                                "patriot" ,
                                "religion", 
                                "samesex", 
                                "security", 
                                "values", 
                                "welfare")) %>% 
          gather(key = "Category", value = "R2", c("D", "C", "S", "CSD"), factor_key=TRUE)
        
# 3) Figures   

plot_figure1 <- function(label,d,xlabels,scale_limits=NULL, fig_color){
  vlabels= c("Demographics", "Cognitive Tasks", "Personality Surveys", "Combined")
  hlines <- c(.05, .10, .15, .20, .25, .30, .35)

  p = ggplot(data = d, aes(x=Targets, y=R2, fill = Category))
  
  for (h in hlines){
      p <- p + geom_hline(yintercept=h, color="gray88")
  }  
    
  p<- p + geom_bar(stat="identity", position = position_dodge()) + 
    labs(x = "Outcome", y = "R2") +
    scale_fill_manual(values = fig_color , 
                    name = "Predictor Type", 
                    labels = vlabels) + 
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  if (!is.null(scale_limits)){
    p = p + scale_x_discrete(limits=scale_limits, 
                   labels = xlabels) 
  }

  print(p)
  dev.copy(pdf,sprintf("R2_%s.pdf",label), width=10, height=8)

  dev.off()
  svg(file = sprintf("R2_%s.svg", label), width=10, height=8)
  print(p)
  dev.off()
  return(p)
}

## color spectrums for Figure 1
# blue shades
Blue <- brewer.pal(11, "RdBu")
blue_shades <- Blue[8:11]


Fig1_Pol <- plot_figure1('politics',R2Values_Politics, fig_color=blue_shades,
            c("Social Conservatism", "Economic Conservatism", "Nationalism", "Patriotism", "Authoritarianism", "Social Dominance", "System Justification", "Extreme Pro-group Values", "Factor: Political Conservatism"),
            c("secs_social",  "secs_econ", "nation", "patriot_scale", "auth", "sdo", "SJ", "progrSum", "EFA_MR1_Pol"))

Fig1_Dog <- plot_figure1('dogmatism',R2Values_Dogmatism,fig_color=blue_shades,
             c("Dogmatism", "IH1: Intellectual Ego", "IH2: Belief Updating", "IH3: Respect for Others", "IH4: Overconfidence", "Factor: Dogmatism"),
             c("dogma",  "IH1", "IH2", "IH3", "IH4", "EFA_MR2_Dog"))

Fig1_Rel <- plot_figure1('religiosity',R2Values_Religiosity,fig_color=blue_shades,
             c("Religious Attendance",  "Religion Importance", "Religious Prayer", "Factor: Religiosity"),
             c("rel_atten",  "rel_imp", "rel_pray", "EFA_MR3_Rel"))

Fig1_DogRel <- plot_figure1('dogmatism_religiosity',R2Values_Dogmatism_Religiosity,fig_color=blue_shades,
             c("Dogmatism", "IH1: Intellectual Ego", "IH2: Belief Updating", "IH3: Respect for Others", "IH4: Overconfidence", "Factor: Dogmatism", "Religious Attendance",  "Religion Importance", "Religious Prayer", "Factor: Religiosity"),
             c("dogma",  "IH1", "IH2", "IH3", "IH4", "EFA_MR2_Dog","rel_atten",  "rel_imp", "rel_pray", "EFA_MR3_Rel"))


Fig1_Policies <- plot_figure1('policies',R2Values_Policies, fig_color=blue_shades, xlabels = c("Abortion", "Business", "Family Unit", "Fiscal Responsibility", "Gun Ownership", "Immigration" , "Limited Gov't", "MAGA", "Traditional Marriage", "Paris Climate Agree.", "Patriotism", "Religion", "Same-Sex Marriage", "National Security", "Traditional Values", "Welfare Benefits"), 
             c("abortion",  "business",  "familyunit",  "fiscalres", "gunown", "immigr",  "limgovt",  "maga", "marriage",  "parisagre", "patriot" ,  "religion", "samesex",  "security",  "values",  "welfare"))


Fig1_Panel <- ggarrange(Fig1_Pol, Fig1_DogRel, Fig1_Policies, ncol=1, nrow=3, common.legend = TRUE, legend = "right")

#save as pdf and svg
print(Fig1_Panel)
  dev.copy(pdf,sprintf("Fig1_Panel.pdf", Fig1_Panel), width=10, height=16)
  dev.off()
  svg(file = "Fig1_Panel.svg", width=10, height=16)
  print(Fig1_Panel)
  dev.off()
  


```



### Figure 4, 5, and supplement - Plot figures from cognitive ontology



```{r mkFiguresCognitiveOnly}

# Figure using cognitive predictors only
# Outcome: 3 factors
StdEstimates_C <-all_results$C[2:6]

StdEstimates_C <- rownames_to_column(StdEstimates_C, 
                                       var="Targets")
        
C_StdEst_Factors <- filter(StdEstimates_C, 
                             Targets %in% c("EFA_MR1_Pol_C", 
                                            "EFA_MR2_Dog_C", 
                                            "EFA_MR3_Rel_C"))

domain_levels_c <- list(
                    "Caution"='stdbeta:Caution', 
                    "Temporal Discounting"='stdbeta:Discount', 
                    "Perceptual Processing Time"='stdbeta:Percep', 
                    "Speed of Evidence Accumulation"='stdbeta:SpeededIP', 
                    "Strategic Information Processing"='stdbeta:StrategicIP'
                   )




domain_levels_reverse_c = names(domain_levels_c)
names(domain_levels_reverse_c)=domain_levels_c

C_StdEst_Factors_tidy <- gather(C_StdEst_Factors, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
    mutate(PsychDomain = factor(PsychDomain,
                                levels=domain_levels_c)) %>%
    mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_c))


p = ggplot(data = C_StdEst_Factors_tidy, 
       aes(x=Targets, y=StdEstimate, fill = PsychDomainNames))

hlines <- c(-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 

p<- p +   geom_bar(stat="identity", position = position_dodge()) + 
  labs(x = "Outcome", y = "Standardized Estimate", fill = "Cognitive Domain") + 
   scale_fill_viridis(discrete= TRUE, option="B", begin = .25, end=.8)  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  scale_x_discrete(limits=c("EFA_MR1_Pol_C", "EFA_MR2_Dog_C", "EFA_MR3_Rel_C"), 
                   labels = c("Political Conservatism", "Dogmatism", "Religiosity")) +
  theme(plot.margin = margin(1.5, 0, 0, .5, "cm")) 

print(p)
dev.copy(pdf,sprintf("stdbeta_Factors_Cognitive.pdf"), width=10, height=8)
dev.off()
svg(file = "stdbeta_Factors_Cognitive.svg", width=10, height=8)
  print(p)
  dev.off()
           


# outcome: political  variables
C_StdEst_Political <- filter(StdEstimates_C, 
                             Targets %in% c("SECS_Social_C", 
                                            "SECS_Econ_C", 
                                            "Nationalism_C",
                                            "Patriotism scale_C",
                                            "Authoritarianism_C",
                                            "System Justification_C",
                                            "SDO_C",
                                            "Progroup_Actions_C"))
                           
C_StdEst_Political_tidy <- gather(C_StdEst_Political, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
    mutate(PsychDomain = factor(PsychDomain,
                                levels=domain_levels_c)) %>%
    mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_c))


p = ggplot(data = C_StdEst_Political_tidy, 
       aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 

hlines <- c( -.3, -.25,-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 

p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
  labs(x = "Outcome", y = "Standardized Estimate", fill = "Cognitive Domain") + 
   scale_fill_viridis(discrete= TRUE, option="B", begin = .25, end=.8)  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  scale_x_discrete(limits=c("SECS_Social_C", 
                                            "SECS_Econ_C", 
                                            "Nationalism_C",
                                            "Patriotism scale_C", "Authoritarianism_C",
                                            "System Justification_C",
                                            "SDO_C",
                                            "Progroup_Actions_C" ), 
                  labels = c("Social Conservatism", "Economic Conservatism", "Nationalism", "Patriotism", "Authoritarianism", "System Justification", "Social Dominance", "Extreme Pro-Group Attitudes" )) +
  theme(plot.margin = margin(1.5, 0, 0, .5, "cm")) 

print(p)
dev.copy(pdf,sprintf("stdbeta_Political_Cognitive.pdf"), width=12, height=8)
dev.off() 
svg(file = "stdbeta_Political_Cognitive.svg", width=12, height=8)
  print(p)
  dev.off()


# dogmatism + religiosity outcome variables



C_StdEst_DogRel <- filter(StdEstimates_C, 
                             Targets %in% c("Dogmatism_C", 
                                            "IH1_C", 
                                            "IH2_C",
                                            "IH3_C",
                                            "IH4_C",
                                            "Religious_Attendance_C",
                                            "Religion_Importance_C",
                                            "Religious_Prayer_C"))
                           
C_StdEst_DogRel_tidy <- gather(C_StdEst_DogRel, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
    mutate(PsychDomain = factor(PsychDomain,
                                levels=domain_levels_c)) %>%
    mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_c))


p = ggplot(data = C_StdEst_DogRel_tidy, 
       aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 
hlines <- c(-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 

p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
  labs(x = "Outcome", y = "Standardized Estimate", fill = "Cognitive Domain") + 
   scale_fill_viridis(discrete= TRUE, option="B", begin = .25, end=.8)  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  scale_x_discrete(limits=c("Dogmatism_C", 
                                            "IH1_C", 
                                            "IH2_C",
                                            "IH3_C",
                                            "IH4_C",
                                            "Religious_Attendance_C",
                                            "Religion_Importance_C",
                                            "Religious_Prayer_C"), 
                  labels =  c("Dogmatism", "IH1: Intellectual Ego", "IH2: Belief Updating", "IH3: Respect for Others", "IH4: Overconfidence", "Religious Attendance",  "Religion Importance", "Religious Prayer")) +
  theme(plot.margin = margin(1.5, 0, 0, .5, "cm")) 

print(p)
dev.copy(pdf,sprintf("stdbeta_DogRel_Cognitive.pdf"), width=12, height=8)
dev.off() 
svg(file = "stdbeta_DogRel_Cognitive.svg", width=12, height=8)
  print(p)
  dev.off()



```

### Figure 4 and supplement - Plot figures from personality survey ontology

                                    
```{r make survey only figures}

StdEstimates_S <-all_results$S[2:13]

StdEstimates_S <- rownames_to_column(StdEstimates_S, 
                                       var="Targets")
        
S_StdEst_Factors <- filter(StdEstimates_S, 
                             Targets %in% c("EFA_MR1_Pol_S", 
                                            "EFA_MR2_Dog_S", 
                                            "EFA_MR3_Rel_S"))

domain_levels_S <- list(
                   "Goal-Directedness"="stdbeta:goal_dir",
                   "Impulsivity"="stdbeta:impuls",

                    "Reward Sensitivity"="stdbeta:rew_sens",
                    "Sensation Seeking"="stdbeta:senseseek",

                    "Emotional Control"="stdbeta:emot_cont",
                    "Agreeableness"="stdbeta:agree", 
                    "Ethical Risk-Taking"="stdbeta:ethic_risk",
                       "Risk Perception"="stdbeta:risk_per", 

                    "Eating Control"="stdbeta:eat_cont",
                    "Mindfulness"="stdbeta:mindful",
                    "Financial Risk-Taking"="stdbeta:fin_risk",
                    "Social Risk-Taking"="stdbeta:soc_risk"
                   )




domain_levels_reverse_S = names(domain_levels_S)
names(domain_levels_reverse_S)=domain_levels_S

S_StdEst_Factors_tidy <- gather(S_StdEst_Factors, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
    mutate(PsychDomain = factor(PsychDomain,
                                levels=domain_levels_S)) %>%
    mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_S))



#Outcome: 3 factors


p = ggplot(data = S_StdEst_Factors_tidy, 
       aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 

hlines <- c(-.4, -.35, -.30, -.25, -.2, -.15, -.10, -.05, 0, .05, .10, .15, .2, .25, .3)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 
p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
   labs(x = "Outcome", y = "Standardized Estimate", fill = "Personality Domain") + 
   scale_fill_viridis(discrete= TRUE, option="D", begin=.1)  +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
   scale_x_discrete(limits=c("EFA_MR1_Pol_S", "EFA_MR2_Dog_S", "EFA_MR3_Rel_S"), 
                   labels = c("Political Conservatism", "Dogmatism", "Religiosity")) +
  theme(plot.margin = margin(1.5, 0, 0, .5, "cm")) 

print(p)
dev.copy(pdf,sprintf("stdbeta_Factors_Survey.pdf"), width=10, height=8)
dev.off()
svg(file = "stdbeta_Factors_Survey.svg", width=10, height=8)
  print(p)
  dev.off()

##Outcome: Political Conservatism - A
S_StdEst_Political <- filter(StdEstimates_S, 
                             Targets %in% c("SECS_Social_S", 
                                            "SECS_Econ_S", 
                                            "Nationalism_S",
                                            "Patriotism scale_S",
                                            "Authoritarianism_S",
                                            "System Justification_S",
                                            "SDO_S",
                                            "Progroup_Actions_S"
                                            ))

S_StdEst_Political_tidy <- gather(S_StdEst_Political, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
              mutate(PsychDomain = factor(PsychDomain,
                                          levels=domain_levels_S)) %>%
              mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_S))



p = ggplot(data = S_StdEst_Political_tidy, 
   aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 

hlines <- c(-.4, -.35, -.3, -.25,-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 

p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
   labs(x = "Outcome", y = "Standardized Estimate", fill = "Personality Domain") + 
   scale_fill_viridis(discrete= TRUE, option="D", begin=.1)  +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
   scale_x_discrete(limits=c("SECS_Social_S", 
                                            "SECS_Econ_S", 
                                            "Nationalism_S",
                                            "Patriotism scale_S" ), 
                  labels = c("Social Conservatism", "Economic Conservatism", "Nationalism", "Patriotism" )) +
   theme(plot.margin = margin(1.5, 0, 0, .5, "cm"), legend.position = "right") 

print(p)
dev.copy(pdf,sprintf("stdbeta_PoliticalA_Survey.pdf"), width=10, height=8)
dev.off()
svg(file = "stdbeta_PoliticalA_Survey.svg", width=10, height=8)
  print(p)
  dev.off()


##Outcome: Political Conservatism - B

p = ggplot(data = S_StdEst_Political_tidy, 
   aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 

hlines <- c(-.4, -.35, -.3, -.25,-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 
p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
   labs(x = "Outcome", y = "Standardized Estimate", fill = "Personality Domain") + 
   scale_fill_viridis(discrete= TRUE, option="D", begin=.1)  +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
   scale_x_discrete(limits=c("Authoritarianism_S",
                                            "System Justification_S",
                                            "SDO_S",
                                            "Progroup_Actions_S"), 
                  labels = c("Authoritarianism", "System Justification", "Social Dominance", "Extreme Pro-Group Attitudes")) +
   theme(plot.margin = margin(1.5, 0, 0, .5, "cm"), legend.position = "right") 

print(p)
dev.copy(pdf,sprintf("stdbeta_PoliticalB_Survey.pdf"), width=10, height=8)
dev.off()
svg(file = "stdbeta_PoliticalB_Survey.svg", width=10, height=8)
  print(p)
  dev.off()



# Outcome: dogmatism 

S_StdEst_Dogmatism <- filter(StdEstimates_S, 
                             Targets %in% c("Dogmatism_S", 
                                            "IH1_S", 
                                            "IH2_S",
                                            "IH3_S",
                                            "IH4_S" ))

S_StdEst_Dogmatism_tidy <- gather(S_StdEst_Dogmatism, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
              mutate(PsychDomain = factor(PsychDomain,
                                          levels=domain_levels_S)) %>%
              mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_S))



p = ggplot(data = S_StdEst_Dogmatism_tidy, 
   aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 

hlines <- c( -.3, -.25,-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2, .25,.3, .35, .4)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5, 5.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 

p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
   labs(x = "Outcome", y = "Standardized Estimate", fill = "Personality Domain") + 
   scale_fill_viridis(discrete= TRUE, option="D", begin=.1)  +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
   scale_x_discrete(limits=c("Dogmatism_S", 
                                            "IH1_S", 
                                            "IH2_S",
                                            "IH3_S",
                                            "IH4_S"  ), 
                  labels = c("Dogmatism", "IH1: Intellectual Ego", "IH2: Belief Updating", "IH3: Respect for Others", "IH4: Overconfidence")) +
   theme(plot.margin = margin(1.5, 0, 0, .5, "cm"), legend.position = "right") 

print(p)
dev.copy(pdf,sprintf("stdbeta_Dogmatism_Survey.pdf"), width=10, height=8)
dev.off()
svg(file = "stdbeta_Dogmatism_Survey.svg", width=10, height=8)
  print(p)
  dev.off()


# Outcome: religiosity


S_StdEst_Religiosity <- filter(StdEstimates_S, 
                             Targets %in% c("Religious_Attendance_S",
                                            "Religion_Importance_S",
                                            "Religious_Prayer_S"))

S_StdEst_Religiosity_tidy <- gather(S_StdEst_Religiosity, -"Targets", 
                             key="PsychDomain", 
                             value="StdEstimate") %>%
              mutate(PsychDomain = factor(PsychDomain,
                                          levels=domain_levels_S)) %>%
              mutate(PsychDomainNames = recode_factor(PsychDomain,!!!domain_levels_reverse_S))



p = ggplot(data = S_StdEst_Religiosity_tidy, 
   aes(x=Targets, y=StdEstimate, fill = PsychDomainNames)) 

hlines <- c( -.3, -.25,-.2, -.15, -.10, -.05, 0, .05, .10, .15, .2, .25)
for (h in hlines){
         p <- p + geom_hline(yintercept=h, color="gray88")
} 

vlines <- c(1.5,2.5, 3.5, 4.5, 5.5)
for (v in vlines){
         p <- p + geom_vline(xintercept=v, color="gray88")
} 

p<- p +  geom_bar(stat="identity", position = position_dodge()) + 
   labs(x = "Outcome", y = "Standardized Estimate", fill = "Personality Domain") + 
   scale_fill_viridis(discrete= TRUE, option="D", begin=.1)  +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
   scale_x_discrete(limits=c( "Religious_Attendance_S",
                                            "Religion_Importance_S",
                                            "Religious_Prayer_S"), 
                  labels =  c("Religious Attendance",  "Religion Importance", "Religious Prayer")) +
   theme(plot.margin = margin(1.5, 0, 0, .5, "cm"), legend.position = "right") 

print(p)
dev.copy(pdf,sprintf("stdbeta_Religiosity_Survey.pdf"), width=10, height=8)
dev.off()
svg(file = "stdbeta_Religiosity_Survey.svg", width=10, height=8)
  print(p)
  dev.off()




```



#### Bayesian regressions


```{r BFsetup}

# Set up Bayesian Regressions
# use listwise deletion because BF computation requires that all analyses include the same dataset

BData <- data.frame(
  EFA_MR3_Rel = combined_data$Ideol_Factor3_Religion,
  EFA_MR2_Dog = combined_data$Ideol_Factor2_Dogmatism,
  EFA_MR1_Pol = combined_data$Ideol_Factor1_Politics,
  
  Age = combined_data$Age,
  Edu = combined_data$Edu,
  Gender = combined_data$Gender,
  Income = combined_data$Income,
  
  agree = combined_data$agree,
  fin_risk = combined_data$fin_risk,
  soc_risk = combined_data$soc_risk,
  ethic_risk = combined_data$ethic_risk,
  eat_cont= combined_data$eat_cont,
  risk_per = combined_data$risk_per,
  goal_dir = combined_data$goal_dir,
  rew_sens = combined_data$rew_sens,
  emot_cont = combined_data$emot_cont,
  impuls = combined_data$impuls,
  mindful = combined_data$mindful,
  senseseek = combined_data$senseseek,
  
  SpeededIP = combined_data$SpeededIP,
  StrategicIP = combined_data$StrategicIP,
  Discount = combined_data$Discount,
  Caution = combined_data$Caution,
  Percep = combined_data$Percep
)

BData <- na.omit(BData)

```



```{r runBayesReg}
targetvars_bayes <- c("EFA_MR1_Pol","EFA_MR2_Dog","EFA_MR3_Rel")
# reuse varsets from earlier regressions

# run BayesFactor on regression models

bfreg_bayes <- list()
models_bayes <- list()

for (v in names(varsets)){
  for (t in targetvars_bayes){
    modelname <- paste(t,v,sep='_')
    models_bayes[[modelname]] <- sprintf('%s ~ %s',t,varsets[v])
    # cat(sprintf('running model: %s\n',modelname))
    bfreg_bayes[[modelname]] <- lmBF(data=BData,
                              formula=formula(models_bayes[[modelname]]))
  }
}



 
```

Perform variable selection for each target using Bayesian model averaging, to get best model for each variable.

```{r bmaVariableSelection}
pthresh <- 75  #posterior probability threshold

bma_models <- list()
bma_vars <- list()

for (t in targetvars_bayes){
  y = BData %>% 
     dplyr::select(!!t)
  ix <- which(names(BData) %in% targetvars_bayes)
  BData_subset <- BData[,-ix]
  #BData_subset <- BData %>% 
  #dplyr::select(-!!targetvars_bayes)
  cat(sprintf('\ntarget variable: %s\n',t))
  bma_models[[t]] <- bic.glm(BData_subset,as.matrix(y),'gaussian')
  cat(sprintf('BMA results (variables with p> %0.3f: %s\n',pthresh/100,t))
  bma_vars[[t]] <- bma_models[[t]]$probne0[bma_models[[t]]$probne0>pthresh]
  print(bma_vars[[t]])
  # compute the model for this set
  modelname <- paste(t,"bma",sep='_')
  models_bayes[[modelname]] <- sprintf('%s ~ %s',t,paste(names(bma_vars[[t]]),collapse=" + "))
  bfreg_bayes[[modelname]] <- lmBF(data=BData,
                              formula=formula(models_bayes[[modelname]]))
  
}
```

Compute BF for all variables vs. demographic only within each target var

```{r bfComputation}
bf_vs_demog = list()
BayesFactorResults=c()
targetvar_labels = list("EFA_MR1_Pol" = "Political",
                        "EFA_MR2_Dog" = "Dogmatism",
                        "EFA_MR3_Rel"="Religiosity")
varset_names = list("D" = "Demographics",
                    "C" = "Cognitive Tasks",
                    "S" = "Personality Surveys",
                    "CS" = "Cognitive Tasks + Personality Surveys",
                    "CSD" = "Combined",
                    "bma" = "Best")

for (t in targetvars_bayes){
  for (v in names(varset_names)){
    if (v != "D"){
      modelname <- paste(t,v,sep='_')
      bf_vs_demog[[modelname]] <- bfreg_bayes[[modelname]]/bfreg_bayes[[paste(t,"D",sep='_')]]
      
    }
  }
}

```

Convert BF results into data frame for plotting.

```{r bfResults}
BayesFactor_Results <- data.frame(Model=character(),
                     Outcome=character(),
                     BF_10=double(),
                     BF_Dem=double())

for (t in targetvars_bayes){
  for (v in names(varset_names)){
    modelname <- paste(t,v,sep='_')
    if (!v=="D"){
      bf_d <- extractBF(bf_vs_demog[[modelname]])$bf
    } else {
      bf_d <- NA
    }
    info <- data.frame(Model=varset_names[[v]],
                       Outcome=targetvar_labels[[t]],
                       BF_10=extractBF(bfreg_bayes[[modelname]])$bf,
                       BF_Dem=bf_d)
    BayesFactor_Results <- rbind(BayesFactor_Results,info)
  }
}


```

#### Figure 2 and Supplement - Bayes Factors

Bayes factor figures

```{r mkBfFigures}


BF_Data <- filter(BayesFactor_Results, Outcome %in% c("Political", "Dogmatism", "Religiosity"))                
BF_Data$Model <- factor(BF_Data$Model, levels = c("Best", "Combined","Cognitive Tasks + Personality Surveys", "Personality Surveys", "Cognitive Tasks", "Demographics"), ordered=TRUE)
    
p <- ggplot(data=BF_Data, aes(y=BF_10, x=Outcome, fill=Model)) + 
  geom_bar(stat="identity", position="dodge") + 
  labs(x = "Outcome", y = "Bayes Factor") + 
  scale_x_discrete(limits=c("Political", "Dogmatism", "Religiosity")) +
  scale_fill_viridis(discrete= TRUE) +
  scale_y_log10()+ 
  coord_flip() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 


hlines <- c(1,3,10,30,100,.3333,.1,.01,.03)

for (h in hlines){
  p <- p + geom_hline(yintercept=h, linetype="dashed", color="grey")
}
print(p)

dev.copy(pdf,sprintf("bf10.pdf"), width=12, height=8)
dev.off()
  svg(file = "bf10.svg", width=12, height=8)
  print(p)
  dev.off()

# supplement figure - BF_Dem
p <- ggplot(data=BF_Data %>% filter(Model !="Demographics"), aes(y=BF_Dem, x=Outcome, fill=Model)) + 
  geom_bar(stat="identity", position="dodge") + 
  labs(x = "Outcome", y = "Bayes Factor") + 
  scale_x_discrete(limits=c("Political", "Dogmatism", "Religiosity")) +
  scale_fill_viridis(discrete=TRUE, end = .85) +
  scale_color_discrete( guide = guide_legend(reverse=TRUE) )+
  scale_y_log10() + 
  coord_flip() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

hlines <- c(1,3,10,30,100,.3333,.1,.01,.03)

for (h in hlines){
  p <- p + geom_hline(yintercept=h, linetype="dashed", color="grey")
}
print(p)

dev.copy(pdf,sprintf("bf_dem.pdf"), width=12, height=8)
dev.off()
 svg(file = "bfdem.svg", width=12, height=8)
  print(p)
  dev.off()


```

### Figure 3 and Supplementary Figures - BMA Model Selection Plots

```{r bmaPlots}

# plot models used in BMA


# outcome: political conservatism factor
PoliticalConservatismBMA <- imageplot.bma(bma_models$EFA_MR1_Pol, order = "probne0", color = c("red", "blue", "white"))
dev.copy(pdf,sprintf("PoliticalConservatismBMA.pdf"), width=10, height=8)
dev.off()

# outcome: dogmatism factor
DogmatismBMA <- imageplot.bma(bma_models$EFA_MR2_Dog, order = "probne0", color = c("red", "blue", "white"))
dev.copy(pdf,sprintf("DogmatismBMA.pdf"), width=10, height=8)
dev.off()

# outcome: religiosity factor
ReligiosityBMA <- imageplot.bma(bma_models$EFA_MR3_Rel, order = "probne0", color = c("red", "blue", "white"))
dev.copy(pdf,sprintf("ReligiosityBMA.pdf"), width=10, height=8)
dev.off()





# plot posterior probabilities of each variable not being included in model
probne0_factors <- data.frame(vars = names(bma_models$EFA_MR1_Pol$probne0),
                        probne0_pol = bma_models$EFA_MR1_Pol$probne0,
                        probne0_dog = bma_models$EFA_MR2_Dog$probne0,
                        probne0_rel = bma_models$EFA_MR3_Rel$probne0)

prob_fact <- probne0_factors %>% gather(key = "Factor", 
                 value = "Prob0", c("probne0_pol", "probne0_dog", "probne0_rel"), 
                 factor_key=TRUE)

variable_order <- c("Caution", "Discount", "Percep", "SpeededIP", "StrategicIP", "goal_dir","impuls","rew_sens", "senseseek", "emot_cont", "agree", "ethic_risk", "risk_per", "eat_cont", "mindful", "fin_risk", "soc_risk", "Age", "Edu", "Gender", "Income")

label_order <-  c("Caution", "Temporal Discounting", "Perceptual Processing", "Speed of Evidence Accum.", "Strategic IP", "Goal-Directedness",  "Impulsivity", "Reward Sensitivity", "Sensation Seeking", "Emotional Control", "Agreeableness", "Ethical Risk-Taking", "Risk Perception", "Eating Control", "Mindfulness", "Financial Risk-Taking", "Social Risk-Taking","Age", "Education", "Gender", "Income")

color_blind_friendly <-  c("#000000", "#E69F00", "#56B4E9")

bma_prob_plot <- ggplot(data=prob_fact, aes(x=vars, y=Prob0, group=Factor, color=Factor)) + 
              geom_segment(aes(xend=vars, yend=0), color="grey50") +
              geom_point(size=2.8) + 
              geom_hline(yintercept=75, linetype="dashed", color="red", alpha=.5) +  
              scale_x_discrete(limits=rev(variable_order), labels = rev(label_order)) + 
              scale_color_manual(values=color_blind_friendly, labels = c("Political Conservatism", "Dogmatism", "Religiosity")) +
              xlab(bquote("Variables (" ~ beta[i] ~ ")")) + ylab(bquote("Pr(" ~ beta[i] != "0" ~ "|" ~ "D)")) +
              ggtitle("Bayesian Model Averaging Posterior Probability") +
              coord_flip() +
              theme(panel.grid.major = element_line(color = "gray90"), panel.grid.minor = element_line(color = "gray90"), legend.position = "right", legend.title=element_blank())

print(bma_prob_plot)
dev.copy(pdf,sprintf("postprob_BMA_Factors.pdf"), width=12, height=8)
dev.off()
svg(file = "postprob_BMA_Factors.svg", width=12, height=8)
print(bma_prob_plot)
dev.off()


```



#### Figure CV:  Crossvalidation

```{r mkFigureCV, fig.width=10, fig.height=8}
       
## import data of cross validation prediction        

## Outcomes: all variables, including factors

CV_Data <- read.csv("predictions_R2.csv")

CV_Data$Outcome <- gsub(".*_", "", CV_Data$Outcome)

outcomes <- c("Social Conservatism",  "Economic Conservatism", "Nationalism", "Patriotism", "Authoritarianism", "Social Dominance Orientation", "System Justification", "Ingroup Attitude", "Political Factor", "Dogmatism", "IH1", "IH2", "IH3", "IH4", "Dogmatism Factor", "Religious Attendance", "Religion Importance", "Religious Prayer", "Religiosity Factor" )


CV_All <- CV_Data %>% filter(Outcome %in% outcomes) %>% filter(prediction_type %in% c( "shuffled95", "CV"))

CVR2_all <- ggplot(data=CV_All, aes(x=Outcome, y= full_ontology_oblimin, fill=prediction_type)) +
  geom_bar(stat="identity", position="identity") + 
  scale_fill_manual(values=c("#2166AC", "#4393C3"), name = "Prediction Type") +  
  scale_x_discrete(limits=outcomes) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("R2 (cross-validated)") +
  theme(plot.margin = margin(0, 0, 0, 1, "cm")) 

print(CVR2_all)        
dev.copy(pdf,sprintf("R2_crossvalidation.pdf",label), width=10, height=8)
dev.off()
svg(file = "R2_crossvalidation.svg", width=10, height=8)
print(CVR2_all)
dev.off()
  
  
## outcomes: only 3 factors

outcomes_factors <- c( "Political Factor", "Dogmatism Factor", "Religiosity Factor" )

CV_Factors <- CV_Data %>% filter(Outcome %in% outcomes_factors) %>% filter(prediction_type %in% c("CV"))

CVR2_factors <- ggplot(data=CV_Factors, aes(x=Outcome, y= full_ontology_oblimin, fill=prediction_type)) +
  geom_bar(stat="identity", position="identity") + 
  scale_fill_manual(values=c("#2166AC")) +  
  scale_x_discrete(limits=outcomes_factors, labels = c("Political Conservatism", "Dogmatism", "Religiosity")) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("R2 (cross-validated)") +
  theme(plot.margin = margin(0, 0, 0, 1, "cm"), legend.position = "none") 

print(CVR2_factors)        
dev.copy(pdf,sprintf("R2_crossvalidation_factors.pdf",label), width=5, height=8)
dev.off()
svg(file = "R2_crossvalidation_factors.svg", width=5, height=8)
print(CVR2_factors)
dev.off()

## Outcomes: all outcomes except factors
outcomes_nofactors <- c("Social Conservatism",  "Economic Conservatism", "Nationalism", "Patriotism", "Authoritarianism", "Social Dominance Orientation", "System Justification", "Ingroup Attitude",  "Dogmatism", "IH1", "IH2", "IH3", "IH4", "Religious Attendance", "Religion Importance", "Religious Prayer")


CV_Orientations <- CV_Data %>% filter(Outcome %in% outcomes_nofactors) %>% filter(prediction_type %in% c("CV"))

CVR2_orient <- ggplot(data=CV_Orientations, aes(x=Outcome, y= full_ontology_oblimin, fill=prediction_type)) +
  geom_bar(stat="identity", position="identity") + 
  scale_fill_manual(values=c("#2166AC", "#4393C3"), name = "Prediction Type") +  
  scale_x_discrete(limits=outcomes_nofactors) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("R2 (cross-validated)") +
  theme(plot.margin = margin(0, 0, 0, 1, "cm"), legend.position = "none") 
print(CVR2_orient)        
dev.copy(pdf,sprintf("R2_crossvalidation_orientations.pdf",label), width=10, height=8)
dev.off()
svg(file = "R2_crossvalidation_orientations.svg", width=10, height=8)
print(CVR2_orient)
dev.off()

```




#### Figure Correlation Matrix  


Print correlogram of ideological orientations x 17 ontology for supplementary figure

```{r correlogramIdeolOntolFigure, fig.width=8, fig.height=10}
## Correlogram of ideological orientations x ontology - Supplementary Figure

data_subset <- combined_data[-c(38)]

Ontols <- data_subset[7:23]
Ideols <- data_subset[24:40]

M2 <- cor(y = Ideols, x=Ontols, use= "pairwise.complete.obs")

corrplot(M2, method="color", 
         addCoef.col = "black", number.cex = .75, 
         tl.col = "black", cl.pos = "b", 
         cl.length=4, insig="blank")

# save figure to file

dev.copy(pdf,"correlogram_ideol_ontol.pdf", width=12, height=12)
dev.off()



```




Print correlogram of 17 ontology for supplementary figure

```{r correlogramOntolFigure, fig.width=8, fig.height=10}
## Correlogram of ideological orientations x ontology - Supplementary Figure

data_subset <- combined_data[-c(38)]

Ontols <- data_subset[7:23]

M2 <- cor(x = Ontols, y=Ontols, use= "pairwise.complete.obs")

corrplot(M2, method="color", 
         addCoef.col = "black", number.cex = .75, 
         tl.col = "black", cl.pos = "b", 
         cl.length=4, insig="blank")

# save figure to file

dev.copy(pdf,"correlogram_ontol.pdf", width=12, height=12)
dev.off()



```

