Summary and Reactions

Summary

Nichols et al. (2020) attempted to replicate a study by Lang et al. to examine the effects of religious music on dishonest behaviour. Specifically, their aim was to replicate the interaction effect between religiosity and condition (i.e., religious music only reduces dishonest behaviour in religious people, and no other individuals).

Lang and colleagues proposed that only religious music affected highly religious people because the association between religion and auditory religious cues subconsciously influences them to be honest. Nichols and colleagues added on to Lang et al.’s experiment by examining individuals from Japan (sometimes religious), the Czech Republic (not really religious) or the USA (super religious), and using the Dots task instead of the Matrix task.

The justification for replicating this study was to see if subliminal religious cues could induce normative behaviour, thus gain a deeper understanding as to how human behaviour can be influenced.

Nichols et al. randomly allocated individuals to one of four conditions: no music (control), religious music, secular music (i.e., not related to religion), or white noise. Depending on their condition, music would play during the Dots task. The Dots task consisted of participants indicating whether the left or right side of the computer screen had more dots by pressing either the M or Z key on a keyboard. Participants were told that there was more monetary incentive if the right side had more dots as opposed to the left side. Thus, this allowed participants to cheat by indicating that the right side had more dots when it did not.

Nichols and colleagues were unable to replicate the interaction effect between religiosity and condition, as Lang et al. did. However, religious music did decrease dishonest behaviour for individuals who were frequently engaged in religious practices or were affiliated with a religious institution.

The results support the theory that religious music only reduced dishonest behaviour in individuals who had a strong association between religion and religious cues.

Reactions

I was surprised to see that… Nichols et al. were unable to replicate the interaction effect between religiosity and condition (i.e., unable to find that higher self-reported religiosity did not reduce dishonest behaviour when religious music was played) but found that “sacred” cues affected those who practice religious rituals/ are affiliated with a religious organisation. It makes sense that the latter finding was observed, since those with a stronger association between religious music and religion due to rituals will feel influenced to act honestly. However, I would have assumed that the prior result would have also been found because higher religiosity would indicate that they think more about religion than those who did not rate as high.

The most interesting part of the paper… was the decision to add Japan as a site, in addition to the Czech Republic and the USA. I was unaware of the religious practices of both Japan and the Czech Republic and it was interesting to read the differences (e.g., Japan is focused on religious practices while the USA is focused on personal beliefs). It was also interesting to read that one of the reasons that Nichols et al. believe they could not replicate the interaction effect between religiosity and condition was because Japan had greater religious diversity. It made me wonder if another East Asian sample would have similar results, or have similar religious practices to Japan.

Future work in this area needs to consider… any possible links between religiosity and morality. Nichols et al. admit that non religious and religious people have differing beliefs on whether religion is necessary to be “a good person”. It could very well be that religious people who do believe that religion is needed for their morals may be more affected by religious cues more than religious or non religious people who do not think that religions dictate their morality. Moreover, I think it would be interesting to ask people who do not believe that religion dictates their morality what they think dictates their morality and see if those can be influenced, similar to religious cues.

Verification

I undertook the verification task with Victor, Jessica and Sam. Against Jenny and Kate’s advice, we decided to divide and conquer the plots and table we were to reproduce.

We had to first recreate the demographic statistics, figures and tables:

Data were collected at universities across three sites: the USA, the Czech Republic, and Japan. A total of 460 (228 females) adults were randomly assigned to one of four conditions: religious, secular, white noise, or control (no sound). Due to self-reported suspicion on the goals of the experiment and previous experience with studies using the Dots Game, a total of 52 participants were excluded from the analysis. In support of this decision, supplementary analyses of the incorrectly claimed higher-paying sides in the Dots game showed that the excluded participants had 5.5 higher odds of incorrectly claiming all money from the higher-paying side (100% claimed) compared to participants in the included sample (see S1 Table G in S1 File for the analysis of the full sample). Five additional participants had missing data on crucial variables such as sex and age. Data included in the analyses therefore comprise: 123 American participants (Mage = 25.5, SD = 9.8), 128 Czech participants (Mage = 24.4, SD = 3.4), and 157 Japanese participants (Mage = 19.8, SD = 0.9). Within the four experimental conditions, there were: 100 participants in the control condition, 103 participants in the white-noise condition, 103 participants in the secular condition, and 102 participants in the religious condition.

Figure 1

Figure 2

Table 1

I took on Figure 1 (with Victor’s assistance), Victor tackled Figure 2, and Sam and Jessica worked on Table 1.

We created a codebook to make sense of the variables in the csv.

include: 0 for included, 1 for excluded
site: 1 for USA, 2 for Czech Republic, 3 for Japan
id: id of participant
cond: 1 for control, 2 for white noise, 3 for secular, 4 for religious
claimpercent: percentage claimed incorrectly 
claimmoney: how much money was made from the task
sex: 0 = female, 1 = male, NA = missing data
age: age
religiosity: 1 - 5, least to most self-rating of religiosity
affil: 0 nonaffiliated, 1 affiliated
religion: religion
ritual: self-rated ritual frequency

This would prove to be helpful, especially when we would find out that some variables (e.g., condition) were inconsistent with those used for the plots.

Before recreating the demographic statistics and plots, we had to load the correct libraries, and read the csv.

library(tidyverse)
library(dplyr) #used for cleaning up data
library(janitor) #used to read csv
library(gt) #used to create the table
library(ggplot2) #used to plot graphs for Fig 1 and Fig 2
library(ggdist) #used for raincloud plot in Fig 1
library(gridExtra) #combines plots for Fig 1 and Fig 2
# Data for demographic statistics
demographics <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")

# Data for Fig 1
data1 <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")

