1. Introduction

This file accompanies manuscript by Nichols et al. titled “Replicating the Effects of Auditory Religious Cues on Dishonest Behavior”, which investigate the effects of religious audiorty cues on dishonest behavior. The present file documents statistical analyses reported in the main text and supplementary materials.

2.0 Data set up

{
###############################################################
### General set up ###
  
# Load libraries
libs <- c("yarrr","dplyr","psych", "glmmADMB", "lsmeans", "nlme", "stargazer")
lapply(libs, require, character.only = TRUE)
# if those packages are not installed, uncomment the line below:
# install.packages(libs)

## Note that glmmADMB may not install using this method. Here is is a code to install glmmADMB 

#install.packages("R2admb")
#install.packages("glmmADMB", 
#    repos=c("http://glmmadmb.r-forge.r-project.org/repos",
#            getOption("repos")),
#    type="source") 


####################################################################################  
## Load data
dA<-read.csv("Nichols_et_al_data.csv", stringsAsFactors=FALSE)

# check missing values
dA[sapply(dA, function(dA) any(is.na(dA)))]

# Center & factor variables
dat <- dA %>% 
  mutate(id = as.factor(id),
         age.c = age - mean(age, na.rm=T),
         sex = as.factor(sex),
         relig = abs(relig-5), # reverse coding
         ritual = abs(ritual - 7), # reverse coding
         affil = as.factor(affil_cong), # factor religious affiliation (yes/no)
         impact = (deep + powerful)/2, # as per Lang et al. 2016
         tempo = (fast + abs(slow-7))/2,  # as per Lang et al. 2016
         positivity = (interesting + pleasant + exciting + relaxing + happy)/5,  # as per Lang et al. 2016
         negativity = (distressing + irritating + boring + sad)/4, # as per Lang et al. 2016
         difficult = abs(difficult_task-6),
         CTdiff = CT - CT_cheat,
         claim.b = (claim*459 + 0.5)/460,
         claim = claim*100
         )



# Re-order conditions to: religous, secular, noise, and control
dat$con[dat$con==4] <- 0 # make religious prime the reference category
dat$con[dat$con==1] <- 4
dat$con[dat$con==3] <- 1
dat$con[dat$con==4] <- 3

# treatment variable
dat$con<-factor(dat$con,levels= c(0,1,2,3),
                  labels = c("Religious", "Secular", "Noise","Control"))


# site variable - USA is the reference category
dat$sit <- dat$site
dat$site<-factor(dat$site,levels= c(1,2,3),
                 labels = c("USA", "CZ","JP"))


## Exclude PPTs ###########################################################################
# Exclude participants who guessed the hypotheses and played the game before
dat.full <- dat
dat <- dat[dat$include!=1,]

# Create site subsamples
USA <- dat[dat$site=="USA",]
CZ <- dat[dat$site=="CZ",]
JP <- dat[dat$site=="JP",]
sites <- c("USA", "CZ", "JP")
}

2.1. Histograms

Histogram of percentage of money claimed dishonestly

Histograms of claim per site

2.2. Distributions of main variables

Pirate plots of various variables by site these plots combine scatter plots (dots), density function (beans), and mean with 95% confidence interval (bar with a box). We recommend to zoom in the result to see the details. These four plots show between-site differences in the percentage of money claimed, age, religiousness, and the frequency of ritual behavior.

Pirate plots of the same variables by condition

Assess task difficulty across sites. First two plots show data for completion times during cheating trials and difficulty for full data set. The second row limits data span only to participants who cheated in less than 20% of the cases. One Czech outlier removed from the completion times.

2.3. Means and SD

