#############################################################
####Religion and Resource Stress -- Analysis Script 01-2020####
#############################################################

library(Hmisc)
library(lme4)
library(lmerTest) 
library(ggplot2)
library(ggthemes)
library(phia) #also loads car package
library(plyr)
library(lm.beta)
lm.beta.lmer <- function(mod) { b <- fixef(mod)[-1]; sd.x <- apply(matrix(getME(mod,"X")[,-1]),2,sd); sd.y <- sd(getME(mod,"y")); b*sd.x/sd.y; }

#G = High gods
#SG = Superior gods
#MS = minor spirits
#Gor_SG = broader supernatural agents

#re-naming data files to shorter names (arbitrary-just for easier typing in R)
mor<-read.csv("DT-MoralHighGods_Long_EP.csv")
ang<-read.csv("DT-GodsAndResourceStress_Long.csv")


#####################################
##Analyses on Newly Coded Variables##
#####################################

##How many overlapping language families are there?
table(ang$language_family)

##Dummy-coding language families
ang$Afroasiatic<-ifelse(as.character(ang$language_family)=="Afro-Asiatic",1,0)
ang$Afroasiatic<-ifelse(is.na(ang$Afroasiatic)==T,0,ang$Afroasiatic)

ang$Austronesian<-ifelse(as.character(ang$language_family)=="Austronesian",1,0)
ang$Austronesian<-ifelse(is.na(ang$Austronesian)==T,0,ang$Austronesian)

ang$AtlanticCongo<-ifelse(as.character(ang$language_family)=="Atlantic-Congo",1,0)
ang$AtlanticCongo<-ifelse(is.na(ang$AtlanticCongo)==T,0,ang$AtlanticCongo)

ang$Indoeuropean<-ifelse(as.character(ang$language_family)=="Indo-European",1,0)
ang$Indoeuropean<-ifelse(is.na(ang$Indoeuropean)==T,0,ang$Indoeuropean)

ang$Austroasiatic<-ifelse(as.character(ang$language_family)=="Austroasiatic",1,0)
ang$Austroasiatic<-ifelse(is.na(ang$Austroasiatic)==T,0,ang$Austroasiatic)

ang$NuclearTransNewGuinea<-ifelse(as.character(ang$language_family)=="Nuclear Trans New Guinea",1,0)
ang$NuclearTransNewGuinea<-ifelse(is.na(ang$NuclearTransNewGuinea)==T,0,ang$NuclearTransNewGuinea)

ang$SinoTibetan<-ifelse(as.character(ang$language_family)=="Sino-Tibetan",1,0)
ang$SinoTibetan<-ifelse(is.na(ang$SinoTibetan)==T,0,ang$SinoTibetan)

ang$Turkic<-ifelse(as.character(ang$language_family)=="Turkic",1,0)
ang$Turkic<-ifelse(is.na(ang$Turkic)==T,0,ang$Turkic)

ang$Uralic<-ifelse(as.character(ang$language_family)=="Uralic",1,0)
ang$Uralic<-ifelse(is.na(ang$Uralic)==T,0,ang$Uralic)


##Dummy variables for separate stressors
#Famine (re-coded during restructuring as 1 in variable "ResourceStress")
ang$famine<-ifelse(as.character(ang$ResourceStress)=="1",1,0)
#Natural Hazards (re-coded during restructuring as 2 in variable "ResourceStress")
ang$naturalhazards<-ifelse(as.character(ang$ResourceStress)=="2",1,0)
#Chronic Scarcity (re-coded during restructuring as 3 in variable "ResourceStress")
ang$chron_scarcity<-ifelse(as.character(ang$ResourceStress)=="3",1,0)
#Abundance (re-coded during restructuring as 4 in variable "ResourceStress")
ang$abundance<-ifelse(as.character(ang$ResourceStress)=="4",1,0)

#Recoding resource stress
ang$stress2<-ifelse(ang$abundance==1,(0-ang$Stress),
                    ifelse(ang$abundance==0,ang$Stress,
                           NA))

#Dichotomized stress variable (primary analyses do not use this variable)
ang$stress3<-ifelse(ang$chron_scarcity==1&ang$stress2<(-.11),(-1),
                    ifelse(ang$chron_scarcity==1&ang$stress2>(-.11),1,
                           ifelse(ang$famine==1&ang$stress2<(-.11),(-1),
                                  ifelse(ang$famine==1&ang$stress2>(-.11),1,
                                         ang$stress2))))


ang_fam<-ang[ang$famine==1,]
ang_naturalhazards<-ang[ang$naturalhazards==1,]
ang_abundance<-ang[ang$abundance==1,]
ang_chron_scarcity<-ang[ang$chron_scarcity==1,]
ang_hurthelp<-ang[is.na(ang$stress2)==F&is.na(ang$G3_Help)==F&is.na(ang$G2_Harm)==F,]

