###SUPPLEMENTARY ANALYSIS FOR PEPPER ET AL
###THE INFLUENCE OF MORTALITY AND SOCIOECONOMIC STATUS ON RISK AND DELAYED REWARDS: A REPLICATION WITH BRITISH PARTICIPANTS

####Setup basics####

#Open packages
library(ggplot2)
library(dplyr)
library(car)
library(psych)
library(metafor)
library(pwr)

#Clear workspace
rm(list=ls())

#Load data
Replication1 = read.csv("PepperReplication1Data.csv")
Replication2 = read.csv("PepperReplication2Data.csv")
Replication3 = read.csv("PepperReplication3Data.csv")

##Set parameters for directed tests
#Set values of alpha, delta and gamma - these parameters must be set for directed tests (see main analyses, below) to work
#To replicate the approach by Griskevicius et al (2011), delta/gamma = 0.01/0.04, and gamma/alpha=0.04/0.05
alpha = 0.05
delta = 0.01
gamma = 0.04

####Replication 1 - main models####

#Risk, interaction model to check for predicted interaction between child SES and condition
RiskInteractionMain1 = lm(RiskTotal~scale(ChildSES)*ConditionString, data=Replication1)
Anova(RiskInteractionMain1, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteractionMain1)
confint(RiskInteractionMain1)

##Extract betas and standar errors for meta-analaysis (see "Meta-analysis of main interactions" section, below) 
#Extract and store betas and standard errors for main effect of condition from the model "RiskInteractionMain1" for meta-analysis (see below) 
RiskConditionEffectb1 = coef(summary(RiskInteractionMain1))["ConditionStringMortality","Estimate"]
RiskConditionEffectse1 = coef(summary(RiskInteractionMain1))["ConditionStringMortality","Std. Error"]

#Extract and store betas and standard errors for main effect of child SES from the model "RiskInteractionMain1" for meta-analysis (see below) 
RiskSESEffectb1 = coef(summary(RiskInteractionMain1))["scale(ChildSES)","Estimate"]
RiskSESEffectse1 = coef(summary(RiskInteractionMain1))["scale(ChildSES)","Std. Error"]

#Extract and store betas and standard errors for interaction between condition and child SES from the model "RiskInteractionMain1" for meta-analysis (see below) 
RiskInteractionb1 = coef(summary(RiskInteractionMain1))["scale(ChildSES):ConditionStringMortality","Estimate"]
RiskInteractionse1 = coef(summary(RiskInteractionMain1))["scale(ChildSES):ConditionStringMortality","Std. Error"]


##Calculate directed p-values for interaction between child SES and condition on risk for replication 1
#To set the parameters for this test, run the code in the head of the script, under startup basics

#Extract t-values from summary of RiskInteractionMain1
TValueRisk1 = coef(summary(RiskInteractionMain1))["scale(ChildSES):ConditionStringMortality","t value"]
TValueRisk1

#Extract degrees of freedom from RiskInteractionMain1
dfRisk1 = df.residual(RiskInteractionMain1)
dfRisk1

#Calculate 1-tailed p-value for RiskInteractionMain1
p1tRisk1 = pt(TValueRisk1, df=dfRisk1, lower.tail=FALSE)
p1tRisk1

##Get directed test p-value
# Check if p1t is greater than or less than gamma/alpha and apply formula from Rice & Gaines (1994), box 1, as appropriate
pDirectedRisk1=ifelse(p1tRisk1 < gamma/alpha, (p1tRisk1*(1+delta/gamma)), ((1-p1tRisk1)*(1+gamma/delta)))  
pDirectedRisk1

#Discounting interaction model to check for predicted interaction between child SES and condition
DiscountingInteractionMain1 = lm(DiscountingTotal~scale(ChildSES)*ConditionString, data=Replication1)
Anova(DiscountingInteractionMain1, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteractionMain1)
confint(DiscountingInteractionMain1)

##Extract betas and standar errors for meta-analaysis (see "Meta-analysis of main interactions" section, below) 
#Extract and store betas and standard errors for main effect of condition from the model "DiscountingInteractionMain1" for meta-analysis (see below) 
DiscountingConditionEffectb1 = coef(summary(DiscountingInteractionMain1))["ConditionStringMortality","Estimate"]
DiscountingConditionEffectse1 = coef(summary(DiscountingInteractionMain1))["ConditionStringMortality","Std. Error"]

#Extract and store betas and standard errors for main effect of child SES from the model "DiscountingInteractionMain1" for meta-analysis (see below) 
DiscountingSESEffectb1 = coef(summary(DiscountingInteractionMain1))["scale(ChildSES)","Estimate"]
DiscountingSESEffectse1 = coef(summary(DiscountingInteractionMain1))["scale(ChildSES)","Std. Error"]

#Extract and store betas and standard errors for interaction between condition and child SES from the model "DiscountingInteractionMain1" for meta-analysis (see below) 
DiscountingInteractionb1 = coef(summary(DiscountingInteractionMain1))["scale(ChildSES):ConditionStringMortality","Estimate"]
DiscountingInteractionse1 = coef(summary(DiscountingInteractionMain1))["scale(ChildSES):ConditionStringMortality","Std. Error"]