# Data for Fig 2
data2 <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")

# Data for Table 1
data3 <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")

Then, we could rename the variables, to align with the code used for the plots/table.

data1 <- data1 %>%  
  filter(include == 0) %>% 
  rename(cond = con,
         claimpercent = claim,
         claimmoney = moneyclaim,
         CT_practice = `completion time (practice included)`,
         CT_payments = `completion time (payments only)`,
         religiosity = relig,
         religion = Religion)

For our individual plots, we will also use the select() function from the dplyr package to use the relevant variables.

Demographic statistics

We needed to calculate the number of males and females, as well as the mean age and SD, and number of each site. We also had to relabel the sex and sites with demoSex$sex[demoSex$sex == 0] <- "Female". gt() was used to visualise the data in a table. Additionally, na.rm = TRUE skips all the NA values, fmt_number() changes the values to 1 decimal place and use_seps() removes any digit separators.

demographics <- demographics %>% 
  select(sex, site, include, age) 

demoSex <- demographics %>% 
  count(sex)

demoSite <- demographics %>% 
  filter(include == 0) 

demoSite <- demoSite %>% 
  group_by(site) %>% 
  summarise(mean = mean(age, na.rm = TRUE),
            sd = sd(age, na.rm = TRUE),
            n = n())

demoSex$sex[demoSex$sex == 0] <- "Female"
demoSex$sex[demoSex$sex == 1] <- "Male"

demoSite$site[demoSite$site == 1] <- "USA"
demoSite$site[demoSite$site == 2] <- "Czech Republic"
demoSite$site[demoSite$site == 3] <- "Japan"

gt(demoSite) %>%
  fmt_number(
    columns = c("mean", "sd"),
    decimals = 1,
    use_seps = FALSE)
site mean sd n
USA 25.5 9.8 123
Czech Republic 24.4 3.4 128
Japan 19.8 0.9 157
 
gt(demoSex)
sex n
Female 228
Male 227
NA 5

Figure 1

Before recreating Figure 1, I used the select() function to create a smaller table with the relevant variables. The mutate function was used to make claimpercent a percentage instead of a decimal.

data1 <- data1 %>% 
  select(site, claimpercent, cond) %>% 
  mutate(claimpercent = claimpercent * 100)

My first step was to figure out what kind of plot Figure 1 was so I could create it using the ggplot2 package.

My first course of action was to check the R Markdown file provided by Nichols and colleagues. I was unable to knit the file because of outdated packages, so I instead tried running the code line by line to see if the code related to Figure 1. However, Victor was able to run the whole R Markdown File and discovered that the file only contained new plots that were not in the paper.

However, we figured out that the columns in Figure 1 was just geom_col() combined with another function, Victor and I tried to find if we could combine geom_area() with with geom_col().

fig1_areacol <- ggplot(data = data1) +
  geom_col(aes(x = claimpercent, y = cond)) +
  geom_area(aes(x = claimpercent, y = cond)) +
  scale_x_continuous(
    name = "Percent Claimed",
    labels = c(0, 20, 40, 60, 80, 100),
    breaks = c(0, 20, 40, 60, 80, 100)) +
  scale_y_discrete(name = NULL) +
  ggtitle(label = "Data By Condition")

plot(fig1_areacol)

Clearly, this was not the right way to go about it. After asking for help on Slack, Jenny S recommended to Google “raincloud plots”.

I came across this blog which showed that ggdist::stat_halfeye() could be used to create the density plot.

fig1_raincloud <- ggplot(data1, aes(x = cond, y = claimpercent)) +
  ggdist::stat_halfeye(
    adjust = .5,
    width = .5,
    .width = 0,
    point_colour = NA
  ) +
  coord_flip()

plot(fig1_raincloud)

Now we have the density plot! I then combined ggdist::stat_halfeye() with geom_col().

fig1_densitycol <- ggplot(data1, aes(x = cond, y = claimpercent)) +
  geom_col(
    aes(x = as.numeric(cond)+.5),
    width = .3
  ) +
  ggdist::stat_halfeye(
    adjust = .5,
    width = .5,
    .width = 0,
    point_colour = NA
  ) +
  coord_flip()

plot(fig1_densitycol)

For some reason, the density plots were “squished” because the claimpercent values were way too high. It was then we realised that geom_col() was summing all of the values rather than taking the means from each condition.

Before calculating the means, we had to rename the conditions so that they better aligned with the original code. We copied the code from their R Markdown file to reorder the conditions, where they made religious (cond = 4) a reference category.

For “Data by Condition”, this made 1 = “Control”, 2 = “White Noise”, 3 = “Secular”, and 4 = “Religious”, and for “Data by Site”, 0 = “USA”, 1 = “Japan”, and 2 = “Czech Republic”.

# data by condition
data1$cond[data1$cond==4] <- 0 
data1$cond[data1$cond==1] <- 4 
data1$cond[data1$cond==3] <- 1 
data1$cond[data1$cond==4] <- 3

# data by site
data1$site[data1$site==1] <- 0 
data1$site[data1$site==3] <- 1

Then, we calculated the means by multiplying each condition/site with the amount of values they had, then dividing claimpercent by the number of values. We had to use the mutate function from dplyr to create a new variable.

# data by condition

data1 <- data1 %>% 
  mutate(numberOf = (cond == 0) * 100 + (cond == 1 | cond == 2) * 103 + (cond == 3) * 102) %>% 
  mutate(claimpercent2 = claimpercent / numberOf)

# data by site

data1 <- data1 %>% 
  mutate(numberOf2 = (site == 0) * 123 + (site == 2) * 127 + (site == 1) * 156) %>% 
  mutate(claimpercent3 = claimpercent / numberOf2)

And our efforts paid off!

