Chapter 5 Appendix

5.1 Simulated Data Creation

Below is the code used to create the simulated data that will be used to create the example SGBA seen throughout the rest of this example analysis.

# load libraries ----------------------------------------------------------
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# set seed for predictability
set.seed(42, kind = "Mersenne-Twister")

# write simulation functions ----------------------------------------------

# function for simulating bimodal gender variables
binom_bounded <- function(
    n, prop, mean1, sd1, mean2, sd2, lim_low = -Inf, lim_up = Inf, round
    ){
  # error checking
  if(prop > 1 || prop< 0){stop("binomial proportion should be between 0 and 1")}
  if(n < 1){stop("n must be greater than or equal to 1")}
  if(sd1 < 0 || sd2 < 0){stop("standard deviations must be non-negative")}
  # create index
  index <- rbinom(n, size = 1, prob = prop)
  # simulate each peak using a bounded `rnorm()`
  binom_sim1 <- index * qnorm(
    runif(
      n, pnorm(lim_low, mean1, sd1), pnorm(lim_up, mean1, sd1)
      ), 
    mean1, sd1)
  binom_sim0 <- (1 - index) * qnorm(
    runif(
      n, pnorm(lim_low, mean2, sd2), pnorm(lim_up, mean2, sd2)
      ), 
    mean2, sd2)
  binom_sim <- binom_sim0 + binom_sim1
  binom_sim <- round(binom_sim, round)
  binom_sim
}


# function for simulating Likert/ordinal outcome variable
ordinal_sim <- function(
    n, mean, sd, lim_low = -Inf, lim_up = Inf
    ){
  ord_sim <- qnorm(
    runif(
      n, pnorm(lim_low, mean, sd), pnorm(lim_up, mean, sd)
      ),
    mean, sd)
  ord_sim <- round(ord_sim, 0)
  ord_sim
}


# simulate SGBA-5 data with n of 30 -------------------------------------

# biological sex: categorical (female, intersex, male)
bio_sex <- factor(
  sample(c('Female', 'Male'), 30, replace=TRUE, prob=c(0.6, 0.4)),
  )

# gender identity: ordered (0 [feminine] to 100 [masculine])
gen_id <- binom_bounded(
  n = 30, prop = 0.6, mean1 = 20, sd1 = 20, mean2 = 85, sd2 = 15, lim_low = 0,
  lim_up = 100, round = 0
)

# gender role: ordered (0 [feminine] to 100 [masculine])
gen_role <- binom_bounded(
  n = 30, prop = 0.6, mean1 = 25, sd1 = 20, mean2 = 80, sd2 = 20, lim_low = 0,
  lim_up = 100, round = 0
)


# simulate example positive outcome data with n of 30 -----------------------

## simulated numerical outcome: --------------------------------------------

# by sex: males (mean = 12, SD = 3), females (mean = 4, sd = 2)
outcome_num_pos_m <- rnorm(n = 16, mean = 12, sd = 3)
outcome_num_pos_f <- rnorm(n = 14, mean = 4, sd = 2)
outcome_num_pos_sex <- c(outcome_num_pos_f, outcome_num_pos_m)

# by g_id: high (mean = 12, SD = 3), low (mean = 4, sd = 2)
outcome_num_pos_gid_high <- rnorm(n = 15, mean = 12, sd = 3)
outcome_num_pos_gid_low <- rnorm(n = 15, mean = 4, sd = 2)
outcome_num_pos_gid <- c(outcome_num_pos_gid_high, outcome_num_pos_gid_low)

# by g_role: high (mean = 12, SD = 3), low (mean = 4, sd = 2)
outcome_num_pos_grol_high <- rnorm(n = 12, mean = 12, sd = 3)
outcome_num_pos_grol_low <- rnorm(n = 18, mean = 7, sd = 2)
outcome_num_pos_grol <- c(outcome_num_pos_grol_high, outcome_num_pos_grol_low)


## simulated Likert outcome: ---------------------------------------------

