Summary

The behavioural immune system helps humans to avoid infection, before it enters the body. Disgust is a key emotion in this system and can be used interpersonally. The literature suggests people use trade-offs to balance the need fo social interaction (and thus direct and indirect contact with people) and the risk of infection (avoiding pathogens). However, this system can be subjective – people tend to view infection threats as being lower in those they’re more familiar with. This study investigated the link between willingness to engage in infection-risky behaviours and the level of value put on a person. (Tybur et al. 2020) conducted three studies. In studies one and two, participants reported how comfortable they were with infection-risky behaviour with someone they knew personally (a romantic partner, a friend, an acquaintance or someone they disliked). They then completed a task in which researchers determined their welfare-trade off ratio (WTR). This ratio reflected people’s willingness to trade off their own welfare for that of another. In study three, participants reported their comfort with these same infection-risky behaviours, with a stranger, who was described as either high or low on honesty-humility and agreeableness which supposedly inform how much we value a person (interpersonal value). It was found that people are less opposed to engaging in infection-risky behaviour with friends compared to those they dislike, and with honest and agreeable strangers relative to dishonest and disagreeable ones. The level at which a person values someone is commensurate with them being comfortable engaging in infection-risky acts with that person (more value, less aversion to risky contact). Motivations to avoid infection risky contact varied across people with no symptoms, mostly accounted for their interpersonal value. Thus, infection transmission within social networks may be worsened by being less worried about contamination with social partners you value highly. Prejudices towards obese, elderly, physically disabled or foreigners may result from perceptions of interpersonal value rather than perceptions of infectiousness. Contact rituals such as handshaking or hugging maybe be used to indicate and support interpersonal valuation and refusals of such rituals may be interpreted as suggesting the target is not valued enough to risk infection, is perceived as having some symptom of contagious illness, or both. Thus, interventions to reduce infection, could focus on reducing infection within social networks. (Tybur et al. 2020) findings might help inform approaches to reducing disease transmission during outbreaks (such as COVID-19) – people need help understanding the risk of engaging in risky behaviours with those they value.

Reaction

1. I was surprised that…humans are willing to put people they like more, at more risk of infection, than those they don’t like. Although this makes sense when you think about the ways in which COVID-19 spread, it seems counterintuitive to Evolutionary Psychology. If we were all in a wolf pack as cavemen, we would want to do anything to protect our pack, and ensure we were not ousted from the pack. Surely risking infection would be a risk to our caveman pack?

2. It seems that the next step in this area of research would be to…test whether this relaxed pathogen avoidance towards offspring and mates results only from their high interpersonal value or whether sexual value and genetic relatedness (which inform interpersonal value) additionally affect pathogen avoidance. Secondly, to clarify whether the behavioural immune system produces prejudice toward the obese, elderly and physically disabled because they’re perceived as not offering the interpersonal benefits that offset the infection risks posed by any social interaction.

3. The take home message that I am left with after reading this paper is…that people don’t need help not interacting with most people, especially in a pandemic. However, when it comes to people they like, a lot of education needs to happen occur related to our biases towards infection risk. Yet again, it appears that our cognitive biases are not as clever as they may seem, and are actually putting us at risk.

Load packages

library(tidyverse)
library(patchwork)
library(dplyr)
library(readr)
library(ggplot2)
library(infer)
library(janitor)
library(here)
library(rstatix)
library(ggeasy)
library(ggpubr)
library(corrplot)
library(car)
library(psych)
library(BayesFactor)
library(gt)
library(jmv)
library(tidyverse)
library(cowplot)
library(extrafont)
library(report)
library(papaja)
library(apaTables)

Read data

The data in the OSF repo is in .csv format, so I need to use the read_csv() function from the reader package.

data_1_raw = read_csv('WTR_Comfort_S1.csv')
data_2 = read_csv('WTR_Comfort_S2.csv')
data_3 = read_csv('WTR_Comfort_S3.csv')

Reproducing Demographic descriptives

nparticipants = nrow(data_1_raw)
print (nparticipants)

## [1] 504

nparticipants = nrow(data_2)
print (nparticipants)

## [1] 430

nparticipants = nrow(data_3)
print (nparticipants)

## [1] 905

data_1_raw$sex <- as.factor(data_1_raw$sex)
data_1_raw %>% count(sex=="1")

## # A tibble: 2 x 2
##   `sex == "1"`     n
##   <lgl>        <int>
## 1 FALSE          226
## 2 TRUE           278

data_2$sex <- as.factor(data_2$sex)
data_2 %>% count(sex == "1")

## # A tibble: 2 x 2
##   `sex == "1"`     n
##   <lgl>        <int>
## 1 FALSE          188
## 2 TRUE           242

data_3$sex <- as.factor(data_3$sex)
data_3 %>% count(sex=="1")

## # A tibble: 2 x 2
##   `sex == "1"`     n
##   <lgl>        <int>
## 1 FALSE          453
## 2 TRUE           452

mean_age = mean(data_1_raw$age)
print(round(mean_age, digits = 2))

## [1] 35.88

sd_age = sd(data_1_raw$age)
print(round(sd_age, digits = 2))

## [1] 10.2

mean_age = mean(data_2$age)
print(round(mean_age, digits = 2))

## [1] 36.29

sd_age = sd(data_2$age)
print(round(sd_age, digits = 2))

## [1] 10.99

mean_age = mean(data_3$age)
print(round(mean_age, digits = 2))

## [1] 38.22

sd_age = sd(data_3$age)
print(round(sd_age, digits = 2))

## [1] 11.64

Exclude_english_s1 <- filter(data_1_raw, English_exclude == "1")
print(nrow(Exclude_english_s1))

## [1] 3

data_2_english_filtered <- data_2 %>%
  filter(English_exclude != "1") 

data_3_english_filtered <- data_3 %>%
  filter(English_Exclude != "1") 