fig1_densitycol2 <- ggplot(data1, aes(x = cond)) +
  geom_col(
    aes(x = as.numeric(cond)+.5, y = claimpercent2),
    width = .3
  ) +
  ggdist::stat_halfeye(
    aes(y = claimpercent),
    adjust = .5,
    width = .5,
    .width = 0,
    point_colour = NA
  ) +
  coord_flip()

plot(fig1_densitycol2)

First of all, we changed both cond and site to a factor with as.numeric(). We had to change it before we plotted but after mutating all the variables we needed because the variables would not calculate a factor. We later use factor() instead, because we realised as.numeric() was the wrong function.

In the ggdist::stat_halfeye section, adjust changed the height so we could leave room for the column. .width and point_colour were a part of the raincloud plot code from the blog, but we removed them since they are not in Fig 1.

We used coord_flip() from ggplot2 to flip the graph.

All we had left to do was:

  • label the x and y axis, as well as the x and y axis values,
  • move the columns below the density plot,
  • add title and subtitle,
  • add the error bars,
  • add the colours,
  • do the same for the “Data by Site” figure.

To calculate the error bars, we “group_cond” and “group_site” so geom_error(). group_by() creates a tibble so we can summarise() to make a new data frame.

The mean was calculated with mean = mean(claimpercent), while n = n() assigns n the number of rows in the dataframe.

Like na.rm = TRUE, drop_na() from tidyverse “drops” rows that contains missing values.

Then we created two new variables, “lowerCI” and “upperCI”, by calculating se and then calculating “lowerCI” and “upperCI”.

We had to change “group_cond” and “group_site” into a factor so it would align with our plot later. Without this, the error bars would not align with the cond and site bars which will also be turned into factors.

group_cond <- data1 %>% 
  group_by(cond) %>% 
  summarise(mean = mean(claimpercent, na.rm = TRUE),
            sd = sd(claimpercent, na.rm = TRUE),
            n = n()) %>% 
  drop_na() %>% 
  mutate(se = sd / sqrt(n),
         lowerCI = mean - qt(1 - (0.05/2), n - 1) * se,
         upperCI = mean + qt(1 - (0.05/2), n - 1) * se) %>% 
  ungroup()

group_cond$cond <- as.factor(group_cond$cond)

# data by site
group_site <- data1 %>% 
  group_by(site) %>% 
  summarise(mean = mean(claimpercent, na.rm = TRUE),
            sd = sd(claimpercent, na.rm = TRUE),
            n = n()) %>% 
  drop_na() %>% 
  mutate(se = sd / sqrt(n),
         lowerCI = mean - qt(1 - (0.05/2), n - 1) * se,
         upperCI = mean + qt(1 - (0.05/2), n - 1) * se) %>% 
  ungroup()

group_site$site <- factor(group_site$site)

The final Figure 1 plot needed the following:

  1. Create error bars

We used geom_errorbar() from ggplot2 to plot the error bars, using the SE and assigning the lower and upper CI from the previous chunk. I adjusted the width to fit the columns.

  1. Rename the x-axis, y-axis, as well as the axis values.

For scale_x_discrete() and scale_y_continuous(), expand made sure the x and y axis so it included zero. Then, we used label() to name each value, and used break() so labels and breaks were the same length. name = NULL removed the x axis label, since the original figure did not have one, and name = relabeled the y axis.

  1. Add the aesthetics of the Figure 1 plot.

Both from ggplot 2, theme(legend.position = "none") removed the legend.

scale_fill_manual allowed us to add the exact colours used in the original plot, and the colours were obtained by using a colour picking site.

ggtitle() from ggplot2 was used to create the “A.” title in the corner of the graph, and facet_grid() allowed us to put the title “Data by _____” in its own box.

The code for fig_cond and fig_site are the same - just with different colours and labels!

# data by condition

data1$cond <- factor(data1$cond)

fig1_cond <- ggplot(data1, aes(x = cond)) +
  geom_col(
    aes(y = claimpercent2, fill = cond),
    width = .3
  ) +
  ggdist::stat_halfeye(
    aes(y = claimpercent, fill = cond),
    adjust = .5,
    width = .4,
    .width = 0,
    point_colour = NA,
    position = position_nudge(x = 0.17, y = 0)
  ) +
  geom_errorbar(data = group_cond,
                aes(ymin = lowerCI, ymax = upperCI),
                show.legend = NA,
                width = .3
                ) +
  coord_flip() +
  expand_limits(x = 0, y = 0) +
  scale_x_discrete(name = NULL,
                     expand = c(0,0),
                     breaks = c(0, 1, 2, 3),
                     labels = c("Control", "Noise", "Secular", "Religious")
  ) +
  scale_y_continuous(name = "Percent Claimed",
                     expand = c(0,0),
                     breaks = c(0, 20, 40, 60, 80, 100)
                     ) +
  ggtitle(
    label = "Data by condition"
  ) +
  theme_classic() +
  theme(legend.position = "none") +
  scale_fill_manual(values = c("#adadad","#336ca3","#77c2ec","#ecb234")) +
  ggtitle(label = "A.") +
  facet_grid(. ~ "Data by condition")

plot(fig1_cond)

Then, all we had to do was use grid.arrange() from gridExtra to combine the two plots!

grid.arrange(fig1_site, fig1_cond, ncol = 2)

We really appreciated the power of Google for Figure 1 - without that raincloud plot blog we would have never figured out how to plot this!

Figure 2

Like we did for Figure 1, we had to select the relevant variables for Figure 2. To match with Nichols’ code, we had to reverse code “ritual” by subtracting 7 (the number of points on the scale) from ritual and taking the absolute value.

Cond was renamed again, and changed to a factor.