# by sex males (mean = 2, SD = 1), females (mean = 5, sd = 2)
outcome_ord_pos_m <- ordinal_sim(
  n = 16, mean = 2, sd = 1, lim_low = 1, lim_up = 7
)
outcome_ord_pos_f <- ordinal_sim(
  n = 14, mean = 5, sd = 1, lim_low = 1, lim_up = 7
)
outcome_ord_pos_sex <- c(outcome_ord_pos_f, outcome_ord_pos_m) %>%
  factor(., levels = c(1,2,3,4,5,6,7), ordered = TRUE)


# by g_id: high (mean = 5, SD = 2), low (mean = 1, sd = 2)
outcome_ord_pos_high <- ordinal_sim(
  n = 15, mean = 5, sd = 2, lim_low = 1, lim_up = 7
)
outcome_ord_pos_low <- ordinal_sim(
  n = 15, mean = 1, sd = 2, lim_low = 1, lim_up = 7
)
outcome_ord_pos_gid <- c(outcome_ord_pos_high, outcome_ord_pos_low) %>%
  factor(., levels = c(1,2,3,4,5,6,7), ordered = TRUE)


# by g_role: high (mean = 5, SD = 2), low (mean = 1, sd = 2)
outcome_ord_pos_high <- ordinal_sim(
  n = 12, mean = 5, sd = 2, lim_low = 1, lim_up = 7
)
outcome_ord_pos_low <- ordinal_sim(
  n = 18, mean = 1, sd = 2, lim_low = 1, lim_up = 7
)
outcome_ord_pos_grol <- c(outcome_ord_pos_high, outcome_ord_pos_low) %>%
  factor(., levels = c(1,2,3,4,5,6,7), ordered = TRUE)


## simulated binary outcome: -----------------------------------------------

# by sex males (yes = .2, no = .8), females (yes = .8, no = .2)
outcome_bin_pos_m <- sample(
  c('yes','no'), 16, replace = TRUE, prob = c(.2, .8)
)
outcome_bin_pos_f <- sample(
  c('yes','no'), 14, replace = TRUE, prob = c(.8, .2)
)
outcome_bin_pos_sex <- append(outcome_bin_pos_f, outcome_bin_pos_m) %>%
  factor()


# by g_id: high (yes = .4, no = .6), low (yes = .8, no = .2)
outcome_bin_pos_high <- sample(
  c('yes','no'), 15, replace = TRUE, prob = c(.4, .6)
)
outcome_bin_pos_low <- sample(
  c('yes','no'), 15, replace = TRUE, prob = c(.8, .2)
)
outcome_bin_pos_gid <- append(outcome_bin_pos_high, outcome_bin_pos_low) %>%
  factor()


# by g_rol: high (yes = .3, no = .5), low (yes = .8, no = .2)
outcome_bin_pos_high <- sample(
  c('yes','no'), 12, replace = TRUE, prob = c(.3, .7)
)
outcome_bin_pos_low <- sample(
  c('yes','no'), 18, replace = TRUE, prob = c(.8, .2)
)
outcome_bin_pos_grol <- append(outcome_bin_pos_high, outcome_bin_pos_low) %>%
  factor()


# simulate example negative outcome data with n of 30 -----------------------

# simulated numerical outcome: continuous (mean = 10, SD = 3)
outcome_num_neg <- rnorm(n = 30, mean = 10, sd = 3)

# simulated Likert outcome: 7-point Likert scale (mean = 4, SD = 2)
outcome_ord_neg <- ordinal_sim(
  n = 30, mean = 4, sd = 2, lim_low = 1, lim_up = 7
)

# simulated binary outcome: categorical (yes, no)
outcome_bin_neg <- factor(
  sample(c('yes','no'), 30, replace = TRUE, prob = c(0.67, 0.33))
)


# create example data frame -----------------------------------------------

# combine into dataframe
sim_data <- tibble(
  bio_sex, gen_id, gen_role, #outcome_bin_pos, 
  outcome_num_neg, outcome_ord_neg, outcome_bin_neg
  ) %>% arrange(., bio_sex) %>% 
  cbind(., outcome_num_pos_sex) %>%
  cbind(., outcome_ord_pos_sex) %>%
  cbind(., outcome_bin_pos_sex) %>%
  arrange(., gen_id) %>% 
  cbind(., outcome_num_pos_gid) %>% 
  cbind(., outcome_ord_pos_gid) %>%
  cbind(., outcome_bin_pos_gid) %>%
  arrange(., gen_role) %>% 
  cbind(., outcome_num_pos_grol) %>% 
  cbind(., outcome_ord_pos_grol) %>%
  cbind(., outcome_bin_pos_grol)


