Part 1: Summary & Reaction

Summary

In today’s media, news readers are often exposed to conflicting headlines on health and nutrition. For example:

  • Fasting diets could reverse diabetes, researchers say.
  • Fasting diets may raise risk of diabetes, researchers say.

Conflicting headlines such as the above lead to confusion and backlash against scientific advice. This is because the term “researchers say” leads readers to think “all researchers”, resulting in news readers thinking scientific advice is unreliable, confusing, and contradictory. An implication of this confusion is news readers disregard credible health advice.

Previous studies have found adding the qualifier “some” can help readers to infer “some but not all”. The phrase “some researchers” implies absence of complete consensus, as researchers are not a homogeneous group. A gap in literature that remains unanswered is whether qualifiers can be used to reduce perceived conflict.

Haigh & Birch (2021) thus aimed to investigate whether hedging news headlines by adding the qualifier “some” before “researchers” can reduce readers’ perceived conflict when reading conflicting headlines. The authors hypothesized participants would perceive greater conflict when reading conflicting news headlines, and the use of qualifiers would reduce perceived conflict.

In two 2 x 2 between-subjects experiments, participants read either a series of conflicting or non-conflicting headlines on health and nutrition, in either a generic (“Researchers say”) or qualified (“Some researchers say”) headline format. After reading the headlines, participants in Experiment 1 self-rated their perceived contradiction, confusion, and scientific advancement using Likert scales. Experiment 2 measured general beliefs on science, where participants self-rated their nutritional confusion, nutritional backlash, mistrust of expertise, confidence in the scientific community, and beliefs about development and certainty of knowledge.

Unfortunately, the authors’ findings did not support their original hypothesis that participants who read qualified, conflicting headlines would rate these headlines as less conflicting. Compared to participants who saw non-conflicting headlines, Experiment 1’s participants who saw conflicting headlines rated them as more contradictory, confusing, and resulted in them knowing less about how to be healthy (i.e., conflict manipulation was effective). Conversely, Experiment 2 found no effect of the conflict manipulation (i.e., participants did not rate conflicting vs. non-conflicting headlines differently). In both experiments, adding the qualifier “some” did not reduce perceived conflict.

While qualifiers did not reduce perceived conflict, finding effective strategies to reduce perceived conflict remains practically important. Strategies resulting in less reader confusion and conflict are beneficial to ensure news readers follow credible health advice. Therefore, the authors recommended future studies examine the effects of other hedging strategies on reducing perceived conflict, such as adding uncertainty terms (i.e., may, might) and reporting effect sizes.

Reactions

  1. I can see that this paper relates to the issue of conflicting headlines that are in abundance on social media, blogs/news websites, and magazines/newspapers/TV. I really enjoyed reading the authors’ paper and found their study on the use of qualifiers in reducing perceived confict practically important. It was unfortunate that the use of qualifiers in hedging headlines - which I thought was a simple and elegant solution - did not yield any effect on reducing perceived conflict. However, I am curious to learn whether future studies on hedging headlines by using uncertainty terms (e.g., may/might/could) and reporting effect size will show more promising effects on reducing news readers’ perceived conflict.

  2. I wonder whether the authors could have also included questions that assessed participants’ attitudes before and after reading the headlines. For example, asking participants “How likely are you to take Vitamin C” before and after headlines on Vitamin C are presented, and assessing the differences in scores pre- and post-manipulation could reveal additional insight on participants’ attitudes and behaviours. Furthermore, pre- and post- manipulation measures can reveal more accurate results, as it takes into account individual differences in ratings (i.e., some participants might generally rate higher - using differences in scores by comparing with participants’ baseline scores will remove this effect). The authors also stated assessing “engagement in healthy behaviour” was one of the rationales for conducting the study. Yet, their questions to participants were worded generally and did not assess specific health behaviours. To me, it seems like there was a missed opportunity to include additional questions framed in the context of participants’ health behaviours.

  3. I was surprised that participants from a wide age range were recruited in this study, which I thought might introduce potential confounding variables. From personal experience, I have noticed differences in attitudes between the older and younger generations when it comes to advice on health and nutrition. Yet, the researchers did not factor in age when analysing their data. Specifically, I would have thought participants in older age groups would be more likely to trust the health advice (lower mistrust) and have more confidence in the scientific community. Meanwhile, I expected younger participants to have higher mistrust due to today’s widespread misinformation and conflicting headlines. However, the authors made no comment on how different demographic variables may influence their measures. Since their findings were not significant, this made me wonder whether the wide age range of 18-73 could have played a factor in their (lack of) findings.

Part 2: Verification

Reproducibility Goals

I will be going into each one (with screenshots in more detail) with each section below. However, for a quick overview we had to reproduce:

  • Descriptive statistics for demographics in Experiment 1 & 2
  • Means for dependent variables measured in Experiment 1
  • Figures 1, 2, 3

Setting Up

First, we load the packages that we will be using throughout this report.

library(tidyverse) 
library(ggeasy)
library(ggbeeswarm)
library(patchwork)
library(gt) 
library(rstatix) 
library(corrplot)
  • tidyverse contains a set of packages that allows for ‘tidy data’. It includes dplyr which we will use for data manipulation, and ggplot which is a plotting package we will use for data visualisation
  • ggeasy allows easy customisation of our ggplot graphics
  • ggbeeswarm helps us to visualise density of data on our plots (i.e., add a swarm plot)
  • patchwork allows us to combine plots and specify the display layout of combined plots
  • gt is used for tables
  • rstatix is used for running statistical tests (i.e., calculating effect size)
  • corrplot is used for producing correlation matrices

Next, we load the data from CSV files using readcsv() which reads the CSV into a tibble. We faced challenge #1 here where the experiment data files were named “Study 7 data” and “Study 8 data” This was not in accordance to the paper which reported findings from Experiment 1 and Experiment 2.

Lesson #1 on reproducibility: Make sure file names are clearly labeled and easily understandable at first glance!

This took awhile to figure out - we had to go into each spreadsheet separately and look at the variables recorded. By matching the variables measured to each experiment, we finally figured out:

Experiment 1 = Study 8 data.csv; Experiment 2 = Study 7 data.csv

Next, I’m creating new objects expone and exptwo using the assignment operator <- for the data in Experiment 1 and Experiment 2, respectively. Using %>% (also called “a pipe operator”), I can forward a value into the next function and pipe the data through a series of steps!

Once I discovered %>%, I knew I would be using it a lot as the pipe operator is really great for the DRY (Don’t repeat yourself) principle!

