rm(list=ls(all=TRUE))  # Clear all variables
graphics.off()  # Clear all graphics

# Load a package you'll need for Bayesian modeling and one with the data
library(brms)

# Load the relevant packages
library(ggplot2)
library(reshape2)
library(scales)
library(dplyr)
library(stringr)
library(blavaan)
library(semPlot)

# Turning off scientific notation will facilitate interpretting the results later
options(scipen=999)


get_prior(formula=t1_t2_change~0+Condition+PerspectiveTaking+Openness,
      data=S1D,family=gaussian(link="identity"))

#Run an initial model of the change in intergroup attitude from t1 to t2 with condition, perspective taking,
# and openness to experience as independent variables
M1<-brm(t1_t2_change~0+Condition+PerspectiveTaking+Openness,
        data = S1D,family = gaussian(link = "identity"),
        prior = c(set_prior("normal(.1,1)", class = "b", coef = "ConditionControl"),
                  set_prior("normal(.2,1)", class = "b", coef = "ConditionExperimental1"),
                  set_prior("normal(.3,1)", class = "b", coef = "ConditionExperimental2"),
                  set_prior("normal(.3,.3)", class = "b", coef = "PerspectiveTaking"),
                  set_prior("normal(.3,.3)", class = "b", coef = "Openness"),
                  set_prior("normal(3,2)", class = "sigma")),
        warmup=1000,iter=5000,chains=3,
        control=list(adapt_delta=0.95))
        
print(M1)


# Assess how well the model did
dev.new()
plot(M1)

# Visualize the Bayesian model
dev.new()
plot(marginal_effects(M1), points = TRUE)

#Extract posterior distributions
M1Post<-posterior_samples(M1)

# Compute posterior probability that participants in condition 2 will have a greater intergroup attitude than participants in condition 1
pM1ConditionExperimental1<-length(M1Post$`b_ConditionExperimental1`[M1Post$`b_ConditionExperimental1` > M1Post$`b_ConditionControl`])/length(M1Post$`b_ConditionExperimental1`)
cat("Model 1 probability that participants in experimental condition 1 show a greater improvement in intergroup attitude than participants in the control condition from time 1 to time 2 = ", pConditionExperimental1, "\n")

pM1ConditionExperimental2<-length(M1Post$`b_ConditionExperimental2`[M1Post$`b_ConditionExperimental2` > M1Post$`b_ConditionControl`])/length(M1Post$`b_ConditionExperimental2`)
cat("Model 1 probability that participants in experimental condition 2 show a greater improvement in intergroup attitude than participants in the control condition from time 1 to time 2 = ", pConditionExperimental2, "\n")

pM1ConditionExperimental2<-length(M1Post$`b_ConditionExperimental2`[M1Post$`b_ConditionExperimental2` > M1Post$`b_ConditionExperimental1`])/length(M1Post$`b_ConditionExperimental2`)
cat("Model 1 probability that participants in experimental condition 2 show a greater improvement in intergroup attitude than participant in experimental condition 1 from time 1 to time 2 = ", pConditionExperimental2, "\n")

pM1Openness<-length(M1Post$`b_Openness`[M1Post$`b_Openness` < 0])/length(M1Post$`b_Openness`)
cat("Model 1 probability that openness is less than 0 in the posterior distribution = ", pOpenness, "\n")

pM1PerspectiveTaking<-length(M1Post$`b_PerspectiveTaking`[M1Post$`b_PerspectiveTaking` > 0])/length(M1Post$`b_PerspectiveTaking`)
cat("Model 1 probability that Perspective Taking is greater than 0 in the posterior distribution = ", pPerspectiveTaking, "\n")

#Try a simplified model that includes only condition as an IV
get_prior(formula=t1_t2_change~0+Condition,
          data=S1D,family=gaussian(link="identity"))

M2<-brm(t1_t2_change~0+Condition,
        data = S1D,family = gaussian(link = "identity"),
        prior = c(set_prior("normal(.1,.5)", class = "b", coef = "ConditionControl"),
                  set_prior("normal(.2,.5)", class = "b", coef = "ConditionExperimental1"),
                  set_prior("normal(.3,.5)", class = "b", coef = "ConditionExperimental2"),
                  set_prior("normal(3,2)", class = "sigma")),
        warmup=1000,iter=5000,chains=3,
        control=list(adapt_delta=0.95))

print(M2)

# Assess how well the model did
dev.new()
plot(M2)