exclude_sex_s1 <- filter(data_1_raw, sex == "3")
print(nrow(exclude_sex_s1))

## [1] 2

data_2_sex_filtered <- data_2 %>% filter(sex != "3") 
print(data_2_sex_filtered)

## # A tibble: 429 x 85
##    participant sex     age relat income poli_soc poli_econ trust_gen   HH1   HH2
##          <dbl> <fct> <dbl> <dbl>  <dbl>    <dbl>     <dbl>     <dbl> <dbl> <dbl>
##  1           1 2        20     2     10        4         4         6     3     3
##  2           2 2        46     1      3        1         1         4     4     2
##  3           3 1        38     1      6        1         1         7     5     1
##  4           4 1        42     1      9        1         2         8     5     2
##  5           5 2        56     1     17        6         6         8     3     4
##  6           6 2        42     1      3        2         2         2     3     1
##  7           7 1        30     1     14        6         6        10     5     1
##  8           8 2        39     2     10        1         2         6     4     4
##  9           9 2        27     1     16        5         5         9     2     4
## 10          10 2        43     1     20        6         6         9     5     5
## # … with 419 more rows, and 75 more variables: HH3 <dbl>, HH4 <dbl>, HH5 <dbl>,
## #   HH6 <dbl>, HH7 <dbl>, HH8 <dbl>, HH9 <dbl>, HH10 <dbl>, AG1 <dbl>,
## #   AG2 <dbl>, AG3 <dbl>, AG4 <dbl>, AG5 <dbl>, AG6 <dbl>, AG7 <dbl>,
## #   AG8 <dbl>, AG9 <dbl>, AG10 <dbl>, DS1 <dbl>, DS2 <dbl>, DS3 <dbl>,
## #   DS4 <dbl>, DS5 <dbl>, DS6 <dbl>, DS7 <dbl>, target_sex <dbl>,
## #   target_category <dbl>, relate_length <dbl>, partner_age <dbl>,
## #   partner_proximity <dbl>, partner_height <dbl>, partner_weight <dbl>,
## #   partner_attractive <dbl>, partner_hygiene <dbl>, comf1 <dbl>, comf2 <dbl>,
## #   comf3 <dbl>, comf4 <dbl>, comf5 <dbl>, comf6 <dbl>, comf7 <dbl>,
## #   comf8 <dbl>, comf9 <dbl>, comf10 <dbl>, 75_109 <dbl>, 75_94 <dbl>,
## #   75_79 <dbl>, 75_64 <dbl>, 75_49 <dbl>, 75_34 <dbl>, 75_19 <dbl>,
## #   75_4 <dbl>, 75_-11 <dbl>, 75_-26 <dbl>, 19_28 <dbl>, 19_24 <dbl>,
## #   19_20 <dbl>, 19_16 <dbl>, 19_12 <dbl>, 19_9 <dbl>, 19_5 <dbl>, 19_1 <dbl>,
## #   19_-3 <dbl>, 19_-7 <dbl>, 46_67 <dbl>, 46_58 <dbl>, 46_48 <dbl>,
## #   46_39 <dbl>, 46_30 <dbl>, 46_21 <dbl>, 46_12 <dbl>, 46_2 <dbl>,
## #   46_-7 <dbl>, 46_-16 <dbl>, English_exclude <dbl>

data_3_sex_filtered <- data_3 %>% filter(sex != "3") 
print(data_3_sex_filtered)

## # A tibble: 897 x 68
##    participant sex     age relat income DS_Inst   DS1   DS2   DS3   DS4   DS5
##          <dbl> <fct> <dbl> <dbl>  <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1         599 2        28     1     10       1     7     6     7     5     7
##  2         542 2        37     1      7       1     5     7     5     6     5
##  3         685 1        37     1      6       1     7     7     7     7     7
##  4         778 2        77     1      4       1     7     7     5     7     7
##  5         484 2        73     2      4       1     4     6     4     4     4
##  6          23 1        70     1      7       1     3     3     2     5     4
##  7         643 1        70     1      5       1     4     5     3     6     5
##  8         276 2        69     2      2       1     7     5     4     7     4
##  9         360 1        69     1     10       1     7     6     4     6     5
## 10          33 2        68     2      4       1     6     7     4     4     7
## # … with 887 more rows, and 57 more variables: DS6 <dbl>, DS7 <dbl>,
## #   Target <dbl>, Value <dbl>, Alex_inst <dbl>, Alex_height <dbl>,
## #   Alex_weight <dbl>, Alex_age <dbl>, Alex_wealth <dbl>,
## #   Alex_intelligent <dbl>, Alex_attractive <dbl>, Alex_kind <dbl>,
## #   Alex_honest <dbl>, Alex_care <dbl>, comfort_inst <dbl>, comf1 <dbl>,
## #   comf2 <dbl>, comf3 <dbl>, comf4 <dbl>, comf5 <dbl>, comf6 <dbl>,
## #   comf7 <dbl>, comf8 <dbl>, comf9 <dbl>, comf10 <dbl>, WTR_INST <dbl>,
## #   75_109 <dbl>, 75_94 <dbl>, 75_79 <dbl>, 75_64 <dbl>, 75_49 <dbl>,
## #   75_34 <dbl>, 75_19 <dbl>, 75_4 <dbl>, 75_-11 <dbl>, 75_-26 <dbl>,
## #   19_28 <dbl>, 19_24 <dbl>, 19_20 <dbl>, 19_16 <dbl>, 19_12 <dbl>,
## #   19_9 <dbl>, 19_5 <dbl>, 19_1 <dbl>, 19_-3 <dbl>, 19_-7 <dbl>, 46_67 <dbl>,
## #   46_58 <dbl>, 46_48 <dbl>, 46_39 <dbl>, 46_30 <dbl>, 46_21 <dbl>,
## #   46_12 <dbl>, 46_2 <dbl>, 46_-7 <dbl>, 46_-16 <dbl>, English_Exclude <dbl>