Specifically, after readcsv(), I am going to use another function rename() which is part of the dplyr package and allows us to change names of individual variables. This was the result of looking at the data file and seeing some variables that do not have human readable names (enter: challenge #2). Upon further inspection, it looks like these variables (SC0 and FL_10_D0/FL_12_D0) correspond to variables headline recall score and participants’ condition, respectively! To avoid later confusion, I know my future self will thank me for renaming them to ‘recall_score’ and ‘condition’ now.

Lesson #2 on reproducibility: Make sure variable names have clear, easily understandable labels!!!

Finally, the first two rows of the data file do not contain data but other irrelevant information due to the survey software the researchers used. As such, I am using slice() which is from dplyr, typically used to select/remove/duplicate rows. To remove the first two rows, I’m specifying the argument as -1:-2 (the minus sign denotes remove). The colon operator : allows us to define a sequence (in this case, 1 to 2).

Some of our group members decided not to use slice() and instead replicate the authors’ code where the qualtRics package was loaded. qualtRics automatically removes the first 2 rows, converts variables to numeric/factor automatically, among other things (more details here: https://cran.r-project.org/web/packages/qualtRics/vignettes/qualtRics.html). I personally decided against this, as I am unfamiliar with qualtRics and wanted greater control over my data. Specifically, as someone who is completely new to R and coding, I wanted to know exactly what was happening to my data at every single time point - for me, I decided to KISS (keep things simple).

I’m also setting the theme using theme_set() which overrides the current default theme, and specifying theme_bw() (a simple black and white theme). I will be using this theme for all the plots I will be creating using ggplot later on, so this line of code will save me adding a line in each plot (aka, good for the DRY (Don’t Repeat Yourself) principle).

expone <- "Study 8 data.csv" %>% 
  read_csv() %>% 
  rename(
    recall_score = SC0, #nicer name for recall score
    condition = FL_10_DO #nicer name for condition
  ) %>% 
  slice(-1:-2) #removes rows 1 to 2

exptwo <- "Study 7 data.csv" %>% 
  read_csv() %>% 
  rename(
    recall_score = SC0, 
    condition = FL_12_DO 
  ) %>% 
  slice(-1:-2)

theme_set(theme_bw()) #set theme

Descriptive Statistics for Demographics - Experiment 1

Reproducibility goal: To reproduce the descriptive statistics for total number of participants, sex (number of males/females), age (mean, SD, range).

The researchers had a mishap with Experiment 1. Due to a technical error, 59 participants completed the survey twice. Luckily, unique profile IDs Prolific_PID were captured and the data file lists completion in chronological order. We use the $ operator to specify a variable in the data frame (i.e., expone$Prolific_PID refers to participant IDs in Experiment 1).

Researchers reported 312 participants signed up. I decided to use n_distinct() from dplyr. It counts unique values, and is similar to length(unique(x)) from baseR (but in my opinion better, as it’s more concise!)

n_distinct(expone$Prolific_PID)
## [1] 312

We now need to exclude the duplicates (second attempts), and only keep the first attempts. Using square brackets [] (also known as extraction operators) allows us to specify a logical vector and extract specific elements. duplicated() from baseR then finds duplicate values in expone$Prolific_PID. The chunk below lists all the duplicate IDs (a total of 59).

# List all duplicate IDs
expone$Prolific_PID[duplicated(expone$Prolific_PID)] 
##  [1] "5cd1836da6f34300017e240c" "5c76f2d92819ab0015b94b4c"
##  [3] "5c72de7b96a9600001870966" "5c3d48955fd1050001a99364"
##  [5] "5eb18525b95d6127da6815ee" "5eaac092d8a78e0172680a0e"
##  [7] "5edff6e7c53e0f33aed588ce" "5d6a55cb78fce300014e078e"
##  [9] "5e8b765e88065404c4d50fe2" "5ed22fc8bc5d7c01191e3f78"
## [11] "5ea4397b9b06662614d6e4f6" "5d371c527a04ed00018ac1c8"
## [13] "5e9c29bb77f86e13d6435f3d" "5ea3ef637e4bd537fe45cf6c"
## [15] "5ebb0f6aff8da30d1f29dc29" "5cd59eea8618af0001bdab4b"
## [17] "5eac3189c4a262061aed16d5" "5e3842b92fdfd2000fc286a4"
## [19] "5ec504eab4fbc3423f29cd97" "5e332f72c43078000cda48dd"
## [21] "5ea5e76a147a4909b3ff0482" "5c28b31a0091e40001ca5030"
## [23] "5eb089604b5931137dd66a19" "5ed0d82c81c0bd1bccc8ceea"
## [25] "5ed94fdd9a2ae04cf21bfa48" "5be1cc598c6a19000137a503"
## [27] "5cd813aca9a8c4001963a0aa" "5e925105a981b55934a34813"
## [29] "5eb30ae9bb223e0e176edecd" "5ec637f5e0bca2000ae6f211"
## [31] "5eb179d6d6fcb726f9275f97" "5dcbd8542bdeaa8740d52630"
## [33] "5b570f4cc146600001b82d8d" "5ed793104268812282fdc90d"
## [35] "5edf50e72ef80a1fe0267aeb" "5d7f598628843a00181eb444"
## [37] "5ec2ec9fdaef0d11e3109f74" "5d1290103b20b0000102e8e3"
## [39] "57eee744e62704000199d5ec" "5e6a839e2fd3b003a0f7f248"
## [41] "5e7365f19674532b961c2bb6" "5eb94e5731298c01178531c0"
## [43] "5e9db2995b38950c6f669b55" "584bb2b8bd873800015531da"
## [45] "5eac351c63858608351866b4" "5ec0260375bf15077a00e645"
## [47] "5cc3289b9e21e200015f0bb3" "5eb16a0602af57258fd3f8e1"
## [49] "5ed7578ef05e671db283844e" "5ebc1ae1ea22c801479541ab"
## [51] "5edfd94cb54d22309545cf06" "559ab96cfdf99b219a612bcf"
## [53] "5eac35df3043c62536d14f14" "5e9ff2a0cf50621a9b17c94f"
## [55] "5e5d7f349238db09c60bcbab" "5ea611dda778214a5e89fbf2"
## [57] "5c62d8e2a34174000187a003" "5ed6a937eb466b1029493c39"
## [59] "5ebfabc7676c2502837188cf"

Next, we will remove the duplicates from subsequent analyses by again specifying a logical vector using [], where exclamation mark ! essentially means ‘not’/‘exclude’. expone[!duplicated(expone$Prolific_PID), ] thus means we will be selecting all UNIQUE Prolific_PID from expone. I’m then overwriting expone with only unique IDs using <-.

# Removing second attempts for 59 duplicate IDs 
expone <- expone[!duplicated(expone$Prolific_PID), ]

We need to apply the pre-registered exclusion criteria: did not complete the task (Finished=0), did not respond seriously (Serious_check=0), failed attention check (recall_score <4). I am going to create this data set as a new object exponefinaldata, by using expone and piping with filter() from the dplyr package to pick cases based on their values stated in the chunk below. I can then specify the arguments based on exclusion criteria, == being the equality operator used to state the value in the data frame that we want using "".

Finally, I use count() to check the total number of participants remaining after exclusion criteria is applied. This should correspond to N = 294.

exponefinaldata <- expone %>% 
  filter(
    Consent == "1", #filter to include those who Consented
    Finished == "1", #filter to include those who Finished
    Serious_check == "1", #filter to include those who answered they passed Serious Check
    recall_score >= "4", #filter to include Recall Score of 4 and above only
    )

exponefinaldata %>% 
  count(Serious_check)
## # A tibble: 1 x 2
##   Serious_check     n
##   <chr>         <int>
## 1 1               294

I am also tidying our final data exponefinaldata for ease of use. The original data frame contains many unnecessary variables - 340 to be exact! This makes getting a glimpse of the data using glimpse() difficult, which we will be using a lot in RStudio’s Console, as we analyse the data and reproduce the plots.

Challenge #3: Data is quite untidy, and could be improved by better workflow!

Lesson on learning R/coding: Taking time to tidy the data is a worthwhile investment of time and effort.

Since only a fraction of variables are interesting and useful, I am going to tidy the data using select() from dplyr to pick out the variables based on their names so we can keep keep all rows and RETAIN only the specified columns. To do this, I used glimpse() in the console to view all variables, allowing me to pick out the important ones, and since some were listed in sequence I could use the operator : to specify sequence.

By default, all the data is in ‘chr’ which means Character. From Jenny S’s Q&A, she mentioned that it is best practice to convert the data from character to factor/numeric/etc, if they will be later used for subsequent analyses. This avoids future problems like when averaging values (a numeric variable is needed, not character). Factors are essentially categories, allowing us to categorize and represent categorical data that have multiple levels. Numerics are numeric values. Based on this, I am going to change variables from Character to Numeric/Factor accordingly, using use mutate() and across() (both from dplyr). In this case, we are using mutate() to overwrite the variables (specified within c() which creates a vector) as numeric via as.numeric or factor via as.factor. Using across() makes it easy to apply the same transformation across multiple columns, and saves time and lines of coding!

Mutate() will soon be our new best friend! While we used it to overwrite existing data in this instance, it is typically used to create NEW variables within the original data set. Examples to follow!

exponefinaldata <- exponefinaldata %>% 
  select(Finished, Gender, Age, Serious_check, recall_score, condition, contradiction_1:advancement)

# Change variables from Chr to Numeric/Factor accordingly
exponefinaldata <- exponefinaldata %>% 
  mutate(
    across(c(Age,contradiction_1:advancement), as.numeric),
    across(c(Gender, condition), as.factor)
         )

I had previously used as.numeric and as.factor to change each variable one by one! Thanks to Jenny S’s guidance and introduction to across() I am now able to do it within 2 lines of code. :) Case in point (screenshot of what I was able to remove below), much more aligned to DRY principle!

Lesson on learning R: If it looks repetitive, it probably is… (Find a way to make it DRY and KISS!)

Glimpsing the data, I notice that another label that isn’t clear is for the Gender variable. Males/females are listed as 1 and 2 - but it is unclear which is which!! Another challenge accepted - I will need to use count() to see how many males/females there are and match them back to the original reproducibility goal (126 males, 168 females).

exponefinaldata %>% 
  count(Gender)
## # A tibble: 2 x 2
##   Gender     n
##   <fct>  <int>
## 1 1        126
## 2 2        168

From the above output, we can work backwards to see that 1 = male; 2 = females. However, I am not going to remember this later, so let’s change them into ‘male’ and ‘female’ now. mutate() will be used to modify/overwrite the column Gender, and case_when() from dplyr allows us to specify that when Gender is 1, ~ (the tilde tells R what value we want in the cell instead) as “Male” instead. When Gender is 2, recode as “Female” instead. Just to check my work, I am using count again to show that instead of 1/2, it is now human readable (and part of our reproducibility goal for descriptive statistics checked off)!!

exponefinaldata <- exponefinaldata %>% 
  mutate(
    Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female"
    )
)

