--- title: "ONLINE SUPPLEMENT - Beyond sex differences in mean: meta-analysis of differences in skewness, kurtosis, and correlation" author: "**Pietro Pollo, Szymon M. Drobniak, Hamed Haselimashhadi, Malgorzata Lagisz, Ayumi Mizuno, Daniel W. A. Noble, Laura A. B. Wilson, Shinichi Nakagawa**" format: html: toc: true toc-location: left toc-depth: 3 toc-title: "**Table of Contents**" output-file: "index.html" theme: simplex embed-resources: true code-fold: show code-tools: true number-sections: true #bibliography: ./bib/ref.bib fontsize: "12" max-width: "10" code-overflow: wrap crossref: fig-title: Figure # (default is "Figure") tbl-title: Table # (default is "Table") title-delim: — # (default is ":") fig-prefix: Fig. # (default is "Figure") tbl-prefix: Tab. # (default is "Table") editor_options: chunk_output_type: console editor: markdown: wrap: sentence --- ```{r setup} #| include: false knitr::opts_chunk$set( collapse = TRUE, message = FALSE, warnings = FALSE, echo = TRUE#, #comment = "#>" ) ``` # Update [ GitHub repository ](https://github.com/pietropollo/new_effect_size_statistics) .# Introduction `metafor` package in `R` .# Content # Prerequisites ## Loading packages `R` statistical software and existing `R` packages, which you will first need to download and install.`install.packages()` to install them.`metafor` , you can execute `install.packages("metafor")` in the console (bottom left pane of `R Studio` ).```{r packages} if (!require("pacman")) {install.packages("pacman")} pacman::p_load(corrr, DT, ggdist, ggtext, here, janitor, metafor, pander, patchwork, tidyverse) options(DT.options = list(rownames = FALSE, dom = "Blfrtip", scrollX = TRUE, pageLength = 5, columnDefs = list(list(targets = '_all', className = 'dt-center')), buttons = c('copy', 'csv', 'excel', 'pdf'))) source("layout.R") ``` ## Custom functions `Enter` to have `R` ‘learn’ these custom functions.`R` code and you will see what we mean).```{r} #| code-fold: true # calculate effect sizes ---- ## skewness ---- <- function (x, output = "est" ) {<- length (x)if (output == "est" ) { # skewness estimate sqrt (n * (n - 1 )) / (n - 2 )) * 1 / n) * sum ((x - mean (x)) ^ 3 )) / 1 / n) * sum ((x - mean (x)) ^ 2 )) ^ (3 / 2 )))else if (output == "var" ) { # skewness sampling variance 6 * n * (n - 1 )) / - 2 ) * (n + 1 ) * (n + 3 ))## kurtosis ---- <- function (x, output = "est" ) {<- length (x)if (output == "est" ) { # kurtosis estimate + 1 ) * n * (n - 1 )) / ((n - 2 ) * (n - 3 ))) * sum ((x - mean (x)) ^ 4 ) / (sum ((x - mean (x)) ^ 2 ) ^ 2 ))) - 3 * ((n - 1 ) ^ 2 ) / ((n - 2 ) * (n - 3 )))else if (output == "var" ) { # kurtosis sampling variance 24 * n * ((n - 1 ) ^ 2 )) / - 3 ) * (n - 2 ) * (n + 3 ) * (n + 5 ))## Zr ---- <- # Zr estimate function (r) { 0.5 * log ((1 + r) / (1 - r))<- # Zr variance function (n) {1 / (n - 3 )## other effect sizes (lnRR and lnVR) ---- <- function (data = raw_data, # calculates other already established effect size statistics escalc (measure = m,m1i = data$ mean_male,m2i = data$ mean_female,sd1i = data$ sd_male,sd2i = data$ sd_female,n1i = data$ n_male,n2i = data$ n_female,var.names = c (paste0 (m,"_est" ),paste0 (m,"_var" )))# processing functions ---- <- function (chosen_trait = "fat_mass" ,measure = "KU_delta" ) {<- %>% filter (trait_name == chosen_trait,%in% c ("CCP-IMG" ,"HMGU" ,"JAX" ,"MRC H" ,"TCP" )) %>% mutate (type = "individual" ) %>% select (phenotyping_center,n = n_total,est = paste0 (measure, "_" , "est" ),var = paste0 (measure, "_" , "var" ),lower = paste0 (measure, "_" , "lower" ),upper = paste0 (measure, "_" , "upper" ))<- rma.mv (data = ind_effects,yi = est,V = var,test = "t" ,random = list (~ 1 | phenotyping_center, ~ 1 | strain_fig))<- data.frame (trait_name = chosen_trait,est = model$ beta[1 ],var = model$ se ^ 2 ,lower = model$ ci.lb,upper = model$ ci.ub,phenotyping_center = "Mean" ,strain_fig = "ES" )%>% bind_rows (df_model) %>% mutate (est_type = measure,centre_and_strain = factor (paste0 (phenotyping_center," \n " ,%>% mutate (centre_and_strain = factor (centre_and_strain,levels = c ("Mean \n ES" ,rev (levels (centre_and_strain)[- 6 ]))))<- function (chosen_trait_1 = "fat_mass" ,chosen_trait_2 = "heart_weight" ) {<- %>% filter (trait_name %in% c (chosen_trait_1,%in% c ("CCP-IMG" ,"HMGU" ,"JAX" ,"MRC H" ,"TCP" )) %>% pivot_wider (id_cols = c (specimen_id,names_from = trait_name) %>% clean_names () %>% drop_na () %>% group_by (strain_fig,%>% group_modify (~ correlate (.x)) %>% drop_na (all_of (chosen_trait_2)) %>% ungroup () %>% left_join (df_raw %>% filter (trait_name %in% c (chosen_trait_1,%in% c ("CCP-IMG" ,"HMGU" ,"JAX" ,"MRC H" ,"TCP" )) %>% pivot_wider (id_cols = c (specimen_id,names_from = trait_name) %>% clean_names () %>% drop_na () %>% group_by (strain_fig,%>% summarise (n = n ())) %>% rename (r_est = chosen_trait_2) %>% mutate (zr_est = r.to.zr (r_est),zr_var = zr.variance (n)) %>% select (- c (4 : 6 )) %>% pivot_wider (names_from = sex,values_from = c (n,%>% mutate (delta_zr_est = zr_est_male - zr_est_female,delta_zr_var = zr_var_male + zr_var_female,delta_zr_upper = delta_zr_est + qt (0.975 , n_male + n_female - 2 ) * sqrt (delta_zr_var),delta_zr_lower = delta_zr_est - qt (0.975 , n_male + n_female - 2 ) * sqrt (delta_zr_var))<- rma.mv (data = df_effects_cor,yi = delta_zr_est,V = delta_zr_var,test = "t" ,random = list (~ 1 | phenotyping_center, ~ 1 | strain_fig))<- data.frame (delta_zr_est = mlma_zr$ beta[1 ],delta_zr_lower = mlma_zr$ ci.lb,delta_zr_upper = mlma_zr$ ci.ub,phenotyping_center = "Mean" ,strain_fig = "ES" )%>% bind_rows (df_model) %>% mutate (centre_and_strain = factor (paste0 (phenotyping_center," \n " ,%>% mutate (centre_and_strain = factor (centre_and_strain,levels = c ("Mean \n ES" ,rev (levels (centre_and_strain)[- 5 ]))))# visualisation functions ---- <- function (chosen_trait = "fat_mass" ,measure = "KU_delta" ) {<- process.ind_effects (chosen_trait = chosen_trait,measure = measure) %>% ggplot (aes (y = centre_and_strain,x = est,xmax = upper,xmin = lower,shape = strain_fig,col = phenotyping_center)) + geom_pointrange () + geom_vline (xintercept = 0 ,linetype = "dotted" ) + theme_classic () + theme (legend.position = "none" ,axis.text.y = element_blank (),axis.title.y = element_blank (),plot.tag.position = c (0.15 , 0.98 ))if (measure == "ROM" ) {+ labs (x = "lnRR" ) + scale_x_continuous (limits = c (- 0.51 , 0.51 ),breaks = c (- 0.5 , 0 , 0.5 )) + theme (axis.title.x = ggtext:: element_markdown (face = "italic" ))else if (measure == "VR" ) {+ labs (x = "lnVR" ) + scale_x_continuous (limits = c (- 1 , 1 ),breaks = c (- 1 , 0 , 1 )) + theme (axis.title.x = ggtext:: element_markdown (face = "italic" ))else if (measure == "SK_delta" ) {+ labs (x = "Δ*sk*" ) + scale_x_continuous (limits = c (- 2.1 , 2.1 ),breaks = c (- 2 , 0 , 2 )) + theme (axis.title.x = ggtext:: element_markdown ())else if (measure == "KU_delta" ) {+ labs (x = "Δ*ku*" ) + scale_x_continuous (limits = c (- 15 , 15 ),breaks = c (- 15 , 0 , 15 )) + theme (axis.title.x = ggtext:: element_markdown ())<- function (chosen_trait = "fat_mass" ) {%>% filter (trait_name == chosen_trait,%in% c ("CCP-IMG" ,"HMGU" ,"JAX" ,"MRC H" ,"TCP" )) %>% add_row (phenotyping_center = "Mean" ,strain_fig = "ES" ) %>% mutate (centre_and_strain = factor (paste0 (phenotyping_center," \n " ,%>% mutate (centre_and_strain = factor (centre_and_strain,levels = c ("Mean \n ES" ,rev (levels (centre_and_strain)[- 5 ]))),value_s = scale (value)) %>% ggplot (aes (x = value_s,y = centre_and_strain,fill = sex,linetype = sex)) + stat_slab (scale = 0.7 , alpha = 0.4 ,linewidth = 0.6 ,col = "black" ) + scale_fill_manual (values = c ("white" ,"black" )) + scale_linetype_manual (values = c ("solid" ,"dashed" )) + labs (x = paste0 (str_to_sentence (str_replace_all (chosen_trait,"_" ," " ))," \n (scaled)" ),y = "Phenotyping centre and mice strain" ) + theme_classic () + theme (legend.position = "none" ,axis.title.x = element_text (size = 12 , margin = margin (t = 0.2 ,unit = "cm" )),axis.title.y = element_text (size = 12 , margin = margin (r = 0.2 ,unit = "cm" )),axis.text.x = element_text (size = 10 ),axis.text.y = element_text (size = 10 ),plot.tag.position = c (0.53 , 0.98 ))<- function (chosen_trait_1 = "fat_mass" ,chosen_trait_2 = "heart_weight" ) {process.cor_effects (chosen_trait_1 = chosen_trait_1,chosen_trait_2 = chosen_trait_2) %>% ggplot (aes (y = centre_and_strain,x = delta_zr_est,xmax = delta_zr_upper,xmin = delta_zr_lower,shape = strain_fig,col = phenotyping_center)) + geom_pointrange () + geom_vline (xintercept = 0 ,linetype = "dotted" ) + labs (y = "Phenotyping centre and mice strain" ,x = "Δ*Zr*" , shape = "Strain" ) + scale_x_continuous (limits = c (- 1 , 1 ),breaks = c (- 1 , 0 , 1 )) + theme_classic () + theme (legend.position = "none" ,axis.title.x = ggtext:: element_markdown (size = 12 , margin = margin (t = 0.2 ,unit = "cm" )),axis.title.y = element_text (size = 12 ,margin = margin (r = - 0.1 ,unit = "cm" )),axis.text.x = element_text (size = 10 ),axis.text.y = element_text (size = 10 ),plot.tag.position = c (0.3 , 0.99 ))<- function (chosen_trait_1 = "fat_mass" ,chosen_trait_2 = "heart_weight" ,chosen_lims = c (- 3 , 5 )) {<- %>% filter (trait_name %in% c (chosen_trait_1,%in% c ("CCP-IMG" ,"HMGU" ,"JAX" ,"MRC H" ,"TCP" )) %>% pivot_wider (id_cols = c (specimen_id,names_from = trait_name) %>% clean_names () %>% drop_na () %>% mutate (centre_and_strain = factor (paste0 (phenotyping_center,%>% mutate (centre_and_strain = factor (centre_and_strain,levels = rev (levels (centre_and_strain))),trait_1_s = scale (get (chosen_trait_1))[,1 ],trait_2_s = scale (get (chosen_trait_2))[,1 ])<- list ()for (i in 1 : length (levels (df_cor$ centre_and_strain))) {<- sort (levels (df_cor$ centre_and_strain))[i]<- %>% filter (centre_and_strain == level_i) %>% ggplot (aes (x = trait_1_s,y = trait_2_s,shape = sex,linetype = sex)) + geom_point (alpha = 0.008 ,+ geom_abline (intercept = 0 ,slope = 1 ,linewidth = 0.5 ,linetype = "dotted" ) + geom_smooth (method = "lm" ,se = F,col = "black" ) + scale_shape_manual (values = c (3 , 4 )) + scale_linetype_manual (values = c ("solid" ,"dashed" )) + scale_x_continuous (limits = chosen_lims) + scale_y_continuous (limits = chosen_lims) + labs (x = paste0 (str_to_sentence (str_replace_all (chosen_trait_1, "_" , " " ))," \n (scaled)" ),y = paste0 (str_to_sentence (str_replace_all (chosen_trait_2, "_" , " " ))," (scaled)" )) + theme_classic () + theme (legend.position = "none" ,plot.tag.position = c (0.05 , 0.91 ),axis.title.x = element_text (size = 12 , margin = margin (t = 0.2 ,unit = "cm" )),axis.title.y = element_text (size = 12 , margin = margin (r = 0.2 ,unit = "cm" )),axis.text.x = element_text (size = 10 ),axis.text.y = element_text (size = 10 ))if (i != 6 ) {<- + theme (axis.title.x = element_blank (),axis.text.x = element_blank (),axis.line.x = element_blank (),axis.ticks.x = element_blank ())<- plotreturn (plot_list)``` # Equations and custom functions to calculate effect sizes ## Skewness [ \frac{1}{n} \sum_{i = 1}^{n}(x - \bar{x}) ^ 2 ] ^ \frac{3}{2}}## Kurtosis [ \sum_{i = 1}^{n}(x_{i} - \bar{x}) ^ 2 ] ^ 2} -## Zr # Data loading and preparation ```{r data} # raw data ---- df_raw <- read_csv("mice_data_sample.csv") %>% # small adjustments to make plots more readable: mutate(phenotyping_center = ifelse(phenotyping_center == "MRC Harwell", "MRC H", phenotyping_center), strain_fig = case_when(strain_accession_id == "MGI:2159965" ~ "N", strain_accession_id == "MGI:2683688" ~ "NCrl", strain_accession_id == "MGI:2164831" ~ "NTac", strain_accession_id == "MGI:3056279" ~ "NJ", strain_accession_id == "MGI:2160139" ~ "NJcl")) df_meta_analysed <- df_raw %>% group_by(sex, trait_name, phenotyping_center, strain_fig) %>% summarize(mean = mean(value, na.rm = T), sd = sd(value, na.rm = T), n = n(), SK_est = calc.skewness(value), SK_var = calc.skewness(value, output = "var"), KU_est = calc.kurtosis(value), KU_var = calc.kurtosis(value, output = "var")) %>% pivot_wider(id_cols = c(trait_name, phenotyping_center, strain_fig), names_from = sex, values_from = c(mean:KU_var)) %>% mutate(SK_delta_est = SK_est_male - SK_est_female, SK_delta_var = SK_var_male + SK_var_female, KU_delta_est = KU_est_male - KU_est_female, KU_delta_var = KU_var_male + KU_var_female) %>% bind_cols(calc.effect(., m = "ROM")) %>% # lnRR bind_cols(calc.effect(., m = "CVR")) %>% # lnCVR bind_cols(calc.effect(., m = "VR")) %>% # lnVR filter(!is.na(CVR_est)) %>% mutate(n_total = n_female + n_male, prop_females = n_female / (n_female + n_male)) %>% select(trait_name, phenotyping_center, strain_fig, n_total, prop_females, ROM_est, ROM_var, CVR_est, CVR_var, VR_est, VR_var, SK_delta_est, SK_delta_var, KU_delta_est, KU_delta_var) %>% mutate(ROM_upper = ROM_est + qt(0.975, n_total - 1) * sqrt(ROM_var), ROM_lower = ROM_est - qt(0.975, n_total - 1) * sqrt(ROM_var), CVR_upper = CVR_est + qt(0.975, n_total - 1) * sqrt(CVR_var), CVR_lower = CVR_est - qt(0.975, n_total - 1) * sqrt(CVR_var), VR_upper = VR_est + qt(0.975, n_total - 1) * sqrt(VR_var), VR_lower = VR_est - qt(0.975, n_total - 1) * sqrt(VR_var), SK_delta_upper = SK_delta_est + qt(0.975, n_total - 1) * sqrt(SK_delta_var), SK_delta_lower = SK_delta_est - qt(0.975, n_total - 1) * sqrt(SK_delta_var), KU_delta_upper = KU_delta_est + qt(0.975, n_total - 1) * sqrt(KU_delta_var), KU_delta_lower = KU_delta_est - qt(0.975, n_total - 1) * sqrt(KU_delta_var)) ``` # Meta-analytical models ## Single variable effect sizes ```{r} map2_dfr (.x = rep (c ("fat_mass" ,"heart_weight" ,"glucose" ,"total_cholesterol" ),each = 4 ),.y = rep (c ("ROM" ,"VR" ,"SK_delta" ,"KU_delta" ), 4 ),.f = process.ind_effects) %>% mutate (est_type = case_when (est_type == "ROM" ~ "lnRR" ,== "VR" ~ "lnVR" ,== "SK_delta" ~ "delta_sk" ,== "KU_delta" ~ "delta_ku" )) %>% datatable (.,extensions = "Buttons" ,rownames = FALSE )``` ## Correlational effect sizes ```{r} map2_dfr (.x = c ("fat_mass" ,"glucose" ),.y = c ("heart_weight" ,"total_cholesterol" ),.f = process.cor_effects) %>% mutate (relationship = rep (c ("fat mass and heart weight" ,"glucose and total cholesterol" ),each = 7 )) %>% relocate (relationship) %>% datatable (.,extensions = "Buttons" ,rownames = FALSE )``` # Visualisations ```{r} ## figure_2 ---- <- list ()1 ]] <- ridgeline.custom ("fat_mass" )2 : 5 ] <- map2 (.x = rep ("fat_mass" , 4 ),.y = c ("ROM" ,"VR" ,"SK_delta" ,"KU_delta" ),.f = caterpillar.custom)6 ]] <- ridgeline.custom ("heart_weight" )7 : 10 ] <- map2 (.x = rep ("heart_weight" , 4 ),.y = c ("ROM" ,"VR" ,"SK_delta" ,"KU_delta" ),.f = caterpillar.custom)<- %>% wrap_plots (ncol = 5 ) + plot_annotation (tag_levels = "A" ))## figure_3 ---- <- list ()1 : 6 ] <- cor.plot.custom (chosen_trait_1 = "fat_mass" ,chosen_trait_2 = "heart_weight" )7 ]] <- cor.caterpillar.custom (chosen_trait_1 = "fat_mass" ,chosen_trait_2 = "heart_weight" )8 : 13 ] <- cor.plot.custom (chosen_trait_1 = "glucose" ,chosen_trait_2 = "total_cholesterol" )14 ]] <- cor.caterpillar.custom (chosen_trait_1 = "glucose" ,chosen_trait_2 = "total_cholesterol" )<- %>% wrap_plots () + plot_layout (design = layout_2,heights = c (rep (1 , 6 ), 0.6 ),widths = c (rep (0.23 , 2 ), 0.02 , rep (0.23 , 2 )),axes = "collect" ,guides = "collect" ) + plot_annotation (tag_levels = list (c ("A" , rep ("" , 5 ), "B" ,"C" ,rep ("" , 5 ),"D" ))))## figure_4 ---- <- list ()1 ]] <- ridgeline.custom ("glucose" )2 : 5 ] <- map2 (.x = rep ("glucose" , 4 ),.y = c ("ROM" ,"VR" ,"SK_delta" ,"KU_delta" ),.f = caterpillar.custom)6 ]] <- ridgeline.custom ("total_cholesterol" )7 : 10 ] <- map2 (.x = rep ("total_cholesterol" , 4 ),.y = c ("ROM" ,"VR" ,"SK_delta" ,"KU_delta" ),.f = caterpillar.custom)<- %>% wrap_plots (ncol = 5 ) + plot_annotation (tag_levels = "A" ))``` # Software and package versions ```{r versions} sessionInfo() %>% pander() ``` 
