1 Setting up

First, I load all libraries that we need for the analysis. The pacman package just makes it easier to load packages. Note that the final version of this page doesn’t evaluate the code chunk below to avoid installing packages on your machine. If you’re fine with them being installed, set eval = TRUE.

If you run the entire project with the renv private library, installing packages should’ve happened with the renv::restore call.

if (!requireNamespace("pacman"))
  install.packages("pacman")

library(pacman)

# load packages
p_load(
  tidyverse,
  lubridate,
  here,
  MBESS,
  ggridges,
  GGally,
  ggalt,
  cowplot,
  brms,
  ggbeeswarm,
  extrafont,
  kableExtra,
  osfr
)

# set seed
set.seed(42)

# set theme
theme_set(theme_cowplot())

Below custom functions that I use throughout the project.

dens_with_points <- 
  function(
    data,
    variable
  ) {
    p <- 
      ggplot(data, aes_string(x = variable, y = 0)) +
      geom_density_ridges(
        jittered_points = TRUE, 
        position = "raincloud", 
        fill = "darkslateblue", 
        point_color = "darkslateblue", 
        color = "darkslateblue", 
        alpha = 0.5
      ) +
      theme_cowplot() +
      theme(
        axis.line=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        legend.position="none",
        panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank()
      )
    
    return(p)
  }

word2num <- function(word){
  
  # added to deal with NA
  if (is.na(word)){
    return(NA)
  }
  else {
    wsplit <- strsplit(tolower(word)," ")[[1]]
    one_digits <- list(zero=0, one=1, two=2, three=3, four=4, five=5,
                       six=6, seven=7, eight=8, nine=9)
    teens <- list(eleven=11, twelve=12, thirteen=13, fourteen=14, fifteen=15,
                  sixteen=16, seventeen=17, eighteen=18, nineteen=19)
    ten_digits <- list(ten=10, twenty=20, thirty=30, forty=40, fifty=50,
                       sixty=60, seventy=70, eighty=80, ninety=90)
    doubles <- c(teens,ten_digits)
    out <- 0
    i <- 1
    while(i <= length(wsplit)){
      j <- 1
      if(i==1 && wsplit[i]=="hundred")
        temp <- 100
      else if(i==1 && wsplit[i]=="thousand")
        temp <- 1000
      else if(wsplit[i] %in% names(one_digits))
        temp <- as.numeric(one_digits[wsplit[i]])
      else if(wsplit[i] %in% names(teens))
        temp <- as.numeric(teens[wsplit[i]])
      else if(wsplit[i] %in% names(ten_digits))
        temp <- (as.numeric(ten_digits[wsplit[i]]))
      if(i < length(wsplit) && wsplit[i+1]=="hundred"){
        if(i>1 && wsplit[i-1] %in% c("hundred","thousand"))
          out <- out + 100*temp
        else
          out <- 100*(out + temp)
          j <- 2
      }
      else if(i < length(wsplit) && wsplit[i+1]=="thousand"){
        if(i>1 && wsplit[i-1] %in% c("hundred","thousand"))
          out <- out + 1000*temp
        else
          out <- 1000*(out + temp)
          j <- 2
      }
      else if(i < length(wsplit) && wsplit[i+1] %in% names(doubles)){
        temp <- temp*100
        out <- out + temp
      }
      else{
        out <- out + temp
      }
      i <- i + j
    }
    return(out)
  }
}

describe <- function(
  dat,
  variable,
  trait = FALSE
){
  # if variable is not repeated-measures, take only one measure per participant
  if (trait == TRUE){
    dat <- 
      dat %>%
      group_by(id) %>% 
      slice(1) %>% 
      ungroup()
  }
  
  # then get descriptives
  descriptives <-
    dat %>%
    filter(!is.na(UQ(sym(variable)))) %>% # remove missing values
    summarise(
      across(
        !! variable,
        list(
          mean = mean,
          sd = sd,
          median = median,
          min = min,
          max = max,
          cilow = ~Rmisc::CI(.x)[[3]], # lower CI
          cihigh = ~Rmisc::CI(.x)[[1]] # upper CI
        )
      )
    )

  descriptives <-
    descriptives %>%

    # only keep measure
    rename_all(
      ~ str_remove(
        .,
        paste0(variable, "_")
      )
    ) %>%
    mutate(
      variable = variable,
      range = max - min
    ) %>%
    relocate(variable) %>%
    relocate(
      range,
      .after = max
    )
  
  return(descriptives)
}

# raincloud plot function from https://github.com/RainCloudPlots/RainCloudPlots/blob/master/tutorial_R/R_rainclouds.R
# Defining the geom_flat_violin function ----
# Note: the below code modifies the
# existing github page by removing a parenthesis in line 50

"%||%" <- function(a, b) {
  if (!is.null(a)) a else b
}