# Count males and females
exponefinaldata %>% 
  count(Gender)
## # A tibble: 2 x 2
##   Gender     n
##   <chr>  <int>
## 1 Female   168
## 2 Male     126

After all this tidying the data, we can now easily use summarise() from dplyr to calculate summary statistics: mean, SD, min and max for Age. Using gt() from the gt package we can display the data in a table, and fmt_number() (also from gt pacakge) further allows us to format the numeric values within the table in 2 decimal places for Mean and SD (using vars() to select those specific columns). And finally using tab_header() from gt package we add a header with title.

Jenny S introduced us to the gt(), which I decided to adopt! Previously I was printing the summary statistics as a data frame using as.data.frame() from baseR. However,the gt package affords better/easier customisation and flexibiility (e.g. I can easily change decimal points, add a title to the table etc.)!

The age range should be 18-69, mean 34.29, SD 12.97.

# Age descriptive statistics
exponefinaldata %>% 
  summarise(
    Mean = mean(Age),
    SD = sd(Age),
    Min = min(Age),
    Max = max(Age)
) %>% 
  gt() %>% 
  fmt_number(
    columns = vars(Mean, SD),     
    decimals = 2
  ) %>% 
  tab_header (
    title = "Descriptive Statistics for Age in Experiment 1"
  )
Descriptive Statistics for Age in Experiment 1
Mean SD Min Max
34.29 12.97 18 69

Means & Figure 1 - Experiment 1

Reproducibility goal: Reproduce the three plots below, along with the mean data.

Again, the first thing to do is get the data tidy and ready for the plots! We need to:

  1. Sum the six contradiction scores into one variable
  2. Separate the condition variable that is currently in the format “block_number_format_conflict”. We are interested in the format and conflict variables for reproducing our plots

rowwise() from dplyr allows us to compute the data frame one row at a time. This is important as we want to sum across the contradiction scores in the same row - row-wise rather than down the column-wise. We will use mutate() to create this sum of scores as a new variable, ‘contradiction’, in the data frame. To separate the condition column into 4 new columns, we use separate() from tidyr package (within tidyverse) and specify the arguments, data= with our data, col= with the Condition variable. Finally, into=c("block", "number", "format", "conflict") allows us to create a vector using c() and separate the old data into 4 new variables block/number/format/conflict. Finally, we convert format/conflict into factor variables using the same process described previously.

# Sum contradiction variables into a new variable 'contradiction'
exponefinaldata <- exponefinaldata %>% 
  rowwise() %>% 
  mutate(
    contradiction = sum(contradiction_1, contradiction_2, contradiction_3, 
                        contradiction_4, contradiction_5, contradiction_6)
  )

# Separate the data into 4 variables
exponefinaldata <- separate(
  data = exponefinaldata, 
  col = condition, 
  into = c("block", "number", "format", "conflict"))

# Convert format and conflict to Factor variables
exponefinaldata <- exponefinaldata %>% 
  mutate(
    across(c(format:conflict), as.factor)
         )

At this point, I noticed that the authors’ plots compare between “Conf.” vs. “Non-Conf” but the data is recorded as “Conflict” vs. “Consistent”. This made me wonder why the authors did not just set the data up as Conf/Non-Conf from the get-go? Either way, I’m renaming them using levels() which provides access to levels attributes of a variable, since conflict is a Factor variable.

# Rename from Conflict/Consistent to Conf./Non-Conf.
levels(exponefinaldata$conflict)[levels(exponefinaldata$conflict)=="Conflict"] <- "Conf."
levels(exponefinaldata$conflict)[levels(exponefinaldata$conflict)=="Consistent"] <- "Non-Conf."

Calculating the Mean for dependent variables

group_by() is part of dplyr and allows us to group by the variable “conflict”. Any future operations will then be performed separately for each group. Then, we can calculate the means.

exponefinaldata %>% 
  group_by(conflict) %>% 
  summarise(mean.contradiction = mean(contradiction),
            mean.confusion = mean(confusion),
            mean.advancement = mean(advancement)) %>% 
  gt() %>% 
  fmt_number(
    columns = vars(mean.contradiction, mean.confusion, mean.advancement),     
    decimals = 3
  ) %>% 
  tab_header (
    title = "Mean for contradiction, confusion, and advanement in Experiment 1"
  )
Mean for contradiction, confusion, and advanement in Experiment 1
conflict mean.contradiction mean.confusion mean.advancement
Conf. 25.259 4.524 −0.245
Non-Conf. 13.449 3.646 0.007

This corresponds to our authors’ reported for conflicting vs. non-conflicting:

Mean Confusion: 4.52 vs. 3.65

Mean Scientific Advancement: -0.25 vs. 0.007

Mean Contradiction: 25.3 vs. 13.4

Finally - to reproduce our plots!!!

Another challenge at this stage was whether to reproduce our plots using the package the authors used, which is pirateplot - or attempt to use ggplot which we had learnt from Danny’s tutorials.

Lesson on troubleshooting AND learning R: As a group, we decided to use ggplot as we wanted to develop a better understanding of how to use ggplot. On the other hand, pirateplot was very foreign to us… For us, this was a matter of KISS (keep is super simple)!

We’ll start by plotting the contradiction plot. We are using ggplot package to visualise our data. We first specify the arguments for ggplot() - by specifying the argument data. Next, aes() is a quoting function that makes it easy to work with variables and lets us specify the aesthetics mapped to variables (such as colour/shape/fill etc.). Here, we want the x-axis to be conflict, y-axis to be contradiction, and fill to be conflict group the participants were in. The authors also used a violin plot.To add a layer of geom (geometric object) to our ggplot, we use the plus sign + to add a layer, and geom_violin() for adding a violin plot. print() from baseR is the most common method to print an output - in this case I am printing the plot we just coded.

Interestingly, you do not actually need to specify data and aes. However, as a newbie to R, I thought it would help me understand better if I explicitly state each argument/function within my ggplot().

Geometric objects can also be boxplots, violin, histogram - essentially, it allows us to visually render the data! This will be really useful as we continue…

Attempt #1

contradiction1 <- ggplot(
  data = exponefinaldata,
  aes(
    x = conflict,
    y = contradiction,
    fill = conflict
  )
)  + 
  geom_violin()

print(contradiction1)