##Calculate directed p-values for interaction between child SES and condition on delay discounting for replication 1
#To set the parameters for this test, run the code in the head of the script, under startup basics

#Extract t-values from summary of DiscountingInteractionMain1
TValueDiscounting1 = coef(summary(DiscountingInteractionMain1))["scale(ChildSES):ConditionStringMortality","t value"]
TValueDiscounting1

#Extract degrees of freedom from RiskInteractionMain1
dfDiscounting1 = df.residual(DiscountingInteractionMain1)
dfDiscounting1

#Calculate 1-tailed p-value for RiskInteractionMain1
p1tDiscounting1 = pt(TValueDiscounting1, df=dfDiscounting1, lower.tail=FALSE)
p1tDiscounting1

##Get directed test p-value
# Check if p1tDiscounting1 is greater than or less than gamma/alpha and apply formula from Rice & Gaines (1994), box 1, as appropriate
pDirectedDiscounting1=ifelse(p1tDiscounting1 < gamma/alpha, (p1tDiscounting1*(1+delta/gamma)), ((1-p1tDiscounting1)*(1+gamma/delta)))  
pDirectedDiscounting1

####Replication 1 - figures####

##Make data frame & models to plot interactions

#Find limits for "High" (>1 sd above mean) and "Low" (<1 sd below mean) SES
plusone1=mean(Replication1$ChildSES)+sd(Replication1$ChildSES)
minusone1=mean(Replication1$ChildSES)-sd(Replication1$ChildSES)
#Also store the mean ChildSES as a variable
MeanChildSES1 = mean(Replication1$ChildSES)

#Create a new data frame for comparing "High" and "Low" SES for risk
HighLow1=data.frame(Condition=c(1, 2, 1, 2), ChildSES=c(plusone1, plusone1, minusone1, minusone1))

#Add row names for "high" and "low" SES
row.names(HighLow1)=c("Condition 1, high", "Condition 2, high", "Condition 1, low", "Condition 2, low")
HighLow1$cSES=c("High", "High", "Low", "Low")

#Build model to get predictions and standard errors for risk
RiskModel1 = lm(RiskTotal~scale(ChildSES)+Condition, data=Replication1)
predict(RiskModel1, newdata=HighLow1, se.fit=TRUE)

#Build model to get predictions and standard errors for delay discounting
DiscountingModel1 = lm(DiscountingTotal~scale(ChildSES)+Condition, data=Replication1)
predict(DiscountingModel1, newdata=HighLow1, se.fit=TRUE)

#Add risk predictions to data frame
HighLow1$RiskPrediction=predict(RiskModel1, newdata=HighLow1, se.fit=TRUE)$fit

#Add risk standard errors to data frame
HighLow1$RiskSE=predict(RiskModel1, newdata=HighLow1, se.fit=TRUE)$se.fit

#Add delay discounting predictions to data frame
HighLow1$DiscountingPrediction=predict(DiscountingModel1, newdata=HighLow1, se.fit=TRUE)$fit

#Add delay discounting standard errors to data frame
HighLow1$DiscountingSE=predict(DiscountingModel1, newdata=HighLow1, se.fit=TRUE)$se.fit

#Add condition labels to data frame
HighLow1$Condition=c("Mortality", "Control", "Mortality", "Control")

#View the data frame
View(HighLow1)

#Bar plot of risk by condition for high and low child SES (Figure 1a)
Fig1a=ggplot(HighLow1, aes(x=Condition, y=RiskPrediction, fill=cSES)) +
  theme_bw() + labs (x = "Experimental condition", y = "Risk acceptance") +
  geom_bar(stat="identity", position="dodge", colour = "black") + ggtitle('(a)') +
  geom_errorbar(width=0.2, position=position_dodge(0.9), aes(ymin=RiskPrediction-RiskSE, ymax=RiskPrediction+RiskSE)) + 
  scale_fill_manual(values=c("black", "darkgrey"))

#Save Figure 1a to working directory
png("Figure1a.png", res=300, width=3, height=3, units="in")
plot(Fig1a)
dev.off()

#Show Figure 1a in the environment
Fig1a

#Bar plot of delay discounting by condition for high and low child SES (Figure 1b) 
Fig1b=ggplot(HighLow1, aes(x=Condition, y=DiscountingPrediction, fill=cSES)) +
  theme_bw() + labs (x = "Experimental condition", y = "Delay discounting") +
  geom_bar(stat="identity", position="dodge", colour = "black") + ggtitle('(b)') +
  geom_errorbar(width=0.2, position=position_dodge(0.9), aes(ymin=DiscountingPrediction-DiscountingSE, ymax=DiscountingPrediction+DiscountingSE))+
  scale_fill_manual(values=c("black", "darkgrey"))

#Save Figure 1b to working directory
png("Figure1b.png", res=300, width=3, height=3, units="in")
plot(Fig1b)
dev.off()

#Show Figure 1b in the environment
Fig1b

####Replication 2 - main models####

#Risk, interaction model to check for predicted interaction between child SES and condition
RiskInteractionMain2 = lm(RiskTotal~scale(ChildSES)*ConditionString, data=Replication2)
Anova(RiskInteractionMain2, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteractionMain2)
confint(RiskInteractionMain2)