## Report means and SD ###########################################################################
# THis approach is fof the sacred variable. You can just replce the word "sacred" 
# with e.g. "claim" to get all the same parameters for the claim variable
{
# Religious condition
mean0 <- mean(dat$claim[dat$con=="Religious"], na.rm = T) # mean
sd0 <- sd(dat$claim[dat$con=="Religious"], na.rm = T) # standard deviation
n0 <- length(dat$con[dat$con=="Religious" & !is.na(dat$con)]) # number of observations
se0 <- sd0/sqrt(n0) # standard error
lCI0 <- mean0 - (1.96*se0) # lower 95% CI
uCI0 <- mean0 + (1.96*se0) # upper 95% Ci

# Secular condition
mean1 <- mean(dat$claim[dat$con=="Secular"], na.rm = T) # mean
sd1 <- sd(dat$claim[dat$con=="Secular"], na.rm = T) # standard deviation
n1 <- length(dat$con[dat$con=="Secular" & !is.na(dat$con)]) # number of observations
se1 <- sd1/sqrt(n1) # standard error
lCI1 <- mean1 - (1.96*se1) # lower 95% CI
uCI1 <- mean1 + (1.96*se1) # upper 95% Ci

# Noise condition
mean2 <- mean(dat$claim[dat$con=="Noise"], na.rm = T) # mean
sd2 <- sd(dat$claim[dat$con=="Noise"], na.rm = T) # standard deviation
n2 <- length(dat$con[dat$con=="Noise" & !is.na(dat$con)]) # number of observations
se2 <- sd2/sqrt(n2) # standard error
lCI2 <- mean2 - (1.96*se2) # lower 95% CI
uCI2 <- mean2 + (1.96*se2) # upper 95% Ci

# Control condition
mean3 <- mean(dat$claim[dat$con=="Control"], na.rm = T) # mean
sd3 <- sd(dat$claim[dat$con=="Control"], na.rm = T) # standard deviation
n3 <- length(dat$con[dat$con=="Control" & !is.na(dat$con)]) # number of observations
se3 <- sd3/sqrt(n3) # standard error
lCI3 <- mean3 - (1.96*se3) # lower 95% CI
uCI3 <- mean3 + (1.96*se3) # upper 95% Ci

# compute Cohen's d

# Relig vs. Secular
d1 <- (mean0-mean1)/sqrt((sd0^2+sd1^2)/2)
# Relig vs. Noise
d2 <- (mean0-mean2)/sqrt((sd0^2+sd2^2)/2)
# Relig vs. Control
d3 <- (mean0-mean3)/sqrt((sd0^2+sd3^2)/2)
}

# Display the results for your table 1
# number of participants; mean value; standard deviation; 95% CI; and Cohen's d
rel <- c(n0,mean0,sd0,lCI0,uCI0)
sec <- c(n1,mean1,sd1,lCI1,uCI1,d1)
noi <- c(n2,mean2,sd2,lCI2,uCI2,d2)
con <- c(n3,mean3,sd3,lCI3,uCI3,d3)

print(rez <- as.data.frame(cbind(rel,sec,noi,con)), digits = 2)
##   rel     sec     noi    con
## 1 102 103.000 103.000 100.00
## 2  27  30.316  29.396  22.38
## 3  30  30.982  31.091  23.89
## 4  22  24.332  23.392  17.69
## 5  33  36.299  35.401  27.06
## 6 102  -0.098  -0.068   0.18
print(dat %>% group_by(site) %>% summarise(M = mean(age, na.rm = T), SD = sd(age, na.rm = T), NROW(id)))
## # A tibble: 3 x 4
##   site      M    SD `NROW(id)`
##   <fct> <dbl> <dbl>      <int>
## 1 USA    25.5 9.84         123
## 2 CZ     24.4 3.37         128
## 3 JP     19.8 0.902        157
print(dat %>% group_by(site) %>% summarise(M = mean(claim, na.rm = T), SD = sd(claim, na.rm = T), NROW(id)))
## # A tibble: 3 x 4
##   site      M    SD `NROW(id)`
##   <fct> <dbl> <dbl>      <int>
## 1 USA    40.6  34.3        123
## 2 CZ     11.7  13.7        128
## 3 JP     29.9  28.4        157

3.0 Main models

We approached the modeling of our dependent variable quite simply: while we could have used linear mixed-models (LMM) with site as the nesting factor, or generalized linear models (GLMs) accounting for the bounded nature of our data, we opted for simplicity and intepretability. So, we proceed with ordinary least squares models with the religious condition as a reference category (all, other conditions are compared to it) and the USA as the reference site (this choice was arbitrary, we are not really interested in the differences across sites). However, we also report Beta regression, which are appropriate for percentage data. Nevertheless, the results are qualitatively similar.

# Maniuplation check
summary(lm0 <- lm(sacred ~ con, data = dat))

summary(lm0.1 <- lm(religious ~ con, data = dat))
Estimates with 95% CI
sacred religious
(1) (2)
conSecular -0.85*** -1.32***
(-1.21, -0.49) (-1.69, -0.95)
conNoise -1.91*** -1.60***
(-2.27, -1.54) (-1.97, -1.23)
Constant 5.17*** 5.13***
(4.91, 5.43) (4.87, 5.39)
Observations 302 303
Note: x p<.1; * p<.05; ** p<.01; *** p<.001 
# First, simple condition effect holding site intercepts constant (no effect in the original study)
summary(lm1 <- lm(claim ~ con + sex + age.c + site, data = dat))

# Add interaction with religiosity (significant in the original study)
summary(lm2 <- lm(claim ~ con*relig + sex + age.c + site, data=dat))

# Add interaction with ritual frequency (trending in the original study)
summary(lm3 <- lm(claim ~ con*ritual + sex + age.c + site, data=dat))