# Visualize the Bayesian model
dev.new()
plot(marginal_effects(M2), points = TRUE)

#Let's compare Models 1 and 2 to see which better fits the data
loo_compare(loo(M1),loo(M2))

#Model 2 looks only slightly better by loo_compare, but model 1 includes more of our variables. M1 is likely better than M2
#Let's include a separate model to identify whether Allport's 4 criteria might also have individual effects

#Model 3 get priors
get_prior(formula=t1_t2_change~0+Condition+PerspectiveTaking+Openness+(Equals+Interests+Goals),
          data=S1D,family=gaussian(link="identity"))

M3<-brm(t1_t2_change~0+Condition+PerspectiveTaking+Openness+Equals+Interests+Goals,
        data = S1D,family = gaussian(link = "identity"),
        prior = c(set_prior("normal(.1,.5)", class = "b", coef = "ConditionControl"),
                  set_prior("normal(.2,.5)", class = "b", coef = "ConditionExperimental1"),
                  set_prior("normal(.3,.5)", class = "b", coef = "ConditionExperimental2"),
                  set_prior("normal(.3,.3)", class = "b", coef = "PerspectiveTaking"),
                  set_prior("normal(.3,.3)", class = "b", coef = "Openness"),
                  set_prior("normal(1,.5)", class = "b", coef = "Equals"),
                  set_prior("normal(1,.5)", class = "b", coef = "Interests"),
                  set_prior("normal(1,.5)", class = "b", coef = "Goals"),
                  set_prior("normal(3,2)", class = "sigma")),
        warmup=1000,iter=5000,chains=3,
        control=list(adapt_delta=0.95))

print(M3)

# Assess how well the model did
dev.new()
plot(M3)

# Visualize the Bayesian model
dev.new()
plot(conditional_effects(M3), points = TRUE)

#Let's compare Models 1 and 3 to see which better fits the data
loo_compare(loo(M1),loo(M3))

#Since Model 3 contains Allport's conditions and comes out slightly better in loo_compare, let's go with model 3 and extract
# some posterior distributions to explore effects
M3Post<-posterior_samples(M3)

# Compute posterior probability that participants in condition 2 will have a greater intergroup attitude than participants in condition 1
pM3ConditionExperimental1versusControl<-length(M3Post$`b_ConditionExperimental1`[M3Post$`b_ConditionExperimental1` > M3Post$`b_ConditionControl`])/length(M3Post$`b_ConditionExperimental1`)
cat("Model 3 probability that participants in experimental condition 1 show a greater improvement in intergroup attitude than participants in the control condition from time 1 to time 2 = ", pM3ConditionExperimental1versusControl, "\n")

pM3ConditionExperimental2versusControl<-length(M3Post$`b_ConditionExperimental2`[M3Post$`b_ConditionExperimental2` > M3Post$`b_ConditionControl`])/length(M3Post$`b_ConditionExperimental2`)
cat("Model 3 probability that participants in experimental condition 2 show a greater improvement in intergroup attitude than participants in the control condition from time 1 to time 2 = ", pM3ConditionExperimental2versusControl, "\n")

pM3ConditionExperimental2versusExperimental1<-length(M3Post$`b_ConditionExperimental2`[M3Post$`b_ConditionExperimental2` > M3Post$`b_ConditionExperimental1`])/length(M1Post$`b_ConditionExperimental2`)
cat("Model 3 probability that participants in experimental condition 2 show a greater improvement in intergroup attitude than participant in experimental condition 1 from time 1 to time 2 = ", pM3ConditionExperimental2versusExperimental1, "\n")

pM3ConditionControlgreater0<-length(M3Post$`b_ConditionControl`[M3Post$`b_ConditionControl` > 0])/length(M3Post$`b_ConditionControl`)
cat("Model 3 probability that the control condition has an effect greater than 0 in the posterior distribution = ", pM3ConditionControlgreater0, "\n")

pM3ConditionControlgreaterprior<-length(M3Post$`b_ConditionExperimental1`[M3Post$`b_ConditionControl` > .1])/length(M3Post$`b_ConditionControl`)
cat("Model 3 probability that the control condition has an effect greater than .1 in the posterior distribution = ", pM3ConditionControlgreaterprior, "\n")

pM3ConditionExperimental1greater0<-length(M3Post$`b_ConditionExperimental1`[M3Post$`b_ConditionExperimental1` > 0])/length(M3Post$`b_ConditionExperimental1`)
cat("Model 3 probability that experimental condition 1 has an effect greater than 0 in the posterior distribution = ", pM3ConditionExperimental1greater0, "\n")