We also need to include participants’ format in the x-axis, so using facet_wrap(), which is part of ggplot2, we can wrap panels into two-dimensions. The vars() function allows us to specify which faceting variable we want - in this case “format” and strip.position argument allows us to position the labels on the bottom (by default strip is on top).

Attempt #2

contradiction2 <- ggplot(
  data = exponefinaldata,
  aes(
    x = conflict,
    y = contradiction,
    fill = conflict
  )
)  + 
  geom_violin() +
  facet_wrap(
    vars(format),
    strip.position = "bottom"
  ) 

print(contradiction2)

Getting closer! Next, the authors include a bee swarm plot (which is a one-dimensional scatter plot) on their plots, so I am using geom_beeswarm from the beeswarm package, and specifying the argument cex as 0.2 which specifies the width for each data point’s spacing. Finally, authors include 95% confidence intervals and mean on their plots - displayed as crossbar. Using stat_summary(), which is in the ggplot2 package, we can add different summary statistics on our plots. fun.data argument allows stat_summary() to operate on a data frame. To add the 95% CI, we use "mean_cl_normal" within the fun.data argument. To add the geometric object, the argument geom = "crossbar" is used for crossbars. The crossbars automatically take the same colour as the plot, but we want them in white and with some transparency added - which we can specify using fill (changes fill colour) and alpha (changes transparency) values.

Attempt #3

contradiction3 <- ggplot(
  data = exponefinaldata,
  aes(
    x = conflict,
    y = contradiction,
    fill = conflict
  )
)  + 
  geom_violin() +
  facet_wrap(
    vars(format),
    strip.position = "bottom"
  ) +
  geom_beeswarm(
    cex = 0.2 
  ) +
  stat_summary(
    fun.data = "mean_cl_normal",  #adds 95% CI
    geom = "crossbar", #specify add geometric object crossbar
    fill = "white",    #changes crossbar fill colour
    alpha = .7 #changes transparency of fill to 70%
  )

print(contradiction3)

The final part is to beautify the plot (add title/remove legend etc.). ggtitle() from the ggplot2 package allows us to modify labels on the plot. We can then center the title by specifying the theme() and changing the plot.title argument with element_text(hjust = 0.5) (0.5 specifies centred). We remove the x axis label using scale_x_discrete() (since x-axis is a discrete variable) and specify the limits argument of the y-axis and its label name using scale_y_continuous() (since y-axis is a continous variable) and the vector c(0,30) for the y min and y max.

We can easily remove the legend using easy_remove_legend() from the ggeasy package and change the colours of the plot to match the authors’ using scale_fill_manual() and specifying the fill colour values manually slategray2 and lightpink1.

ggeasy was a discovery from Jenny’s workshop! I love the ease of use, and how easy/readable the functions are. The function names literally tell you what it is doing!

Attempt #4 (aka Final Attempt)

# Plot Contradiction
contradiction <- ggplot(
  data = exponefinaldata,
  aes(
    x = conflict,
    y = contradiction,
    fill = conflict
  )
)  + 
  geom_violin() +
  facet_wrap(
    vars(format),
    strip.position = "bottom"
  ) +
  geom_beeswarm(
    cex = 0.2
    ) +
  stat_summary( 
    fun.data = "mean_cl_normal",  
    geom = "crossbar", 
    fill = "white",   
    alpha = .7 
  ) +
  ggtitle(
    label = "Contradiction"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5)
  ) +
  scale_x_discrete(
    name = NULL
  ) +
  scale_y_continuous(
    name = "Perceived Contradiction",
    limits = c(0,30)
  ) +
  easy_remove_legend() +
  scale_fill_manual(
    values = c("slategray2", "lightpink1") 
  )
            
print(contradiction)

The result is a plot that almost identically matches the authors’ plots! A great achievement considering we used ggplot rather than pirateplot!

Doing the same above process for confusion & advancement next. We only need to change the y-axis variable, titles and upper/lower limits to match the ones in the paper.

# Plot Confusion
confusion <- ggplot(
  data = exponefinaldata,
  aes(
    x = conflict,
    y = confusion,
    fill = conflict
  )
)  + 
  geom_violin() +
  ggtitle(
    label = "Confusion"
  ) +
  facet_wrap(
    vars(format),
    strip.position = "bottom"
  ) +
  stat_summary( 
    fun.data = "mean_cl_normal",  
    geom = "crossbar", 
    fill = "white",    
    alpha = .7 
  ) +
  theme(
    plot.title = element_text(hjust = 0.5) 
  ) +
  scale_x_discrete(
    name = NULL
  ) +
  scale_y_continuous(
    name = "Perceived Confusion"
  ) +
  easy_remove_legend() +
  geom_beeswarm(
    cex = 0.2  
    ) + 
  scale_fill_manual(
    values = c("slategray2", "lightpink1")
  )

print(confusion)

# Plot Advancement
advancement <- ggplot(
  data = exponefinaldata,
  aes(
    x = conflict,
    y = advancement,
    fill = conflict
  )
)  + 
  geom_violin() +
  facet_wrap(
    vars(format), 
    strip.position = "bottom"
  ) +
  geom_beeswarm(
    cex = 0.2 
    ) +
  stat_summary(
    fun.data = "mean_cl_normal",  
    geom = "crossbar",
    fill = "white",    
    alpha = .7 
  ) +
  ggtitle(
    label = "Advancement" 
  ) +
  theme(
    plot.title = element_text(hjust = 0.5) 
  ) +
  scale_x_discrete(
    name = NULL 
  ) +
  scale_y_continuous(
    name = "Perceived Scientific Advancement" 
  ) +
  easy_remove_legend()+
  scale_fill_manual(
    values = c("slategray2", "lightpink1") 
  )

print(advancement)

To display it together, we can use the following code:

contradiction + advancement + confusion

Challenge: This is not the layout we want for the plots. We want them to be displayed in two columns… So, we started troubleshooting and finding packages/functions on Google that allows us to do this.

We found the function plot_layout from the patchwork package, which allows us to do just that. We’ll name this object combinedplot and specify the argument ncol=2 for 2 columns.

combinedplot1 <- contradiction + advancement + confusion + plot_layout(ncol = 2)

plot (combinedplot1)

Challenge: The codes for our 3 plots are VERY repetitive! This is not in alignment with the DRY (Don’t Repeat Yourself) principle… This is where we discovered the world of functions! Thanks to Jenny S for walking us through the creation of functions at Q&A…

Firstly, we make our function’s code within the curly brackets {} - typically curly brackets are used to group statements. The previously used code for a single plot can be placed within {}, with the exception of certain arguments that change across the three plots above. These are: the y variable y_var, plot title plot_title, y-axis title y_title, and y-axes limits lim_1 lim_2. We can later specify these when we create our plots. We’ll also call this function figure.1.fun.

# Making our function 
figure.1.fun <- function(y_var, plot_title, y_title, lim_1, lim_2)  {
  
  ggplot(exponefinaldata,aes(x = conflict, y = y_var, fill = conflict)) +
  geom_violin() +
  ggtitle(label = plot_title) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 8), axis.title.y = element_text(size = 6)) +
  scale_x_discrete(name = NULL) +
  scale_y_continuous(name = y_title, limits = c(lim_1, lim_2)) +
  facet_wrap(vars(format), strip.position = "bottom") +
  stat_summary(fun.data = "mean_cl_normal", geom = "crossbar", fill = "white", alpha = .7) +
  easy_remove_legend() +
  geom_beeswarm(cex = 0.2) + 
  scale_fill_manual(values = c("slategray2", "lightpink1"))
             
}

Now to use the function! We simply need to specify arguments for y_var, plot_title, y_title, lim_1 and lim_2 for each corresponding plot: contradiction, advancement, and confusion. Then, combine using the patchwork package again as mentioned above.