# save simulated df
write_csv(sim_data, file = "sim-data.csv")

5.1.1 Session Info for Data Creation

sessioninfo::session_info(pkgs = "loaded")
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.5.2 (2025-10-31)
##  os       macOS Tahoe 26.4
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/Toronto
##  date     2026-04-14
##  pandoc   3.6.3 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
##  quarto   1.8.25 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────
##  package      * version date (UTC) lib source
##  bit            4.6.0   2025-03-06 [1] CRAN (R 4.5.0)
##  bit64          4.6.0-1 2025-01-16 [1] CRAN (R 4.5.0)
##  bookdown       0.46    2025-12-05 [1] CRAN (R 4.5.2)
##  bslib          0.10.0  2026-01-26 [1] CRAN (R 4.5.2)
##  cachem         1.1.0   2024-05-16 [1] CRAN (R 4.5.0)
##  cli            3.6.5   2025-04-23 [1] CRAN (R 4.5.0)
##  crayon         1.5.3   2024-06-20 [1] CRAN (R 4.5.0)
##  digest         0.6.39  2025-11-19 [1] CRAN (R 4.5.2)
##  dplyr        * 1.2.0   2026-02-03 [1] CRAN (R 4.5.2)
##  evaluate       1.0.5   2025-08-27 [1] CRAN (R 4.5.0)
##  farver         2.1.2   2024-05-13 [1] CRAN (R 4.5.0)
##  fastmap        1.2.0   2024-05-15 [1] CRAN (R 4.5.0)
##  forcats      * 1.0.1   2025-09-25 [1] CRAN (R 4.5.0)
##  generics       0.1.4   2025-05-09 [1] CRAN (R 4.5.0)
##  ggplot2      * 4.0.2   2026-02-03 [1] CRAN (R 4.5.2)
##  glue           1.8.0   2024-09-30 [1] CRAN (R 4.5.0)
##  gtable         0.3.6   2024-10-25 [1] CRAN (R 4.5.0)
##  hms            1.1.4   2025-10-17 [1] CRAN (R 4.5.0)
##  htmltools      0.5.9   2025-12-04 [1] CRAN (R 4.5.2)
##  jquerylib      0.1.4   2021-04-26 [1] CRAN (R 4.5.0)
##  jsonlite       2.0.0   2025-03-27 [1] CRAN (R 4.5.0)
##  knitr          1.51    2025-12-20 [1] CRAN (R 4.5.2)
##  lifecycle      1.0.5   2026-01-08 [1] CRAN (R 4.5.2)
##  lubridate    * 1.9.5   2026-02-04 [1] CRAN (R 4.5.2)
##  magrittr       2.0.4   2025-09-12 [1] CRAN (R 4.5.0)
##  otel           0.2.0   2025-08-29 [1] CRAN (R 4.5.0)
##  pillar         1.11.1  2025-09-17 [1] CRAN (R 4.5.0)
##  pkgconfig      2.0.3   2019-09-22 [1] CRAN (R 4.5.0)
##  purrr        * 1.2.1   2026-01-09 [1] CRAN (R 4.5.2)
##  R6             2.6.1   2025-02-15 [1] CRAN (R 4.5.0)
##  RColorBrewer   1.1-3   2022-04-03 [1] CRAN (R 4.5.0)
##  readr        * 2.2.0   2026-02-19 [1] CRAN (R 4.5.2)
##  rlang          1.1.7   2026-01-09 [1] CRAN (R 4.5.2)
##  rmarkdown      2.30    2025-09-28 [1] CRAN (R 4.5.0)
##  rstudioapi     0.18.0  2026-01-16 [1] CRAN (R 4.5.2)
##  S7             0.2.1   2025-11-14 [1] CRAN (R 4.5.2)
##  sass           0.4.10  2025-04-11 [1] CRAN (R 4.5.0)
##  scales         1.4.0   2025-04-24 [1] CRAN (R 4.5.0)
##  sessioninfo    1.2.3   2025-02-05 [1] CRAN (R 4.5.0)
##  stringi        1.8.7   2025-03-27 [1] CRAN (R 4.5.0)
##  stringr      * 1.6.0   2025-11-04 [1] CRAN (R 4.5.0)
##  tibble       * 3.3.1   2026-01-11 [1] CRAN (R 4.5.2)
##  tidyr        * 1.3.2   2025-12-19 [1] CRAN (R 4.5.2)
##  tidyselect     1.2.1   2024-03-11 [1] CRAN (R 4.5.0)
##  tidyverse    * 2.0.0   2023-02-22 [1] CRAN (R 4.5.0)
##  timechange     0.4.0   2026-01-29 [1] CRAN (R 4.5.2)
##  tzdb           0.5.0   2025-03-15 [1] CRAN (R 4.5.0)
##  vctrs          0.7.1   2026-01-23 [1] CRAN (R 4.5.2)
##  vroom          1.7.0   2026-01-27 [1] CRAN (R 4.5.2)
##  withr          3.0.2   2024-10-28 [1] CRAN (R 4.5.0)
##  xfun           0.56    2026-01-18 [1] CRAN (R 4.5.2)
##  yaml           2.3.12  2025-12-10 [1] CRAN (R 4.5.2)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
##  * ── Packages attached to the search path.
## 
## ──────────────────────────────────────────────────────────────────────────────────────────────────