##Extract betas and standar errors for meta-analaysis (see "Meta-analysis of main interactions" section, below) 
#Extract and store betas and standard errors for main effect of condition from the model "RiskInteractionMain2" for meta-analysis (see below) 
RiskConditionEffectb2 = coef(summary(RiskInteractionMain2))["ConditionStringMortality","Estimate"]
RiskConditionEffectse2 = coef(summary(RiskInteractionMain2))["ConditionStringMortality","Std. Error"]

#Extract and store betas and standard errors for main effect of child SES from the model "RiskInteractionMain2" for meta-analysis (see below) 
RiskSESEffectb2 = coef(summary(RiskInteractionMain2))["scale(ChildSES)","Estimate"]
RiskSESEffectse2 = coef(summary(RiskInteractionMain2))["scale(ChildSES)","Std. Error"]

#Extract and store betas and standard errors for interaction between condition and child SES from the model "RiskInteractionMain2" for meta-analysis (see below) 
RiskInteractionb2 = coef(summary(RiskInteractionMain2))["scale(ChildSES):ConditionStringMortality","Estimate"]
RiskInteractionse2 = coef(summary(RiskInteractionMain2))["scale(ChildSES):ConditionStringMortality","Std. Error"]

#Calculate directed p-values for interaction between child SES and condition on risk for replication 2
#To set the parameters for this test, run the code in the head of the script, under startup basics

#Extract t-values from summary of RiskInteractionMain2
TValueRisk2 = coef(summary(RiskInteractionMain2))["scale(ChildSES):ConditionStringMortality","t value"]
TValueRisk2

#Extract degrees of freedom from RiskInteractionMain2
dfRisk2 = df.residual(RiskInteractionMain2)
dfRisk2

#Calculate 1-tailed p-value for RiskInteractionMain2
p1tRisk2 = pt(TValueRisk2, df=dfRisk2, lower.tail=FALSE)
p1tRisk2

##Get directed test p-value
# Check if p1t is greater than or less than gamma/alpha and apply formula from Rice & Gaines (1994), box 1, as appropriate
pDirectedRisk2=ifelse(p1tRisk2 < gamma/alpha, (p1tRisk2*(1+delta/gamma)), ((1-p1tRisk2)*(1+gamma/delta)))  
pDirectedRisk2

#Discounting interaction model to check for predicted interaction between child SES and condition
DiscountingInteractionMain2 = lm(DiscountingTotal~scale(ChildSES)*ConditionString, data=Replication2)
Anova(DiscountingInteractionMain2, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteractionMain2)
confint(DiscountingInteractionMain2)

##Extract betas and standar errors for meta-analaysis (see "Meta-analysis of main interactions" section, below) 
#Extract and store betas and standard errors for main effect of condition from the model "DiscountingInteractionMain2" for meta-analysis (see below) 
DiscountingConditionEffectb2 = coef(summary(DiscountingInteractionMain2))["ConditionStringMortality","Estimate"]
DiscountingConditionEffectse2 = coef(summary(DiscountingInteractionMain2))["ConditionStringMortality","Std. Error"]

#Extract and store betas and standard errors for main effect of child SES from the model "DiscountingInteractionMain2" for meta-analysis (see below) 
DiscountingSESEffectb2 = coef(summary(DiscountingInteractionMain2))["scale(ChildSES)","Estimate"]
DiscountingSESEffectse2 = coef(summary(DiscountingInteractionMain2))["scale(ChildSES)","Std. Error"]

#Extract and store betas and standard errors for interaction between condition and child SES from the model "DiscountingInteractionMain2" for meta-analysis (see below) 
DiscountingInteractionb2 = coef(summary(DiscountingInteractionMain2))["scale(ChildSES):ConditionStringMortality","Estimate"]
DiscountingInteractionse2 = coef(summary(DiscountingInteractionMain2))["scale(ChildSES):ConditionStringMortality","Std. Error"]


#Calculate directed p-values for interaction between child SES and condition on delay discounting for replication 2
#To set the parameters for this test, run the code in the head of the script, under startup basics

#Extract t-values from summary of DiscountingInteractionMain2
TValueDiscounting2 = coef(summary(DiscountingInteractionMain2))["scale(ChildSES):ConditionStringMortality","t value"]
TValueDiscounting2

#Extract degrees of freedom from RiskInteractionMain2
dfDiscounting2 = df.residual(DiscountingInteractionMain2)
dfDiscounting2

#Calculate 1-tailed p-value for RiskInteractionMain2
p1tDiscounting2 = pt(TValueDiscounting2, df=dfDiscounting2, lower.tail=FALSE)
p1tDiscounting2

##Get directed test p-value
# Check if p1tDiscounting2 is greater than or less than gamma/alpha and apply formula from Rice & Gaines (1994), box 1, as appropriate
pDirectedDiscounting2=ifelse(p1tDiscounting2 < gamma/alpha, (p1tDiscounting2*(1+delta/gamma)), ((1-p1tDiscounting2)*(1+gamma/delta)))  
pDirectedDiscounting2

####Replication 2 - figures####

##Make data frame & models to plot interactions