##Creating datasets for each of the god variables (this will allow us to calculate R2 later)
ang2<-ang[is.na(ang$G1_Associated)==F,]
ang3<-ang[is.na(ang$G2_Harm)==F,]
ang4<-ang[is.na(ang$G3_Help)==F,]
ang5<-ang[is.na(ang$G4_Punitive)==F,]
ang6<-ang[is.na(ang$SG1_Associated)==F,]
ang7<-ang[is.na(ang$SG2_Harm)==F,]
ang8<-ang[is.na(ang$SG3_Help)==F,]
ang9<-ang[is.na(ang$SG4_Punitive)==F,]
ang10<-ang[is.na(ang$MS1_Associated)==F,]
ang11<-ang[is.na(ang$MS2_Harm)==F,]
ang12<-ang[is.na(ang$MS3_Help)==F,]
ang13<-ang[is.na(ang$MS4_Punitive)==F,]

#Empty model controlling for language family (anova) (this will allow us to calculate R2 later)
summary(lmer(stress2~(1|SCCS_ID),data=ang)) 
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang2)) 
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang3)) 
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang4)) 
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang5))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang_hurthelp))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+SinoTibetan+Uralic+(1|SCCS_ID),data=ang6))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang7))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+SinoTibetan+(1|SCCS_ID),data=ang8))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang9))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Turkic+NuclearTransNewGuinea+Austroasiatic+SinoTibetan+(1|SCCS_ID),data=ang10))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang11))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang12))
summary(lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang13))

#################################################################
###PART A: Central Analyses on Stressors and Religious Beliefs###
#################################################################

##For the analyses below, you can get the R squared by comparing...
#...the variance explained at the society level with the variance explained by language family alone