5.2 Ch. 3 Descriptive Analysis

5.2.1 Ch. 3 Set Up

# load libraries and simulated data ---------------------------------------
library(tidyverse)
library(ggpubr)

# set seed for reproducibility
set.seed(42)

# load data set (from RDS if using previous script to simulate the data)
sim_data <- read_csv(file = "sim-data.csv")

# default density plots and scale function ----------------------------------
plot_dens_2D <- function(data, x, y){
  plot <- ggplot(data = data, aes(x = x, y = y)) + 
    geom_density_2d_filled(aes(fill = after_stat(level)),
                           bins = 9, contour_var = 'ndensity') +
    scale_fill_brewer(palette = 16, direction = 1, 
                      labels = c('Low', '', '', '','','','','', 'High')) +
    guides(fill = guide_legend(title = 'Density of Occurences', 
                               direction = 'horizontal', nrow = 1, 
                               label.position = 'bottom'))
  plot
}

5.2.2 Figure 3.1

# biological sex bar plot
ggplot(sim_data, aes(x = bio_sex)) +
  geom_bar(aes(fill = bio_sex), show.legend = F) + 
  labs(title = 'Barplot of Simulated Biological Sex') +
  ylab('Count') + xlab('Biological Sex') +
  theme_classic()

5.2.3 Figure 3.2

# gender identity density plot
ggplot(sim_data, aes(x = gen_id)) +
  geom_density(alpha=.5, fill = 'forestgreen', colour = 'black') +
  labs(title = 'Density plot of Simulated\nGender Identity Responses') +
  ylab('Density') + xlab('Feminine - Masculine Continuum') +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

# gender role density plot
ggplot(sim_data, aes(x = gen_role)) +
  geom_density(alpha=.5, fill = 'purple2', colour = 'black') +
  labs(title = 'Densityplot of Simulated\nGender Role Responses') +
  ylab('Density') + xlab('Feminine - Masculine Continuum') +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

5.3 Ch. 4 SGBA of a Continuous Variable

5.3.1 Ch. 4 Set Up

# load libraries and simulated data ---------------------------------------
library(tidyverse)
library(ggpubr)

# set seed for reproducibility
set.seed(42)

# load data set (from RDS if using previous script to simulate the data)
sim_data <- read_csv(file = "sim-data.csv")

# default density plots and scale function ----------------------------------
plot_dens_2D <- function(data, x, y){
  plot <- ggplot(data = data, aes(x = x, y = y)) + 
    geom_density_2d_filled(aes(fill = after_stat(level)),
                           bins = 9, contour_var = 'ndensity') +
    scale_fill_brewer(palette = 16, direction = 1, 
                      labels = c('Low', '', '', '','','','','', 'High')) +
    guides(fill = guide_legend(title = 'Density of Occurences', 
                               direction = 'horizontal', nrow = 1, 
                               label.position = 'bottom'))
  plot
}