# Plotting Contradiction, Advancement and Confusion plots using function 
contradiction.plot <- figure.1.fun(
  y_var = exponefinaldata$contradiction, 
  plot_title = "Contradiction", 
  y_title = "Perceived Contradiction", 
  lim_1 = 1, lim_2= 30
  )

advancement.plot <- figure.1.fun(
  y_var = exponefinaldata$advancement, 
  plot_title = "Advancement", 
  y_title = "Perceived Scientific Advancement", 
  lim_1 = -1, lim_2 = 1
  )

confusion.plot <- figure.1.fun(
  y_var = exponefinaldata$confusion, 
  plot_title = "Confusion", 
  y_title = "Confusion", 
  lim_1 = 1, lim_2 = 5
  )

# Combining plots 
combineplots1 <- contradiction.plot + advancement.plot + confusion.plot + plot_layout(ncol = 2)

print(combineplots1)

And voila! We get the same plots, but more succinctly!! We have saved many lines of unnecessary repeated coding in the process of using the function! Something to keep in mind for Experiment 2’s plots.

Figure 2 - Experiment 1

Next, we reproduce the histogram for Experiment 1 displaying number of participants in each condition who felt the body of research reported in the headlines resulted in us knowing more than before, less than before, or the same as before. Reproducibility goal:

This follows on from the Advancement plot, whereby participants rated -1 (less), 0 (same), or 1 (more). The first step is, using the baseR function ordered(), we can specify the levels of advancement (-1, 0, 1) and label them accordingly (less, same, more).

Next, we plot the histogram using ggplot2’s geom_bar which produces a bar plot. We specify the x-axis variable as exponefinaldata$advancement, and group/fill using both the conditions conflict and format by stating as conflict:format. Specifying the argument position="dodge" places the histogram bars side by side (the default is set to “stack”).

Finally, we format the histogram. We use the function scale_fill_grey() to specify the grey colour scheme as grey, then change the title of the legend to “Condition” and rename the short-formed labels (e.g. Conf.) to long-form (e.g. Conflict) by stating each variable using the labels argument. The x-axis title is “Advancement”, and we change the labels to capitalise (Less, Same, More); y-axis title is named “Number of Participants”.

exponefinaldata$advancement <- ordered(
  exponefinaldata$advancement,
  levels = c(-1, 0, 1),
  labels = c("less", "same", "more")
)

advancement_histogram <- ggplot(
  data = exponefinaldata,
  aes(
    x = exponefinaldata$advancement,
    group = conflict:format,
    fill = conflict:format
    )
) +
  geom_bar (
    position = "dodge" #places bars side by side
) +
  scale_fill_grey( #change colour scheme of legend to grayscale
    name = "Condition", #specify legend title
    labels = c("Conflicting/Generic",
               "Conflicting/Qualified",
               "Non-conflicting/Generic",
               "Non-conflicting/Qualified"
               )
  ) + 
  scale_x_discrete(
    name = "Advancement",
    labels = c("Less",
               "Same",
               "More"
               )
  ) +
  scale_y_continuous(
    name = "Number of Participants",
  ) 

print(advancement_histogram)

Descriptive Statistics - Experiment 2

Reproducibility goal:

Tidying up the data. The same process as Experiment 1 was used here to tidy the data for easy use later: reducing the number of variables displayed to only relevant ones, remove duplicates (3 duplicate IDs found), and changing variables from Character to Factor/Numeric.

# Select only relevant data for analyses
exptwo <- exptwo %>% 
  select(Consent, Finished, Gender, Age, NC_1:condition)

# Remove Duplicates
exptwo$Prolific_PID[duplicated(exptwo$Prolific_PID)]
## [1] "5b8800b14fe6450001f07d71" "5ae8d7b596845f0001c7948d"
## [3] "5e68d2b5576136000ba058a2"
exptwo <- exptwo[!duplicated(exptwo$Prolific_PID), ]

# Change variables from Chr to Numeric/Factor accordingly
exptwo <- exptwo %>% 
  mutate(
    across(c(Age:Development_sci_know_6), as.numeric),
    across(c(Consent, Gender, condition), as.factor)
         )

The authors reported 412 consented to take part. Let’s try to count this:

# Count number of Consent 
exptwo %>% 
  count(Consent)
## # A tibble: 1 x 2
##   Consent     n
##   <fct>   <int>
## 1 1         409

Oh no! This doesn’t match up with the reported. I wonder if it is because they included the 3 duplicates in their consent count? This is really odd.. Let me try to reverse the process to BEFORE I did any data tidying (using the original data file)…

exptwooriginal <- "Study 7 data.csv" %>% 
  read_csv()

exptwoconsent <- exptwooriginal %>%  
  filter(Consent == "1") 

count(exptwoconsent)
## # A tibble: 1 x 1
##       n
##   <int>
## 1   412

There we go! N = 412. So it seems that when reporting Number of total consented participants = 412, the authors included their 3 duplicates.

Lesson on troubleshooting: Always trace your steps back one at a time, and run the line of code line-by-line when things aren’t quite adding up…

Now, applying pre-registered exclusion criteria (same as Experiment 1).

exptwofinaldata <- exptwo %>% 
  filter(
    Consent == "1", #filter to include those who consented
    Finished == "1", #filter to include those who Finished
    Serious_check == "1", #filter to include those who answered they passed Serious Check
    recall_score >= "4", #filter to include recalls core 4 and above only
    )

The authors did not actually remove duplicates in their codebook, as the duplicates did not complete their survey (their data was thus incomplete and wouldn’t have been in the data analysis once the filter() above was applied). However, I thought it might be good to check for duplicates anyway as it felt like a good practice to do so before data analysis.

count() allows us to check how many final participants are included in our data set, and cross reference with the data reported in our paper (N = 400).

exptwofinaldata %>% 
  count (Consent)
## # A tibble: 1 x 2
##   Consent     n
##   <fct>   <int>
## 1 1         400

Based on the code for Experiment 1, I am assuming that 1 and 2 corresponds to Male and Female, respectively. Using the functions described previously, I’ll rename accordingly and count() to check whether the count corresponds with our paper’s reporting of 150 identifying as male, 248 females, and 2 neither of those categories.

exptwofinaldata <- exptwofinaldata %>% 
  mutate(
    Gender = case_when(
      Gender == 1 ~ "Male",
      Gender == 2 ~ "Female",
      Gender == 3 ~ "Neither"
    )
  )

exptwofinaldata %>% 
  count(Gender)
## # A tibble: 3 x 2
##   Gender      n
##   <chr>   <int>
## 1 Female    248
## 2 Male      150
## 3 Neither     2

Finally, to calculate the mean, SD and range for age, we use summarise() and display in a table using the gt package.

Interestingly, for Experiment 2, the authors reported their descriptive statistics to 1 decimal place, in contrast to 2 decimal places used for Experiment 1… This made me wonder why, and whether it would have been better to be more consistent throughout!

This should correspond to - Mean: 33.5, SD: 12, Range: 18-73.

exptwofinaldata %>% 
  summarise(
    Mean = mean(Age),
    SD = sd(Age),
    Min = min(Age),
    Max = max(Age)
) %>% 
  gt() %>% 
  fmt_number(
    columns = vars(Mean, SD),    
    decimals = 1
  ) %>% 
  tab_header (
    title = "Descriptive Statistics for Age in Experiment 2"
  )
Descriptive Statistics for Age in Experiment 2
Mean SD Min Max
33.5 12.0 18 73

Figure 3 - Experiment 2

Reproducibility goal:

Figure 3 plots the six variable measured and is similar to Figure 1 (Nutritional Confusion, Nutritional Backlash, Mistrust of Expertise, Confidence in Scientific Community, Certainty of Knowledge, & Development of Knowledge). As such, we can use the previous function with a few tweaks!

Setting up the data is similar to Experiment 1: split participants’ conditions into 4 variables using separate(), changing the new variables ‘Format’ and ‘Conflict’ from Character to Factor, and renaming Conflict/Consistent to Conf./Non-Conf.