# Rename variables for consistency
colnames(data_1_raw)[16] = 'relationship'
colnames(data_2)[37] = 'relationship'
colnames(data_3)[15] = 'relationship'
colnames(data_3)[68] = 'English_exclude'

# Relevel relationship values for study 3
data_3$relationship = factor(as.factor(data_3$relationship), levels = c(1, 0))

# Function to transform data for each study
transform_data = function(study, names_list) {
  
  # Choose appropriate data
  if (study == '1')       data = data_1_raw             
  else if (study == '2')  data = data_2
  else                    data = data_3
  
  # Filter out exclusions
  # data = filter(data, English_exclude == 0)   
  
  # Create average contact comfort
  data = data %>% mutate(avg_cc = rowMeans(select(., comf1:comf10) - 4)) 
  
  # Create dataframe for WTR
  if (study == '1')       WTR = data[,c(49:40,59:50,69:60,79:70,89:80,99:90)]           
  else if (study == '2')  WTR = data[,c(64:55,74:65,84:75)]
  else                    WTR = data[,c(47:38,57:48,67:58)]
  
  # Adjust WTR values to 0 or 1
  WTR[WTR == -0.45] = 0
  WTR[WTR != 0] = 1
  
  # Create a matrix to store WTR information
  if (study == '1')       caulsum = matrix(,,ncol = 7)            
  else                    caulsum = matrix(,,ncol = 4)
  caulsum = caulsum[-1,]
  colnames(caulsum) = names_list
  
  # For each participant calculate WTR
  for (i in 1:nrow(data)) {
    individual_WTR = matrix(WTR[i,], ncol = 10, byrow = TRUE)
    WTR_anchor = apply(individual_WTR, 1, calculate_WTR)
    if (NA%in%WTR_anchor) WTR_total = NA
    else WTR_total = mean(WTR_anchor)
    if (study == '1')   caulperson = matrix(c(WTR_anchor,WTR_total), ncol = 7)            
    else                caulperson = matrix(c(WTR_anchor,WTR_total), ncol = 4)
    caulsum = rbind(caulsum,caulperson)
  }
  
  # Create final dataframe
  caulsum = as.data.frame(caulsum)
  caulsum = round(caulsum, digits = 3)
  
  final_data = as.data.frame(filter(mutate(caulsum, relationship = data$relationship, avg_cc = data$avg_cc, data$sex != 3, english_exclude = data$English_exclude)))  
  return(final_data)
}

# Function to calculate WTR 
calculate_WTR = function(raw_WTR) {
  raw_WTR = as.numeric(raw_WTR)
  m = seq(-35, 145, by = 20)
  if (sum(raw_WTR) == 10) anchor_WTR = 1.55
  else if (sum(raw_WTR) == 0) anchor_WTR = -0.45
  else {
    shift_count = 0
    shift_point = 0
    for (i in 1:9) {
      if ((raw_WTR[i] - raw_WTR[i + 1]) == 1){
        shift_count = shift_count + 1
        shift_point = (m[i] + m[i + 1]) / 2 + shift_point
      }
    }
    if (shift_count > 2) anchor_WTR = NA
    else anchor_WTR = 0.01 * shift_point / shift_count
  }
}

# Transform data for each study
final_data_1_raw = transform_data(1, c("WTR_37","WTR_23","WTR_75","WTR_19","WTR_46","WTR_68","WTR_total"))
final_data_2 = transform_data(2, c("WTR_75","WTR_19","WTR_46","WTR_total"))
final_data_3 = transform_data(3, c("WTR_75","WTR_19","WTR_46","WTR_total"))

final_data_2 = filter(final_data_2, english_exclude == 0)
final_data_3 = filter(final_data_3, english_exclude == 0)

#Counts number of NAs in data
excl_1 = sum(is.na(final_data_1_raw$WTR_total))
excl_2 = sum(is.na(final_data_2$WTR_total))
excl_3 = sum(is.na(final_data_3$WTR_total))

print(excl_1)

## [1] 35

print(excl_2)

## [1] 19

print(excl_3)

## [1] 30

Reproducing Figure 1 (Violin Plots)

Figure 1 is three violin plots with box plots patched together. Each study has relationship type on the x-axis and contact comfort on the y-axis. The figure includes: median (box plot), interquartile range (box plot), density (shaded violin), whiskers (box plot - extend 1.5 * IQ range), and outliers (dots).

First I’ll use colnames() [] = to rename the variables for consistency.

colnames(data_1_raw)[16] = 'relationship'
colnames(data_2)[37] = 'relationship'
colnames(data_3)[15] = 'relationship'
colnames(data_3)[68] = 'English_exclude'

I’m using levels = c() to relevel relationship values for study 3. They’re the opposite way around.

data_3$relationship = factor(as.factor(data_3$relationship), levels = c(1, 0))

I’m now using a function to transform the data and generate a plot for each study.