# Add interaction with affiliation to the religion associated with our religious prime (not explored in the original study)
summary(lm4 <- lm(claim ~ con*affil + sex + age.c + site, data = dat))
Estimates with 95% CI
claim claim claim claim
(1) (2) (3) (4)
conSecular 3.45 -1.34 -2.70 -0.78
(-4.02, 10.91) (-14.26, 11.57) (-13.84, 8.44) (-9.31, 7.76)
conNoise 3.21 5.81 -5.82 -0.76
(-4.21, 10.62) (-7.31, 18.93) (-16.85, 5.22) (-9.04, 7.52)
conControl -4.91 -6.12 -11.41* -9.52*
(-12.37, 2.55) (-18.88, 6.63) (-22.13, -0.68) (-17.88, -1.17)
relig -2.12
(-6.77, 2.53)
ritual -4.19**
(-7.34, -1.03)
affil1 -26.09***
(-40.12, -12.06)
sex1 2.18 1.93 2.51 2.65
(-3.24, 7.60) (-3.62, 7.48) (-3.04, 8.07) (-2.82, 8.13)
age.c 0.49* 0.47x 0.52* 0.48*
(0.03, 0.96) (-0.004, 0.94) (0.06, 0.99) (0.02, 0.94)
siteCZ -28.56*** -28.90*** -29.74*** -30.86***
(-35.27, -21.85) (-35.84, -21.97) (-36.74, -22.73) (-37.84, -23.88)
siteJP -8.14* -9.23* -10.26* -11.32**
(-15.23, -1.06) (-16.67, -1.79) (-18.03, -2.49) (-18.72, -3.92)
conSecular:relig 2.97
(-3.40, 9.35)
conNoise:relig -1.40
(-7.94, 5.14)
conControl:relig 0.82
(-5.49, 7.13)
conSecular:ritual 3.46
(-1.01, 7.94)
conNoise:ritual 4.49x
(-0.10, 9.08)
conControl:ritual 3.89x
(-0.48, 8.26)
conSecular:affil1 23.83*
(5.64, 42.01)
conNoise:affil1 21.84*
(2.76, 40.92)
conControl:affil1 24.94*
(5.97, 43.91)
Constant 38.06*** 42.08*** 46.61*** 44.61***
(31.09, 45.02) (31.37, 52.80) (36.86, 56.35) (36.75, 52.47)
Observations 403 395 384 397
Note: x p<.1; * p<.05; ** p<.01; *** p<.001

Model diagnostics

4.0 Supplementary Material

4.1 Musical ratings