exptwofinaldata <- separate(
  data = exptwofinaldata, 
  col = condition, 
  into = c("block", "number", "format", "conflict")
  )

exptwofinaldata <- exptwofinaldata %>% 
  mutate(
    across(c(format:conflict), as.factor)
         )

levels(exptwofinaldata$conflict)[levels(exptwofinaldata$conflict)=="Conflict"] <- "Conf."
levels(exptwofinaldata$conflict)[levels(exptwofinaldata$conflict)=="Consistent"] <- "Non-Conf."

Since there were multiple questions asked for each variable, we need to calculate the participants’ average scores. I am going to do this by using mutate() to create a new variable to represent the average scores of each (i.e., confusion, backlash, mistrust, confidence, certainty, development), then I’ll specify the formula to calculate the average of each dependent variable.

# Calculate participants' average scores as a new variable for each of the 6 variables studied
exptwofinaldata <- exptwofinaldata %>% 
  mutate (
    confusion = ((NC_1 + NC_2 + NC_3 + NC_4 + NC_5 + NC_6)/6),
    backlash = ((NBS_1 + NBS_2 + NBS_3 + NBS_4 + NBS_5 + NBS_6)/6),
    mistrust = ((Mistrust_expertise_1 + Mistrust_expertise_2 + Mistrust_expertise_3)/3),
    confidence = exptwofinaldata$GSS,
    certainty = ((Certainty_sci_know_1 + Certainty_sci_know_2 + Certainty_sci_know_3 + 
                    Certainty_sci_know_4 + Certainty_sci_know_5 + Certainty_sci_know_6)/6),
    development = ((Development_sci_know_1 + Development_sci_know_2 + Development_sci_know_3 + 
                      Development_sci_know_4 + Development_sci_know_5 + Development_sci_know_6)/6)
    )

Creating the Function

Since the plots are highly similar to the ones from Experiment 1, we copied over our function and made some slight modifications. Since most plots have a y-axis upper limit of 5, we state this argument within our function() and later for our “confidence” plot we can override this by with the argument lim_2 = 3.

# Figure 3 function
figure.3.fun <- function(y_var, plot_title, y_title, lim_2 = 5) {
  
  ggplot(exptwofinaldata, aes(x= conflict, y = y_var, fill = conflict)) +
  geom_violin() +
  facet_wrap(vars(format), strip.position = "bottom") +
  stat_summary(fun.data = "mean_cl_normal",geom = "crossbar", fill = "white", alpha = .7) +
  geom_beeswarm(cex = 0.2) +
  ggtitle(label = plot_title) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 8), axis.title.y = element_text(size = 6)) +
  easy_center_title() +
  easy_remove_legend() +
  easy_remove_x_axis(what = c("title")) +
  scale_y_continuous(name = y_title, limits= c(1,lim_2)) +
  scale_fill_manual(values = c("slategray2", "lightpink1")) 
  
} 

# Use function to make plots
confusion <- figure.3.fun(
  y_var = exptwofinaldata$confusion, 
  plot_title = "Nutritional Confusion", 
  y_title = "Nutritional Confusion"
  )

backlash <- figure.3.fun(
  y_var = exptwofinaldata$backlash, 
  plot_title = "Nutritional Backlash", 
  y_title = "Nutritional Backlash"
  )

mistrust <- figure.3.fun(
  y_var = exptwofinaldata$mistrust, 
  plot_title = "Mistrust of Expertise", 
  y_title = "Mistrust of Expertise"
  )

confidence <- figure.3.fun(
  y_var = exptwofinaldata$confidence, 
  plot_title = "Confidence in Scientific Community", 
  y_title = "Confidence in Scientific Community", 
  lim_2 = 3
  )

certainty <- figure.3.fun(
  y_var = exptwofinaldata$certainty, 
  plot_title = "Certainty of Knowledge", 
  y_title = "Certainty of Knowledge"
  )

development <- figure.3.fun(
  y_var = exptwofinaldata$development, 
  plot_title = "Development of Knowledge", 
  y_title = "Development of Knowledge"
  )

# Combine plots
combinedplots2 <- confusion + backlash + mistrust + confidence + certainty + 
  development + plot_layout(ncol = 2)

print(combinedplots2)

Hooray! We were able to reproduce ALL the descriptive statistics and figures!!

Part 3: Exploratory Analyses

Q1: Is there a relationship between participants’ recall of headlines and their feelings of mistrust in science?

Enzenbach et al. (2019) found participants’ interest in the research topic can result in bias. In our paper, the headline recall test was used as an attention check. However, it is also possible that participants’ implicit biases play a role in affecting their attention to scientific headlines - and thus affecting their recall score. Those who have more mistrust in health and nutrition advice may pay less attention to related headlines (and vice versa). Therefore, I am interested to see whether there is a relationship between participants’ recall scores and their self-reported feelings of mistrust in science.

Step 1: Descriptive Statistics

To start, I am going to count the number for each recall_score variable using count() to make sure there are enough participants in each group I am interested in. Since the study excluded participants that scored 3 and below, data displayed is from 4-7. From this, we can see that the majority of participants scored 7, and only 7 participants scored a recall score of 4.

exptwofinaldata %>% 
  count(recall_score)
## # A tibble: 4 x 2
##   recall_score     n
##   <chr>        <int>
## 1 4                7
## 2 5               14
## 3 6               61
## 4 7              318

For the purpose of this analysis, I am going to exclude participants who scored 6 from subsequent analyses to allow maximal differentiation between low vs. high recall groups. I am then going to mutate() a new variable called recallers. Using the function case_when(), I specify recall_score 4 - 5 as “Low” and recall_score 7 as “High”. I also use na.omit() which is part of the stats package, and allows me to deal with NAs. Specifically, since those that have recall_score of 6 will result in NA under the new variable recallers, I’ll be able to exclude them from my analysis.

The stats package automatically loads in an RStudio session, hence we do not need to load it at the start.

exptwofinaldata <- exptwofinaldata %>% 
  mutate(
    recallers = case_when(
      recall_score <= 5 ~ "Low",
      recall_score == 7 ~ "High"
    )
  ) %>% 
  na.omit() #removes NAs

exptwofinaldata %>% 
  count(recallers)
## # A tibble: 2 x 2
##   recallers     n
##   <chr>     <int>
## 1 High        318
## 2 Low          21

The vertical bar | functions similarly as the pipe function, so we can filter for recallers “Low” and “High”. I first group by the variable ‘recallers’. I can then find the Mean, SE, and SD for the dependent variable “mistrust of expertise” grouped by recallers using summarise(). Using gt(), we display the descriptive statistics in a table.

byrecall <- exptwofinaldata %>% 
  filter(
    recallers == "Low" | recallers == "High"
    ) %>% 
  group_by(recallers) %>% 
  summarise(mean = mean(mistrust),
            sd = sd(mistrust),
            n = n(),
            se = sd/sqrt(n))  

byrecall %>% 
  gt() %>% 
  fmt_number(
    columns = vars(mean, sd, se),
    decimals = 2
  ) %>% 
  tab_header (
    title = "Perceived mistrust of expertise in high vs. low recall score groups"
  )
Perceived mistrust of expertise in high vs. low recall score groups
recallers mean sd n se
High 1.92 0.64 318 0.04
Low 2.37 0.63 21 0.14

It looks like low recallers (those who have low recall scores/poor memory of the headlines) have on average higher mistrust of expertise compared to high recallers. This is in line with my hypothesis that participants who score higher on recall (perhaps as a result of paying more attention towards headlines) also have on average more trust in science and expert opinion.

Step 2: Data visualisation

Using ggplot and its functions geom_jitter() (adds a small amount of variation to the location of each point, while showing all individual data points) and geom_boxplot, we can visualise the data nicely.

I am using geom_jitter() instead of geom_point() as the geom_jitter() function creates the scatterplot with random noise added to minimise overlapping and visualise the data a lot better. If I used geom_point() it would have just plotted each data point one on top of another and wouldn’t be as informative.