geom_flat_violin <- function(mapping = NULL, data = NULL, stat = "ydensity",
                             position = "dodge", trim = TRUE, scale = "area",
                             show.legend = NA, inherit.aes = TRUE, ...) {
  layer(
    data = data,
    mapping = mapping,
    stat = stat,
    geom = GeomFlatViolin,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      trim = trim,
      scale = scale,
      ...
    )
  )
}

#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
GeomFlatViolin <-
  ggproto("GeomFlatViolin", Geom,
    setup_data = function(data, params) {
      data$width <- data$width %||%
        params$width %||% (resolution(data$x, FALSE) * 0.9)

      # ymin, ymax, xmin, and xmax define the bounding rectangle for each group
      data %>%
        group_by(group) %>%
        mutate(
          ymin = min(y),
          ymax = max(y),
          xmin = x,
          xmax = x + width / 2
        )
    },

    draw_group = function(data, panel_scales, coord) {
      # Find the points for the line to go all the way around
      data <- transform(data,
        xminv = x,
        xmaxv = x + violinwidth * (xmax - x)
      )

      # Make sure it's sorted properly to draw the outline
      newdata <- rbind(
        plyr::arrange(transform(data, x = xminv), y),
        plyr::arrange(transform(data, x = xmaxv), -y)
      )

      # Close the polygon: set first and last point the same
      # Needed for coord_polar and such
      newdata <- rbind(newdata, newdata[1, ])

      ggplot2:::ggname("geom_flat_violin", GeomPolygon$draw_panel(newdata, panel_scales, coord))
    },

    draw_key = draw_key_polygon,

    default_aes = aes(
      weight = 1, colour = "grey20", fill = "white", size = 0.5,
      alpha = NA, linetype = "solid"
    ),

    required_aes = c("x", "y")
  )

single_cloud <- 
  function(
    raw_data,
    summary_data,
    variable,
    color,
    title,
    trait = FALSE
  ){
    
    # take only one row per person if it's a trait variable
    if (trait == TRUE){
      raw_data <-
        raw_data %>% 
        group_by(id) %>% 
        slice(1) %>% 
        ungroup()
    }
    
    # the plot
    p <- 
      ggplot(
        raw_data %>%
          mutate(Density = 1),
        aes(
          x = Density,
          y = get(variable)
        )
      ) +
      geom_flat_violin( # the "cloud"
        position = position_nudge(x = .2, y = 0),
        adjust = 2,
        color = NA,
        fill = color,
        alpha = 0.5
      ) +
      geom_point( # the "rain"
        position = position_jitter(width = .15),
        size = 1,
        color = color,
        alpha = 0.5
      ) +
      geom_point( # the mean from the summary stats
        data = summary_data %>%
          filter(variable == !! variable) %>%
          mutate(Density = 1),
        aes(
          x = Density + 0.175,
          y = mean
        ),
        color = color,
        size = 2.5
      ) +
      geom_errorbar( # error bars
        data = summary_data %>%
          filter(variable == !! variable) %>%
          mutate(Density = 1),
        aes(
          x = Density + 0.175,
          y = mean,
          ymin = cilow,
          ymax = cihigh
        ),
        width = 0,
        size = 0.8,
        color = color
      ) +
      ylab(title) +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_blank(),
        axis.line = element_blank()
      ) +
      guides(
        color = FALSE,
        fill = FALSE
      ) +
      coord_flip()
    
    return(p)
  }

lm_function <- 
  function(
    data, 
    mapping, 
    ...
    ){
  p <- 
    ggplot(
      data = data, 
      mapping = mapping
      ) + 
    geom_point(
      color = "#56B4E9",
      alpha = 0.5
    ) + 
    geom_smooth(
      method=lm, 
      fill="#0072B2", 
      color="#0072B2", 
      ...)
  p
}

dens_function <-
  function(
    data,
    mapping,
    ...
  ){
    p <- 
      ggplot(
        data = data,
        mapping = mapping
      ) +
      geom_density(fill = "#009E73", color = NA, alpha = 0.5)
  }

model_diagnostics <- 
  function(
    model
  ){
  plot_grid(
    pp_check(
      model,
      type = "dens_overlay",
      nsamples = 100
    ),
    pp_check(
      model,
      type = "loo_pit_qq",
      nsamples = 100
    ),
    pp_check(
      model,
      type = "loo_pit_overlay",
      nsamples = 100
    ),
    pp_check(
      model,
      type = "stat",
      stat = "median",
      nsamples = 100
    ),
    pp_check(
      model,
      type = "stat",
      stat = "mean",
      nsamples = 100
    ),
    labels = c("Density overlay", "LOO-PIT QQ", "LOO-PIT Uniform", "Predicted medians", "Predicted means"),
    ncol = 2,
    label_size = 8,
    hjust = 0,
    vjust = 0,
    label_x = 0,
    label_y = 0.93
  )
  }