pM3ConditionExperimental1greaterprior<-length(M3Post$`b_ConditionExperimental1`[M3Post$`b_ConditionExperimental1` > .2])/length(M3Post$`b_ConditionExperimental1`)
cat("Model 3 probability that experimental condition 1 has an effect greater than .2 in the posterior distribution = ", pM3ConditionExperimental1greaterprior, "\n")

pM3ConditionExperimental2greater0<-length(M3Post$`b_ConditionExperimental2`[M3Post$`b_ConditionExperimental2` > 0])/length(M3Post$`b_ConditionExperimental2`)
cat("Model 3 probability that experimental condition 2 has an effect greater than 0 in the posterior distribution = ", pM3ConditionExperimental2greater0, "\n")

pM3ConditionExperimental2greaterprior<-length(M3Post$`b_ConditionExperimental2`[M3Post$`b_ConditionExperimental2` > .3])/length(M3Post$`b_ConditionExperimental2`)
cat("Model 3 probability that experimental condition 2 has an effect greater than .3 in the posterior distribution = ", pM3ConditionExperimental2greaterprior, "\n")

pM3Opennessgreater0<-length(M3Post$`b_Openness`[M3Post$`b_Openness` > 0])/length(M3Post$`b_Openness`)
cat("Model 3 probability that openness is greater than 0 in the posterior distribution = ", pM3Opennessgreater0, "\n")

pM3Opennessgreaterprior<-length(M3Post$`b_Openness`[M3Post$`b_Openness` > .3])/length(M3Post$`b_Openness`)
cat("Model 3 probability that openness is greater than .3 in the posterior distribution = ", pM3Opennessgreaterprior, "\n")

pM3PerspectiveTakinggreater0<-length(M3Post$`b_PerspectiveTaking`[M3Post$`b_PerspectiveTaking` > 0])/length(M3Post$`b_PerspectiveTaking`)
cat("Model 3 probability that Perspective Taking is greater than 0 in the posterior distribution = ", pM3PerspectiveTakinggreater0, "\n")

pM3PerspectiveTakinggreaterprior<-length(M3Post$`b_PerspectiveTaking`[M3Post$`b_PerspectiveTaking` > .3])/length(M3Post$`b_PerspectiveTaking`)
cat("Model 3 probability that Perspective Taking is greater than .3 in the posterior distribution = ", pM3PerspectiveTakinggreaterprior, "\n")

pM3Equalsgreater0<-length(M3Post$`b_Equals`[M3Post$`b_Equals` > 0])/length(M3Post$`b_Equals`)
cat("Model 3 probability that Equals is greater than 0 in the posterior distribution = ", pM3Equalsgreater0, "\n")

pM3Equalsgreaterprior<-length(M3Post$`b_Equals`[M3Post$`b_Equals` > 1])/length(M3Post$`b_Equals`)
cat("Model 3 probability that Equals is greater than 1 in the posterior distribution = ", pM3Equalsgreaterprior, "\n")

pM3Interestsgreater0<-length(M3Post$`b_Interests`[M3Post$`b_Interests` > 0])/length(M3Post$`b_Interests`)
cat("Model 3 probability that Interests is greater than 0 in the posterior distribution = ", pM3Interestsgreater0, "\n")

pM3Interestsgreaterprior<-length(M3Post$`b_Interests`[M3Post$`b_Interests` > 1])/length(M3Post$`b_Interests`)
cat("Model 3 probability that Interests is greater than 1 in the posterior distribution = ", pM3Interestsgreaterprior, "\n")

pM3Goalsgreater0<-length(M3Post$`b_Goals`[M3Post$`b_Goals` > 0])/length(M3Post$`b_Goals`)
cat("Model 3 probability that Goals is greater than 0 in the posterior distribution = ", pM3Goalsgreater0, "\n")

pM3Goalsgreaterprior<-length(M3Post$`b_Goals`[M3Post$`b_Goals` > 1])/length(M3Post$`b_Goals`)
cat("Model 3 probability that Goals is greater than 1 in the posterior distribution = ", pM3Goalsgreaterprior, "\n")






#Create a new model for the change from time 1 to time 3 using condition, perspective taking, openness, and Allport's 4 criteria
get_prior(formula=t1_t3_change~0+Condition+PerspectiveTaking+Openness+Equals+Interests+Goals,
          data=S1D,family=gaussian(link="identity"))