# Maniuplation check
summary(lm1 <- lm(positivity ~ con, data = dat))
## 
## Call:
## lm(formula = positivity ~ con, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7151 -0.3773 -0.1367  0.4633  2.8227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.53673    0.07632  33.237   <2e-16 ***
## conSecular   0.17842    0.10766   1.657   0.0986 .  
## conNoise    -1.15942    0.10822 -10.714   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7556 on 291 degrees of freedom
##   (114 observations deleted due to missingness)
## Multiple R-squared:  0.3827, Adjusted R-squared:  0.3784 
## F-statistic: 90.19 on 2 and 291 DF,  p-value: < 2.2e-16
lsmeans(lm1, pairwise~con)
## $lsmeans
##  con       lsmean     SE  df lower.CL upper.CL
##  Religious   2.54 0.0763 291     2.39     2.69
##  Secular     2.72 0.0759 291     2.57     2.86
##  Noise       1.38 0.0767 291     1.23     1.53
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast            estimate    SE  df t.ratio p.value
##  Religious - Secular   -0.178 0.108 291 -1.657  0.2236 
##  Religious - Noise      1.159 0.108 291 10.714  <.0001 
##  Secular - Noise        1.338 0.108 291 12.394  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(lm2 <- lm(negativity ~ con, data = dat))
## 
## Call:
## lm(formula = negativity ~ con, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5791 -0.3936 -0.1436  0.4209  2.1709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.82071    0.07303  24.932  < 2e-16 ***
## conSecular  -0.17714    0.10276  -1.724   0.0858 .  
## conNoise     0.75837    0.10354   7.325  2.3e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7266 on 295 degrees of freedom
##   (110 observations deleted due to missingness)
## Multiple R-squared:  0.2387, Adjusted R-squared:  0.2336 
## F-statistic: 46.26 on 2 and 295 DF,  p-value: < 2.2e-16
lsmeans(lm2, pairwise~con)
## $lsmeans
##  con       lsmean     SE  df lower.CL upper.CL
##  Religious   1.82 0.0730 295     1.68     1.96
##  Secular     1.64 0.0723 295     1.50     1.79
##  Noise       2.58 0.0734 295     2.43     2.72
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast            estimate    SE  df t.ratio p.value
##  Religious - Secular    0.177 0.103 295  1.724  0.1980 
##  Religious - Noise     -0.758 0.104 295 -7.325  <.0001 
##  Secular - Noise       -0.936 0.103 295 -9.080  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(lm3 <- lm(tempo ~ con, data = dat))
## 
## Call:
## lm(formula = tempo ~ con, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9899 -0.6050 -0.1050  0.5101  2.3950 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.60500    0.07868  33.110   <2e-16 ***
## conSecular   0.17718    0.11099   1.596    0.111    
## conNoise     1.38490    0.11155  12.415   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7868 on 297 degrees of freedom
##   (108 observations deleted due to missingness)
## Multiple R-squared:  0.3806, Adjusted R-squared:  0.3764 
## F-statistic: 91.25 on 2 and 297 DF,  p-value: < 2.2e-16
lsmeans(lm3, pairwise~con)
## $lsmeans
##  con       lsmean     SE  df lower.CL upper.CL
##  Religious   2.60 0.0787 297     2.45     2.76
##  Secular     2.78 0.0783 297     2.63     2.94
##  Noise       3.99 0.0791 297     3.83     4.15
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast            estimate    SE  df t.ratio p.value
##  Religious - Secular   -0.177 0.111 297  -1.596 0.2488 
##  Religious - Noise     -1.385 0.112 297 -12.415 <.0001 
##  Secular - Noise       -1.208 0.111 297 -10.854 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(lm4 <- lm(impact ~ con, data = dat))
## 
## Call:
## lm(formula = impact ~ con, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9100 -0.7626 -0.2300  0.7374  2.7374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.9100     0.1056  27.557  < 2e-16 ***
## conSecular   -0.1800     0.1493  -1.205    0.229    
## conNoise     -1.1474     0.1497  -7.664  2.6e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.056 on 296 degrees of freedom
##   (109 observations deleted due to missingness)
## Multiple R-squared:  0.1865, Adjusted R-squared:  0.181 
## F-statistic: 33.92 on 2 and 296 DF,  p-value: 5.433e-14
lsmeans(lm4, pairwise~con)
## $lsmeans
##  con       lsmean    SE  df lower.CL upper.CL
##  Religious   2.91 0.106 296     2.70     3.12
##  Secular     2.73 0.106 296     2.52     2.94
##  Noise       1.76 0.106 296     1.55     1.97
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast            estimate    SE  df t.ratio p.value
##  Religious - Secular    0.180 0.149 296 1.205   0.4510 
##  Religious - Noise      1.147 0.150 296 7.664   <.0001 
##  Secular - Noise        0.967 0.150 296 6.461   <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

4.2 Site-specific models

In this section, we replicate the site-specific models from Lang et al. (2016) and add site-specific models also for the other moderating variables.

# USA
summary(lm1 <- lm(claim ~ con*relig + sex + age.c, data=USA))

# CZ
summary(lm2 <- lm(claim ~ con*relig + sex + age.c, data=CZ))

# JP
summary(lm3 <- lm(claim ~ con*relig + sex + age.c, data=JP))
Estimates with 95% CI
claim claim claim
(1) (2) (3)
conSecular -19.26 -7.72 7.32
(-52.81, 14.29) (-22.57, 7.13) (-12.09, 26.74)
conNoise -4.00 5.86 12.50
(-38.26, 30.27) (-7.11, 18.84) (-9.51, 34.50)
conControl -24.96 6.25 0.63
(-56.63, 6.70) (-7.79, 20.28) (-19.06, 20.33)
relig -6.30 -1.09 0.52
(-15.81, 3.21) (-6.05, 3.87) (-8.98, 10.03)
sex1 5.41 -0.83 2.67
(-7.82, 18.64) (-6.13, 4.47) (-6.82, 12.16)
age.c 0.52 -0.25 1.73
(-0.13, 1.16) (-1.07, 0.56) (-3.52, 6.98)
conSecular:relig 7.88 7.30* -0.79
(-5.76, 21.51) (0.41, 14.19) (-14.58, 13.01)
conNoise:relig 1.50 -0.63 -4.55
(-13.05, 16.05) (-6.88, 5.62) (-18.40, 9.30)
conControl:relig 6.71 -0.44 -6.03
(-6.27, 19.69) (-7.10, 6.23) (-18.88, 6.82)
Constant 53.92*** 9.98x 32.10**
(30.37, 77.47) (-0.14, 20.10) (10.19, 54.01)
Observations 122 118 155
Note: x p<.1; * p<.05; ** p<.01; *** p<.001 
# USA
summary(lm1 <- lm(claim ~ con*ritual + sex + age.c, data=USA))