dataA$cond[dataA$cond==4] <- 0 
dataA$cond[dataA$cond==1] <- 4 
dataA$cond[dataA$cond==3] <- 1 
dataA$cond[dataA$cond==4] <- 3

dataA$cond <- factor(dataA$cond, levels= c(0,1,2,3),
                labels = c("Religious", "Secular", "Noise","Control"))

dataB was created to take the data from all participants in dataA that are in the religious condition, which is needed to plot the SE.

dataB <- dataA %>% 
  filter(cond == "Religious")

We also had to calculate the CI limits by calculating the mean, sd and n for each grouped variable, then applying the CI limit formula.

groupA <- dataA %>% 
  group_by(cond, affil) %>% 
  summarise(mean = mean(claimpercent, na.rm = TRUE),
            sd = sd(claimpercent, na.rm = TRUE),
            n = n()) %>% 
  drop_na() %>% 
  mutate(se = sd / sqrt(n),
         lowerCI = mean - qt(1 - (0.05/2), n - 1) * se,
         upperCI = mean + qt(1 - (0.05/2), n - 1) * se) %>% 
  ungroup()
  
# labelling the levels of religious affiliation 
groupA$affil <- factor(groupA$affil, levels = c(0, 1), 
                       labels = c("Non-affiliated", "Affiliated"))

For the first plot, we used method = "lm" to draw a line of best fit. geom_smooth also comes from the ggplot2 package and is used to plot figA and figB, and coord_cartesian() allowed us to set the ylimits without losing any data points.

figA <- ggplot(dataA, aes(religiosity, claimpercent, color = cond)) +
  geom_smooth(method = "lm", se = FALSE, size = 1.5) +
  geom_smooth(data = dataB, size = 0, se = TRUE, method = "lm", fill = "#fbf1d8") + 
  theme_light() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") + 
  coord_cartesian(ylim = c(0, 50)) + 
  labs(x = "Religiosity", y = "Percentage claimed") + 
  facet_grid(. ~ "Condition*Religiosity") + 
  scale_color_manual(values = c("#dea520", "#5ab3e4", "#094689", "#a4a4a4"))

The second plot had the same code as the first plot. As a note, figA and figB needed two geom_smooth() to plot the SE line for only the religious condition.

figC was different from figA and figB as we had to plot the CI limits. geom_line() was used to create those connecting lines in figC, and used the function position_dodge to move each line to stop them from overlapping.

figC <- ggplot(groupA, aes(affil, mean, color = cond)) +
  geom_line(aes(group = cond), position = position_dodge(width = 0.5), size = 1.5) + 
  geom_errorbar(data = groupA, aes(ymin = lowerCI, ymax = upperCI), width = 0.2, position = position_dodge(width = 0.5), size = 1.3) +
  theme_light() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ylim(0, 50) +
  scale_x_discrete() +
  labs(x = "Religious affiliation", y = "Percentage claimed") +
  facet_grid(. ~ "Condition*Religious affiliation") +
  scale_color_manual(values = c("#dea520", "#5ab3e4", "#094689", "#a4a4a4"))

Then, we combined all three plots.

grid.arrange(figA, figB, figC, ncol = 3)

Caterpillar coding was especially helpful for Figure 2. Victor had trouble with geom_line() but Sam was able to point out a mistake he might have made, which helped Victor solve the problem!

Table 1

Again, we selected the relevant variables.

data3 <- data3 %>% 
  select(-CT_practice, -CT_payments, -CT, -CT_cheat, -Religion_Text, -affil_cong)

We now had to recreate the new variables (negativity, positivity, tempo and impact). We first consulted the codebook to find out what the variables were, what they measured, and how they were calculated. It was then we found that they were the averages of variables in the csv, though this was not clear at first.

data3 <- data3 %>% 
  mutate(
    negativity = (distressing + irritating + boring + sad)/4,
    positivity = (interesting + pleasant + exciting + relaxing + happy)/5,
    impact = (deep + powerful)/2,
    tempo = (fast + abs(slow-7))/2
  )

After creating the new variables, we needed to produce the descriptive statistics.

To do this, we tried using the group_by() and summarise() function, but this did not work because we needed new values we could put into a table.

data3 %>% 
  group_by(cond) %>%
  summarise(mean_claim = mean(claimpercent, na.rm = TRUE),
            mean_sac = mean(sacred, na.rm = TRUE),
            mean_neg = mean(negativity, na.rm = TRUE),
            mean_pos = mean(positivity, na.rm = TRUE),
            mean_imp = mean(impact, na.rm = TRUE),
            mean_temp = mean(tempo, na.rm = TRUE),
            SD_claim = sd(claimpercent, na.rm = TRUE),
            SD_sac = sd(sacred, na.rm = TRUE),
            SD_neg = sd(negativity, na.rm = TRUE),
            SD_pos = sd(positivity, na.rm = TRUE),
            SD_imp = sd(impact, na.rm = TRUE),
            SD_temp = sd(tempo, na.rm = TRUE)
  )

By consulting the codebook again, we ended up using a long piece of code to obtain the descriptives (means, SDs and CIs) for the six variables (percent claimed, scaredness, negativity, positivity, tempo, impact) for each condition.

The code is separated into six separate parts for each variable, and is further partitionaed into four conditions, which also contains code that computes the six descriptive statistics (mean, SD, number of responses, SE, upper and lower CI limits).

First, we will explain how we labelled each value. The format of each label are as follows:

  • descriptive.condition number.variable, where the:
    • control group is (1)
    • white noise group is (2)
    • secular music group is (3)
    • and religious music group is (4).

For example,‘sd3s’ contains the standard deviation (sd) of the sacredness variable (s) in the secular music group (3).