First I need to choose appropriate data using if, else if and else. Then I need to filter out the exclusions using filter(). Then, I need to create an average for contact comfort using avg_cc = rowMeans(select(., comf1:comf10)-4 and move the values from 1 to 7 to -3 to 3. Then to generate the plot, first we’ll map relationship on x-axis and contact comfort on y-axis. Next, we’ll add violin plot (geom_violin) and boxplot (geom_boxplot). We then need to rename x-axis labels using labels= and specify our colours using scale_fill_manual(values = colour_list) and scale_colour_manual(values = colour_list). Next we’ll adjust the theme using theme(), add our title using ggtitle and finally, adjust x-axis of plot depending on study.

generate_plot = function(study, relationship_labels, colour_list) {
  # Choose appropriate data
  if (study == '1') data = data_1_raw             
  else if (study == '2') data = data_2
  else data = data_3
  # Filter out exclusions
  data = filter(data, English_exclude == 0)   
  # Create average contact comfort
  data = data %>% mutate(avg_cc = rowMeans(select(., comf1:comf10)-4))    
# Generate a plot
  # Relationship on x-axis, contact comfort on y-axis
  plot = ggplot(data, aes(as.factor(relationship), avg_cc)) +     
    # Add violin plot and boxplot
    geom_violin(aes(fill = as.factor(relationship), colour = as.factor(relationship), alpha = .7)) +
    geom_boxplot(aes(colour = as.factor(relationship)), width = 0.1, notch = FALSE, 
                 outlier.color = NULL, outlier.shape = 16, outlier.size = 1) +
    # Rename x-axis labels
    scale_x_discrete(name = NULL, labels = relationship_labels) +
    # Specify colours
    scale_fill_manual(values = colour_list) +
    scale_colour_manual(values = colour_list) +
    # Adjust theme
    theme(
      legend.position = 'none',
      panel.background = element_rect(fill = NA),
      plot.background = element_rect(fill = NA, color = NA),
      panel.border = element_rect(fill = "transparent", colour = NA),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      plot.title = element_text(hjust = 0.5),
      axis.line.x.bottom = element_line(color = 'black'),
      plot.margin = margin(0, 0, 0, 0, 'cm')
    ) +
       # Add title
 ggtitle(paste('Study', study))
  # Adjust x-axis of plot depending on study
  if (study == '1') {
    plot = plot + theme(axis.line.y.left = element_line(color = 'black')) +
      scale_y_continuous(name = 'Contact Comfort', n.breaks = 7)
  } else {
    plot = plot + theme(
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      axis.title.y = element_blank()
    )
  }
  return(plot)
}

Now we can generate a plot for each study using generate_plot

plot_1 = generate_plot(1, c('Romantic\nPartner', 'Friend', 'Acquaintance', 'Enemy'), c('#3f6485', '#b17276', '#829966', '#e89b39'))
plot_2 = generate_plot(2, c('Friend', 'Acquaintance', 'Enemy'), c('#b17276', '#829966', '#e89b39'))
plot_3 = generate_plot(3, c('High-Value\nStranger', 'Low-Value\nStranger'), c('#86699a', '#6494a3'))

and combine plots together using master_plot = plot_1 + plot_2 + plot_3

master_plot = plot_1 + plot_2 + plot_3 + plot_layout(widths = c(3, 2.3, 1.9)) +
  theme(plot.caption = element_text(hjust = 1, family = 'Times New Roman')) +
  labs(caption = "Fig. 1. Mean comfort with infection-risky contact with each target type, separately for Studies 1 through 3. In each violin
    plot, the horizontal line indicates the median, the white box indicates the interquartile range of the data, the shaded area
    indicates the density of the data, and the whiskers extend 1.5 times the interquartile range. Outliers are indicated by dots.")

And now to add the dotted lines between each plot using ggdraw() + draw_line

ggdraw() +
  draw_plot(master_plot) +
  draw_line(
    x = c(0.445, 0.445),
    y = c(0.15, 0.9),
    colour = "dark grey", size = 0.5, linetype = 2
  ) +
  draw_line(
    x = c(0.741, 0.741),
    y = c(0.15, 0.9),
    colour = "dark grey", size = 0.5, linetype = 2
  )

Now this is our final code…I unfortunately don’t have all of our ‘failed’ attempts but I can explain some of the struggles we had.

1. How do we make the outliers the same colour as the violin plot?
2. How do we get rid of the y-axis for all of the plots so they all shre contact comfort?
3. How do we add dotted lines between the plots?
4. How to have the circle in each box of the legend that matches the colour of the box and also for the dashed lines to match the background of the box?
5. How do we change the font?
6. How do we add the x- and y-axis lines back?
7. How do we adjust the width of the individual plots when patching them together?
8. How were WTR switchpoints defined?

Reproducing Figure 2 (Scatter Plots with regression lines)

This figure is three scatterplots with regression lines patched together. Each study has WTR on the x-axis and contact comfort on the y-axis. Points and regression lines are grouped by relationship type. It includes points, regression line between WTR and contact comfort across categories (solid line), regression lines between WTR and contact comfort within categories (dashed line), and 95% confidence intervals (shaded areas around each regression line).

# Function to transform data for each study
transform_data = function(study, names_list) {
  # Choose appropriate data
  if (study == '1')       data = data_1_raw             
  else if (study == '2')  data = data_2
  else                    data = data_3
  # Filter out exclusions
  data = filter(data, English_exclude == 0)   
  # Create average contact comfort
  data = data %>% mutate(avg_cc = rowMeans(select(., comf1:comf10) - 4)) 
  # Create dataframe for WTR
  if (study == '1')       WTR = data[,c(49:40,59:50,69:60,79:70,89:80,99:90)]           
  else if (study == '2')  WTR = data[,c(64:55,74:65,84:75)]
  else                    WTR = data[,c(47:38,57:48,67:58)]
  # Adjust WTR values to 0 or 1
  WTR[WTR == -0.45] = 0
  WTR[WTR != 0] = 1
  # Create a matrix to store WTR information
  if (study == '1')       caulsum = matrix(,,ncol = 7)            
  else                    caulsum = matrix(,,ncol = 4)
  caulsum = caulsum[-1,]
  colnames(caulsum) = names_list
  # For each participant calculate WTR
  for (i in 1:nrow(data)) {
    individual_WTR = matrix(WTR[i,], ncol = 10, byrow = TRUE)
    WTR_anchor = apply(individual_WTR, 1, calculate_WTR)
    if (NA%in%WTR_anchor) WTR_total = NA
    else WTR_total = mean(WTR_anchor)
    if (study == '1')   caulperson = matrix(c(WTR_anchor,WTR_total), ncol = 7)            
    else                caulperson = matrix(c(WTR_anchor,WTR_total), ncol = 4)
    caulsum = rbind(caulsum,caulperson)
  }
  # Create final dataframe
  caulsum = as.data.frame(caulsum)
  caulsum = round(caulsum, digits = 3)
  final_data = as.data.frame(filter(mutate(caulsum, relationship = data$relationship, avg_cc = data$avg_cc, data$sex != 3)))
  return(final_data)
}
# Function to calculate WTR 
calculate_WTR = function(raw_WTR) {
  raw_WTR = as.numeric(raw_WTR)
  m = seq(-35, 145, by = 20)
  if (sum(raw_WTR) == 10) anchor_WTR = 1.55
  else if (sum(raw_WTR) == 0) anchor_WTR = -0.45
  else {
    shift_count = 0
    shift_point = 0
    for (i in 1:9) {
      if ((raw_WTR[i] - raw_WTR[i + 1]) == 1){
        shift_count = shift_count + 1
        shift_point = (m[i] + m[i + 1]) / 2 + shift_point
      }
    }
    if (shift_count > 2) anchor_WTR = NA
    else anchor_WTR = 0.01 * shift_point / shift_count
  }
}
# Function to generate a plot
generate_plot = function(study, relationship_labels, colour_list, fill_list) {
  # Choose appropriate data
  if (study == '1')       data = final_data_1_raw             
  else if (study == '2')  data = final_data_2
  else                    data = final_data_3
  # Generate a plot
  # WTR on x-axis, contact comfort on y-axis
  plot = ggplot(data, aes(WTR_total, avg_cc))
  # Adjust size of points for scatterplot depending on study
  if (study == '3') {
    plot = plot + geom_point(aes(group = relationship, colour = as.factor(relationship), fill = as.factor(relationship)), alpha = 0.5, shape = 16, size = 0.7, show.legend = TRUE)
  } else {
    plot = plot + geom_point(aes(group = relationship, colour = as.factor(relationship), fill = as.factor(relationship)), alpha = 0.5, shape = 16, size = 1.2, show.legend = TRUE)
  }
  plot = plot +
    # Add regression lines for individual relationship groups and overall
    geom_smooth(aes(group = relationship, colour = as.factor(relationship), fill = as.factor(relationship)), size = 0.6, alpha = 0.2, method = 'lm', linetype = 'dashed') +
    geom_smooth(aes(colour = 'black'), size = 0.6, alpha = 0.7, method = 'lm') +
    # Rename and relabel x- and y-axis
    scale_y_continuous(
      name = 'Contact Comfort',
      limits = c(-3, 3),
      n.breaks = 7
    ) +
    scale_x_continuous(
      name = 'WTR',
      breaks = c(-0.45, 0.05, 0.55, 1.05, 1.55)
    ) +
    # Specify colours
    scale_colour_manual(values = colour_list) +
    scale_fill_manual(labels = relationship_labels, values = fill_list) +
    # Adjust theme
    theme(
      panel.background = element_rect(fill = NA),
      plot.background = element_rect(fill = NA, color = NA),
      panel.border = element_rect(fill = "transparent", colour = NA),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      plot.title = element_text(hjust = 0.5),
      axis.line.x.bottom = element_line(color = 'black'),
      axis.line.y.left = element_line(color = 'black'),
      legend.margin = margin(-0.15, -0.15, -0.15, -0.15, unit = 'cm')
    )
  # Add legend and legend subtitle
  if (study == '1') plot = plot + guides(fill = guide_legend('Target\nStudy 1'))
  else plot = plot + guides(fill = guide_legend(paste('Study', study)))
  plot = plot +
    guides(
      # colour = guide_legend('Study 1'),
      colour = FALSE
      # fill = guide_legend()
    ) +
    labs(
      # legend
      # fill = paste('Study', study),
      # colour = paste('Study', study)
    ) +
    # Add title
    ggtitle(paste('Study', study))
  return(plot)
}
# Transform data for each study
final_data_1_raw = transform_data(1, c("WTR_37","WTR_23","WTR_75","WTR_19","WTR_46","WTR_68","WTR_total"))
final_data_2 = transform_data(2, c("WTR_75","WTR_19","WTR_46","WTR_total"))
final_data_3 = transform_data(3, c("WTR_75","WTR_19","WTR_46","WTR_total"))
# Generate plot for each study
plot_1 = generate_plot(1, c('Romantic Partner', 'Friend', 'Acquaintance', 'Enemy'), c('#3f6485', '#b17276', '#829966', '#e89b39', 'black'), c('#3f6485', '#b17276', '#829966', '#e89b39', 'black'))
plot_2 = generate_plot(2, c('Friend', 'Acquaintance', 'Enemy'), c('#b17276', '#829966', '#e89b39', 'black'), c('#b17276', '#829966', '#e89b39', 'black'))
plot_3 = generate_plot(3, c('High-Value Stranger', 'Low-Value Stranger'), c('#6494a3', '#86699a', 'black'), c('#86699a', '#6494a3', 'black'))
# Combine plots together
plot_1 + plot_2 + plot_3 + guide_area() + plot_layout(nrow = 2, guides = 'collect') +
  plot_annotation(caption = 'Fig. 2. Relations between welfare-trade-off ratio (WTR) and comfort with infection risky contact, separately for each target type in Studies 1 through 3.
        In each scatterplot, the solid line indicates the best-fitting relationship between WTR and contact comfort across categories, and the dashed lines indicate the 
        best-fitting relationship within each target category. Shaded areas around the regression lines indicate 95% confidence intervals.',
                  theme = theme(plot.caption = element_text(hjust = 1, family = 'Times New Roman'))) 

3a. Are there gender differences in the threshold for disgust?

My experience is that women tend to feel higher levels of disgust compared to men. This is supported in the literature where women have continually been found to have noticably higher levels of disgust compared to men, across many different types of disgust, such as insects, open sores, dirty clothing, and faeces (Curtis, Aunger, and Rabie 2004; Prokop and Jančovičová 2013), and animals, death, hygiene, food, and sex (Haidt, McCauley, and Rozin 1994). In (Tybur et al. 2020), disgust was measured on a likert scale and disgust statements related to dog faeces, red sores, sweaty palms, mouldy food, body odor, cockroaches, and bloody cuts.

According to the codebook, the gender variable is sex (where 1 = male, 2 = female, 3 = other). The disgust scale is DS1 to DS7. First, I checked that the data had these variables - it does. However, it’s spread across study 1-3.

First, I made a subset of the data to include only the variables I was interested in: sex and DS1-DS7, using select(). I kept participant number there too because I think it’s helpful. I also removed sex = 3 using filter() because the authors’ suggested this data wasn’t very accurate. Finally I made a new column called exp to keep track of which study each participant had come from.

# Make a subset of the original study 1 dataset to only include participant, sex, and the disgust variables

data_1_DSxSex <- data_1_raw %>%
  select(participant, sex, DS1:DS7) %>% 
  filter(sex != '3') %>% # Remove 'other' as a gender category
  mutate(exp = 1) # Make a new variable based to identify the study it came from

# Make a subset of the original study 2 dataset to only include participant, sex, and the disgust variables

data_2_DSxSex <- data_2 %>%
  select(participant, sex, DS1:DS7) %>% 
  filter(sex != '3') %>% # Remove 'other' as a gender category
  mutate(exp = 2) # Make a new variable based to identify the study it came from

# Make a subset of the original study 3 dataset to only include participant, sex, and the disgust variables

data_3_DSxSex <- data_3 %>% 
  select(participant, sex, DS1:DS7) %>% 
  filter(sex != '3') %>% # Remove 'other' as a gender category
  mutate(exp = 3) # Make a new variable based to identify the study it came from

I then combined these subsets into one dataset using rbind() There are 7 disgust variables and it was impossible (for me at least) to work with them all separately, but I thought their mean would give a more holistic account of the data than just one variable. I used mutate() to make a new variable avg_DS which merged the 7 DS columns into one by taking their row means (rowMeans()). I then piped (%>%) this onto a new dataset with the same variables before except replacing DS1:DS7 with avg_DS and called this avgdisgust. Working out how to do this took me so long, but then as I got to know my data, I worked it out! I then made a summary table of the descriptive statistics. Working out the best way to do this took me hours and many meltdowns. I couldn’t work out how Jenny had done it despite being at the Q&A and re-watching her code it. Finally came up with what is quite a pretty descriptives table! Thanks to summarise() and gt().

Total_DSxSex <- rbind(data_1_DSxSex, data_2_DSxSex, data_3_DSxSex) # Combine these 3 new datasets to a total dataset to analyse

# Make a new variable that is a combined average of all disgust variables
avgdisgust <- Total_DSxSex %>% 
  mutate(avg_DS = rowMeans(select(., DS1:DS7)))  %>%
  select(participant, avg_DS, sex) # Use a subset of the total dataset to include only participant, the average disgust variable, and sex

Total_DSxSexSummary <- avgdisgust %>% # Make a summary table of descriptive statistics from the total data
  group_by(sex) %>%
  summarise (mean = mean(avg_DS),
             sd = sd(avg_DS),
             n = n(),
             se = sd/sqrt(n)
            )

 Total_DSxSexSummary %>%
  gt() %>% 
  fmt_number(
    columns = 2:5, 
    decimals = 2)

Well… according to this table, females feel higher levels of disgust on average (4.81) compared to males (4.46).

I then plotted this to have a visual comparison using ggplot() + geom_col(), before doing some inferential statistical testing.

# Plot sex (males and females) as a function of disgust in a box plot

avgdisgust %>%
  ggplot(aes(x = sex, y = avg_DS, fill = sex)) +
  geom_boxplot(alpha = 0.5) +
  stat_summary(fun = "mean") +
    scale_x_discrete(name = "Sex", 
                     labels = c("1" = "Males", "2" = "Females")) +
    scale_y_continuous(name = 'Average Disgust') +
    (easy_text_size(15))

Yep.. females are definitely still higher! Now let’s see if this is statistically significant.

# Complete a t-test to compare the means between average disgust and sex

question1_ttest <- t.test(formula = avg_DS ~ sex, data = avgdisgust)

# Write the t-test findings in a sentence for a report

report(question1_ttest)

Effect sizes were labelled following Cohen’s (1988) recommendations.

The Welch Two Sample t-test testing the difference of avg_DS by sex (mean in group 1 = 4.46, mean in group 2 = 4.81) suggests that the effect is positive, statistically significant, and small (difference = 0.35, 95% CI [-0.45, -0.25], t(1784.16) = -7.02, p < .001; Cohen’s d = -0.33, 95% CI [-0.43, -0.24])

Nice! p is very small and certainly < .05. I’ll take this as a significant result. It turns out that there are sex differences in average disgust felt, with females feeling more disgusted than males (p<.05).

3b. Are their age differences in the personality trait Honesty-Humility?

(Tybur et al. 2020) spoke about honesty-humility being a personality trait that contributed to a higher interpersonal value (i.e. people who score high on honesty-humility are more likable). They also spoke about the elderly being valued less (along with obese or physically disabled people, or immigrants). I wondered if there was any age differences in honesty-humility. The authors don’t specify how old is elderly, so I’ll just look generally at age differences and see what I find. A study by (Kawamoto, Shimotsukasa, and Oshio 2020), studying Japanese also found that the level of Honesty-Humility was positively associated with age. I.e. Honesty-Humility levels increase with age. This doesn’t quite align with (Tybur et al. 2020) hypothesis…Let’s have a look!

According to the codebook, the age variable is age (where 1 = male, 2 = female, 3 = other). The Honesty-Humility scale is HH1 to HH10. First, I checked that the data had these variables - it does. However, it’s only found in study 1 and 2.

First, I made a subset of the data to include only the variables I was interested in: age and HH1-HH10, using select(). I kept participant number there too because I think it’s helpful. I also removed any missing values using na.omit(). Finally I made a new column called exp to keep track of which study each participant had come from.

# Making a new dataset from data 1 with participants, the Honesty-Humility variables and age as the columns
  
data_1_HHxAge <- data_1_raw %>%
  select(participant, HH1:HH10, age) %>%
  na.omit() %>% # Remove any missing values
  mutate(exp = 1) # Make a new variable based to identify the study it came from

# Making a new dataset from data 2 with participants, the Honesty-Humility variables and age as the columns

data_2_HHxAge <- data_2 %>%
  select(participant, HH1:HH10, age) %>%
  na.omit() %>% # Remove any missing values
  mutate(exp = 2) # Make a new variable based to identify the study it came from

I then combined these subsets into one dataset using rbind()

# Combine these 3 new datasets to a total dataset to analyse

Total_HHxAge <- rbind(data_1_HHxAge, data_2_HHxAge) 

There are 10 Honesty-Humility variables and it was impossible (for me at least) to work with them all separately, but I thought their mean would give a more holistic account of the data than just one variable. I used mutate() to make a new variable avg_HH which merged the 10 HH columns into one by taking their row means (rowMeans()). I then piped (%>%) this onto a new dataset with the same variables before except replacing HH1:HH10 with avg_HH and called this avghonestyhumility.

# Make a new dataset that combines two newly created variables: a combined average of all honesty-humility variables and an age variable that separates age into 5 groups

avghonestyhumility <- Total_HHxAge %>% 
  mutate(avg_HH = rowMeans(select(., HH1:HH10)))  %>%
  mutate(AgeGroup = case_when(age >= 18 & age <=29 ~ '18-29',
                              age >= 30 & age <= 39 ~ '30-39',
                              age >= 40 & age <= 49 ~ '40-49',
                              age >= 50 & age <= 59 ~ '50-59',
                              age >= 60 & age <= 79 ~ '60-73')
  ) %>%
  select(participant, avg_HH, AgeGroup) # Use a subset of the total dataset to include only participant, the average contact comfort variable, and sex

I then made a summary table of the descriptive statistics using summarise() and gt().

# Make a summary table of descriptive statistics from the total data

Total_HHxAgeSummary <- avghonestyhumility %>%
         group_by(AgeGroup) %>%
           summarise (mean = mean(avg_HH),
                      sd = sd(avg_HH),
                      n = n(),
                      se = sd/sqrt(n)
                      )

 Total_HHxAgeSummary %>%
  gt() %>% 
  fmt_number(
    columns = 2:5, 
    decimals = 2) 

According to this table, ‘the elderly’ which I would say would be somewhere in the 60-73 age group category, actually had the highest Honesty-Humility on average (3.55) with people aged 18-29 having the lowest (3.12). It seems this trait actually increases with age!

I then plotted this to have a visual comparison using ggplot() + geom_boxplot(), before doing some inferential statistical testing.

# Plot age as a function of honesty-humility in a box plot

  avghonestyhumility %>%
  ggplot(aes(x = AgeGroup, y = avg_HH, fill = AgeGroup)) +
  geom_boxplot(alpha = 0.5) +
  stat_summary(fun = "mean") +
  scale_x_discrete(name = "Age Group") +
  scale_y_continuous(name = 'Average Honesty-Humility') +
  (easy_text_size(15))

Wow! Look at that little trend…Honesty-Humility definitely seems to be increasing across age and the ‘elderly’ category is highest! Now let’s see if this is statistically significant.

# Complete an ANOVA test to compare the means between average honesty-humility and age group
  
question2_anova <- aov(avg_HH ~ AgeGroup, data = avghonestyhumility)

# Write the t-test findings in a sentence for a report

report(question2_anova)  

## The ANOVA (formula: avg_HH ~ AgeGroup) suggests that:
## 
##   - The main effect of AgeGroup is statistically significant and small (F(4, 925) = 3.70, p = 0.005; Eta2 = 0.02, 90% CI [2.75e-03, 0.03])
## 
## Effect sizes were labelled following Field's (2013) recommendations.

Well..It’s a small effect but it’s certainly significant! Looks like there are age differences but they’re more in line with (Kawamoto, Shimotsukasa, and Oshio 2020) findings than (Tybur et al. 2020).

3c. Do people in a relationship have higher levels of contact comfort than those in platonic relationships?

I would assume that people in a relationship would have higher levels of contact comfort than those who aren’t. (Tybur et al. 2020) measured comfort from 1 = very uncomfortable, to 7 = very comfortable. This ranged from using the same water bottle, deodrant stick, comb, towel, hat, and earbuds as the target, as well as touching their used hankerchief or socks, and sitting next to them while they were sneezing or coughing, or eating the same sandwich as them. I would suggest that a person who is in a relationship with someone would be far more comfortable doing many of these things, than someone who is single.

According to the codebook, the relationship variable is relat (where 1 = yes currently in a romantic relationship, 2 = no not currently in a romantic relationship, 3 = other). The contact comfort scale is comf1 to comf10. First, I checked that the data had these variables - it does. However, it’s spread across study 1-3.

First, I made a subset of the data to include only the variables I was interested in: relat and comf1-comf10, using select(). I kept participant number there too because I think it’s helpful. I also removed missing variables using na.omit(). Finally I made a new column called exp to keep track of which study each participant had come from.

# Making a new dataset from data 1 with participants, the contact comfort variables and relat as the columns

data_1_ComfxRelationship <- data_1_raw %>%
  select(participant, comf1:comf10, relat) %>%
  na.omit() %>% # Remove any missing values
  mutate(exp = 1) # Make a new variable based to identify the study it came from

# Making a new dataset from data 2 with participants, the contact comfort variables and relat as the columns

data_2_ComfxRelationship <- data_2 %>%
  select(participant, comf1:comf10, relat) %>%
  na.omit() %>% # Remove any missing values
  mutate(exp = 2) # Make a new variable based to identify the study it came from

# Making a new dataset from data 3 with participants, the contact comfort variables and relat as the columns

data_3_ComfxRelationship <- data_3 %>%
  select(participant, comf1:comf10, relat) %>%
  na.omit() %>% # Remove any missing values
  mutate(exp = 3) # Make a new variable based to identify the study it came from

I then combined these subsets into one dataset using rbind() There are 10 Contact Comfort variables and it was impossible (for me at least) to work with them all separately, but I thought their mean would give a more holistic account of the data than just one variable.

# Combine these 3 new datasets to a total dataset to analyse

Total_ComfxRelationship <- rbind(data_1_ComfxRelationship, data_2_ComfxRelationship, data_3_ComfxRelationship) 

I used mutate() to make a new variable avg_cc which merged the 10 comf columns into one by taking their row means (rowMeans()). I then piped (%>%) this onto a new dataset with the same variables before except replacing comf1:comf10 with avg_cc and called this avgcomf.

# Make a new variable that is a combined average of all contact comfort variables

avgcomf <- Total_ComfxRelationship %>% 
  mutate(avg_cc = rowMeans(select(., comf1:comf10)))  %>%
  select(participant, avg_cc, relat) # Use a subset of the total dataset to include only participant, the average contact comfort variable, and sex

I then made a summary table of the descriptive statistics using summarise() and gt().

# Make a summary table of descriptive statistics from the total data

Total_ComfxRelationshipSummary <- avgcomf %>%
  group_by(relat) %>%
  summarise (mean = mean(avg_cc),
             sd = sd(avg_cc),
             n = n(),
             se = sd/sqrt(n)
             ) 

 Total_ComfxRelationshipSummary %>%
  gt() %>% 
  fmt_number(
    columns = 2:5, 
    decimals = 2)

Well… according to this table, people in romantic relationships do have higher levels of contact comfort on average (3.39) compared to males (3.03).

I then plotted this to have a visual comparison using ggplot() + geom_col(), before doing some inferential statistical testing.

# Plot sex (males and females) as a function of relationship type in a bar plot with error bars

Total_ComfxRelationshipSummary %>%
  ggplot(aes(x=relat,
             y=mean,
             fill=relat)
  ) +
  geom_col() +
  geom_errorbar(aes(
    ymin = mean - se,
    ymax = mean + se),
    width = .2) +
  scale_x_continuous(
    breaks = c(1:2),
    labels = c('In a romantic relationship', 'Not in a romantic relationship'),
    name = "Relationship Status"
  ) +
  scale_y_continuous(
    expand = c(0,0),
    limits = c(0, 6),
    name = "Average Contact Comfort") +
  easy_text_size(15) +
  theme(legend.position = "none")

Oh yeah, definitely some sort of difference! Being in a romantic relationship definitely seems to produce higher average contact comfort levels compared to not being in one! Wonder if it’s statistically significant?

# Complete a t-test to compare the means between comf2 and relationship

question3_ttest <- t.test(formula = avg_cc ~ relat, data = avgcomf)

# Write the t-test findings in a sentence for a report

report(question3_ttest) 

## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of avg_cc by relat (mean in group 1 = 3.39, mean in group 2 = 3.03) suggests that the effect is negative, statistically significant, and small (difference = -0.36, 95% CI [0.20, 0.52], t(1223.07) = 4.42, p < .001; Cohen's d = 0.25, 95% CI [0.14, 0.37])

Wow! Looks like it’s significant too, only small but certainly significant! I actually had no research backing for this, it just seemed like something I’ve found anecdotally. Of course, there’s some chance this is not actually significant so more testing would need to be done, but an interesting find nonetheless!

1. Include a analysis code that produces everything

Although it seemed lovely (especially compared to the other groups) that (Tybur et al. 2020) had included their analysis code, it didn’t have everything. In an ideal world, we would reproduce an article that gave us all the data in the right form, all the variable names, and all the analysis script. I mean ideally we could just download their files of OSF, import them into our RStudio, load their csv files, and run the script, and reproduce EVERYTHING! This didn’t happen. Although their was analysis script code, it only included script for the demographic descriptive statistics. Nothing for the figures! This meant that actually, we were reproducing alot of the code we needed, from scratch, just like the other groups! Indeed, Sharla Gelfand in her talk, “Don’t repeat yourself, talk to yourself,” she suggested going so far as making README set up files and explaining everything you ran and every file you used. Indeed, explaining what you’re doing, why you’re doing it, and what to do next. This would have been very helpful for reproducibility and may allow for more reproducibility to be done, due to massive time-saving.

2. Include an accurate codebook with clear variable descriptions

Again, although we had a codebook, it wasn’t very good. It wasn’t accurate, and many of the variable names were mixed up. This is an easy mistake but can get very confusing and messy for someone who hasn’t worked with your data, and can’t call you to chat and ask you about it or what it means. One example of this is in the study 2 section of our codebook, target_sex is described as Target Relationship, target_category is described as how long they’ve known the target and relate_length is described as whether the target is a man or a woman. This mad it very difficult for us to understand what each variable meant. Of course, these were eventually pretty self-explanatory to match up, but it may not be for less self-explanatory variables.

3. Use consistent code style, formatting etc.

I have absolutely no prior experience in coding. I learned what Danielle Navarro taught us (library(tidyverse)) and that is it. Their coding style was clearly very different and to be honest, I couldn’t understand any of it. Their commenting was almost non-existent so I couldn’t read that to understand it. Finally, it was very lengthy. They repeated code for very similar purposes, and it could have been put into some functions to simplify it considerably. I say this without actually being able to do this myself, however, I’m also not publishing my work to Psychological Science! Apart from this, their code was just inconsistent. Brackets were on different lines and both single and double quotation marks were used interchangably. These things are very much nit-picking, however, it did lead to difficulties in reproducibility and a significantly longer time to reproduce, because the code was so complicated.