# CZ
summary(lm2 <- lm(claim ~ con*ritual + sex + age.c, data=CZ))

# JP
summary(lm3 <- lm(claim ~ con*ritual + sex + age.c, data=JP))
Estimates with 95% CI
claim claim claim
(1) (2) (3)
conSecular -17.38 -3.69 2.94
(-48.57, 13.81) (-14.03, 6.65) (-15.01, 20.89)
conNoise -22.25 2.70 -3.69
(-54.63, 10.14) (-7.49, 12.88) (-21.52, 14.13)
conControl -32.52* 4.45 -10.45
(-60.79, -4.25) (-6.29, 15.19) (-27.35, 6.45)
ritual -7.30* -0.84 -5.65
(-13.35, -1.26) (-4.15, 2.48) (-12.94, 1.63)
sex1 7.28 1.13 2.59
(-5.71, 20.28) (-4.01, 6.27) (-7.24, 12.42)
age.c 0.56x -0.32 3.31
(-0.08, 1.20) (-1.13, 0.48) (-2.90, 9.51)
conSecular:ritual 5.39 6.48** 3.38
(-3.95, 14.74) (1.74, 11.21) (-7.21, 13.96)
conNoise:ritual 8.11 0.95 7.03
(-2.06, 18.29) (-3.54, 5.43) (-4.67, 18.74)
conControl:ritual 8.02x 0.42 3.17
(-0.50, 16.54) (-4.29, 5.12) (-8.09, 14.43)
Constant 59.75*** 8.92* 44.42***
(39.10, 80.40) (1.26, 16.58) (20.24, 68.59)
Observations 122 121 141
Note: x p<.1; * p<.05; ** p<.01; *** p<.001 
# USA
summary(lm1 <- lm(claim ~ con*affil + sex + age.c, data=USA))

# CZ
summary(lm2 <- lm(claim ~ con*affil + sex + age.c, data=CZ))

# JP
summary(lm3 <- lm(claim ~ con*affil + sex + age.c, data=JP))
Estimates with 95% CI
claim claim claim
(1) (2) (3)
conSecular -19.76x 1.09 7.64
(-41.59, 2.06) (-6.39, 8.57) (-5.95, 21.23)
conNoise -22.05* 4.45 6.93
(-43.07, -1.03) (-2.94, 11.83) (-6.25, 20.12)
conControl -32.62** 5.55 -7.42
(-53.76, -11.47) (-2.00, 13.10) (-20.48, 5.64)
affil1 -48.78*** -4.04 -13.37
(-71.92, -25.63) (-20.25, 12.17) (-54.03, 27.29)
sex1 8.07 -0.32 1.90
(-4.31, 20.44) (-5.31, 4.66) (-7.62, 11.42)
age.c 0.53x -0.38 1.85
(-0.08, 1.14) (-1.14, 0.38) (-3.22, 6.92)
conSecular:affil1 41.15* 30.93** -1.40
(8.44, 73.87) (10.91, 50.96) (-50.23, 47.44)
conNoise:affil1 52.75** 2.12 -6.61
(18.22, 87.29) (-17.60, 21.84) (-59.47, 46.25)
conControl:affil1 52.08** -1.65 11.36
(19.36, 84.79) (-23.12, 19.82) (-41.38, 64.10)
Constant 60.13*** 8.52** 34.32***
(44.63, 75.63) (2.73, 14.32) (14.73, 53.91)
Observations 122 120 155
Note: x p<.1; * p<.05; ** p<.01; *** p<.001

4.3 Beta regression

# transform data to accomodate them for the need of the Beta regression (Smithson & Verkulien, 2006)

# create separata data sets, exluding missing values pairwise
dat1 <- na.omit(subset(dat, select = c(claim.b, con, relig, sex, age.c, site)))
dat2 <- na.omit(subset(dat, select = c(claim.b, con, ritual, sex, age.c, site)))
dat3 <- na.omit(subset(dat, select = c(claim.b, con, affil, sex, age.c, site)))


summary(glm1 <- glmmadmb(claim.b ~ con*relig + sex + age.c + site, family = 'beta', data=dat1))
summary(glm2 <- glmmadmb(claim.b ~ con*ritual + sex + age.c + site, family = 'beta', data=dat2))
summary(glm3 <- glmmadmb(claim.b ~ con*affil + sex + age.c + site, family = 'beta', data=dat3))