5.3.2 Figure 4.1

# biological sex bar plot - interaction
ggplot(sim_data, aes(x = outcome_num_pos_sex, fill = bio_sex)) + 
  geom_density(alpha = 0.7) +
  scale_fill_manual(values = c('forestgreen', 'skyblue')) +
  labs(
    title = 'Interaction example',
    subtitle = 'Density plot of continuous \noutcome and biological sex (n=30)'
    ) +
  xlab('Continuous Outcome') + 
  ylab('Response Density') +
  theme_classic() +
  theme(legend.position = 'bottom')

# biological sex bar plot - no interaction
ggplot(sim_data, aes(x = outcome_num_neg, fill = bio_sex)) + 
  geom_density(alpha = 0.7) +
  scale_fill_manual(values = c('forestgreen', 'skyblue')) +
  labs(
    title = 'No interaction example',
    subtitle = 'Density plot of continuous \noutcome and biological sex (n=30)'
    ) +
  xlab('Continuous Outcome') + 
  ylab('Response Density') +
  theme_classic() +
  theme(legend.position = 'bottom')

5.3.3 Table 4.1

# summary stats
pos_cont_sex_sum <- sim_data %>% group_by(bio_sex) %>% 
  summarise(
    n = n(), 
    `mean continuous` = round(mean(outcome_num_pos_sex),1),
    `SD continuous` = round(sd(outcome_num_pos_sex),2),
    `median continuous` = round(median(outcome_num_pos_sex),0),
    `IQR continuous` = round(IQR(outcome_num_pos_sex),)
  ) %>% 
  rename(`biological sex` = bio_sex) 

knitr::kable(pos_cont_sex_sum, booktabs = TRUE, caption = "Interation example summary statisitics: Continuous outcome and biological sex")

5.3.4 Table 4.2

neg_cont_sex_sum <- sim_data %>% group_by(bio_sex) %>% 
  summarise(
    n = n(), 
    `mean continuous` = round(mean(outcome_num_neg),1),
    `SD continuous` = round(sd(outcome_num_neg),2),
    `median continuous` = round(median(outcome_num_neg),0),
    `IQR continuous` = round(IQR(outcome_num_neg),)
  ) %>% 
  rename(`biological sex` = bio_sex) 

knitr::kable(neg_cont_sex_sum, booktabs = TRUE, caption = "No interation example summary statisitics: Continuous outcome and biological sex")

5.3.5 Table 4.3

# t-tests
con_sex_ttest_neg <- t.test(outcome_num_neg ~ bio_sex, data = sim_data)
con_sex_ttest_pos <- t.test(outcome_num_pos_sex ~ bio_sex, data = sim_data)

# example of collected data
con_sex_tab <- tibble(
  "Example" = c("Interaction", "No interaction"),
  "Test" = c("Welch's t-test", "Welch's t-test"),
  "T-score" = c(
    round(con_sex_ttest_pos$statistic,2),
    round(con_sex_ttest_neg$statistic,2)
  ),
  "95% CI" = c(
    paste0(
      "(", 
      round(con_sex_ttest_pos$conf.int[[1]],2), 
      ", ",  
      round(con_sex_ttest_pos$conf.int[[2]],2), 
      ")", sep = ""
    ),
    paste0(
      "(", 
      round(con_sex_ttest_neg$conf.int[[1]],2),
      ", ", 
      round(con_sex_ttest_neg$conf.int[[2]],2), 
      ")", sep = ""
    )
  ),
  "df" = c(
    round(con_sex_ttest_pos$parameter,1), 
    round(con_sex_ttest_neg$parameter,1)
  ),
  "p-value" = c(
    round(con_sex_ttest_pos$p.value, 3),
    round(con_sex_ttest_neg$p.value,3)
  )
)

knitr::kable(con_sex_tab, booktabs = TRUE, caption = "Statistical test of difference: Continuous outcome and biological sex.")

5.3.6 Figure 4.2