#Find limits for "High" (>1 sd above mean) and "Low" (<1 sd below mean) SES
plusone2=mean(Replication2$ChildSES)+sd(Replication2$ChildSES)
minusone2=mean(Replication2$ChildSES)-sd(Replication2$ChildSES)
#Also store the mean ChildSES as a variable
MeanChildSES2 = mean(Replication2$ChildSES)

#Create a new data frame for comparing "High" and "Low" SES for risk
HighLow2=data.frame(Condition=c(1, 2, 1, 2), ChildSES=c(plusone2, plusone2, minusone2, minusone2))

#Add row names for "high" and "low" SES
row.names(HighLow2)=c("Condition 1, high", "Condition 2, high", "Condition 1, low", "Condition 2, low")
HighLow2$cSES=c("High", "High", "Low", "Low")

#Build model to get predictions and standard errors for risk
RiskModel2 = lm(RiskTotal~scale(ChildSES)+Condition, data=Replication2)
predict(RiskModel2, newdata=HighLow2, se.fit=TRUE)

#Build model to get predictions and standard errors for delay discounting
DiscountingModel2 = lm(DiscountingTotal~scale(ChildSES)+Condition, data=Replication2)
predict(DiscountingModel2, newdata=HighLow2, se.fit=TRUE)

#Add risk predictions to data frame
HighLow2$RiskPrediction=predict(RiskModel2, newdata=HighLow2, se.fit=TRUE)$fit

#Add risk standard errors to data frame
HighLow2$RiskSE=predict(RiskModel2, newdata=HighLow2, se.fit=TRUE)$se.fit

#Add delay discounting predictions to data frame
HighLow2$DiscountingPrediction=predict(DiscountingModel2, newdata=HighLow2, se.fit=TRUE)$fit

#Add delay discounting standard errors to data frame
HighLow2$DiscountingSE=predict(DiscountingModel2, newdata=HighLow2, se.fit=TRUE)$se.fit

#Add condition labels to data frame
HighLow2$Condition=c("Mortality", "Control", "Mortality", "Control")

#View the data frame
View(HighLow2)

#Bar plot of risk by condition for high and low child SES (Figure 2a)
Fig2a=ggplot(HighLow2, aes(x=Condition, y=RiskPrediction, fill=cSES)) +
  theme_bw() + labs (x = "Experimental condition", y = "Risk acceptance") +
  geom_bar(stat="identity", position="dodge", colour = "black") + ggtitle('(a)') +
  geom_errorbar(width=0.2, position=position_dodge(0.9), aes(ymin=RiskPrediction-RiskSE, ymax=RiskPrediction+RiskSE)) + 
  scale_fill_manual(values=c("black", "darkgrey"))

#Save Figure 2a to working directory
png("Figure2a.png", res=300, width=3, height=3, units="in")
plot(Fig2a)
dev.off()

#Show Figure 2a in the environment
Fig2a

#Bar plot of delay discounting by condition for high and low child SES (Figure 2b)
Fig2b=ggplot(HighLow2, aes(x=Condition, y=DiscountingPrediction, fill=cSES)) +
  theme_bw() + labs (x = "Experimental condition", y = "Delay discounting") +
  geom_bar(stat="identity", position="dodge", colour = "black") + ggtitle('(b)') +
  geom_errorbar(width=0.2, position=position_dodge(0.9), aes(ymin=DiscountingPrediction-DiscountingSE, ymax=DiscountingPrediction+DiscountingSE))+
  scale_fill_manual(values=c("black", "darkgrey"))

#Save Figure 2b to working directory
png("Figure2b.png", res=300, width=3, height=3, units="in")
plot(Fig2b)
dev.off()

#Show Figure 2b in the environment
Fig2b

####Replication 3 - main models####

#Risk, interaction model to check for predicted interaction between child SES and condition
RiskInteractionMain3 = lm(RiskTotal~scale(ChildSES)*ConditionString, data=Replication3)
Anova(RiskInteractionMain3, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteractionMain3)
confint(RiskInteractionMain3)

##Extract betas and standar errors for meta-analaysis (see "Meta-analysis of main interactions" section, below) 
#Extract and store betas and standard errors for main effect of condition from the model "RiskInteractionMain3" for meta-analysis (see below) 
RiskConditionEffectb3 = coef(summary(RiskInteractionMain3))["ConditionStringMortality","Estimate"]
RiskConditionEffectse3 = coef(summary(RiskInteractionMain3))["ConditionStringMortality","Std. Error"]

#Extract and store betas and standard errors for main effect of child SES from the model "RiskInteractionMain3" for meta-analysis (see below) 
RiskSESEffectb3 = coef(summary(RiskInteractionMain3))["scale(ChildSES)","Estimate"]
RiskSESEffectse3 = coef(summary(RiskInteractionMain3))["scale(ChildSES)","Std. Error"]

#Extract and store betas and standard errors for interaction between condition and child SES from the model "RiskInteractionMain3" for meta-analysis (see below) 
RiskInteractionb3 = coef(summary(RiskInteractionMain3))["scale(ChildSES):ConditionStringMortality","Estimate"]
RiskInteractionse3 = coef(summary(RiskInteractionMain3))["scale(ChildSES):ConditionStringMortality","Std. Error"]


#Calculate directed p-values for interaction between child SES and condition on risk for replication 3
#To set the parameters for this test, run the code in the head of the script, under startup basics