cf <- summary(glm1)$coefficients[,1]
se <- sqrt(diag(vcov(glm1)))
p <- summary(glm1)$coefficients[,4]
ci <- (cbind(est = (plogis(cf)-0.5), LL = (plogis(cf - 1.96 * se)-0.5),
             UL = (plogis(cf + 1.96 * se)-0.5),p))
ci[1,] <- (ci[1,]+0.5)
ci1 <- ci

cf <- summary(glm2)$coefficients[,1]
se <- sqrt(diag(vcov(glm2)))
p <- summary(glm2)$coefficients[,4]
ci <- (cbind(est = (plogis(cf)-0.5), LL = (plogis(cf - 1.96 * se)-0.5),
             UL = (plogis(cf + 1.96 * se)-0.5),p))
ci[1,] <- (ci[1,]+0.5)
ci2 <- ci

cf <- summary(glm3)$coefficients[,1]
se <- sqrt(diag(vcov(glm3)))
p <- summary(glm3)$coefficients[,4]
ci <- (cbind(est = (plogis(cf)-0.5), LL = (plogis(cf - 1.96 * se)-0.5),
             UL = (plogis(cf + 1.96 * se)-0.5),p))
ci[1,] <- (ci[1,]+0.5)
ci3 <- ci
Estimates with 95% CI
est LL UL p
(Intercept) 0.48 0.37 0.60 1.28
conSecular 0.02 -0.12 0.15 0.82
conNoise 0.07 -0.07 0.20 0.32
conControl -0.06 -0.19 0.07 0.38
relig -0.03 -0.08 0.02 0.27
sex1 0.02 -0.04 0.08 0.44
age.c 0.01 0.003 0.01 0.003
siteCZ -0.25 -0.30 -0.19 0
siteJP -0.08 -0.15 0.002 0.05
conSecular:relig 0.03 -0.04 0.10 0.37
conNoise:relig -0.004 -0.07 0.06 0.91
conControl:relig 0.02 -0.05 0.09 0.59
Estimates with 95% CI
est LL UL p
(Intercept) 0.53 0.43 0.64 1.02
conSecular 0.02 -0.10 0.13 0.80
conNoise -0.05 -0.16 0.07 0.44
conControl -0.10 -0.20 0.01 0.08
ritual -0.05 -0.08 -0.01 0.01
sex1 0.03 -0.03 0.09 0.31
age.c 0.01 0.003 0.01 0.001
siteCZ -0.26 -0.31 -0.20 0
siteJP -0.10 -0.17 -0.01 0.02
conSecular:ritual 0.03 -0.02 0.08 0.27
conNoise:ritual 0.05 0.003 0.10 0.04
conControl:ritual 0.04 -0.01 0.09 0.10
Estimates with 95% CI
est LL UL p
(Intercept) 0.51 0.43 0.60 1.25
conSecular 0.03 -0.06 0.12 0.49
conNoise 0.002 -0.08 0.09 0.96
conControl -0.08 -0.17 0.01 0.07
affil1 -0.27 -0.36 -0.14 0.0002
sex1 0.03 -0.03 0.09 0.28
age.c 0.01 0.003 0.01 0.002
siteCZ -0.28 -0.33 -0.22 0
siteJP -0.11 -0.18 -0.03 0.01
conSecular:affil1 0.20 0.01 0.34 0.04
conNoise:affil1 0.28 0.10 0.39 0.003
conControl:affil1 0.24 0.05 0.37 0.02

4.4 Linear mixed models

summary(lm1 <- lme(claim ~ con*relig + sex + age.c, random = ~1|site, data = dat, na.action = na.omit))

summary(lm2 <- lme(claim ~ con*ritual + sex + age.c, random = ~1|site, data = dat, na.action = na.omit))

summary(lm3 <- lme(claim ~ con*affil + sex + age.c, random = ~1|site, data = dat, na.action = na.omit))
Estimates with 95% CI
claim claim claim
(1) (2) (3)
conSecular -1.31 -2.73 -0.80
(-14.22, 11.61) (-13.87, 8.41) (-9.33, 7.74)
conNoise 5.79 -5.85 -0.77
(-7.33, 18.91) (-16.89, 5.19) (-9.05, 7.51)
conControl -6.10 -11.38* -9.53*
(-18.85, 6.66) (-22.11, -0.66) (-17.88, -1.18)
relig -2.09
(-6.74, 2.55)
ritual -4.14*
(-7.30, -0.99)
affil1 -25.86***
(-39.89, -11.84)
sex1 1.91 2.50 2.62
(-3.64, 7.45) (-3.06, 8.05) (-2.85, 8.10)
age.c 0.47x 0.52* 0.48*
(-0.003, 0.94) (0.06, 0.99) (0.02, 0.94)
conSecular:relig 2.93
(-3.45, 9.30)
conNoise:relig -1.41
(-7.96, 5.13)
conControl:relig 0.79
(-5.51, 7.10)
conSecular:ritual 3.47
(-1.01, 7.94)
conNoise:ritual 4.49x
(-0.10, 9.08)
conControl:ritual 3.88x
(-0.49, 8.25)
conSecular:affil1 23.73*
(5.55, 41.92)
conNoise:affil1 21.69*
(2.61, 40.77)
conControl:affil1 24.90*
(5.93, 43.87)
Constant 29.38** 33.23*** 30.55**
(10.26, 48.50) (14.34, 52.12) (11.88, 49.22)
Observations 395 384 397
Note: x p<.1; * p<.05; ** p<.01; *** p<.001