We first calculated the mean and SD for each variable and condition.

We then had to calculate the SE to calculate the CI limits. The original paper used the length() function, which takes the list of responses in the group, though n() performs the same function. !is.na also removed the NA values.

Then, once we calculated SE, we put it in the CI limits formula.

The code below only shows the calculations for the Control Group and percent claimed variable, but we used the same code for the rest of the variables.

# Percent claimed variable
# Control Group (1)
mean1c <- mean(data3$claimpercent[data3$cond=="1"], na.rm = TRUE)*100 # mean
sd1c <- sd(data3$claimpercent[data3$cond=="1"], na.rm = TRUE)*100 # SD
n1c <- length(data3$cond[data3$cond=="1" & !is.na(data3$cond)]) # number of responses
se1c <- sd1c/sqrt(n1c) # standard error
lCI1c <- mean1c - (1.96*se1c) # lower 95% CI
uCI1c <- mean1c + (1.96*se1c) # upper 95% CI

In hindsight, we should have created a function to make the code more concise.

Finally, we had to calculate the Cohen’s d variable, in which the formula is:

The following values are labelled in a similar fashion to the previous format - d stands for Cohen’s d, the number refers to the effect comparison between conditions this variable is measuring, then the final letter represents which of the six variables is being measured.

The following code will only show the first few lines to save space, but the rest of the Cohen’s d values are calculated the same!

# computing Cohen's d
# % claimed
#relig vs sec
d1c <- abs((mean4c-mean3c)/sqrt((sd4c^2+sd3c^2)/2))
#relig vs noise
d2c <- abs((mean4c-mean2c)/sqrt((sd4c^2+sd2c^2)/2))

Now to finally put all these values into a table.

Firstly, we placed the data into a tibble.

The characteristics the table included were listed, then the mutate_if() function rounded the values to 2 decimal places. Using is.numeric() specifies that the functions only apply to columns with numeric values.

Finally, we used labelled the columns and plotted the table.

# Religious table
table1 <- tibble(
  characteristics = c("% claimed", "Sacredness", "Negativity", "Positivity", "Tempo", "Impact"),
  M = c(mean4c, mean4s, mean4n, mean4p, mean4t, mean4i),
  SD = c(sd4c, sd4s, sd4n, sd4p, sd4t, sd4i),
  lCI = c(lCI4c, lCI4s, lCI4n, lCI4p, lCI4t, lCI4i),
  uCI = c(uCI4c, uCI4s, uCI4n, uCI4p, uCI4t, uCI4i),
  d = c("-", "-", "-", "-", "-", "-")
) 
table1 %>% mutate_if(is.numeric, ~round(., 2)) %>%
  gt() %>%
  cols_label(characteristics = "Religious")
Religious M SD lCI uCI d
% claimed 27.33 29.92 21.52 33.14 -
Sacredness 5.17 1.39 4.90 5.44 -
Negativity 1.82 0.65 1.70 1.95 -
Positivity 2.54 0.85 2.37 2.70 -
Tempo 2.60 0.81 2.45 2.76 -
Impact 2.91 1.09 2.70 3.12 -

The same code was used for the secular table, white noise table, and control group table.

Secular M SD lCI uCI d
% claimed 30.32 30.98 24.33 36.30 0.10
Sacredness 4.32 1.35 4.06 4.58 0.62
Negativity 1.64 0.51 1.55 1.74 0.31
Positivity 2.72 0.78 2.56 2.87 0.22
Tempo 2.78 0.82 2.62 2.94 0.22
Impact 2.73 1.13 2.51 2.95 0.16
White Noise M SD lCI uCI d
% claimed 29.40 31.09 23.39 35.40 0.07
Sacredness 3.26 1.19 3.04 3.49 1.48
Negativity 2.58 0.96 2.39 2.76 0.93
Positivity 1.38 0.62 1.26 1.50 1.56
Tempo 3.99 0.72 3.85 4.13 1.80
Impact 1.76 0.93 1.58 1.94 1.13
Control Group M SD lCI uCI d
% claimed 22.38 23.89 17.69 27.06 0.18
Sacredness NA NA NA NA NA
Negativity NA NA NA NA NA
Positivity NA NA NA NA NA
Tempo NA NA NA NA NA
Impact NA NA NA NA NA

We used data.frame to combine the tables.

data_tables <- data.frame(Religious = table1, 
                          Secular = table2,
                          WhiteNoise = table3,
                          Control = table4)

gt(data_tables) let R understand what data frame we wanted to use. tab_source_note(source_note = ...), placed the text at the bottom of the table instead. The function tab_spanner(label = ...,columns = c()) allowed us to place a label that spans across a certain amount of columns, thus labelling each column by each group.

gt_tbl <- gt(data_tables) %>% 
  tab_source_note(
    source_note = "M = Mean; SD = Standard Deviation; CI = 95% Confidence Intervals. Cohen's d represents the effect size of comparisons between musical conditions.") %>%
  tab_spanner(
    label = "Religious (n = 102)",
    columns = c(Religious.M, Religious.SD, Religious.lCI, Religious.uCI, Religious.d)
  ) %>%
  tab_spanner(
    label = "Secular (n = 103)",
    columns = c(Secular.M, Secular.SD, Secular.lCI, Secular.uCI, Secular.d)
  ) %>%
  tab_spanner(
    label = "White Noise (n = 103)",
    columns = c(WhiteNoise.M, WhiteNoise.SD, WhiteNoise.lCI, WhiteNoise.uCI, WhiteNoise.d)
  ) %>%
  tab_spanner(
    label = "Control (n = 100)",
    columns = c(Control.M, Control.SD, Control.lCI, Control.uCI, Control.d)
  )