#Due to the time between measurements (30 days), any effects of the intervention should lessen (e.g., condition, Allport's 4 criteria)
#Model 3 suggested smaller effects from stable traits like openness and perspective-taking.
#But as time from the intervention increases, effects from these more stable traits should become more prominent. 
#As a result, I've reduced the priors for all intervention-related variables and set priors for the variables related to
#more stable traits at their initially anticipated settings
M4<-brm(t1_t3_change~0+Condition+PerspectiveTaking+Openness+Equals+Interests+Goals,
        data = S1D,family = gaussian(link = "identity"),
        prior = c(set_prior("normal(.05,.25)", class = "b", coef = "ConditionControl"),
                  set_prior("normal(.1,.25)", class = "b", coef = "ConditionExperimental1"),
                  set_prior("normal(.15,.25)", class = "b", coef = "ConditionExperimental2"),
                  set_prior("normal(.3,.3)", class = "b", coef = "PerspectiveTaking"),
                  set_prior("normal(.3,.3)", class = "b", coef = "Openness"),
                  set_prior("normal(.2,.25)", class = "b", coef = "Equals"),
                  set_prior("normal(.2,.25)", class = "b", coef = "Interests"),
                  set_prior("normal(.2,.25)", class = "b", coef = "Goals"),
                  set_prior("normal(3,2)", class = "sigma")),
        warmup=1000,iter=5000,chains=3,
        control=list(adapt_delta=0.95))

print(M4)

# Assess how well the model did
dev.new()
plot(M4)

# Visualize the Bayesian model
dev.new()
plot(marginal_effects(M4), points = TRUE)

#For comparison's sake, let's run a more simplified model that doesn't include Allport's criteria
get_prior(formula=t1_t3_change~0+Condition+PerspectiveTaking+Openness,
          data=S1D,family=gaussian(link="identity"))

M5<-brm(t1_t3_change~0+Condition+PerspectiveTaking+Openness,
        data = S1D,family = gaussian(link = "identity"),
        prior = c(set_prior("normal(.05,.25)", class = "b", coef = "ConditionControl"),
                  set_prior("normal(.1,.25)", class = "b", coef = "ConditionExperimental1"),
                  set_prior("normal(.15,.25)", class = "b", coef = "ConditionExperimental2"),
                  set_prior("normal(.3,.25)", class = "b", coef = "PerspectiveTaking"),
                  set_prior("normal(.3,.25)", class = "b", coef = "Openness"),
                  set_prior("normal(3,2)", class = "sigma")),
        warmup=1000,iter=5000,chains=3,
        control=list(adapt_delta=0.95))

print(M5)

# Assess how well the model did
dev.new()
plot(M5)

# Visualize the Bayesian model
dev.new()
plot(marginal_effects(M5), points = TRUE)

# Since marginal_effects is deprecated, use conditional_effects to visualize the Bayesian model instead
dev.new()
plot(conditional_effects(M5), points = TRUE)

#Let's compare Models 4 and 5 to see which better fits the data
loo_compare(loo(M4),loo(M5))

#Since the loo_compare difference isn't too substantial here and Model 4 allows us to compare 
#longitudinal differences across all variables, #we'll go with model 4

#Extract posterior samples
M4Post<-posterior_samples(M4)


pM4ConditionControlgreater0<-length(M4Post$`b_ConditionControl`[M4Post$`b_ConditionControl` > 0])/length(M4Post$`b_ConditionControl`)
cat("Model 4 probability that the control condition has an effect greater than 0 in the posterior distribution = ", pM4ConditionControlgreater0, "\n")

pM4ConditionControlgreaterprior<-length(M4Post$`b_ConditionExperimental1`[M4Post$`b_ConditionControl` > .05])/length(M4Post$`b_ConditionControl`)
cat("Model 4 probability that the control condition has an effect greater than .05 in the posterior distribution = ", pM4ConditionControlgreaterprior, "\n")

pM4ConditionExperimental1greater0<-length(M4Post$`b_ConditionExperimental1`[M4Post$`b_ConditionExperimental1` > 0])/length(M4Post$`b_ConditionExperimental1`)
cat("Model 4 probability that experimental condition 1 has an effect greater than 0 in the posterior distribution = ", pM4ConditionExperimental1greater0, "\n")

pM4ConditionExperimental1greaterprior<-length(M4Post$`b_ConditionExperimental1`[M4Post$`b_ConditionExperimental1` > .1])/length(M4Post$`b_ConditionExperimental1`)
cat("Model 4 probability that experimental condition 1 has an effect greater than .1 in the posterior distribution = ", pM4ConditionExperimental1greaterprior, "\n")