4.5 Adding control variables

summary(lmx <- lm(difficult~site + age.c + sex, data = dat))
lsmeans(lmx, pairwise~site)


# Add interaction with religiosity (significant in the original study)
summary(lm1 <- lm(claim ~ con*relig + sex + age.c + difficult + CTdiff + site, data=dat))
summary(lm2 <- lm(claim ~ con*relig + sex + age.c + impact + tempo + positivity + negativity + difficult + CTdiff + site, data=dat))

# Add interaction with ritual frequency (trending in the original study)
summary(lm3 <- lm(claim ~ con*ritual + sex + age.c + difficult + CTdiff+ site, data=dat))
summary(lm4 <- lm(claim ~ con*ritual + sex + age.c + impact + tempo + positivity + negativity + difficult + CTdiff + site, data=dat))

# Add interaction with affiliation to the religion associated with our religious prime (not explored in the original study)
summary(lm5 <- lm(claim ~ con*affil + sex + age.c + difficult + CTdiff + site, data=dat))
summary(lm6 <- lm(claim ~ con*affil + sex + age.c + impact + tempo + positivity + negativity + difficult + CTdiff + site, data=dat))
Estimates with 95% CI
claim claim claim claim claim claim
(1) (2) (3) (4) (5) (6)
conSecular -1.06 1.93 -3.53 -1.65 -0.24 0.63
(-13.37, 11.26) (-10.99, 14.85) (-14.11, 7.04) (-12.69, 9.38) (-8.43, 7.94) (-7.99, 9.25)
conNoise 3.39 3.40 -5.56 -5.40 -0.05 -0.98
(-9.14, 15.92) (-12.57, 19.37) (-16.08, 4.96) (-19.59, 8.79) (-8.01, 7.91) (-12.70, 10.74)
conControl -5.98 -11.75* -8.22*
(-18.17, 6.20) (-21.93, -1.57) (-16.32, -0.13)
relig -2.54 -2.74
(-7.01, 1.93) (-7.45, 1.98)
ritual -4.53** -4.57**
(-7.56, -1.49) (-7.72, -1.42)
affil1 -22.35** -27.15***
(-36.37, -8.33) (-41.74, -12.56)
sex1 3.55 4.46 3.97 4.96 3.99 5.63x
(-1.82, 8.92) (-2.11, 11.03) (-1.36, 9.30) (-1.52, 11.44) (-1.32, 9.31) (-0.83, 12.09)
age.c 0.52* 0.62* 0.57* 0.65* 0.53* 0.61*
(0.07, 0.98) (0.09, 1.15) (0.13, 1.01) (0.14, 1.16) (0.09, 0.97) (0.10, 1.12)
impact 0.57 0.73 -0.25
(-3.45, 4.59) (-3.20, 4.66) (-4.19, 3.70)
tempo 0.47 0.54 0.90
(-3.70, 4.64) (-3.59, 4.67) (-3.18, 4.98)
positivity 2.99 2.68 4.36
(-2.49, 8.47) (-2.68, 8.03) (-0.99, 9.70)
negativity 5.26* 4.97* 6.09**
(0.69, 9.83) (0.45, 9.48) (1.58, 10.61)
difficult -7.71*** -6.16*** -8.16*** -6.69*** -7.58*** -5.76***
(-10.45, -4.98) (-9.58, -2.73) (-10.89, -5.43) (-10.11, -3.26) (-10.28, -4.89) (-9.11, -2.40)
CTdiff 5.39** 22.39*** 5.36** 21.71*** 5.54*** 23.24***
(2.08, 8.71) (13.62, 31.17) (2.10, 8.61) (13.12, 30.30) (2.27, 8.82) (14.70, 31.78)
siteCZ -23.53*** -24.24*** -24.43*** -25.67*** -25.15*** -26.12***
(-30.46, -16.61) (-34.39, -14.09) (-31.35, -17.52) (-35.68, -15.67) (-32.10, -18.21) (-36.06, -16.19)
siteJP -5.28 -5.37 -6.60x -7.34 -6.72x -6.70
(-12.66, 2.10) (-15.12, 4.39) (-14.16, 0.97) (-17.30, 2.61) (-14.06, 0.62) (-16.29, 2.89)
conSecular:relig 2.43 1.54
(-3.66, 8.53) (-4.78, 7.86)
conNoise:relig -0.23 0.88
(-6.51, 6.04) (-5.83, 7.60)
conControl:relig 0.56
(-5.55, 6.66)
conSecular:ritual 3.49 3.17
(-0.78, 7.76) (-1.21, 7.56)
conNoise:ritual 3.99x 4.79*
(-0.41, 8.38) (0.12, 9.47)
conControl:ritual 3.71x
(-0.48, 7.89)
conSecular:affil1 19.20* 22.01*
(1.32, 37.07) (3.59, 40.44)
conNoise:affil1 17.09x 29.13**
(-1.66, 35.84) (8.96, 49.29)
conControl:affil1 18.46x
(-0.18, 37.10)
Constant 61.30*** 40.20*** 67.23*** 47.10*** 61.65*** 36.40**
(49.09, 73.51) (17.16, 63.23) (55.67, 78.79) (23.88, 70.33) (51.88, 71.42) (14.07, 58.73)
Observations 389 280 379 273 391 282
Note: x p<.1; * p<.05; ** p<.01; *** p<.001