Then, we renamed the columns with the cols_label() function.

We used fmt_number on columns with numeric values to round the values to 2 decimal places. This was difficult for values that were marked with a “-”, as they were not calculated by the researchers. By coming back to the original table visualisation and changing the “-” into NA, the function finally worked, so we applied it to each column.

Next, we needed to hide the columns that repeated the %%claimed, sacredness, negativity, positivity, tempo and impact characteristics by using the function cols_hide(columns = c(...)).

Finally, we used the missing_text function replaced the NA’s with “-”!

gt_tbl %>% cols_label(Religious.characteristics = NULL, 
             Religious.M = "M",
             Religious.SD = "SD",
             Religious.lCI = "lCI",
             Religious.uCI = "uCI",
             Religious.d = "d",
             Secular.characteristics = "Secular",
             Secular.M = "M",
             Secular.SD = "SD",
             Secular.lCI = "lCI",
             Secular.uCI = "uCI",
             Secular.d = "d",
             WhiteNoise.characteristics = "White Noise",
             WhiteNoise.M = "M",
             WhiteNoise.SD = "SD",
             WhiteNoise.lCI = "lCI",
             WhiteNoise.uCI = "uCI",
             WhiteNoise.d = "d",
             Control.characteristics = "Control",
             Control.M = "M",
             Control.SD = "SD",
             Control.lCI = "lCI",
             Control.uCI = "uCI",
             Control.d = "d")%>% 
  fmt_number(
    columns = c(Religious.M, Religious.SD, Religious.lCI, Religious.uCI, Secular.M, Secular.SD, Secular.lCI, Secular.uCI, Secular.d, WhiteNoise.M, WhiteNoise.SD, WhiteNoise.lCI, WhiteNoise.uCI, WhiteNoise.d, Control.M, Control.SD, Control.lCI, Control.uCI, Control.d),
    decimals = 2,
    use_seps = FALSE) %>% 
  cols_hide(
    columns = c(Secular.characteristics, WhiteNoise.characteristics, Control.characteristics)) %>% 
  fmt_missing(
  columns = everything())
M Religious (n = 102) Secular (n = 103) White Noise (n = 103) Control (n = 100)
M SD lCI uCI d M SD lCI uCI d M SD lCI uCI d M SD lCI uCI d
% claimed 27.33 29.92 21.52 33.14 - 30.32 30.98 24.33 36.30 0.10 29.40 31.09 23.39 35.40 0.07 22.38 23.89 17.69 27.06 0.18
Sacredness 5.17 1.39 4.90 5.44 - 4.32 1.35 4.06 4.58 0.62 3.26 1.19 3.04 3.49 1.48 — — — — —
Negativity 1.82 0.65 1.70 1.95 - 1.64 0.51 1.55 1.74 0.31 2.58 0.96 2.39 2.76 0.93 — — — — —
Positivity 2.54 0.85 2.37 2.70 - 2.72 0.78 2.56 2.87 0.22 1.38 0.62 1.26 1.50 1.56 — — — — —
Tempo 2.60 0.81 2.45 2.76 - 2.78 0.82 2.62 2.94 0.22 3.99 0.72 3.85 4.13 1.80 — — — — —
Impact 2.91 1.09 2.70 3.12 - 2.73 1.13 2.51 2.95 0.16 1.76 0.93 1.58 1.94 1.13 — — — — —
M = Mean; SD = Standard Deviation; CI = 95% Confidence Intervals. Cohen's d represents the effect size of comparisons between musical conditions. 
  missing_text = "-"

The lack of a codebook/comprehensive R Markdown file made Table 1 especially hard! This made us really appreciate the power of rubber duck coding.

Exploratory Analyses

Although Nichols et al. raised some interesting questions in their Discussion, I was unable to assess those questions with the current dataset. This led to me coming up with some fun questions that I thought would be relevant to the real world instead!

For our statistical analyses, I used the following packages:

library(ggpubr)
library(report)

Is there a correlation between age and dishonest behaviour?

This research question was based on my curiosity on whether younger or older people were more likely to be dishonest. I found this to be an intriguing question because people tend to think teenagers and university students are more likely to misbehave. Empirical research has also looked into cheating behaviours in university students (Bretag et al., 2019).

The dataset contains ages from 15 to 68 so we can see if young adults are more likely to cheat compared to older people.

exp1 <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")

We renamed and selected the variables for this exploratory analysis.

exp1 <- exp1 %>%  
  filter(include == 0) %>% 
  rename(cond = con,
         claimpercent = claim,
         claimmoney = moneyclaim,
         CT_practice = `completion time (practice included)`,
         CT_payments = `completion time (payments only)`,
         religiosity = relig,
         religion = Religion) %>% 
  select(age, claimpercent) %>% 
  as_tibble()

To visualise the correlation, I created a scatterplot using geom_point() and used geom_smooth() to produce a line of best fit. A scatterplot proved useful in showing me all the data points and where they lie.

fig_correlational <- ggplot(data = exp1, aes(x = age, y = claimpercent), na.rm=TRUE) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_x_continuous(name = "Age") +
  scale_y_continuous(name = "Proportion of Dishonest Behaviour")

plot(fig_correlational)

In this graph, data points on the higher end of the y-axis indicate more dishonest behaviour. The line of best fit appears to show that dishonest behaviour is increasing with age, but it is important to note that this may be because there are fewer data points above the age of 40, hence high outliers will drag up the average.

As such, we should conduct a statistical analysis to see whether or not this is a significant effect.

cor.test(exp1$age,exp1$claimpercent)
## 
##  Pearson's product-moment correlation
## 
## data:  exp1$age and exp1$claimpercent
## t = 2.0056, df = 401, p-value = 0.04557
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.001992444 0.195441749
## sample estimates:
##        cor 
## 0.09965873