#Extract t-values from summary of RiskInteractionMain3
TValueRisk3 = coef(summary(RiskInteractionMain3))["scale(ChildSES):ConditionStringMortality","t value"]
TValueRisk3

#Extract degrees of freedom from RiskInteractionMain3
dfRisk3 = df.residual(RiskInteractionMain3)
dfRisk3

#Calculate 1-tailed p-value for RiskInteractionMain3
p1tRisk3 = pt(TValueRisk3, df=dfRisk3, lower.tail=FALSE)
p1tRisk3

##Get directed test p-value
#Check if p1t is greater than or less than gamma/alpha and apply formula from Rice & Gaines (1994), box 1, as appropriate
pDirectedRisk3=ifelse(p1tRisk3 < gamma/alpha, (p1tRisk3*(1+delta/gamma)), ((1-p1tRisk3)*(1+gamma/delta)))  
pDirectedRisk3

#Discounting interaction model to check for predicted interaction between child SES and condition
DiscountingInteractionMain3 = lm(DiscountingTotal~scale(ChildSES)*ConditionString, data=Replication3)
Anova(DiscountingInteractionMain3, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteractionMain3)
confint(DiscountingInteractionMain3)

##Extract betas and standar errors for meta-analaysis (see "Meta-analysis of main interactions" section, below) 
#Extract and store betas and standard errors for main effect of condition from the model "DiscountingInteractionMain3" for meta-analysis (see below) 
DiscountingConditionEffectb3 = coef(summary(DiscountingInteractionMain3))["ConditionStringMortality","Estimate"]
DiscountingConditionEffectse3 = coef(summary(DiscountingInteractionMain3))["ConditionStringMortality","Std. Error"]

#Extract and store betas and standard errors for main effect of child SES from the model "DiscountingInteractionMain3" for meta-analysis (see below) 
DiscountingSESEffectb3 = coef(summary(DiscountingInteractionMain3))["scale(ChildSES)","Estimate"]
DiscountingSESEffectse3 = coef(summary(DiscountingInteractionMain3))["scale(ChildSES)","Std. Error"]

#Extract and store betas and standard errors for interaction between condition and child SES from the model "DiscountingInteractionMain3" for meta-analysis (see below) 
DiscountingInteractionb3 = coef(summary(DiscountingInteractionMain3))["scale(ChildSES):ConditionStringMortality","Estimate"]
DiscountingInteractionse3 = coef(summary(DiscountingInteractionMain3))["scale(ChildSES):ConditionStringMortality","Std. Error"]


#Calculate directed p-values for interaction between child SES and condition on delay discounting for replication 3
#To set the parameters for this test, run the code in the head of the script, under startup basics

#Extract t-values from summary of DiscountingInteractionMain3
TValueDiscounting3 = coef(summary(DiscountingInteractionMain3))["scale(ChildSES):ConditionStringMortality","t value"]
TValueDiscounting3

#Extract degrees of freedom from RiskInteractionMain3
dfDiscounting3 = df.residual(DiscountingInteractionMain3)
dfDiscounting3

#Calculate 1-tailed p-value for RiskInteractionMain3
p1tDiscounting3 = pt(TValueDiscounting3, df=dfDiscounting3, lower.tail=FALSE)
p1tDiscounting3

##Get directed test p-value
# Check if p1tDiscounting3 is greater than or less than gamma/alpha and apply formula from Rice & Gaines (1994), box 1, as appropriate
pDirectedDiscounting3=ifelse(p1tDiscounting3 < gamma/alpha, (p1tDiscounting3*(1+delta/gamma)), ((1-p1tDiscounting3)*(1+gamma/delta)))  
pDirectedDiscounting3

####Replication 3 - figures####

##Make data frame & models to plot interactions

#Find limits for "High" (>1 sd above mean) and "Low" (<1 sd below mean) SES
plusone3=mean(Replication3$ChildSES)+sd(Replication3$ChildSES)
minusone3=mean(Replication3$ChildSES)-sd(Replication3$ChildSES)
#Also store the mean ChildSES as a variable
MeanChildSES3 = mean(Replication3$ChildSES)

#Create a new data frame for comparing "High" and "Low" SES for risk
HighLow3=data.frame(Condition=c(1, 2, 1, 2), ChildSES=c(plusone3, plusone3, minusone3, minusone3))

#Add row names for "high" and "low" SES
row.names(HighLow3)=c("Condition 1, high", "Condition 2, high", "Condition 1, low", "Condition 2, low")
HighLow3$cSES=c("High", "High", "Low", "Low")

#Build model to get predictions and standard errors for risk
RiskModel3 = lm(RiskTotal~scale(ChildSES)+Condition, data=Replication3)
predict(RiskModel3, newdata=HighLow3, se.fit=TRUE)

#Build model to get predictions and standard errors for delay discounting
DiscountingModel3 = lm(DiscountingTotal~scale(ChildSES)+Condition, data=Replication3)
predict(DiscountingModel3, newdata=HighLow3, se.fit=TRUE)

#Add risk predictions to data frame
HighLow3$RiskPrediction=predict(RiskModel3, newdata=HighLow3, se.fit=TRUE)$fit

#Add risk standard errors to data frame
HighLow3$RiskSE=predict(RiskModel3, newdata=HighLow3, se.fit=TRUE)$se.fit