# scatter plot - interaction
ggplot(data = sim_data, aes(x = gen_id, y = outcome_num_pos_gid)) +
  coord_cartesian(ylim = c(0, 20)) + # ensures that the same y-axis length is used for both plots
  geom_jitter(size = 3, colour = "forestgreen") +
  labs(
    title = 'Interaction Example',
    subtitle = 'Scatter plot of continuous outcome \nand gender identity (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

# scatter plot - no interaction
ggplot(data = sim_data, aes(x = gen_id, y = outcome_num_neg)) +
  coord_cartesian(ylim = c(0, 20)) + # ensures that the same y-axis length is used for both plots
  geom_jitter(size = 3, colour = "forestgreen") +
  labs(
    title = 'No Interaction Example',
    subtitle = 'Scatter plot of continuous outcome \nand gender identity (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

5.3.7 Figure 4.3

# 2d density plot - interaction
ggplot(data = sim_data, aes(x = gen_id, y = outcome_num_pos_gid)) + 
  geom_density_2d_filled(
    aes(fill = after_stat(level)), bins = 9, contour_var = 'ndensity'
    ) +
  scale_fill_brewer(
    palette = 16, direction = 1, 
    labels = c('Low', '', '', '','','','','', 'High')
    ) +
  guides(fill = guide_legend(
    title = 'Density of Occurences', 
    direction = 'horizontal', nrow = 1, 
    label.position = 'bottom')
    ) +
  labs(
    title = 'Interaction Example',
    subtitle = '2-D Density plot of continuous outcome \nand gender identity (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  ) +
  geom_blank(aes(y = 20)) + # extends y axis to ensure plots are comparable
  geom_blank(aes(y = 0))

# 2d density plot - no interaction
ggplot(data = sim_data, aes(x = gen_id, y = outcome_num_neg)) + 
  geom_density_2d_filled(
    aes(fill = after_stat(level)), bins = 9, contour_var = 'ndensity'
    ) +
  scale_fill_brewer(
    palette = 16, direction = 1, 
    labels = c('Low', '', '', '','','','','', 'High')
    ) +
  guides(fill = guide_legend(
    title = 'Density of Occurences', 
    direction = 'horizontal', nrow = 1, 
    label.position = 'bottom')
    ) +
  labs(
    title = 'No Interaction Example',
    subtitle = '2-D Density plot of continuous outcome \nand gender identity (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  ) +
  geom_blank(aes(y = 20)) + # extends y axis to ensure plots are comparable
  geom_blank(aes(y = 0))

5.3.8 Table 4.5

# calculate correlations
con_gi_r_neg <- cor.test(sim_data$outcome_num_neg, sim_data$gen_id)
con_gi_r_pos <- cor.test(sim_data$outcome_num_pos_gid, sim_data$gen_id)

# output table of results
con_gi_tab <- tibble(
  "Example" = c("Interaction", "No interaction"),
  "Correlation Type" = c("Pearson's r", "Pearson's r"),
  "r" = c(
    round(con_gi_r_pos$estimate,2),
    round(con_gi_r_neg$estimate,2)
  ),
  "95% CI" = c(
    paste0(
      "(",
      round(con_gi_r_pos$conf.int[[1]],2),
      ", ",
      round(con_gi_r_pos$conf.int[[2]],2),
      ")", sep = ""
    ),
    paste0(
      "(",
      round(con_gi_r_neg$conf.int[[1]],2),
      ", ",
      round(con_gi_r_neg$conf.int[[2]],2),
      ")", sep = ""
    )
  ),
  "df" = c(
    round(con_gi_r_pos$parameter,1),
    round(con_gi_r_neg$parameter,1)
  ),
  "Significance Test" = c("t-test", "t-test"),
  "T-score" = c(
    round(con_gi_r_pos$statistic,2),
    round(con_gi_r_neg$statistic,2)
  ),
  "p-value" = c(
    round(con_gi_r_pos$p.value, 3),
    round(con_gi_r_neg$p.value,3)
  )
)

knitr::kable(con_gi_tab, booktabs = TRUE, caption = "Statistical test of difference: Continuous outcome and Gender Idenity Item from the SGBA-5.")

5.3.9 Figure 4.4

# scatter plot - interaction
ggplot(data = sim_data, aes(x = gen_id, y = outcome_num_pos_gid)) +
  coord_cartesian(ylim = c(0, 20)) + # ensures that the same y-axis length is used for both plots
  geom_jitter(size = 3, colour = "forestgreen") +
  labs(
    title = 'Gender Identity',
    subtitle = 'Scatter plot of continuous outcome \nand gender identity (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

# scatter plot - no interaction
ggplot(data = sim_data, aes(x = gen_role, y = outcome_num_pos_gid)) +
  coord_cartesian(ylim = c(0, 20)) + # ensures that the same y-axis length is used for both plots
  geom_jitter(size = 3, colour = "forestgreen") +
  labs(
    title = 'Gendered Roles',
    subtitle = 'Scatter plot of continuous outcome \nand gendered roles (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  )

5.3.10 Figure 4.5

# 2d density plot - interaction
ggplot(data = sim_data, aes(x = gen_id, y = outcome_num_pos_gid)) + 
  geom_density_2d_filled(
    aes(fill = after_stat(level)), bins = 9, contour_var = 'ndensity'
    ) +
  scale_fill_brewer(
    palette = 16, direction = 1, 
    labels = c('Low', '', '', '','','','','', 'High')
    ) +
  guides(fill = guide_legend(
    title = 'Density of Occurences', 
    direction = 'horizontal', nrow = 1, 
    label.position = 'bottom')
    ) +
  labs(
    title = 'Gender Identity',
    subtitle = '2-D Density plot of continuous outcome \nand gender identity (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  ) +
  geom_blank(aes(y = 20)) + # extends y axis to ensure plots are comparable
  geom_blank(aes(y = 0))

# 2d density plot - no interaction
ggplot(data = sim_data, aes(x = gen_role, y = outcome_num_pos_gid)) + 
  geom_density_2d_filled(
    aes(fill = after_stat(level)), bins = 9, contour_var = 'ndensity'
    ) +
  scale_fill_brewer(
    palette = 16, direction = 1, 
    labels = c('Low', '', '', '','','','','', 'High')
    ) +
  guides(fill = guide_legend(
    title = 'Density of Occurences', 
    direction = 'horizontal', nrow = 1, 
    label.position = 'bottom')
    ) +
  labs(
    title = 'Gendered Roles',
    subtitle = '2-D Density plot of continuous outcome \nand gendered roles (n=30)'
    ) +
  ylab("Continuous Outcome") + xlab("Feminine - Masculine Continuum") +
  theme_classic() +
  theme(
    legend.position = 'none',
    axis.ticks.x = element_blank(),
    axis.text.x = element_blank()
  ) +
  geom_blank(aes(y = 20)) + # extends y axis to ensure plots are comparable
  geom_blank(aes(y = 0))

5.3.11 Table 4.6

# calculate correlations
con_gi_r_pos <- cor.test(sim_data$outcome_num_pos_gid, sim_data$gen_id)
con_gr_r_pos <- cor.test(sim_data$outcome_num_pos_gid, sim_data$gen_role)

# output table of results
con_2gah_tab <- tibble(
  "Example" = c("Gender Identity", "Gendered Roles"),
  "Correlation Type" = c("Pearson's r", "Pearson's r"),
  "r" = c(
    round(con_gi_r_pos$estimate,2),
    round(con_gr_r_pos$estimate,2)
  ),
  "95% CI" = c(
    paste0(
      "(",
      round(con_gi_r_pos$conf.int[[1]],2),
      ", ",
      round(con_gi_r_pos$conf.int[[2]],2),
      ")", sep = ""
    ),
    paste0(
      "(",
      round(con_gr_r_pos$conf.int[[1]],2),
      ", ",
      round(con_gr_r_pos$conf.int[[2]],2),
      ")", sep = ""
    )
  ),
  "df" = c(
    round(con_gi_r_pos$parameter,1),
    round(con_gr_r_pos$parameter,1)
  ),
  "Significance Test" = c("t-test", "t-test"),
  "T-score" = c(
    round(con_gi_r_pos$statistic,2),
    round(con_gr_r_pos$statistic,2)
  ),
  "p-value" = c(
    round(con_gi_r_pos$p.value, 3),
    round(con_gr_r_pos$p.value,3)
  )
)

knitr::kable(con_2gah_tab, booktabs = TRUE, caption = "Statistical test of difference: Continuous outcome by Gender Idenity and Gendered Roles Items from the SGBA-5.")