Pearson’s correlation indicates that are significant (p < .05), and that there is a small positive correlation between age and dishonest behaviour (r = .01, p < .05). Hence, the Pearson’s correlation implies that older people are slightly more likely to engage in dishonest behaviour than younger people.

However, this effect may still be because of the smaller sample sizes for older age groups, so we conducted a linear regression to determine if this effect is still significant.

lm_claimage <- lm(age ~ claimpercent, data = exp1)

summary(lm_claimage)
## 
## Call:
## lm(formula = age ~ claimpercent, data = exp1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.777 -3.437 -1.955  1.419 45.383 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    22.385      0.428  52.303   <2e-16 ***
## claimpercent    2.138      1.066   2.006   0.0456 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.263 on 401 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.009932,   Adjusted R-squared:  0.007463 
## F-statistic: 4.023 on 1 and 401 DF,  p-value: 0.04557

The linear regression model supports the Pearson’s correlation with a significant main effect (p < .05) of age and dishonest behaviour, with a significant regression equation (F(1, 401) = 4.023) with an R^2 of .001.

To conclude, dishonest behaviour is slightly correlated with age, where older ages are more likely to present dishonest behaviour than younger ages.

Are ages between affiliated and unaffiliated groups significantly different?

This question assesses the difference between younger or older age groups in regards to religious affiliation. From statistics of a USA population in 2017, it appears that younger age groups (18 to 29 years of age) are much more likely to be unaffiliated than older age groups (30 to 49, 60 to 64, 65+) (Statista Research Department, 2017). Therefore, I would expect that younger age groups will be more likely to be unaffiliated than older age groups.

exp2 <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")
exp2 <- exp2 %>%  
  filter(include == 0) %>% 
  rename(cond = con,
         claimpercent = claim,
         claimmoney = moneyclaim,
         CT_practice = `completion time (practice included)`,
         CT_payments = `completion time (payments only)`,
         religiosity = relig,
         religion = Religion) %>% 
  select(affil, id, age) %>% 
  drop_na() %>% 
  as_tibble()

exp2$affil <- factor(exp2$affil, labels = c("Unaffiliated","Affiliated"))

I also created a descriptive statistics table to visualise which ages are more likely to be affiliated with a religious organisation.

affil_summary <- exp2 %>% 
  group_by(affil) %>% 
  summarise(mean = mean(age),
            sd = sd(age),
            n = n(),
            se = sd/sqrt(n))

affil_summary %>% 
  gt() %>% 
  cols_label(affil = "Affiliation",
             mean = "Mean",
             sd = "SD",
             se = "SE") %>% 
  fmt_number(
    columns = c(mean, sd, se),
    decimals = 2)
Affiliation Mean SD n SE
Unaffiliated 22.46 5.37 304 0.31
Affiliated 24.63 8.67 91 0.91

The table shows that the ages for “Affiliated” are slightly higher (M = 24.63) than "Unaffiliated (M = 22.46). It is also interesting to note that the majority of the sample are unaffliated (n = 304).

To better visualise the data, I created a violin plot. My main reason for using a violin plot is that I saw my classmates using it, so I wanted to see if I preferred them to boxplots. Since violin plots also plot all the data points, it’s better to use for smaller sample sizes (according to this website). This is perfect since the “Affliated” group is much smaller than “Unaffiliated” group.

fig_affilage <- ggplot(data = exp2) +
  geom_violin(aes(x = affil, y = age, fill = affil)) +
  stat_summary(aes(x = affil, y = age), fun = "mean") +
  scale_y_continuous(name = "Age") +
  scale_x_discrete(name = "Affiiation") +
  theme_classic() +
  theme(legend.position = "none") +
  ggtitle(label = "Affiliation Across Ages",
          subtitle = "Nichols et al., 2020")

plot(fig_affilage)

This violin plot shows that there are a lot more younger people (y = 20) in the “Unaffiliated” group than the “Affiliated” group. Additionally, it appears that those over the age of 70 are unaffiliated, rather than affiliated. Interestingly, the affiliated group has more people aged 30-40.

To determine if this effect is significant, I conducted a t test to compare the mean ages of “Affiliated” and “Unaffiliated” groups.

affiliated <- exp2 %>% 
  filter(affil == "Affiliated")

unaffiliated <- exp2 %>% 
  filter(affil == "Unaffiliated")

t.test(affiliated$age, unaffiliated$age)
## 
##  Welch Two Sample t-test
## 
## data:  affiliated$age and unaffiliated$age
## t = 2.26, df = 111.42, p-value = 0.02576
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2673378 4.0709358
## sample estimates:
## mean of x mean of y 
##  24.62637  22.45724

The t test reveals that there is a significant effect between Affiliated and Unaffiliated groups (p < .05). The 95% confidence interval [0.27, 4.07] shows us that the average ages of the unaffiliated group are at most 4.07 years older, and at least 0.27 years older.

Ultimately, those in the affiliated group have a slightly higher average age than those in the unaffiliated group.

Are there significant differences of dishonest behaviour between different religions? (Christian, Buddhism, Judaism, Shinto)

Although the paper assessed religiosity in general, they did not assess whether differing religions would have differing levels of dishonest behaviour. I thought that would be interesting to look at because religions may have different moral code or religious practices that may have a unique influence on dishonest behaviour.

exp3 <- read_csv("~/RFiles/verification_report/Nichols_et_al_dataset_V2.0.csv")

The code below is the same as the previous analyses, except that religion is made a factor and religion is filtered to exclude “Other”, “Atheism” and “Agnosticism”.