#Add delay discounting predictions to data frame
HighLow3$DiscountingPrediction=predict(DiscountingModel3, newdata=HighLow3, se.fit=TRUE)$fit

#Add delay discounting standard errors to data frame
HighLow3$DiscountingSE=predict(DiscountingModel3, newdata=HighLow3, se.fit=TRUE)$se.fit

#Add condition labels to data frame
HighLow3$Condition=c("Mortality", "Control", "Mortality", "Control")

#View the data frame
View(HighLow3)

#Bar plot of risk by condition for high and low child SES (Figure 3a)
Fig3a=ggplot(HighLow3, aes(x=Condition, y=RiskPrediction, fill=cSES)) +
  theme_bw() + labs (x = "Experimental condition", y = "Risk acceptance") +
  geom_bar(stat="identity", position="dodge", colour = "black") + ggtitle('(a)') +
  geom_errorbar(width=0.2, position=position_dodge(0.9), aes(ymin=RiskPrediction-RiskSE, ymax=RiskPrediction+RiskSE)) + 
  scale_fill_manual(values=c("black", "darkgrey"))

#Save Figure 3a to working directory
png("Figure3a.png", res=300, width=3, height=3, units="in")
plot(Fig3a)
dev.off()

#Show Figure 3a in the environment
Fig3a

#Bar plot of delay discounting by condition for high and low child SES (Figure 3b)
Fig3b=ggplot(HighLow3, aes(x=Condition, y=DiscountingPrediction, fill=cSES)) +
  theme_bw() + labs (x = "Experimental condition", y = "Delay discounting") +
  geom_bar(stat="identity", position="dodge", colour = "black") + ggtitle('(b)') +
  geom_errorbar(width=0.2, position=position_dodge(0.9), aes(ymin=DiscountingPrediction-DiscountingSE, ymax=DiscountingPrediction+DiscountingSE))+
  scale_fill_manual(values=c("black", "darkgrey"))

#Save Figure 3b to working directory
png("Figure3b.png", res=300, width=3, height=3, units="in")
plot(Fig3b)
dev.off()

#Show Figure 3b in the environment
Fig3b

####Prime perception effects in replication 3####