A few aesthetic changes previously explained in the verification report are also used here. Finally, using ggeasy’s functions we can center title with easy_center_title(), add legend title with easy_add_legend_title(), and remove the x-axis title using easy_remove_x_axis() then specifying the argument what = c("title").

recallplot <- ggplot(
  data = exptwofinaldata,
  aes(
    x = conflict,
    y = mistrust,
    fill = recallers
  ) 
) +
  geom_jitter(
    alpha = .4,
    aes(
      color = recallers
      )
  ) +
  geom_boxplot() +
  ggtitle(
    label = "Mistrust of Expertise in Low vs. High Recall Score Groups"
  ) +
  facet_wrap(
    vars(format),
    strip.position = "bottom" 
  ) +
  easy_center_title() +
  easy_remove_x_axis(what = c("title")) +
  easy_add_legend_title("Recall Score") +
  scale_y_continuous(
    name = "Mistrust of expertise"
  )

print(recallplot)

Visually inspecting these plots, it looks like those who score low on recall also have on average higher mistrust, when compared to participants who scored high on recall.

I decided to use geom_boxplot() as it clearly shows the distribution of mistrust of expertise ratings in high vs. low recall score groups. I can visually inspect the shape of the distribution, central value, and variability.

Step 3: Statistics - t-test

We can run a t-test to find out if there is a statistically significant difference. Since high/low recallers are a categorical variable and mistrust of expertise is a continuous variable - I am going to use a t-test. To do this, I am going to create new objects based on participants’ recall scores (i.e., high vs. low recallers). Then, using t.test() in the stats package we test whether the mistrust between low vs. high recalls are equal.

lowrecall <- exptwofinaldata %>% 
  filter (
    recallers == "Low"
    )

highrecall <- exptwofinaldata %>% 
  filter (
    recallers == "High"
  ) 

t.test(lowrecall$mistrust, highrecall$mistrust)
## 
##  Welch Two Sample t-test
## 
## data:  lowrecall$mistrust and highrecall$mistrust
## t = 3.1594, df = 22.797, p-value = 0.004415
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.15521 0.74476
## sample estimates:
## mean of x mean of y 
##  2.365079  1.915094

Participants who score low scores on the memory recall test have higher mistrust of expertise compared to participants who score high scores, t(22.797) = 3.1594, p = 0.004415.

Step 4: Calculate effect size

I am going to use cohens_d() from the rstatix package to calculate effect size of the above finding. var.equal = TRUE treats the two variances as being equal.

exptwofinaldata %>% 
  cohens_d(mistrust ~ recallers, var.equal = TRUE)
## # A tibble: 1 x 7
##   .y.      group1 group2 effsize    n1    n2 magnitude
## * <chr>    <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 mistrust High   Low     -0.704   318    21 moderate

This effect size of is moderate (d > 0.5), indicating that the means between high vs. low recall participants differ by approximately 0.7 SDs.

Therefore, we can say that participants who have high headline recall scores show on average greater trust in science and expert opinion. Conversly, participants who have low headline recall scores show on average greater mistrust and skepticism. This effect is at least moderate in size.

Q2: Is nutritional confusion correlated with backlash?

Based on previous studies (Clark et al., 2019; Vijaykumar et al., 2021), confusing headlines can result backlash (i.e. negative beliefs about recommendations and research). This analysis thus aims to look at whether there is a correlation between these two variables. We expect that more confusing headlines will be correlated with higher negative beliefs about research (i.e. more backlash).

Step 1: Data Visualisation

cor() from the stats package allows us to calculate the correlation coefficient between variables. Then, using corrplot() from the corrplot package we can create a correlation matrix to visually show correlation between dependent variables.