exp3 <- exp3 %>%  
  filter(include == 0) %>% 
  rename(cond = con,
         claimpercent = claim,
         claimmoney = moneyclaim,
         CT_practice = `completion time (practice included)`,
         CT_payments = `completion time (payments only)`,
         religiosity = relig,
         religion = Religion) %>%
  select(claimpercent, id, religion) %>% 
  drop_na()  %>% 
  as_tibble()

exp3 <- filter(exp3, religion == "Christian" | religion == "Buddhism" | religion == "Judaism" | religion == "Shinto")

exp3$religion <- factor(exp3$religion)

I created a table to easily visualise the differences in means and n.

religion_summary <- exp3 %>% 
  group_by(religion) %>% 
  summarise(mean = mean(claimpercent),
            sd = sd(claimpercent),
            n = n(),
            se = sd/sqrt(n))

religion_summary %>% 
  gt() %>% 
  cols_label(religion = "Religion",
             mean = "Mean",
             sd = "SD",
             se = "SE") %>% 
  fmt_number(
    columns = c(mean, sd, se),
    decimals = 2)
Religion Mean SD n SE
Buddhism 0.26 0.30 37 0.05
Christian 0.28 0.29 109 0.03
Judaism 0.14 0.11 4 0.06
Shinto 0.19 0.12 13 0.03

The table shows us that the Christian group had the highest proportion of dishonest behaviour (M = 0.28), with Buddhism in second (M = 0.26). Shinto (M = 0.19) and Judaism (M = 0.14) had the lowest proportions of dishonest behaviour, but this maybe because they had extremely small sample sizes. We can better visualise this with a violin plot.

For this plot, I also decided to add a boxplot to see if it would provide more information than just a violin plot.

fig_religion <- ggplot(data = exp3, aes(x = religion, y = claimpercent)) +
  geom_violin(aes(fill = religion),
              show.legend = FALSE) +
  geom_boxplot(width= .2) +
  stat_summary(aes(x = religion, y = claimpercent),
               fun = "mean") +
  scale_y_continuous(name = "Percent Claimed") +
  scale_x_discrete(name = "Types of Religion") +
  ggtitle(label = "Percent Claimed Differences Between  Religions",
    subtitle = "Nichols et al., 2020")

plot(fig_religion)

This plot shows similar trends to the table, where Christian and Buddhism had the higher proportions of dishonest behaviour. We can see that Judaism does have the smallest proportion of dishonest behaviour.

An ANOVA test was conducted to see if this effect is significant.

aov_claimrelig <- aov(claimpercent ~ religion, data = exp3)

summary(aov_claimrelig)
##              Df Sum Sq Mean Sq F value Pr(>F)
## religion      3  0.147 0.04894   0.634  0.594
## Residuals   159 12.265 0.07714
report(aov_claimrelig)
## The ANOVA (formula: claimpercent ~ religion) suggests that:
## 
##   - The main effect of religion is statistically not significant and small (F(3, 159) = 0.63, p = 0.594; Eta2 = 0.01, 90% CI [0.00, 0.04])
## 
## Effect sizes were labelled following Field's (2013) recommendations.

It appears that the difference of average dishonest behaviour between Christian, Buddhism, Judaism and Shinto groups is not of practical importance (p > .05), and should be researched further.

These results indicate that there may or may not be a significant effect of dishonest behaviour between differing religions.

Recommendations

A Comprehensive R Markdown File

My main recommendation is to have a more comprehensive R Markdown file. The R Markdown file provided did a poor job of explaining their own code.

As aforementioned, the plots used in the R Markdown file were independent of the ones in the paper, so we could not use it to help us reproduce the plots. This proved especially difficult when I was trying to recreate Figure 1 and had no idea which functions to use.

A more comprehensive R Markdown file would have also provided explanations for how and why code was renamed/organised. One prime example was when I had to reorder “cond” and “site”. Another example includes the reverse coding of the slow variable. We were unsure as to why they reverse coded, which left us guessing but the point of an open and reproducible code means it should be clear why each line of code was done.

By creating a comprehensive R Markdown file, this means that the plots and statistics would be easier to reproduce, without needing to contact the original authors for clarification.

Rubber duck coding

Rubber duck coding would be highly recommended to make their code more understandable. Specifically, rubber duck coding allows for a beginner to R to understand what functions and packages you used, especially if the functions or packages inevitably become outdated. This blog also highlights that it helps people, even if they are proficient in R, understand your problem-solving process.

Although this blog is centered around solving problems in a workplace environment, it has some useful suggestions that would help with reproducibility:

  • Let go of assumptions: assume your audience knows little to nothing about your situation
  • Explain everything: give the full context and background of your code

This recommendation is especially relevant to Nichols et al. because even though they had a R Markdown file, I unable to knit the file because of outdated packages, so there was no way of knowing what the code referred to.

Include a codebook

One of our main struggles when starting the verification part of the report was trying to figure out what the variables meant! This meant we had to go deep into the R Markdown file and reread the paper to make our own codebook. This could have been made much simpler if there was already a codebook provided.

If needed, this guide is a good and simple introduction to how to write a codebook.

Overall, Nichols et al. did take the right steps by providing a R Markdown file and providing their dataset, but more steps need to be taken in order to make the process of computational reproducibility easier.

References

Bretag, T., Harper, R., Rundle, K., Newton, P. M., Ellis, C., Saddiqui, S., & van Haeringen, K. (2019). Contract cheating in Australian higher education: a comparison of non-university higher education providers and universities. Assessment & Evaluation in Higher Education. doi: 10.1080/02602938.2019.1614146

Statista Research Department (2017). Religious affiliation of the population in the United States in 2017, by age. Retrieved from https://www.statista.com/statistics/245453/religious-affiliation-in-the-united-states-by-age/.