#Resource Stress and Gods Associated with Weather# Table 2, Model 1
summary(model1<-lmer(stress2~G1_Associated+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
lm.beta.lmer(model1)

##Secondary analyses looking at specific types of resource stress# Table S4, Model S1
#looking at chronic scarcity
summary(model1b<-lmer(stress2~G1_Associated*famine+G1_Associated*naturalhazards+G1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
#looking at natural hazards
summary(model1c<-lmer(stress2~G1_Associated*famine+G1_Associated*chron_scarcity+G1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
#Looking at abundance
summary(model1d<-lmer(stress2~G1_Associated*famine+G1_Associated*naturalhazards+G1_Associated*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
#Looking at famine 
summary(model1e<-lmer(stress2~G1_Associated*chron_scarcity+G1_Associated*naturalhazards+G1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))

#Resource Stress and Gods that Harm with Weather#
summary(model2<-lmer(stress2~G2_Harm+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
lm.beta.lmer(model2)

##Secondary analyses looking at specific types of resource stress # Table S4, Model S2 
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model2b<-lmer(stress2~G2_Harm*famine+G2_Harm*naturalhazards+G2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model2c<-lmer(stress2~G2_Harm*famine+G2_Harm*chron_scarcity+G2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model2d<-lmer(stress2~G2_Harm*famine+G2_Harm*naturalhazards+G2_Harm*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model2e<-lmer(stress2~G2_Harm*chron_scarcity+G2_Harm*naturalhazards+G2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))

#Resource Stress and Gods that Help with Weather# Table 2, Model 3
summary(model3<-lmer(stress2~G3_Help+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
lm.beta.lmer(model3)

##Secondary analyses looking at specific types of resource stress# Table S4, Model S3
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model3b<-lmer(stress2~G3_Help*famine+G3_Help*naturalhazards+G3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model3c<-lmer(stress2~G3_Help*famine+G3_Help*chron_scarcity+G3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model3d<-lmer(stress2~G3_Help*famine+G3_Help*naturalhazards+G3_Help*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model3e<-lmer(stress2~G3_Help*chron_scarcity+G3_Help*naturalhazards+G3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))

#Resource Stress and Gods that Punish with Weather# Table 2, Model 4
summary(model4<-lmer(stress2~G4_Punitive+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
lm.beta.lmer(model4)

##Secondary analyses looking at specific types of resource stress# Table S4, Model S4
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model4b<-lmer(stress2~G4_Punitive*famine+G4_Punitive*naturalhazards+G4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model4c<-lmer(stress2~G4_Punitive*famine+G4_Punitive*chron_scarcity+G4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model4d<-lmer(stress2~G4_Punitive*famine+G4_Punitive*naturalhazards+G4_Punitive*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model4e<-lmer(stress2~G4_Punitive*chron_scarcity+G4_Punitive*naturalhazards+G4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))

#Modeling Helping and Harming Together# Table 2, Model 5
summary(model5<-lmer(stress2~G3_Help+G2_Harm+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang_hurthelp)) 
vif(model5) ##Look at Multicollinearity, since the predictors are highly correlated

#Resource Stress and Superior Gods Associated With Weather# Table 3, Model 6
summary(model6<-lmer(stress2~SG1_Associated+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+SinoTibetan+Uralic+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S5, Model S5
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model6b<-lmer(stress2~SG1_Associated*famine+SG1_Associated*naturalhazards+SG1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Indoeuropean+NuclearTransNewGuinea+SinoTibetan+Uralic+(1|SCCS_ID),data=ang))
summary(model6c<-lmer(stress2~SG1_Associated*famine+SG1_Associated*chron_scarcity+SG1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Indoeuropean+NuclearTransNewGuinea+SinoTibetan+Uralic+(1|SCCS_ID),data=ang))
summary(model6d<-lmer(stress2~SG1_Associated*famine+SG1_Associated*naturalhazards+SG1_Associated*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Indoeuropean+NuclearTransNewGuinea+SinoTibetan+Uralic+(1|SCCS_ID),data=ang))
summary(model6e<-lmer(stress2~SG1_Associated*chron_scarcity+SG1_Associated*naturalhazards+SG1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Indoeuropean+NuclearTransNewGuinea+SinoTibetan+Uralic+(1|SCCS_ID),data=ang))

#Resource Stress and Superior Gods that Harm With Weather# Table 3, Model 7
summary(model7<-lmer(stress2~SG2_Harm+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S5, Model S6
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model7b<-lmer(stress2~SG2_Harm*famine+SG2_Harm*naturalhazards+SG2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model7c<-lmer(stress2~SG2_Harm*famine+SG2_Harm*chron_scarcity+SG2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model7d<-lmer(stress2~SG2_Harm*famine+SG2_Harm*naturalhazards+SG2_Harm*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))
summary(model7e<-lmer(stress2~SG2_Harm*chron_scarcity+SG2_Harm*naturalhazards+SG2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=ang))

#Resource Stress and Superior Gods that Help With Weather# Table 3, Model 8
summary(model8<-lmer(stress2~SG3_Help+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+SinoTibetan+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S5, Model S7
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model8b<-lmer(stress2~SG3_Help*famine+SG3_Help*naturalhazards+SG3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+SinoTibetan+(1|SCCS_ID),data=ang))
summary(model8c<-lmer(stress2~SG3_Help*famine+SG3_Help*chron_scarcity+SG3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+SinoTibetan+(1|SCCS_ID),data=ang))
summary(model8d<-lmer(stress2~SG3_Help*famine+SG3_Help*naturalhazards+SG3_Help*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+SinoTibetan+(1|SCCS_ID),data=ang))
summary(model8e<-lmer(stress2~SG3_Help*chron_scarcity+SG3_Help*naturalhazards+SG3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+SinoTibetan+(1|SCCS_ID),data=ang))

#Resource Stress and Superior Gods that Punish With Weather# Table 3, Model 9
summary(model9<-lmer(stress2~SG4_Punitive+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S5, Model S8
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model9b<-lmer(stress2~SG4_Punitive*famine+SG4_Punitive*naturalhazards+SG4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model9c<-lmer(stress2~SG4_Punitive*famine+SG4_Punitive*chron_scarcity+SG4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model9d<-lmer(stress2~SG4_Punitive*famine+SG4_Punitive*naturalhazards+SG4_Punitive*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model9e<-lmer(stress2~SG4_Punitive*chron_scarcity+SG4_Punitive*naturalhazards+SG4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

#Resource Stress and Minor Spirits Associated With Weather# Table 4, Model 10
summary(model10<-lmer(stress2~MS1_Associated+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Turkic+NuclearTransNewGuinea+Austroasiatic+SinoTibetan+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S6, Model S9
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model0b<-lmer(stress2~MS1_Associated*famine+MS1_Associated*naturalhazards+MS1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Turkic+NuclearTransNewGuinea+Austroasiatic+SinoTibetan+(1|SCCS_ID),data=ang))
summary(model0c<-lmer(stress2~MS1_Associated*famine+MS1_Associated*chron_scarcity+MS1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Turkic+NuclearTransNewGuinea+Austroasiatic+SinoTibetan+(1|SCCS_ID),data=ang))
summary(model0d<-lmer(stress2~MS1_Associated*famine+MS1_Associated*naturalhazards+MS1_Associated*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Turkic+NuclearTransNewGuinea+Austroasiatic+SinoTibetan+(1|SCCS_ID),data=ang))
summary(model0e<-lmer(stress2~MS1_Associated*chron_scarcity+MS1_Associated*naturalhazards+MS1_Associated*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Turkic+NuclearTransNewGuinea+Austroasiatic+SinoTibetan+(1|SCCS_ID),data=ang))

#Resource Stress and Minor Spirits That Harm With Weather# Table 4, Model 11
summary(model11<-lmer(stress2~MS2_Harm+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S6, Model S10
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model11b<-lmer(stress2~MS2_Harm*famine+MS2_Harm*naturalhazards+MS2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model11c<-lmer(stress2~MS2_Harm*famine+MS2_Harm*chron_scarcity+MS2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model11d<-lmer(stress2~MS2_Harm*famine+MS2_Harm*naturalhazards+MS2_Harm*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model11e<-lmer(stress2~MS2_Harm*chron_scarcity+MS2_Harm*naturalhazards+MS2_Harm*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))


#Resource Stress and Minor Spirits That Help With Weather# Table 4, Model 12
summary(model12<-lmer(stress2~MS3_Help+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S6, Model S11
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model12b<-lmer(stress2~MS3_Help*famine+MS3_Help*naturalhazards+MS3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model12c<-lmer(stress2~MS3_Help*famine+MS3_Help*chron_scarcity+MS3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model12d<-lmer(stress2~MS3_Help*famine+MS3_Help*naturalhazards+MS3_Help*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model12e<-lmer(stress2~MS3_Help*chron_scarcity+MS3_Help*naturalhazards+MS3_Help*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

#Resource Stress and Minor Spirits That Punish With Weather# Table 4, Model 13
summary(model13<-lmer(stress2~MS4_Punitive+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

##Secondary analyses looking at specific types of resource stress# Table S6, Model S12
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model13b<-lmer(stress2~MS4_Punitive*famine+MS4_Punitive*naturalhazards+MS4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model13c<-lmer(stress2~MS4_Punitive*famine+MS4_Punitive*chron_scarcity+MS4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model13d<-lmer(stress2~MS4_Punitive*famine+MS4_Punitive*naturalhazards+MS4_Punitive*chron_scarcity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))
summary(model13e<-lmer(stress2~MS4_Punitive*chron_scarcity+MS4_Punitive*naturalhazards+MS4_Punitive*abundance+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+NuclearTransNewGuinea+(1|SCCS_ID),data=ang))

#################################################
##For the Existing SCCS Moralizing God Variable##
#################################################

##Dummy-coding language family (there are more language families represented in these data)
mor$Afroasiatic<-ifelse(as.character(mor$language_family)=="Afro-Asiatic",1,0)
mor$Afroasiatic<-ifelse(is.na(mor$Afroasiatic)==T,0,mor$Afroasiatic)

mor$Austronesian<-ifelse(as.character(mor$language_family)=="Austronesian",1,0)
mor$Austronesian<-ifelse(is.na(mor$Austronesian)==T,0,mor$Austronesian)

mor$AtlanticCongo<-ifelse(as.character(mor$language_family)=="Atlantic-Congo",1,0)
mor$AtlanticCongo<-ifelse(is.na(mor$AtlanticCongo)==T,0,mor$AtlanticCongo)

mor$Indoeuropean<-ifelse(as.character(mor$language_family)=="Indo-European",1,0)
mor$Indoeuropean<-ifelse(is.na(mor$Indoeuropean)==T,0,mor$Indoeuropean)

mor$Algic<-ifelse(as.character(mor$language_family)=="Algic",1,0)
mor$Algic<-ifelse(is.na(mor$Algic)==T,0,mor$Algic)

mor$Arawakan<-ifelse(as.character(mor$language_family)=="Arawakan",1,0)
mor$Arawakan<-ifelse(is.na(mor$Arawakan)==T,0,mor$Arawakan)

mor$AthabaskanEyakTlingit<-ifelse(as.character(mor$language_family)=="Athabaskan-Eyak-Tlingit",1,0)
mor$AthabaskanEyakTlingit<-ifelse(is.na(mor$AthabaskanEyakTlingit)==T,0,mor$AthabaskanEyakTlingit)

mor$Chibchan<-ifelse(as.character(mor$language_family)=="Chibchan",1,0)
mor$Chibchan<-ifelse(is.na(mor$Chibchan)==T,0,mor$Chibchan)

mor$SinoTibetan<-ifelse(as.character(mor$language_family)=="Sino-Tibetan",1,0)
mor$SinoTibetan<-ifelse(is.na(mor$SinoTibetan)==T,0,mor$SinoTibetan)

mor$Tupian<-ifelse(as.character(mor$language_family)=="Tupian",1,0)
mor$Tupian<-ifelse(is.na(mor$Tupian)==T,0,mor$Tupian)

mor$Turkic<-ifelse(as.character(mor$language_family)=="Turkic",1,0)
mor$Turkic<-ifelse(is.na(mor$Turkic)==T,0,mor$Turkic)

mor$Uralic<-ifelse(as.character(mor$language_family)=="Uralic",1,0)
mor$Uralic<-ifelse(is.na(mor$Uralic)==T,0,mor$Uralic)

mor$UtoAztecan<-ifelse(as.character(mor$language_family)=="Uto-Aztecan",1,0)
mor$UtoAztecan<-ifelse(is.na(mor$UtoAztecan)==T,0,mor$UtoAztecan)

mor$Dravidian<-ifelse(as.character(mor$language_family)=="Dravidian",1,0)
mor$Dravidian<-ifelse(is.na(mor$Dravidian)==T,0,mor$Dravidian)

mor$EskimoAleut<-ifelse(as.character(mor$language_family)=="Eskimo-Aleut",1,0)
mor$EskimoAleut<-ifelse(is.na(mor$EskimoAleut)==T,0,mor$EskimoAleut)

mor$Mande<-ifelse(as.character(mor$language_family)=="Mande",1,0)
mor$Mande<-ifelse(is.na(mor$Mande)==T,0,mor$Mande)

mor$Salishan<-ifelse(as.character(mor$language_family)=="Salishan",1,0)
mor$Salishan<-ifelse(is.na(mor$Salishan)==T,0,mor$Salishan)

table(mor$language_family)

##Dummy coding stressors
#Famine=1 under ResourceStress Variable
mor$famine<-ifelse(as.character(mor$ResourceStress)=="1",1,0)
#Natural Hazards=2 under ResourceStress Variable
mor$naturalhazards<-ifelse(as.character(mor$ResourceStress)=="2",1,0)
#Chronic Scarcity=3 under ResourceStress Variable
mor$chron_scarcity<-ifelse(as.character(mor$ResourceStress)=="3",1,0)
#Abundance=4 under ResourceStress Variable 
mor$abundance<-ifelse(as.character(mor$ResourceStress)=="4",1,0)

#Reverse-coding abundance to be able to merge it with other stressors
mor$stress2<-ifelse(mor$abundance==1,(0-mor$Stress),
                    ifelse(mor$abundance==0,mor$Stress,
                           NA))

#Dichotomized stress variable (not necessary for this analysis)
mor$stress3<-ifelse(mor$chron_scarcity==1&mor$stress2<(-.11),(-1),
                    ifelse(mor$chron_scarcity==1&mor$stress2>(-.11),1,
                           ifelse(mor$famine==1&mor$stress2<(-.11),(-1),
                                  ifelse(mor$famine==1&mor$stress2>(-.11),1,
                                         mor$stress2))))

##Creating datasets for each of the god variables (this can be used to compute R2 later)
mor2<-mor[is.na(mor$MHG_PresAbs)==F,]
mor3<-mor[is.na(mor$HG_Pres_MoralNoYes)==F,]

#Empty model (anova) (this can be used to compute R2 later)
summary(fit1<-lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=mor)) 
summary(fit1<-lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=mor2)) 
summary(fit1<-lmer(stress2~Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+(1|SCCS_ID),data=mor3))


##For the analyses below, you can get the R squared by comparing...
#...the variance explained at the society level with the variance explained by language family alone

#Association between presence of moralizing high gods and resource stress# Table 5, Model 14
summary(model14<-lmer(stress2~MHG_PresAbs+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+Dravidian+EskimoAleut+Mande+Salishan+(1|SCCS_ID),data=mor))
summary(model14<-lmer(scale(stress2)~scale(MHG_PresAbs)+scale(Social_Complexity)+scale(Afroasiatic)+scale(Austronesian)+scale(AtlanticCongo)+scale(Indoeuropean)+scale(Algic)+scale(Arawakan)+scale(AthabaskanEyakTlingit)+scale(Chibchan)+scale(SinoTibetan)+scale(Tupian)+scale(Turkic)+scale(Uralic)+scale(UtoAztecan)+scale(Dravidian)+scale(EskimoAleut)+scale(Mande)+scale(Salishan)+(1|SCCS_ID),data=mor))
##Extract effect size from model on line 355, since beta.lmer is having trouble with this model

##Secondary analyses looking at specific types of resource stress# Table S7, Model S13
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model14b<-lmer(stress2~MHG_PresAbs*famine+MHG_PresAbs*naturalhazards+MHG_PresAbs*abundance+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+Dravidian+EskimoAleut+Mande+Salishan+(1|SCCS_ID),data=mor))
summary(model14c<-lmer(stress2~MHG_PresAbs*famine+MHG_PresAbs*chron_scarcity+MHG_PresAbs*abundance+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+Dravidian+EskimoAleut+Mande+Salishan+(1|SCCS_ID),data=mor))
summary(model14d<-lmer(stress2~MHG_PresAbs*famine+MHG_PresAbs*naturalhazards+MHG_PresAbs*chron_scarcity+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+Dravidian+EskimoAleut+Mande+Salishan+(1|SCCS_ID),data=mor))
summary(model14e<-lmer(stress2~MHG_PresAbs*chron_scarcity+MHG_PresAbs*naturalhazards+MHG_PresAbs*abundance+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+Dravidian+EskimoAleut+Mande+Salishan+(1|SCCS_ID),data=mor))


#Association between presence of moralizing high gods amongst societies that have high gods and resource stress# Table 5, Model 15
summary(model15<-lmer(stress2~HG_Pres_MoralNoYes+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+(1|SCCS_ID),data=mor))
summary(model15<-lmer(scale(stress2)~scale(HG_Pres_MoralNoYes)+scale(Social_Complexity)+scale(Afroasiatic)+scale(Austronesian)+
                        scale(AtlanticCongo)+scale(Indoeuropean)+scale(Algic)+scale(Arawakan)+scale(AthabaskanEyakTlingit)+
                        scale(Chibchan)+scale(SinoTibetan)+scale(Tupian)+scale(Turkic)+scale(Uralic)+scale(UtoAztecan)+(1|SCCS_ID),data=mor))
##Extract effect size from model on line 368, since beta.lmer is having trouble with this model

##Secondary analyses looking at specific types of resource stress# Table S7, Model S14
#Looking at Chronic scarcity, then natural hazards, then abundance, then famine
summary(model15b<-lmer(stress2~HG_Pres_MoralNoYes*famine+HG_Pres_MoralNoYes*naturalhazards+HG_Pres_MoralNoYes*abundance+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+(1|SCCS_ID),data=mor))
summary(model15c<-lmer(stress2~HG_Pres_MoralNoYes*famine+HG_Pres_MoralNoYes*chron_scarcity+HG_Pres_MoralNoYes*abundance+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+(1|SCCS_ID),data=mor))
summary(model15d<-lmer(stress2~HG_Pres_MoralNoYes*famine+HG_Pres_MoralNoYes*naturalhazards+HG_Pres_MoralNoYes*chron_scarcity+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+(1|SCCS_ID),data=mor))
summary(model15e<-lmer(stress2~HG_Pres_MoralNoYes*chron_scarcity+HG_Pres_MoralNoYes*naturalhazards+HG_Pres_MoralNoYes*abundance+Social_Complexity+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan+(1|SCCS_ID),data=mor))

###########################################################################
###PART B: Central Analyses on Stressors, Religious Beliefs, and Sharing###
###########################################################################


mor<-read.csv("DT-MoralHighGods_EP.csv")
ang<-read.csv("DT-GodsAndResourceStress.csv")

library(Kendall)

##Dummy-coding language family information
mor$Afroasiatic<-ifelse(as.character(mor$language_family)=="Afro-Asiatic",1,0)
mor$Afroasiatic<-ifelse(is.na(mor$Afroasiatic)==T,0,mor$Afroasiatic)
mor$Austronesian<-ifelse(as.character(mor$language_family)=="Austronesian",1,0)
mor$Austronesian<-ifelse(is.na(mor$Austronesian)==T,0,mor$Austronesian)
mor$AtlanticCongo<-ifelse(as.character(mor$language_family)=="Atlantic-Congo",1,0)
mor$AtlanticCongo<-ifelse(is.na(mor$AtlanticCongo)==T,0,mor$AtlanticCongo)
mor$Indoeuropean<-ifelse(as.character(mor$language_family)=="Indo-European",1,0)
mor$Indoeuropean<-ifelse(is.na(mor$Indoeuropean)==T,0,mor$Indoeuropean)
mor$Afroasiatic<-ifelse(as.character(mor$language_family)=="Afro-Asiatic",1,0)
mor$Afroasiatic<-ifelse(is.na(mor$Afroasiatic)==T,0,mor$Afroasiatic)
mor$Austronesian<-ifelse(as.character(mor$language_family)=="Austronesian",1,0)
mor$Austronesian<-ifelse(is.na(mor$Austronesian)==T,0,mor$Austronesian)
mor$AtlanticCongo<-ifelse(as.character(mor$language_family)=="Atlantic-Congo",1,0)
mor$AtlanticCongo<-ifelse(is.na(mor$AtlanticCongo)==T,0,mor$AtlanticCongo)
mor$Indoeuropean<-ifelse(as.character(mor$language_family)=="Indo-European",1,0)
mor$Indoeuropean<-ifelse(is.na(mor$Indoeuropean)==T,0,mor$Indoeuropean)
mor$Algic<-ifelse(as.character(mor$language_family)=="Algic",1,0)
mor$Algic<-ifelse(is.na(mor$Algic)==T,0,mor$Algic)
mor$Arawakan<-ifelse(as.character(mor$language_family)=="Arawakan",1,0)
mor$Arawakan<-ifelse(is.na(mor$Arawakan)==T,0,mor$Arawakan)
mor$AthabaskanEyakTlingit<-ifelse(as.character(mor$language_family)=="Athabaskan-Eyak-Tlingit",1,0)
mor$AthabaskanEyakTlingit<-ifelse(is.na(mor$AthabaskanEyakTlingit)==T,0,mor$AthabaskanEyakTlingit)
mor$Chibchan<-ifelse(as.character(mor$language_family)=="Chibchan",1,0)
mor$Chibchan<-ifelse(is.na(mor$Chibchan)==T,0,mor$Chibchan)
mor$SinoTibetan<-ifelse(as.character(mor$language_family)=="Sino-Tibetan",1,0)
mor$SinoTibetan<-ifelse(is.na(mor$SinoTibetan)==T,0,mor$SinoTibetan)
mor$Tupian<-ifelse(as.character(mor$language_family)=="Tupian",1,0)
mor$Tupian<-ifelse(is.na(mor$Tupian)==T,0,mor$Tupian)
mor$Turkic<-ifelse(as.character(mor$language_family)=="Turkic",1,0)
mor$Turkic<-ifelse(is.na(mor$Turkic)==T,0,mor$Turkic)
mor$Uralic<-ifelse(as.character(mor$language_family)=="Uralic",1,0)
mor$Uralic<-ifelse(is.na(mor$Uralic)==T,0,mor$Uralic)
mor$UtoAztecan<-ifelse(as.character(mor$language_family)=="Uto-Aztecan",1,0)
mor$UtoAztecan<-ifelse(is.na(mor$UtoAztecan)==T,0,mor$UtoAztecan)
mor$Dravidian<-ifelse(as.character(mor$language_family)=="Dravidian",1,0)
mor$Dravidian<-ifelse(is.na(mor$Dravidian)==T,0,mor$Dravidian)
mor$EskimoAleut<-ifelse(as.character(mor$language_family)=="Eskimo-Aleut",1,0)
mor$EskimoAleut<-ifelse(is.na(mor$EskimoAleut)==T,0,mor$EskimoAleut)
mor$Mande<-ifelse(as.character(mor$language_family)=="Mande",1,0)
mor$Mande<-ifelse(is.na(mor$Mande)==T,0,mor$Mande)
mor$Salishan<-ifelse(as.character(mor$language_family)=="Salishan",1,0)
mor$Salishan<-ifelse(is.na(mor$Salishan)==T,0,mor$Salishan)

ang$Afroasiatic<-ifelse(as.character(ang$language_family)=="Afro-Asiatic",1,0)
ang$Afroasiatic<-ifelse(is.na(ang$Afroasiatic)==T,0,ang$Afroasiatic)
ang$Austronesian<-ifelse(as.character(ang$language_family)=="Austronesian",1,0)
ang$Austronesian<-ifelse(is.na(ang$Austronesian)==T,0,ang$Austronesian)
ang$AtlanticCongo<-ifelse(as.character(ang$language_family)=="Atlantic-Congo",1,0)
ang$AtlanticCongo<-ifelse(is.na(ang$AtlanticCongo)==T,0,ang$AtlanticCongo)
ang$Indoeuropean<-ifelse(as.character(ang$language_family)=="Indo-European",1,0)
ang$Indoeuropean<-ifelse(is.na(ang$Indoeuropean)==T,0,ang$Indoeuropean)
ang$NuclearTransNewGuinea<-ifelse(as.character(ang$language_family)=="Nuclear Trans New Guinea",1,0)
ang$NuclearTransNewGuinea<-ifelse(is.na(ang$NuclearTransNewGuinea)==T,0,ang$NuclearTransNewGuinea)
ang$SinoTibetan<-ifelse(as.character(ang$language_family)=="Sino-Tibetan",1,0)
ang$SinoTibetan<-ifelse(is.na(ang$SinoTibetan)==T,0,ang$SinoTibetan)
ang$Uralic<-ifelse(as.character(ang$language_family)=="Uralic",1,0)
ang$Uralic<-ifelse(is.na(ang$Uralic)==T,0,ang$Uralic)
ang$Turkic<-ifelse(as.character(ang$language_family)=="Turkic",1,0)
ang$Turkic<-ifelse(is.na(ang$Turkic)==T,0,ang$Turkic)
ang$Austroasiatic<-ifelse(as.character(ang$language_family)=="Austroasiatic",1,0)
ang$Austroasiatic<-ifelse(is.na(ang$Austroasiatic)==T,0,ang$Austroasiatic)


##Creating single stressor variable
ang$ZSCCS_Abundance_imp_rev<-(0-ang$ZSCCS_Abundance_imp)
mor$ZSCCS_Abundance_imp_rev<-(0-mor$ZSCCS_Abundance_imp)
ang_stress<-data.frame(ang$ZFamine,ang$ZNatural_Hazards,ang$ZChronic_Scarcity,ang$ZSCCS_Abundance_imp_rev)
ang$stress<-rowMeans(ang_stress,na.rm=T)
mor_stress<-data.frame(mor$ZFamine,mor$ZNatural_Hazards,mor$ZChronic_Scarcity,mor$ZSCCS_Abundance_imp_rev)
mor$stress<-rowMeans(mor_stress,na.rm=T)

##Regressions on Food Sharing
summary(model16<-lm(S3_Resolved_IA_TR~stress+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model16)
confint(lm.beta(model16),level = .95)
summary(model17<-lm(S3_Resolved_IA_TR~stress+G3_Help+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model17)
confint(lm.beta(model17),level = .95)
summary(model18<-lm(S3_Resolved_IA_TR~stress+G2_Harm+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model18)
confint(lm.beta(model18),level = .95)
summary(model19<-lm(S3_Resolved_IA_TR~stress+G4_Punitive+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model19)
confint(lm.beta(model19),level = .95)
summary(model20<-lm(S3_Resolved_IA_TR~HG_Pres_MoralNoYes+Social_Complexity+stress+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Arawakan+SinoTibetan+Uralic+UtoAztecan,data=mor))
lm.beta(model20)
confint(lm.beta(model20),level = .95)

##Just regressions for labor sharing
summary(model21<-lm(LS4_Res_IA_2~stress+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model21)
confint(lm.beta(model21),level = .95)
summary(model22<-lm(LS4_Res_IA_2~stress+G3_Help+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model22)
confint(lm.beta(model22),level = .95)
summary(model23<-lm(LS4_Res_IA_2~stress+G2_Harm+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model23)
confint(lm.beta(model23),level = .95)
summary(model24<-lm(LS4_Res_IA_2~stress+G4_Punitive+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean,data=ang))
lm.beta(model24)
confint(lm.beta(model24),level = .95)
summary(model25<-lm(LS4_Res_IA_2~HG_Pres_MoralNoYes+Social_Complexity+stress+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Arawakan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan,data=mor))
lm.beta(model25)
confint(lm.beta(model25),level = .95)

#making social complexity and moral high god interaction term
mor$HGPres_x_Complex=mor$HG_Pres_MoralNoYes * mor$Social_Complexity

summary(model20a<-lm(S3_Resolved_IA_TR~HG_Pres_MoralNoYes+Social_Complexity+stress+ HGPres_x_Complex+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan,data=mor))
summary(model20a)
summary(lm.beta(model20a))
summary.lm.beta(model20a)
confint(lm.beta(model20a),level = .95)

summary(model25a<-lm(LS4_Res_IA_2~HG_Pres_MoralNoYes+Social_Complexity+stress+ HGPres_x_Complex+Afroasiatic+Austronesian+AtlanticCongo+Indoeuropean+Algic+Arawakan+AthabaskanEyakTlingit+Chibchan+SinoTibetan+Tupian+Turkic+Uralic+UtoAztecan,data=mor))
summary(model25a)
lm.beta(model25a)
summary.lm.beta(model25a)
confint(lm.beta(model25a),level = .95)

###########################
##Monte Carlo Simulations##
##########################
##Model involving Resource Stress, Helping Gods, and Food Sharing
require(MASS)
a=0.9401
b=-0.2995
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.07160976, 0,
  0, 0.04109729
),2,2)
acov
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Harming Gods, and Food Sharing
require(MASS)
a=1.194
b=-0.69148
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.11048976, 0,
  0, 0.0482022
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Punitive Gods, and Food Sharing
require(MASS)
a=0.2929
b=-0.3853
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.08988004, 0,
  0, 0.05175625
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Moralizing Gods, and Food Sharing
require(MASS)
a=0.468225
b=-0.286069
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.03375671, 0,
  0, 0.03394769
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Helping Gods, and Labor Sharing
require(MASS)
a=0.9401
b=0.08712
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.07160976, 0,
  0, 0.076016
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Harming Gods, and Labor Sharing
require(MASS)
a=1.194
b=0.1561
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.11048976, 0,
  0, 0.18198756
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Punitive Gods, and Labor Sharing
require(MASS)
a=0.2929
b=0.09117
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.08988004, 0,
  0, 0.05847208
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')

##Model involving Resource Stress, Moralizing Gods, and Labor Sharing
require(MASS)
a=0.468225
b=0.188035
rep=20000
conf=95
pest=c(a,b)
acov <- matrix(c(
  0.03375671, 0,
  0, 0.02658661
),2,2)
mcmc <- mvrnorm(rep,pest,acov,empirical=FALSE)
ab <- mcmc[,1]*mcmc[,2]
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
     main='Distribution of Indirect Effect')