I decided to visualise the data before jumping straight into plotting, as the correlation plot allows us to see from a birds-eye view where other correlations might exist (and might helpful to further analyses!

corr_data <- exptwofinaldata %>% 
  select(
    confusion,
    backlash,
    mistrust,
    confidence,
    certainty,
    development
  ) %>% 
  cor() 

corrplot(
  corr_data,
  method = "pie", #shows correlation as pie chart
  type = "lower" #displays only lower of triangular matrix to remove duplicates
)

From the correlation plot, it seems that there is a positive correlation between nutrition confusion and backlash.

We can also see another large correlation is between mistrust and confidence. While I won’t be looking at this specifically, being able to see it all laid out in the correlation plot is quite insightful!

Next, I am going to plot a graph with nutrition confusion vs. backlash using geom_jitter() which allows me to see whether there is a pattern in the relationship between both variables. Then, using geom_smooth() to create a smooth line, and specifying the argument for method="lm" allows us to plot the line of best fit using the linear regression model y = mx + c.

ggplot(
  exptwofinaldata, 
  aes(
    x = confusion, 
    y = backlash
  )
) +
geom_jitter() +
geom_smooth(
  method = "lm"
) +
scale_x_continuous(
  name = "Nutritional Confusion"
) +
scale_y_continuous(
  name = "Nutritional Backlash"
) +
ggtitle(
  label = "Nutritional backlash as a function of nutritional confusion"
)

As expected, there is a positive linear correlation between confusion and backlash. Next, let’s look at statistical analysis.

Step 2: Statistics - Correlation Test

I am using cor.test() from stats package to look at the Pearson’s correlation between confusion and backlash.

cor.test (exptwofinaldata$confusion, exptwofinaldata$backlash)
## 
##  Pearson's product-moment correlation
## 
## data:  exptwofinaldata$confusion and exptwofinaldata$backlash
## t = 10.606, df = 337, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4158959 0.5760755
## sample estimates:
##       cor 
## 0.5002534

Correlation is 0.466 - a strong positive relationship.

There is a significant strong positive relationship between nutrition confusion and backlash, t(337) = 10.606, p < 0.05.

Q3: Does age influence participants’ ratings of scientific advancement?

Related to one of my initial reactions, I wondered whether participants’ age influences their ratings. This is supported by Klaczynski and Robinson (2000), which found that reasoning biases like personal theories and epistemological beliefs differs with age. Haigh and Birch used a 3-point scale (know less, know the same, know more) to assess how participants perceive scientific advancement. I am interested to see whether age influences participants’ ratings of scientific advancement (i.e. do older adults rate scientific advancement differently to younger adults in this study)? If so, this may be a confounding variable that could influence the results, and further studies should account for this variable when interpreting results.

Step 1: Descriptive Statistics

First, I am going to group the age variables in categories to allow easier interpretation. From the verification report, Experiment 1’s age range is 18-69. This allows us to group as below:

  • Gen Z (1997 – 2012): 9 – 24
  • Millennials (1981 – 1996): 25 – 40
  • Gen X (1965 – 1980): 41 – 56
  • Boomers (1946 – 1964): 57 - 75

I am thus going to use the mutate() and case_when to create a new variable called generation where I will specify using %in% the age ranges and generation names. The %in% operator allows us to identify an element (e.g. numbers) within the data frame.

Then, I am converting generations to a factor variable using baseR function factor(), and ordering it chronologically (R defaults to alphabetical) using the levels argument. Then, I am going to create a frequency table with both advancement and generations. print() function from baseR print/display the output of our table.

exponefinaldata <- exponefinaldata %>% 
  mutate(
      generations = case_when(
        Age %in% 18:24 ~ "GenZ",
        Age %in% 25:40 ~ "Millennials",
        Age %in% 41:56 ~ "GenX",
        Age %in% 57:69 ~ "Boomers"
    )
  ) 

exponefinaldata$generations <- factor(
  exponefinaldata$generations, 
  levels = c("GenZ", "Millennials", "GenX", "Boomers")
  ) #change to factor variable and specify levels chronologically

generations_table <- table(exponefinaldata$advancement, exponefinaldata$generations)

print(generations_table)
##       
##        GenZ Millennials GenX Boomers
##   less   28          33   19       5
##   same   41          69   32      17
##   more   14          25    9       2

Unfortunately this isn’t too informative, as the numbers of each group are different. Displaying frequency table as % of each generation would be be more insightful!

We can use prop.table() from baseR for this, specifying we want the proportion calculated across columns (not rows) by specifying the argument margin = 2. Since I want it displayed as a % (out of 100), I then use *100 to multiply the proportion data to convert it into percentages.

prop.table(
  generations_table,
  margin = 2 #indicates we want to calculate proportion across columns of the table
)*100
##       
##             GenZ Millennials      GenX   Boomers
##   less 33.734940   25.984252 31.666667 20.833333
##   same 49.397590   54.330709 53.333333 70.833333
##   more 16.867470   19.685039 15.000000  8.333333

From the table above, it looks like significantly more Boomers say they know the same before/after reading the headlines (70%), contrasted with approximately 50% of this response from all other generations. This could mean they are more resistant to conflicting headlines, and that manipulation of conflict was not as strong on Boomers as for the other age groups!

Step 2: Data Visualisation

I am creating a frequency plot (in percentage, after my lesson above from the frequency table) using ggplot() and geom_bar(), I’m going to facet_wrap() by generation and display it in one row for easy comparison using the argument nrow. easy_center_title() and easy_remove_legend() from ggeasy once again allows me to easily center my plot title and remove legend. We need to specify the argument group = 1 in geom_bar() so that the bar plot considers the three advancement scores as 1 group for each generation. ..prop.. allows us to display y as proportions, and the argument labels=scales::percent changes this to % on the y-axis.

ggplot(
  exponefinaldata, aes(advancement)) +
  geom_bar(aes(y=..prop.., group=1, fill=generations), position="dodge") + 
  facet_wrap(vars(generations),strip.position = "bottom", nrow = 1) +
  scale_y_continuous(labels=scales::percent, name = "% of ratings") + 
  scale_x_discrete(name = "Perceived Scientific Advancement", labels = c("Less", "Same", "More")) +
  ggtitle(label = "Frequency plot for Perceived Scientific Advancement") +
  easy_center_title() +
  easy_remove_legend()

This is in line with what we observed from the frequency table - Boomers show a higher percentage for “know the same” compared to other age groups. To determine whether this observed pattern is significantly different, we need to run a Chi-square test.

Step 3: Statistics - Chi-square test

A chi-square test is selected as the two variables are of a categorial nature.

chisq.test(generations_table) 
## 
##  Pearson's Chi-squared test
## 
## data:  generations_table
## X-squared = 5.0736, df = 6, p-value = 0.5344

We fail to reject the null hypothesis that both age and perceived advancement variables are independent of each other, X^2(6, 294) = 5.0736, p = 0.5344.

We have insufficient evidence to conclude that a participant’s generation and their perceived scientific advance are related.

Part 4: Recommendations

The authors did a LOT well, which helped us reproduce their descriptive statistics and plots! For example, they included their R script (codebook), and even had checkpoint package that would run their code using the version of each package they used so we could easily run their code and reproduce the same results (something that would be helpful in years to come, when packages are updated or become retired). That said, a few recommendations on what they can do to improve reproducibility:

  1. Improve reproducibility workflow by adding a Read Me file and using human readable names/clear labels.
  • The authors included two CSV files in their OSF - titled Experiment 7.csv and Experiment 8.csv. Yet, there was no mention to Experiments 7/8 anywhere - only Experiments 1 and 2. This had our group puzzled at first, wondering which was which! We had to look in the CSV files to check the variables measured (as the authors had a fancy line of code that downloaded it from a link in their code book, whereas we wanted to use the provided CSV files to load the data) - before we were able to load our data on RStudio. I felt like this could have been improved by including a Read Me file in the OSF folder, where the authors could have explained what was what (it seems like they ran a series of experiments and these two might have been 7 and 8), and how a third party looking to reproduce their results could get started.
  • Also, dependent variables “condition” and “recall scores” were not human readable. At first glance, as we were completely new to the data, we had no idea what we were looking at when we saw SC0, FL_12_DO, and FL_10_DO. We had seen that the author renames these in their codebook, and that helped us to figure it out eventually. However, as a team, we tried not to rely on the authors’ codebook! It would have been better if the authors had used score/condition instead (or any other human readable friendly names). To improve, the authors can choose human readable names that on first glance makes it clear to an external party/reader what the variable is a measure of!
  • Resources:
  1. Improve reproducibility by implementing rubber duck technique.
  • Firstly, the authors installed 17 packages, but next to these packages (and in their codebook) it was not clear what these packages will be used for… This created certain problems for our group. Namely, we tried not to copy the authors’ code and work things out for ourselves. But when we got stuck, we would have a quick look at their code and try the same functions used (to no success, since the function would be from a package we did not install)! To follow the KISS principle, we tried to limit our use of unnecessary packages so this posed an issue for our group. It would be better for the author to help others reproduce their results by explaining why they used each package.
  • The authors also had very simple commenting on their chunks of code (e.g., ####pirate plots / descriptives #### was followed by all the codes for plots/descriptives). Their code was thus quite difficult to follow along, especially for us newbies to R with no knowledge on most of the packages/codes they used. To implement rubber ducking, the authors could try to comment on their code as if they were explaining to a 3rd party that had very basic knowledge on coding.
  • Resources:
  1. Tidy the data by implementing (consistent) data wrangling.
  • The authors’ data was quite untidy and overwhelming. To give you an idea, Experiment 1’s data had 340 columns. This made glimpsing the data using glimpse() time consuming, as there were too many variables to look at. An option is the authors could provide a tidied data file for BOTH experiments, on top of the raw data files. We tried to power through but eventually realised that taking time to tidy the data was worthwhile. Once we tidied our data using select() to only display the relevant variables, we were able to condense the 340 variables to less than 20!! The authors could easily implement something similar, with the help of tidyr package which includes dplyr that we found had really helpful functions for tidying our data!

In our presentation, we mentioned tidying our data makes it much easier to work with and allows us to use the glimpse function to check our data without having to sift through all the irrelevant variables.

  • The authors also had a few inconsistencies in their paper. They reported to 2 decimal places in descriptive statistics for Experiment 1, yet reported to 1 decimal place for Experiment 2. They included a ‘tidied’ version of their data file in CSV for Experiment 2, but not Experiment 1. These were minor issues and didn’t hinder our process too much, but it did have us a little confused when we encountered them. Ensuring a consistent data wrangling process is implemented throughout would be more helpful and make reproducing their results a smoother journey! If multiple people are working on the data, ensuring communication amongst team members on which functions are used for what purpose might also help with this!
  • Resources:

References

  1. Clark, D., Nagler, R., & Niederdeppe, J. (2019). Confusion and nutritional backlash from news media exposure to contradictory information about carbohydrates and dietary fats. Public Health Nutrition, 22(18), 3336-3348. https://doi.org/10.1017/S1368980019002866
  2. Enzenbach, C., Wicklein, B., Wirkner, K., & Loeffler, M. (2019). Evaluating selection bias in a population-based cohort study with low baseline participation: The LIFE-Adult-Study. BMC medical research methodology, 19(1), 135. https://doi.org/10.1186/s12874-019-0779-8
  3. Haigh, M., & Birch, H. A. (2021). When ‘scientists say’ coffee is good for you one day and bad for you the next: Do generic attributions to ‘scientists’ and ‘experts’ amplify perceived conflict? Collabra: Psychology, 7(1), 23447. https://doi.org/10.1525/collabra.23447
  4. Klaczynski, P. A., & Robinson, B. (2000). Personal theories, intellectual ability, and epistemological beliefs: Adult age differences in everyday reasoning biases. Psychology and Aging, 15(3), 400–416. https://doi.org/10.1037/0882-7974.15.3.400
  5. Vijaykumar, S., McNeill, A., & Simpson, J. (2021). Associations between conflicting nutrition information, nutrition confusion and backlash among consumers in the UK. Public Health Nutrition, 24(5), 914-923. https://doi.org/10.1017/S1368980021000124