#Check whether participant demographics influence their tendency to find the prime convincing
ConvincingModel = lm(Convincing~scale(ChildSES)+scale(AdultSES)+GenderString, data=Replication3)
Anova(ConvincingModel, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(ConvincingModel)

#Check whether participants' belief in the prime (when controlling for child SES) influenced their risk response
ConvincingRisk = lm(RiskTotal~scale(ChildSES)+Convincing, data=Replication3)
Anova(ConvincingRisk, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(ConvincingRisk)

#Check whether participants' belief in the prime (when controlling for child SES) influenced their discounting response
ConvincingDiscounting = lm(DiscountingTotal~scale(ChildSES)+Convincing, data=Replication3)
Anova(ConvincingDiscounting, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(ConvincingDiscounting)

#Check whether participants' responses to the prime (when controlling for child SES) influenced their risk response
PrimePerceptionsRisk = lm(RiskTotal~scale(ChildSES)+ Dangerous + Unsafe + Unpredictable + Uncertain + Aroused, data=Replication3)
Anova(PrimePerceptionsRisk, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(PrimePerceptionsRisk)

#Check whether participants' responses to the prime (when controlling for child SES) influenced their discounting response
PrimePerceptionsDiscounting = lm(DiscountingTotal~scale(ChildSES)+ Dangerous + Unsafe + Unpredictable + Uncertain + Aroused, data=Replication3)
Anova(PrimePerceptionsDiscounting, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(PrimePerceptionsDiscounting)

#Check whether participants' belief (after reading the prime) that the world is unsafe has an effect on risk in interaction with childhood SES
UnsafePerceptionsRisk = lm(RiskTotal~scale(ChildSES)*Unsafe, data=Replication3)
Anova(UnsafePerceptionsRisk, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(UnsafePerceptionsRisk)

#Check whether participants' belief (after reading the prime) that the world would become dangerous has an effect on discounting in interaction with childhood SES
DangerPerceptionsDiscounting = lm(DiscountingTotal~scale(ChildSES)*Dangerous, data=Replication3)
Anova(DangerPerceptionsDiscounting, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DangerPerceptionsDiscounting)

####Meta-analysis of main interactions####

##Meta-analysis (reported in main paper)##
slab=c("Replication 1", "Replication 2", "Replication 3")

# Main effect of condition
b.cond.risk=c(RiskConditionEffectb1,RiskConditionEffectb2, RiskConditionEffectb3)
se.cond.risk=c(RiskConditionEffectse1, RiskConditionEffectse2, RiskConditionEffectse3)
b.cond.delay=c(DiscountingConditionEffectb1, DiscountingConditionEffectb2, DiscountingConditionEffectb3)
se.cond.delay=c(DiscountingConditionEffectse1, DiscountingConditionEffectse2, DiscountingConditionEffectse3)

# Main effect of childhood SES
b.cses.risk=c(RiskSESEffectb1, RiskSESEffectb2, RiskSESEffectb3)
se.cses.risk=c(RiskSESEffectse1, RiskSESEffectse2, RiskSESEffectse3)
b.cses.delay=c(DiscountingSESEffectb1, DiscountingSESEffectb2, DiscountingSESEffectb3)
se.cses.delay=c(DiscountingSESEffectse1, DiscountingSESEffectse2, DiscountingSESEffectse3)

# Interaction
b.int.risk=c(RiskInteractionb1, RiskInteractionb2, RiskInteractionb3)
se.int.risk=c(RiskInteractionse1, RiskInteractionse2, RiskInteractionse3)
b.int.delay=c(DiscountingInteractionb1, DiscountingInteractionb2, DiscountingInteractionb3)
se.int.delay=c(DiscountingInteractionse1, DiscountingInteractionse2, DiscountingInteractionse3)

#Make Figure 4 and save it to the working directory
png("Figure4.png", res=300, width=6, height=6, units="in")
par(mfrow=c(3, 2))
par(mar=c(4.1, 2.1, 4.1, 2.1))
forest(rma(yi=b.cond.risk, sei=se.cond.risk), slab=slab, xlab="Parameter estimate")
text(-0.7, 5.5, "(a) Main effect of condition, risk acceptance",     pos=1)
forest(rma(yi=b.cond.delay, sei=se.cond.delay), slab=slab, xlab="Parameter estimate")
text(-0.8, 5.5, "(b) Main effect of condition, delay discounting",     pos=1)
forest(rma(yi=b.cses.risk, sei=se.cses.risk), slab=slab, xlab="Parameter estimate")
text(-0.35, 5.5, "(c) Main effect of childhood SES, risk acceptance",     pos=1)
forest(rma(yi=b.cses.delay, sei=se.cses.delay), slab=slab, xlab="Parameter estimate")
text(-0.15, 5.5, "(d) Main effect of childhood SES, delay discounting",     pos=1)
forest(rma(yi=b.int.risk, sei=se.int.risk), slab=slab, xlab="Parameter estimate")
text(-1.2, 5.5, "(e) Interaction, risk acceptance",     pos=1)
forest(rma(yi=b.int.delay, sei=se.int.delay), slab=slab, xlab="Parameter estimate")
text(-1.8, 5.5, "(f) Interaction, delay discounting",     pos=1)
dev.off()
rma(yi=b.cses.risk, sei=se.cses.risk)
rma(yi=b.cses.delay, sei=se.cses.delay)

####Supplementary analyses - descriptive stats####

## Descriptive statistics for replication 1
# Make data frame containing key summary variables for replication 1
keydata1 = select(Replication1, Age, RiskTotal, DiscountingTotal, ChildSES, AdultSES)
# Describe data by condition for replication 1 - for supplement table S1.
describeBy(keydata1, Replication1$ConditionString)

## Descriptive statistics for replication 2
# Make data frame containing key summary variables for replication 2
keydata2 = select(Replication2, Age, RiskTotal, DiscountingTotal, ChildSES, AdultSES)
# Describe data by condition for replication 2 - for supplement table S2.
describeBy(keydata2, Replication2$ConditionString)

## Descriptive statistics for replication 3
# Make data frame containing key summary variables for replication 3
keydata3 = select(Replication3, Age, RiskTotal, DiscountingTotal, ChildSES, AdultSES)
# Describe data by condition for replication 3 - for supplement table S3.
describeBy(keydata3, Replication3$ConditionString)

####Supplementary analyses - deprivation scores####

## Correlations between child and adult subjective SES scores and IMD deprivation scores from postcodes
##Replication 1
cor.test(Replication1$DeprivationScore, Replication1$ChildSES)
cor.test(Replication1$DeprivationScore, Replication1$AdultSES)

##Replication 2
cor.test(Replication2$DeprivationScore, Replication2$ChildSES)
cor.test(Replication2$DeprivationScore, Replication2$AdultSES)

##Replication 3
cor.test(Replication3$DeprivationScore, Replication3$ChildSES)
cor.test(Replication3$DeprivationScore, Replication3$AdultSES)
         
##Model for interaction between childhood postcode deprivation score and condition for risk in replication 1
# Supplement table S5
RiskInteractionIMD = lm(RiskTotal~DeprivationScore*ConditionString, data=Replication1)
Anova(RiskInteractionIMD, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteractionIMD)
confint(RiskInteractionIMD)

##Model for interaction between childhood postcode deprivation score and condition for delay discounting in replication 1
# Supplement table S6
DiscountingInteractionIMD = lm(DiscountingTotal~DeprivationScore*ConditionString, data=Replication1)
Anova(DiscountingInteractionIMD, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteractionIMD)
confint(DiscountingInteractionIMD)

####Supplementary analyses - interaction models####

##Replication 1

#Risk interaction model to check for 3-way interactions as in original study
#Table S7
RiskInteraction1 = lm(RiskTotal~(scale(ChildSES)+scale(AdultSES))*ConditionString*GenderString, data=Replication1)
Anova(RiskInteraction1, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteraction1)

#Delay discounting interaction model to check for 3-way interactions as in original study
#Table S8
DiscountingInteraction1 = lm(DiscountingTotal~(scale(ChildSES)+scale(AdultSES))*ConditionString*GenderString, data=Replication1)
Anova(DiscountingInteraction1, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteraction1)
#Extract estimate for meta-analysis
DiscountingInteractionb1 = coef(summary(DiscountingInteraction1))["ConditionStringMortality:GenderStringMale","Estimate"]
#Extract standard error
DiscountingInteractionse1 = coef(summary(DiscountingInteraction1))["ConditionStringMortality:GenderStringMale","Std. Error"]

##Replication 2

#Risk interaction model to check for 3-way interactions as in original study
#Table S9
RiskInteraction2 = lm(RiskTotal~(scale(ChildSES)+scale(AdultSES))*ConditionString*GenderString, data=Replication2)
Anova(RiskInteraction2, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteraction2)

#Delay discounting interaction model to check for 3-way interactions as in original study
#Table S10
DiscountingInteraction2 = lm(DiscountingTotal~(scale(ChildSES)+scale(AdultSES))*ConditionString*GenderString, data=Replication2)
Anova(DiscountingInteraction2, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteraction2)
#Extract estimate for meta-analysis
DiscountingInteractionb2 = coef(summary(DiscountingInteraction2))["ConditionStringMortality:GenderStringMale","Estimate"]
#Extract standard error
DiscountingInteractionse2 = coef(summary(DiscountingInteraction2))["ConditionStringMortality:GenderStringMale","Std. Error"]

##Replication 3

#Risk interaction model to check for 3-way interactions as in original study
#Table S11
RiskInteraction3 = lm(RiskTotal~(scale(ChildSES)+scale(AdultSES))*ConditionString*GenderString, data=Replication3)
Anova(RiskInteraction3, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(RiskInteraction3)

#Delay discounting interaction model to check for 3-way interactions as in original study
#Table S12
DiscountingInteraction3 = lm(DiscountingTotal~(scale(ChildSES)+scale(AdultSES))*ConditionString*GenderString, data=Replication3)
Anova(DiscountingInteraction3, contrasts=list(topic=contr.sum, sys=contr.sum), type =3)
summary(DiscountingInteraction3)
#Extract estimate for meta-analysis
DiscountingInteractionb3 = coef(summary(DiscountingInteraction3))["ConditionStringMortality:GenderStringMale","Estimate"]
#Extract standard error
DiscountingInteractionse3 = coef(summary(DiscountingInteraction3))["ConditionStringMortality:GenderStringMale","Std. Error"]

####Supplementary analyses - Meta-analysis#####

# Meta-analysis for interaction between sex and experimental condition for delay discounting
#Figure S1
b.sex.cond.delay=c(DiscountingInteractionb1, DiscountingInteractionb2, DiscountingInteractionb3)
se.sex.cond.delay=c(DiscountingInteractionse1, DiscountingInteractionse2, DiscountingInteractionse3)
forest(rma(yi=b.sex.cond.delay, sei=se.sex.cond.delay), slab=slab, xlab="Parameter estimate")
rma(yi=b.sex.cond.delay, sei=se.sex.cond.delay)
#Save FigureS1 to working directory and show it in environment
png("FigureS1.png", res=300, width=6, height=3, units="in")
forest(rma(yi=b.sex.cond.delay, sei=se.sex.cond.delay), slab=slab, xlab="Parameter estimate")
dev.off()

# Make Figure S2 (sex by prime interaction)
colnames(Replication1)
summary.data=summarise(group_by(Replication1, GenderString, ConditionString), mean=mean(DiscountingTotal), sd=sd(DiscountingTotal), n=length(DiscountingTotal))
summary.data=rbind(summary.data, summarise(group_by(Replication2, GenderString, ConditionString), mean=mean(DiscountingTotal), sd=sd(DiscountingTotal), n=length(DiscountingTotal)))
summary.data=rbind(summary.data, summarise(group_by(Replication3, GenderString, ConditionString), mean=mean(DiscountingTotal), sd=sd(DiscountingTotal), n=length(DiscountingTotal)))
summary.data$replication=c(rep("Replication 1", times=4), rep("Replication 2", times=4), rep("Replication 3", times=4))
summary.data$se=summary.data$sd/sqrt(summary.data$n)
summary.data
figS2=ggplot(summary.data, aes(x=replication, y=mean, fill=paste(GenderString, ConditionString))) + 
  geom_bar(stat="identity", position=position_dodge(), colour="black") + 
  geom_errorbar(aes(ymax=mean+se, ymin=mean-se), width=0.2, position=position_dodge(0.9)) + 
  scale_fill_manual(values=c("white", "black", "lightgrey", "darkgrey")) +
  theme_bw() + 
  theme(legend.title=element_blank()) + 
  ylab("Mean discounting score") + 
  xlab(" ") 

# Save Figure S2 to working directory
png("FigureS2.png", res=300, units="in", width=6, height=4)
figS2
dev.off()

#Show Figure S2 in environment
figS2

####Minimum detectable effect size calculations####

#Replication 1 - minimum detectable effect size
pwr.f2.test(u=1, v=68, sig.level = 0.05, power=0.8)

#Replication 2 - minimum detectable effect size
pwr.f2.test(u=1, v=155, sig.level = 0.05, power=0.8)

#Replication 3 - minimum detectable effect size
pwr.f2.test(u=1, v=158, sig.level = 0.05, power=0.8)