pM4ConditionExperimental2greater0<-length(M4Post$`b_ConditionExperimental2`[M4Post$`b_ConditionExperimental2` > 0])/length(M4Post$`b_ConditionExperimental2`)
cat("Model 4 probability that experimental condition 2 has an effect greater than 0 in the posterior distribution = ", pM4ConditionExperimental2greater0, "\n")

pM4ConditionExperimental2greaterprior<-length(M4Post$`b_ConditionExperimental2`[M4Post$`b_ConditionExperimental2` > .15])/length(M4Post$`b_ConditionExperimental2`)
cat("Model 4 probability that experimental condition 2 has an effect greater than .15 in the posterior distribution = ", pM4ConditionExperimental2greaterprior, "\n")

pM4Opennessgreater0<-length(M4Post$`b_Openness`[M4Post$`b_Openness` > 0])/length(M4Post$`b_Openness`)
cat("Model 4 probability that openness is greater than 0 in the posterior distribution = ", pM4Opennessgreater0, "\n")

pM4Opennessgreaterprior<-length(M4Post$`b_Openness`[M4Post$`b_Openness` > .3])/length(M4Post$`b_Openness`)
cat("Model 4 probability that openness is greater than .3 in the posterior distribution = ", pM4Opennessgreaterprior, "\n")

pM4PerspectiveTakinggreater0<-length(M4Post$`b_PerspectiveTaking`[M4Post$`b_PerspectiveTaking` > 0])/length(M4Post$`b_PerspectiveTaking`)
cat("Model 4 probability that Perspective Taking is greater than 0 in the posterior distribution = ", pM4PerspectiveTakinggreater0, "\n")

pM4PerspectiveTakinggreaterprior<-length(M4Post$`b_PerspectiveTaking`[M4Post$`b_PerspectiveTaking` > .3])/length(M4Post$`b_PerspectiveTaking`)
cat("Model 4 probability that Perspective Taking is greater than .3 in the posterior distribution = ", pM4PerspectiveTakinggreaterprior, "\n")

pM4Equalsgreater0<-length(M4Post$`b_Equals`[M4Post$`b_Equals` > 0])/length(M4Post$`b_Equals`)
cat("Model 4 probability that Equals is greater than 0 in the posterior distribution = ", pM4Equalsgreater0, "\n")

pM4Equalsgreaterprior<-length(M4Post$`b_Equals`[M4Post$`b_Equals` > .2])/length(M4Post$`b_Equals`)
cat("Model 4 probability that Equals is greater than .2 in the posterior distribution = ", pM4Equalsgreaterprior, "\n")

pM4Interestsgreater0<-length(M4Post$`b_Interests`[M4Post$`b_Interests` > 0])/length(M4Post$`b_Interests`)
cat("Model 4 probability that Interests is greater than 0 in the posterior distribution = ", pM4Interestsgreater0, "\n")

pM4Interestsgreaterprior<-length(M4Post$`b_Interests`[M4Post$`b_Interests` > .2])/length(M4Post$`b_Interests`)
cat("Model 4 probability that Interests is greater than .2 in the posterior distribution = ", pM4Interestsgreaterprior, "\n")

pM4Goalsgreater0<-length(M4Post$`b_Goals`[M4Post$`b_Goals` > 0])/length(M4Post$`b_Goals`)
cat("Model 4 probability that Goals is greater than 0 in the posterior distribution = ", pM4Goalsgreater0, "\n")

pM4Goalsgreaterprior<-length(M4Post$`b_Goals`[M4Post$`b_Goals` > .2])/length(M4Post$`b_Goals`)
cat("Model 4 probability that Goals is greater than .2 in the posterior distribution = ", pM4Goalsgreaterprior, "\n")

pM4Supportivegreater0<-length(M4Post$`b_Supportive`[M4Post$`b_Supportive` > 0])/length(M4Post$`b_Supportive`)
cat("Model 4 probability that Supportive is greater than 0 in the posterior distribution = ", pM4Supportivegreater0, "\n")

pM4Supportivegreaterprior<-length(M4Post$`b_Supportive`[M4Post$`b_Supportive` > .2])/length(M4Post$`b_Supportive`)
cat("Model 3 probability that Supportive is greater than .2 in the posterior distribution = ", pM4Supportivegreaterprior, "\n")