4.6 Full-data models

# First, explore the differences in the incorrectly claimed higher-paying sides between the excluded and included particiapants

summary(lm1 <- lm(claim ~ include + sex + age.c + site, data=dat.full))

## look only at 100% incorrect claims
dat.full$claim.bin <- dat.full$claim
dat.full$claim.bin[dat.full$claim.bin < 100] <- 0
dat.full$claim.bin[dat.full$claim.bin == 100] <- 1

summary(lm2 <- glm(claim.bin ~ include + sex + age.c + site, data=dat.full, family = "binomial"))
exp(lm2$coefficients[2])


# run basic models on full sample

summary(lm3 <- lm(claim ~ con + sex + age.c + site, data = dat.full))

summary(lm4 <- lm(claim ~ con*relig + sex + age.c + site, data=dat.full))

summary(lm5 <- lm(claim ~ con*ritual + sex + age.c + site, data=dat.full))

summary(lm6 <- lm(claim ~ con*affil + sex + age.c + site, data=dat.full))
Estimates with 95% CI
claim claim claim claim
(1) (2) (3) (4)
conSecular 3.04 -1.79 -2.89 -0.34
(-4.06, 10.13) (-14.09, 10.51) (-13.49, 7.71) (-8.47, 7.80)
conNoise 1.26 3.16 -6.28 -1.90
(-5.78, 8.31) (-9.17, 15.49) (-17.03, 4.47) (-9.76, 5.96)
conControl -5.80 -5.81 -9.72x -9.08*
(-12.81, 1.20) (-17.68, 6.05) (-19.86, 0.42) (-16.93, -1.24)
relig -1.31
(-5.67, 3.04)
ritual -2.86x
(-5.85, 0.14)
affil1 -19.98**
(-33.23, -6.73)
sex1 4.31 3.91 4.65x 4.63x
(-0.82, 9.44) (-1.35, 9.17) (-0.63, 9.92) (-0.59, 9.85)
age.c 0.55* 0.54* 0.57* 0.55*
(0.11, 0.99) (0.09, 0.98) (0.13, 1.02) (0.10, 0.99)
siteCZ -29.75*** -30.05*** -30.52*** -31.99***
(-35.93, -23.57) (-36.46, -23.63) (-37.03, -24.01) (-38.52, -25.47)
siteJP -10.64** -11.11** -12.07** -13.37***
(-17.45, -3.83) (-18.26, -3.96) (-19.54, -4.60) (-20.50, -6.23)
conSecular:relig 2.90
(-3.16, 8.96)
conNoise:relig -1.09
(-7.28, 5.10)
conControl:relig 0.03
(-5.84, 5.90)
conSecular:ritual 3.15
(-1.07, 7.37)
conNoise:ritual 3.49
(-0.90, 7.88)
conControl:ritual 2.15
(-1.95, 6.26)
conSecular:affil1 18.13*
(0.78, 35.49)
conNoise:affil1 16.10x
(-2.42, 34.62)
conControl:affil1 16.88x
(-1.29, 35.04)
Constant 40.24*** 42.85*** 46.35*** 45.66***
(33.69, 46.80) (32.72, 52.98) (36.86, 55.84) (38.16, 53.17)
Observations 455 447 436 449
Note: x p<.1; * p<.05; ** p<.01; *** p<.001