# Sexual dimorphism in trait variability and its eco-evolutionary and statistical implications

### Electronic Supporting File (doi:10.5281/zenodo.4146948)

#### Susanne Zajitschek, Felix Zajitschek, Russell Bonduriansky, Robert Brooks, Will Cornwell, Daniel Falster, Malgorzata Lagisz, Jeremy Mason, Alistair Senior, Daniel Noble, & Shinichi Nakagawa

#### 2020-10-28

# 1 General Introduction

This supplementary material follows the general structure described in our methods and is outlined here in Figure 1.1. Generally speaking, the first few sections describe the processes by which the data were cleaned and re-organised for analysis. This is followed by how the first order and second order meta-analyses were conducted. Supplemental tables and figures that are referenced within the main manuscript can be found at the end of the document. Various parts of the workflow can be conveniently referenced using the table of contents panel (located along the side). We will refer back to this figure as we describe the various aspects of the code and data processing in relevant places.

Figure 1.1: Workflow of data processing and meta-analysis [Figure 3 in manuscipt]

# 2 Set-up

## 2.1 Packages

    library(pacman)
    pacman::p_load(readr, dplyr, metafor, devtools, purrr, tidyverse, tidyr, tibble,
        kableExtra, robumeta, ggpubr, ggplot2, png, grid, here, pander)

## 2.2 Functions

Functions needed for preparing the data for meta analyses

### 2.2.1 Loading and cleaning data

These functions will load the raw and and do some cleaning to make it ready for further processing and analysis.

    # loads the raw data, setting some default types for various columns

    load_raw <- function(filename) {
        read_csv(filename, col_types = cols(.default = col_character(), project_id = col_character(),
            id = col_character(), parameter_id = col_character(), age_in_days = col_integer(),
            date_of_experiment = col_datetime(format = ""), weight = col_double(), phenotyping_center_id = col_character(),
            production_center_id = col_character(), weight_date = col_datetime(format = ""),
            date_of_birth = col_datetime(format = ""), procedure_id = col_character(),
            pipeline_id = col_character(), biological_sample_id = col_character(), biological_model_id = col_character(),
            weight_days_old = col_integer(), datasource_id = col_character(), experiment_id = col_character(),
            data_point = col_double(), age_in_weeks = col_integer(), `_version_` = col_character()))
    }

    # Apply some standard cleaning to the data
    clean_raw_data <- function(mydata) {

        group <- read_csv(here("data", "ParameterGrouping.csv"))

        tmp <- mydata %>%
        # Filter to IMPC source (recommend by Jeremey in email to Susi on 20 Aug 2018)
        filter(datasource_name == "IMPC") %>%
        # standardise trait names
        mutate(parameter_name = tolower(parameter_name)) %>%
        # remove extreme ages
        filter(age_in_days > 0 & age_in_days < 500) %>%
        # remove NAs
        filter(!is.na(data_point)) %>%
        # subset to reasonable set of variables, date_of_experiment used as an indicator
        # of batch-level effects
        select(production_center, strain_name, strain_accession_id, biological_sample_id,
            pipeline_stable_id, procedure_group, procedure_name, sex, date_of_experiment,
            age_in_days, weight, parameter_name, data_point) %>%
        # sort
        arrange(production_center, biological_sample_id, age_in_days)

        # filter to groups with > 1 centre
        merge(tmp, tmp %>% group_by(parameter_name) %>% summarise(center_per_trait = length(unique(production_center,
            na.rm = TRUE)))) %>% filter(center_per_trait >= 2) %>%
        # Define population variable
        mutate(population = sprintf("%s-%s", production_center, strain_name)) %>%
        # add grouping variable: these were decided based on functional groups and
        # procedures
        mutate(parameter_group = group$parameter[match(parameter_name, group$parameter_name)]) %>%

        # Assign unique IDs (per trait) each unique parameter_name (=trait,use trait
        # variable) gets a unique number ('id')

        # We add a new variable, where redundant traits are combined [note however, at
        # this stage the dataset still contains nonsensical traits, i.e. traits that may
        # not contain any information on variance]
        mutate(id = match(parameter_name, unique(parameter_name))) %>% as_tibble()
    }

### 2.2.2 Sub-setting data

Create a function for sub-setting the data to choose only one data point per individual per trait: “data_subset_parameterid_individual_by_age”

    data_subset_parameterid_individual_by_age <- function(mydata, parameter, age_min = 0,
        age_center = 100) {
        tmp <- mydata %>% filter(age_in_days >= age_min, id == parameter) %>% # Take results for single individual closest to age_center
        mutate(age_diff = abs(age_center - age_in_days)) %>% group_by(biological_sample_id) %>%
            filter(age_diff == min(age_diff)) %>% select(-age_diff)

        # Still some individuals with multiple records (because same individual appears
        # under different procedures, so filter to one record)
        j <- match(unique(tmp$biological_sample_id), tmp$biological_sample_id)
        tmp[j, ]
    }

### 2.2.3 Population statistics

Create a function called: “calculate_population_stats” This function groups animals from the same strain and institution together. This is done for each trait separately, and only for traits that have been measured in both sexes. Any group containing fewer than 6 individuals is excluded.

    # Function re-organises the data so that the male and female data are side-by-side as columns to make effect size calculations easier.
    multi_spread <- function(df, key, value) {
        # quote key
          keyq <- rlang::enquo(key)
        # break value vector into quotes
        valueq <- rlang::enquo(value)
             s <- rlang::quos(!!valueq)

             df                           %>%
            gather(variable, value, !!!s) %>%
            unite(temp, !!keyq, variable) %>%
            spread(temp, value)
    }

    # Function will calculate population stats and exclude data that is irrelevant and not useful, it will also just convert the data straight away to wide format making it ready for effect size calculations

    calculate_population_stats <- function(mydata, min_individuals = 5){
        mydata %>%
        group_by(population,
                 production_center,
                 strain_name,
                 sex) %>%
        summarise(trait  = parameter_name[1],
                  x_bar  = mean(data_point, na.rm = TRUE),
                  x_sd   =   sd(data_point, na.rm = TRUE),
                  n_ind  = n()) %>%
        ungroup() %>%
        filter(n_ind > min_individuals) %>%
        group_by(population) %>%
        filter(all(c("male", "female") %in% sex)) %>% # Removes any entries where only 1 sex is present
        multi_spread(sex, c(x_bar, x_sd, n_ind))      # Spreads the data out into wide format, so sexes are beside each other which is more typical of meta-analytic data
    }

### 2.2.4 Calculate effect sizes and sample variances

Function: “create_meta_analysis_effect_sizes”

    # This function takes the data and computes lnCVR, lnVR and ROM (lnRR) using male
    # and female data to generate effect size statistics for meta-analyses
    create_meta_analysis_effect_sizes <- function(mydata, measure = c("CVR", "VR", "ROM")) {
        for (i in 1:length(measure)) {
            mydata <- metafor::escalc(m1i = male_x_bar, m2i = female_x_bar, sd1i = male_x_sd,
                sd2i = female_x_sd, n1i = male_n_ind, n2i = female_n_ind, data = mydata,
                measure = measure[i], var.names = c(paste0("effect_size_", measure[i]),
                    paste0("sample_variance_", measure[i])), append = TRUE)
        }
        return(mydata)
    }

## 2.3 Load & clean data

We have already done this step and provide a cleaned up file which is less computationally intensive to deal with. The file has been saved in a folder called `export`. However, the `.csv` is provided in case this is preferred and can be imported and processed using following the steps below: [This code is currently excluded from showing in this `.html` file, but can be viewed and enabled / run in the down-loadable .rmd file. However, it requires the .gz file containing the raw data in the “export” folder. Unfortunately, that file can not be provided via Github, as it is too large (274MB). However, it is freely accessible and uploaded to zenodo.org (https://zenodo.org/deposit/3759701)

    # Load raw data - save cleaned dataset as RDS for reuse
    data_raw <- load_raw(here("data", "dr7.0_all_control_data.csv.gz"))
    dir.create("export", F, F)

    data <- data_raw %>% clean_raw_data()
    saveRDS(data, "export/data_clean.rds")

For analysis we load the RDS created above and other datasets

    data <- readRDS(here("export", "data_clean.rds"))

    procedures <- read_csv(here("data", "procedures.csv"))

## 2.4 Mean-variance relationships

We can also have a look at the mean-variance relationships in these data (Figure ??). Note that these data have not been filtered to the same extent that the data is that is used in the meta-analyses, but this is simply here to demonstrate that mean-variance relationships exist and are strong providing strong justification for using the log coefficient of variation ratio (lnCVR).

    # First, clean up the data

     results <- data %>%
                group_by(population,
                 production_center,
                 strain_name,
                 parameter_name,
                 sex) %>%
        summarise(trait  = parameter_name[1],
                  x_bar  = mean(data_point, na.rm = TRUE),
                  x_sd   =   sd(data_point, na.rm = TRUE),
                  n_ind  = n()) %>%
        ungroup() %>%
        filter(n_ind > 5) %>%   # NOTE here that you are excluding data with 5 exactly....
        group_by(population) %>%
        filter(all(c("male", "female") %in% sex)) %>% # Removes any entries where only 1 sex is present
        multi_spread(sex, c(x_bar, x_sd, n_ind)) %>%
        mutate(lnMean_m = log(male_x_bar),
               logSD_m  = log(male_x_sd),
               lnMean_f = log(female_x_bar),
               lnSD_f   = log(female_x_sd)) %>%
        filter(!is.na(logSD_m) & !is.inf(logSD_m)) %>%
        filter(!is.na(lnSD_f) & !is.inf(lnSD_f))

    male <- ggplot(results, aes(x = lnMean_m, y = logSD_m)) +
        geom_point(color = "#66C2A5", size = 2) +
       geom_smooth(method = "lm", se = FALSE, color="black") +
        labs(x = "log (Mean)",
             y = "log (SD)") +
        theme_classic()

    female <- ggplot(results, aes(x = lnMean_f, y = lnSD_f)) +
        geom_point(color = "#FC8D62", size = 2) +
      geom_smooth(method = "lm", se = FALSE, color="black")+
        labs(x = "log (Mean)",
             y = "log (SD)") +
        theme_classic()

    ggarrange(male, female, ncol = 2, nrow = 1, labels = c("A", "B"))

Mean-variance relationships (log(Mean) vs log(SD, standard deviation)) across all traits for males (A) and females (B).[Figure 1-figure supplement 1 in manuscript]

We can see that there are strong relationships between the log transformed means and standard deviations for both males (Figure ??A) and females (Figure ??B), with the Person’s correlation coefficient being 0.949 and 0.95, respectively.

# 3 Meta-analyses

## 3.1 Population as analysis unit

In this section, we cover the workflow described in Figure 1.1C and D.

### 3.1.1 Loop through meta-analyses on all traits

This section covers Figure 1.1 C, by conducting a meta-analysis for all traits and producing the meta-analysed effect sizes (lnVR, lnCVR and lnRR).

- The loop combines the functions mentioned above and fills the data matrix with results from our meta analysis.
- Error messages indicate traits that either did not reach convergence, or that did not return meaningful results in the meta-analysis, due to absence of variance. Those traits will be removed in later steps, outlined below.

    n <- length(unique(data$id))

    # Create dataframe to store results
    results_alltraits_grouping <- data.frame(tibble(id = 1:n, lnCVR = 0, lnCVR_lower = 0,
        lnCVR_upper = 0, lnCVR_se = 0, lnCVR_I2 = 0, lnVR = 0, lnVR_lower = 0, lnVR_upper = 0,
        lnVR_se = 0, lnVR_I2 = 0, lnRR = 0, lnRR_lower = 0, lnRR_upper = 0, lnRR_se = 0,
        lnRR_I2 = 0, k = 0, trait = 0, male_N = 0, female_N = 0, min_male_N = 0, max_male_N = 0,
        mean_male_N = 0, min_female_N = 0, max_female_N = 0, mean_female_N = 0))

    for (t in 1:n) {
        tryCatch({

            results <- data %>% data_subset_parameterid_individual_by_age(t) %>% calculate_population_stats() %>%
                create_meta_analysis_effect_sizes() %>% mutate(err = seq_len(n()))

            # lnCVR, log repsonse-ratio of the coefficient of variance
            cvr <- metafor::rma.mv(yi = effect_size_CVR, V = sample_variance_CVR, random = list(~1 |
                strain_name, ~1 | production_center, ~1 | err), control = list(optimizer = "optim",
                optmethod = "Nelder-Mead", maxit = 1000), verbose = F, data = results)

            # lnVR, comparison of standard deviations
            cv <- metafor::rma.mv(yi = effect_size_VR, V = sample_variance_VR, random = list(~1 |
                strain_name, ~1 | production_center, ~1 | err), control = list(optimizer = "optim",
                optmethod = "Nelder-Mead", maxit = 1000), verbose = F, data = results)

            # for means, lnRR
            means <- metafor::rma.mv(yi = effect_size_ROM, V = sample_variance_ROM, random = list(~1 |
                strain_name, ~1 | production_center, ~1 | err), control = list(optimizer = "optim",
                optmethod = "Nelder-Mead", maxit = 1000), verbose = F, data = results)

            f <- function(x) unlist(x[c("b", "ci.lb", "ci.ub", "se")])

            extract_I2 <- function(mod) {
                sigma2 <- mod$sigma2
                w <- mod$vi
                k <- mod$k

                sigmaM <- sum(w * (k - 1))/(sum(w)^2 - sum(w^2))

                I2_tot <- round(sum(sigma2)/sum(sigma2 + sigmaM), digits = 3) * 100
                return(I2_tot)
            }
            results_alltraits_grouping[t, c(2:17, 19:26)] <- c(f(cvr), lnCVR_I2 = extract_I2(cvr),
                f(cv), lnVR_I2 = extract_I2(cv), f(means), lnRR_I2 = extract_I2(means),
                k = means$k, male_N = sum(results$male_n_ind), female_N = sum(results$female_n_ind),
                min_male_N = min(results$male_n_ind), max_male_N = max(results$male_n_ind),
                mean_male_N = mean(results$male_n_ind), min_female_N = min(results$female_n_ind),
                max_female_N = max(results$female_n_ind), mean_female_N = mean(results$female_n_ind))
            results_alltraits_grouping[t, 18] <- unique(results$trait)
        }, error = function(e) {
            cat("ERROR :", t, conditionMessage(e), "\n")
        })
    }

    ## ERROR : 84 Optimizer (optim) did not achieve convergence (convergence = 10).
    ## ERROR : 160 Processing terminated since k <= 1.
    ## ERROR : 168 Processing terminated since k <= 1.

In the above function, we use ‘tryCatch’ and ‘conditionMessage’ to prevent the loop from aborting when the first error at row 84 is produced. As convergence in the two listed non-converging cases can’t be achieved by sensibly tweaking (other optim etc.), and we only learn about non-convergence in the loop, it is not possible to exclude the traits (N = 2) beforehand.Similarly, there are 8 traits with very low variation, which can not be excluded prior to running the loop.

Any produced “Warnings” indicate cases where variance components are set to zero during likelihood optimization.

### 3.1.2 Merging datasets & removal of non-converged traits

Here we cover Figure 1.1D, removing non-coverging model results and non-informative traits. Procedure names, grouping variables and trait names (“parameter_names”) are the merged back together with the results from the metafor analysis above.

    results_alltraits_grouping2 <- results_alltraits_grouping %>%
    # Join data with results_alltraits_grouping. We filter duplicated id's to get
    # only one unique row per id (and there is one id per parameter_name)
    left_join(by = "id", data %>% select(id, parameter_group, procedure = procedure_name,
        procedure_name, parameter_name) %>% filter(!duplicated(id))) %>%
    # Below we add 'procedure' (from the previously loaded 'procedures.csv') as a
    # variable; n <- length(unique(results_alltraits_grouping2$parameter_name)))
    # should equal 232
    left_join(by = "procedure", procedures %>% distinct())

    # We exclude 14 parameter names for which metafor models didn't converge ('dp t
    # cells', 'mzb (cd21/35 high)'), and of parameters that don't harbour enough
    # variation
    meta_clean <- results_alltraits_grouping2 %>% filter(!parameter_name %in% c("dp t cells",
        "mzb (cd21/35 high)", "number of caudal vertebrae", "number of cervical vertebrae",
        "number of digits", "number of lumbar vertebrae", "number of pelvic vertebrae",
        "number of ribs left", "number of ribs right", "number of signals", "number of thoracic vertebrae",
        "total number of acquired events in panel a", "total number of acquired events in panel b",
        "whole arena permanence"))

    # Summary of the cleaned data for 218 traits
    sd(meta_clean$k)
    range(meta_clean$male_N)
    range(meta_clean$female_N)
    range(meta_clean$lnCVR_I2)
    range(meta_clean$lnVR_I2)
    range(meta_clean$lnRR_I2)
    range(meta_clean$k)
    mean(meta_clean$mean_male_N)
    range(meta_clean$min_male_N)
    median(meta_clean$mean_male_N)
    sd(meta_clean$mean_male_N)
    sd(meta_clean$mean_male_N)/sqrt(length(meta_clean$mean_male_N))
    mean(meta_clean$mean_female_N)
    median(meta_clean$mean_female_N)
    range(meta_clean$min_female_N)
    sd(meta_clean$mean_female_N)
    sd(meta_clean$mean_female_N)/sqrt(length(meta_clean$mean_female_N))

A total of 14 traits from the original 232 that had been included are removed because they either did not achieve convergence or are nonsensical for analysis of variance (such as traits that show no variation, see list below).

The traits where models did not converge include: “dp t cells”, “mzb (cd21/35 high)”

The traits where there was not enough variation include: “number of caudal vertebrae”, “number of cervical vertebrae”, “number of digits”, “number of lumbar vertebrae”, “number of pelvic vertebrae”, “number of ribs left”,“number of ribs right”, “number of signals”, “number of thoracic vertebrae”, “total number of acquired events in panel a”,“total number of acquired events in panel b”, “whole arena permanence”.

## 3.2 Meta-analysis: condensing non-independent traits

This dataset contained a number of highly correlated traits, such as different kinds of cell counts (for example hierarchical parameterization within immunological assays). As those data-points are not independent of each other, we conducted meta analyses on these correlated parameters to collapse the number of levels. This section describes the workflow depicted in Figure 1.1E & F.

### 3.2.1 Collapsing and merging correlated parameters

Count the number of parameter names (correlated sub-traits) in each parameter group (par_group_size). This serves to identify and separate the traits that are correlated from the full dataset that can be processed as is. If the sample size (n) for a given “parameter group” equals 1, the trait is unique and uncorrelated. All instances, where there are 2 or more traits associated with the same parameter group (90 cases), are selected for a “mini-meta analysis”, which removes the issue of correlation. The workflow here relates to Figure 1.1E.

    # Here we double check numbers of trait parameters in the dataset
    meta1 <- meta_clean
    # length(unique(meta1$procedure)) #18 length(unique(meta1$GroupingTerm)) #9
    # length(unique(meta1$parameter_group)) # 148 levels. To be used as grouping
    # factor for meta-meta-analysis / collapsing down based on things that are
    # classified identically in 'parameter_group' but have different 'parameter_name'
    # length(unique(meta1$parameter_name)) #218 and prepare the dataset for
    # second-order meta-analysis
    meta1_sub <- meta1 %>% # Add summary of number of parameter names in each parameter group
    group_by(parameter_group) %>% mutate(par_group_size = length(unique(parameter_name)),
        sampleSize = as.numeric(k)) %>% ungroup() %>% # Create subsets with > 1 count (par_group_size > 1)
    filter(par_group_size > 1)  # 90 observations

### 3.2.2 Meta-analyses on correlated (sub-)traits, using `robumeta`

Here we pepare the subset of the data (using nest()), and in this first step the model of the meta-analysis effect sizes are calculated. This section specifically undertakes the workflow described in Figure 1.1 F.

    # Create summary of number of parameter names in each parameter group, and merge
    # back together

    meta1b <- meta1 %>% group_by(parameter_group) %>% summarize(par_group_size = length(unique(parameter_name,
        na.rm = TRUE)))

    meta1$par_group_size <- meta1b$par_group_size[match(meta1$parameter_group, meta1b$parameter_group)]

    # Create subsets with > 1 count (par_group_size > 1)

    meta1_sub <- subset(meta1, par_group_size > 1)  # 90 observations
    meta1_sub$sampleSize <- as.numeric(meta1_sub$k)

    # Nesting and meta-analyses on correlated traits, using robumeta

    n_count <- meta1_sub %>% group_by(parameter_group) %>% mutate(raw_N = sum(sampleSize)) %>%
        nest() %>% ungroup()

    model_count <- n_count %>% mutate(model_lnRR = map(data, ~robu(.x$lnRR ~ 1, data = .x,
        studynum = .x$id, modelweights = c("CORR"), rho = 0.8, small = TRUE, var.eff.size = (.x$lnRR_se)^2)),
        model_lnVR = map(data, ~robu(.x$lnVR ~ 1, data = .x, studynum = .x$id, modelweights = c("CORR"),
            rho = 0.8, small = TRUE, var.eff.size = (.x$lnVR_se)^2)), model_lnCVR = map(data,
            ~robu(.x$lnCVR ~ 1, data = .x, studynum = .x$id, modelweights = c("CORR"),
                rho = 0.8, small = TRUE, var.eff.size = (.x$lnCVR_se)^2)))

    #### Extracting and save parameter estimates Here we apply an additional function to
    #### collect the outcomes of the 'mini-meta-analysis' that has condensed our
    #### non-independent traits. Values from our second-order meta-analysis using
    #### robu-meta are then extracted

    count_fun <- function(mod_sub) {
        return(c(mod_sub$reg_table$b.r, mod_sub$reg_table$CI.L, mod_sub$reg_table$CI.U,
            mod_sub$reg_table$SE))
    }  # estimate, lower ci, upper ci, SE

    # Extraction of values created during meta-analyses using robumeta

    robusub_RR <- model_count %>% transmute(parameter_group, estimatelnRR = map(model_lnRR,
        count_fun)) %>% mutate(r = map(estimatelnRR, ~data.frame(t(.)))) %>% unnest(r) %>%
        select(-estimatelnRR) %>% purrr::set_names(c("parameter_group", "lnRR", "lnRR_lower",
        "lnRR_upper", "lnRR_se"))

    robusub_CVR <- model_count %>% transmute(parameter_group, estimatelnCVR = map(model_lnCVR,
        count_fun)) %>% mutate(r = map(estimatelnCVR, ~data.frame(t(.)))) %>% unnest(r) %>%
        select(-estimatelnCVR) %>% purrr::set_names(c("parameter_group", "lnCVR", "lnCVR_lower",
        "lnCVR_upper", "lnCVR_se"))

    robusub_VR <- model_count %>% transmute(parameter_group, estimatelnVR = map(model_lnVR,
        count_fun)) %>% mutate(r = map(estimatelnVR, ~data.frame(t(.)))) %>% unnest(r) %>%
        select(-estimatelnVR) %>% purrr::set_names(c("parameter_group", "lnVR", "lnVR_lower",
        "lnVR_upper", "lnVR_se"))

    robu_all <- full_join(robusub_CVR, robusub_VR) %>% full_join(., robusub_RR)

    #### Combining data Merge the two data sets (the new [robu_all] and the initial
    #### [uncorrelated sub-traits with count = 1])
    meta_all <- meta1 %>% filter(par_group_size == 1) %>% as_tibble()
    # glimpse(meta_all) glimpse(robu_all)

    # Step 1: Columns are matched by name (in our case, 'parameter_group'), and any
    # missing columns will be filled with NA
    combinedmeta <- bind_rows(robu_all, meta_all)
    # glimpse(combinedmeta)

    # Steps 2&3: Add information about number of traits in a parameter group,
    # procedure, and grouping term
    metacombo <- combinedmeta
    metacombo$counts <- meta1$par_group_size[match(metacombo$parameter_group, meta1$parameter_group)]
    metacombo$procedure2 <- meta1$procedure[match(metacombo$parameter_group, meta1$parameter_group)]
    metacombo$GroupingTerm2 <- meta1$GroupingTerm[match(metacombo$parameter_group, meta1$parameter_group)]

    # Clean-up, reorder, and rename

    metacombo <- metacombo[c("parameter_group", "counts", "procedure2", "GroupingTerm2",
        "lnCVR", "lnCVR_lower", "lnCVR_upper", "lnCVR_se", "lnVR", "lnVR_lower", "lnVR_upper",
        "lnVR_se", "lnRR", "lnRR_lower", "lnRR_upper", "lnRR_se")]

    names(metacombo)[names(metacombo) == "procedure2"] <- "procedure"
    names(metacombo)[names(metacombo) == "GroupingTerm2"] <- "GroupingTerm"

## 3.3 Second-order meta-analysis for functional groups

Here, we describe the workflow depicted in Figure 1.1H. We now use our uncorrelated traits from our first order meta-analysis to conduct a secondary meta-analysis by analysing the 9 functional trait groups. We estimate an overall mean lnCVR, lnVR and lnRR for each functional trait group.

### 3.3.1 Preparing the data

Nesting, calculating the number of parameters within each grouping term, and running the meta-analyses.

    metacombo_final <- metacombo %>% group_by(GroupingTerm) %>% nest()

    # **Calculate number of parameters per grouping term

    metacombo_final <- metacombo_final %>% mutate(para_per_GroupingTerm = map_dbl(data,
        nrow))

    # For all grouping terms
    metacombo_final_all <- metacombo %>% nest(data = everything())

    # Final random effects meta-analyses within grouping terms, with SE of the
    # estimate

    overall1 <- metacombo_final %>%
    mutate(model_lnCVR = map(data, ~metafor::rma.uni(yi = .x$lnCVR, sei = (.x$lnCVR_upper -
        .x$lnCVR_lower)/(2 * 1.96), control = list(optimizer = "optim", optmethod = "Nelder-Mead",
        maxit = 1000), verbose = F)), model_lnVR = map(data, ~metafor::rma.uni(yi = .x$lnVR,
        sei = (.x$lnVR_upper - .x$lnVR_lower)/(2 * 1.96), control = list(optimizer = "optim",
            optmethod = "Nelder-Mead", maxit = 1000), verbose = F)), model_lnRR = map(data,
        ~metafor::rma.uni(yi = .x$lnRR, sei = (.x$lnRR_upper - .x$lnRR_lower)/(2 * 1.96),
            control = list(optimizer = "optim", optmethod = "Nelder-Mead", maxit = 1000),
            verbose = F)))

    # Final random effects meta-analyses ACROSS grouping terms, with SE of the
    # estimate

    overall_all1 <- metacombo_final_all %>%
    mutate(model_lnCVR = map(data, ~metafor::rma.uni(yi = .x$lnCVR, sei = (.x$lnCVR_upper -
        .x$lnCVR_lower)/(2 * 1.96), control = list(optimizer = "optim", optmethod = "Nelder-Mead",
        maxit = 1000), verbose = F)), model_lnVR = map(data, ~metafor::rma.uni(yi = .x$lnVR,
        sei = (.x$lnVR_upper - .x$lnVR_lower)/(2 * 1.96), control = list(optimizer = "optim",
            optmethod = "Nelder-Mead", maxit = 1000), verbose = F)), model_lnRR = map(data,
        ~metafor::rma.uni(yi = .x$lnRR, sei = (.x$lnRR_upper - .x$lnRR_lower)/(2 * 1.96),
            control = list(optimizer = "optim", optmethod = "Nelder-Mead", maxit = 1000),
            verbose = F)))

### 3.3.2 Re-structuring the data for each grouping term

Here we delete unused variables, and select the respective effect sizes. Please note - the referencing of the cells does NOT depend on previous ordering of the data. This would only be affected if the output structure from metafor::rma.uni changes.

    # Function will take the averall results and extract the relevant trait of
    # interest
    extract_trait <- function(data, trait) {
        tmp <- data %>% filter(., GroupingTerm == trait) %>% mutate(lnCVR = .[[4]][[1]]$b,
            lnCVR_lower = .[[4]][[1]]$ci.lb, lnCVR_upper = .[[4]][[1]]$ci.ub, lnCVR_se = .[[4]][[1]]$se,
            lnCVR_I2 = .[[4]][[1]]$I2, lnVR = .[[5]][[1]]$b, lnVR_lower = .[[5]][[1]]$ci.lb,
            lnVR_upper = .[[5]][[1]]$ci.ub, lnVR_se = .[[5]][[1]]$se, lnVR_I2 = .[[5]][[1]]$I2,
            lnRR = .[[6]][[1]]$b, lnRR_lower = .[[6]][[1]]$ci.lb, lnRR_upper = .[[6]][[1]]$ci.ub,
            lnRR_se = .[[6]][[1]]$se, lnRR_I2 = .[[6]][[1]]$I2) %>% select(., GroupingTerm,
            lnCVR:lnRR_I2)

        return(tmp)
    }

    All <- overall_all1 %>% mutate(lnCVR = .[[2]][[1]]$b, lnCVR_lower = .[[2]][[1]]$ci.lb,
        lnCVR_upper = .[[2]][[1]]$ci.ub, lnCVR_se = .[[2]][[1]]$se, lnCVR_I2 = .[[2]][[1]]$I2,
        lnVR = .[[3]][[1]]$b, lnVR_lower = .[[3]][[1]]$ci.lb, lnVR_upper = .[[3]][[1]]$ci.ub,
        lnVR_se = .[[3]][[1]]$se, lnVR_I2 = .[[3]][[1]]$I2, lnRR = .[[4]][[1]]$b, lnRR_lower = .[[4]][[1]]$ci.lb,
        lnRR_upper = .[[4]][[1]]$ci.ub, lnRR_se = .[[4]][[1]]$se, lnRR_I2 = .[[4]][[1]]$I2,
        ) %>% select(., lnCVR:lnRR_I2)

    All <- All %>% mutate(GroupingTerm = "All")

    overall2 <- bind_rows(extract_trait(overall1, "Behaviour"), extract_trait(overall1,
        "Morphology"), extract_trait(overall1, "Metabolism"), extract_trait(overall1,
        "Physiology"), extract_trait(overall1, "Immunology"), extract_trait(overall1,
        "Hematology"), extract_trait(overall1, "Heart"), extract_trait(overall1, "Hearing"),
        extract_trait(overall1, "Eye"), All)

## 3.4 Second-order meta-analysis of heterogenity

To see how much heterogeneity there is for each trait across different populations on average (for different functional groups and overall), we also conducted secondary meta-analyses by combining the total heterogeneities. More specifically, we combined the logarithm of the sum of the three variance components (Strain, Location and Unit; see above) with 1/2(N[effect size] - 1) as the sampling variance via the random effect model. That is, we conducted 3 sets of 10 secondary meta-analyses as with those of mean effect sizes for each trait, i.e. meta-analyzing log(total heterogeneities) for 9 functional groups and all the groups combined (Figure 1.1I). The analysis for heterogeneity follows the workflow of the above steps for the different meta-analyses and specifically Figure 1.1I & J. However, in the initial meta-analysis we extract sigma^2 and errors for mouse strains and centers (Institutions). As above a Loop is run and parameters extracted from metafor (sigma2’s, s.nlevels)

    ### Preparing heterogeneity

    # Create dataframe to store results
    results.allhetero.grouping <- as.data.frame(cbind(c(1:n), matrix(rep(0, n * 30),
        ncol = 30)))
    names(results.allhetero.grouping) <- c("id", "sigma2_strain.CVR", "sigma2_center.CVR",
        "sigma2_error.CVR", "s.nlevels.strain.CVR", "s.nlevels.center.CVR", "s.nlevels.error.CVR",
        "sigma2_strain.VR", "sigma2_center.VR", "sigma2_error.VR", "s.nlevels.strain.VR",
        "s.nlevels.center.VR", "s.nlevels.error.VR", "sigma2_strain.RR", "sigma2_center.RR",
        "sigma2_error.RR", "s.nlevels.strain.RR", "s.nlevels.center.RR", "s.nlevels.error.RR",
        "lnCVR", "lnCVR_lower", "lnCVR_upper", "lnCVR_se", "lnVR", "lnVR_lower", "lnVR_upper",
        "lnVR_se", "lnRR", "lnRR_lower", "lnRR_upper", "lnRR_se")

    # Loop

    for (t in 1:n) {
        tryCatch({
            results <- data %>% data_subset_parameterid_individual_by_age(t) %>% calculate_population_stats() %>%
                create_meta_analysis_effect_sizes() %>% mutate(err = seq_len(n()))

            # lnCVR, logaritm of the ratio of male and female coefficients of variance

            cvr. <- metafor::rma.mv(yi = effect_size_CVR, V = sample_variance_CVR, random = list(~1 |
                strain_name, ~1 | production_center, ~1 | err), control = list(optimizer = "optim",
                optmethod = "Nelder-Mead", maxit = 1000), data = results)

            results.allhetero.grouping[t, 2] <- cvr.$sigma2[1]
            results.allhetero.grouping[t, 3] <- cvr.$sigma2[2]
            results.allhetero.grouping[t, 4] <- cvr.$sigma2[3]
            results.allhetero.grouping[t, 5] <- cvr.$s.nlevels[1]
            results.allhetero.grouping[t, 6] <- cvr.$s.nlevels[2]
            results.allhetero.grouping[t, 7] <- cvr.$s.nlevels[3]
            results.allhetero.grouping[t, 20] <- cvr.$b
            results.allhetero.grouping[t, 21] <- cvr.$ci.lb
            results.allhetero.grouping[t, 22] <- cvr.$ci.ub
            results.allhetero.grouping[t, 23] <- cvr.$se

            # lnVR, male to female variability ratio (logarithm of male and female standard
            # deviations)

            vr. <- metafor::rma.mv(yi = effect_size_VR, V = sample_variance_VR, random = list(~1 |
                strain_name, ~1 | production_center, ~1 | err), control = list(optimizer = "optim",
                optmethod = "Nelder-Mead", maxit = 1000), data = results)

            results.allhetero.grouping[t, 8] <- vr.$sigma2[1]
            results.allhetero.grouping[t, 9] <- vr.$sigma2[2]
            results.allhetero.grouping[t, 10] <- vr.$sigma2[3]
            results.allhetero.grouping[t, 11] <- vr.$s.nlevels[1]
            results.allhetero.grouping[t, 12] <- vr.$s.nlevels[2]
            results.allhetero.grouping[t, 13] <- vr.$s.nlevels[3]
            results.allhetero.grouping[t, 24] <- vr.$b
            results.allhetero.grouping[t, 25] <- vr.$ci.lb
            results.allhetero.grouping[t, 26] <- vr.$ci.ub
            results.allhetero.grouping[t, 27] <- vr.$se

            # lnRR, response ratio (logarithm of male and female means)

            rr. <- metafor::rma.mv(yi = effect_size_ROM, V = sample_variance_ROM, random = list(~1 |
                strain_name, ~1 | production_center, ~1 | err), control = list(optimizer = "optim",
                optmethod = "Nelder-Mead", maxit = 1000), data = results)

            results.allhetero.grouping[t, 14] <- rr.$sigma2[1]
            results.allhetero.grouping[t, 15] <- rr.$sigma2[2]
            results.allhetero.grouping[t, 16] <- rr.$sigma2[3]
            results.allhetero.grouping[t, 17] <- rr.$s.nlevels[1]
            results.allhetero.grouping[t, 18] <- rr.$s.nlevels[2]
            results.allhetero.grouping[t, 19] <- rr.$s.nlevels[3]
            results.allhetero.grouping[t, 28] <- rr.$b
            results.allhetero.grouping[t, 29] <- rr.$ci.lb
            results.allhetero.grouping[t, 30] <- rr.$ci.ub
            results.allhetero.grouping[t, 31] <- rr.$se
        }, error = function(e) {
            cat("ERROR :", conditionMessage(e), "\n")
        })
    }

    ## ERROR : Optimizer (optim) did not achieve convergence (convergence = 10).
    ## ERROR : Processing terminated since k <= 1.
    ## ERROR : Processing terminated since k <= 1.

Here we exclude traits without variation between mouse strains and merge data sets containing metafor results with procedure etc. names. Dealing with the correlated parameters, and running the second -order meta-analysis for heterogeneity data (Figure 1.1I)

                  results.allhetero.grouping2 <- results.allhetero.grouping[results.allhetero.grouping$s.nlevels.strain.VR != 0, ]
    # nrow(results.allhetero.grouping) #232

    # Merging dataset with correct procedure names
    # procedures <- read.csv(here("export", "procedures.csv"))

    results.allhetero.grouping2$parameter_group <- data$parameter_group[match(results.allhetero.grouping2$id, data$id)]
          results.allhetero.grouping2$procedure <- data$procedure_name[match(results.allhetero.grouping2$id, data$id)]

       results.allhetero.grouping2$GroupingTerm <- procedures$GroupingTerm[match(results.allhetero.grouping2$procedure, procedures$procedure)]
     results.allhetero.grouping2$parameter_name <- data$parameter_name[match(results.allhetero.grouping2$id, data$id)]

    #We deal with correlated parameters following the procedures described above (extraction of effects sizes).

    metahetero1 <- results.allhetero.grouping2
    # length(unique(metahetero1$procedure)) #19
    # length(unique(metahetero1$GroupingTerm)) #9
    # length(unique(metahetero1$parameter_group)) #152
    # length(unique(metahetero1$parameter_name)) #223

    # Count of number of parameter names (correlated sub-traits) in each parameter group (par_group_size)

    metahetero1b <- metahetero1                    %>%
                    group_by(parameter_group)      %>%
                    mutate(par_group_size = n_distinct(parameter_name))

    metahetero1$par_group_size <- metahetero1b$par_group_size[match(metahetero1$parameter_group, metahetero1b$parameter_group)]

    # Create subsets with > 1 count (par_group_size > 1)

    metahetero1_sub <- subset(metahetero1, par_group_size > 1) # 92 observations

    # Nest data

    n_count. <- metahetero1_sub %>%
      group_by(parameter_group) %>%
      nest()

    # meta-analysis preparation

    model_count. <- n_count. %>%
      mutate(
        model_lnRR = map(data, ~ robu(.x$lnRR ~ 1,
          data = .x, studynum = .x$id, modelweights = c("CORR"), rho = 0.8,
          small = TRUE, var.eff.size = (.x$lnRR_se)^2
        )),
        model_lnVR = map(data, ~ robu(.x$lnVR ~ 1,
          data = .x, studynum = .x$id, modelweights = c("CORR"), rho = 0.8,
          small = TRUE, var.eff.size = (.x$lnVR_se)^2
        )),
        model_lnCVR = map(data, ~ robu(.x$lnCVR ~ 1,
          data = .x, studynum = .x$id, modelweights = c("CORR"), rho = 0.8,
          small = TRUE, var.eff.size = (.x$lnCVR_se)^2
        ))
      )

    # Robumeta object details:  str(model_count.$model_lnCVR[[1]])

    # Perform meta-analyses on correlated sub-traits, using robumeta
    # Extract and save parameter estimates

    count_fun. <- function(mod_sub) {
      return(c(as.numeric(mod_sub$mod_info$term1), mod_sub$N))
    }

    robusub_RR. <- model_count. %>%
      transmute(estimatelnRR = map(model_lnRR, count_fun.)) %>%
      mutate(r = map(estimatelnRR, ~ data.frame(t(.)))) %>%
      unnest(r) %>%
      select(-estimatelnRR) %>%
      purrr::set_names(c("parameter_group", "var.RR", "N.RR"))

    robusub_CVR. <- model_count. %>%
      transmute(estimatelnCVR = map(model_lnCVR, count_fun.)) %>%
      mutate(r = map(estimatelnCVR, ~ data.frame(t(.)))) %>%
      unnest(r) %>%
      select(-estimatelnCVR) %>%
      purrr::set_names(c("parameter_group", "var.CVR", "N.CVR"))

    robusub_VR. <- model_count. %>%
      transmute(estimatelnVR = map(model_lnVR, count_fun.)) %>%
      mutate(r = map(estimatelnVR, ~ data.frame(t(.)))) %>%
      unnest(r) %>%
      select(-estimatelnVR) %>%
      purrr::set_names(c("parameter_group", "var.VR", "N.VR"))

    robu_all. <- full_join(robusub_CVR., robusub_VR.) %>% full_join(., robusub_RR.)

    #Merge the two data sets (the new [robu_all.] and the initial [uncorrelated sub-traits with count = 1])

    #In this step, we   1) merge the N from robumeta and the  N from metafor (s.nlevels.error) together into the same columns (N.RR, N.VR, N.CVR)

    # 2) calculate the total variance for metafor models as the sum of random effect variances and the residual error, then add in the same columns together with the residual variances from robumeta

    metahetero_all <- metahetero1 %>%
      filter(par_group_size == 1) %>%
      as_tibble()
    metahetero_all$N.RR <- metahetero_all$s.nlevels.error.RR
    metahetero_all$N.CVR <- metahetero_all$s.nlevels.error.CVR
    metahetero_all$N.VR <- metahetero_all$s.nlevels.error.VR
    metahetero_all$var.RR <- log(sqrt(metahetero_all$sigma2_strain.RR + metahetero_all$sigma2_center.RR + metahetero_all$sigma2_error.RR))
    metahetero_all$var.VR <- log(sqrt(metahetero_all$sigma2_strain.VR + metahetero_all$sigma2_center.VR + metahetero_all$sigma2_error.VR))
    metahetero_all$var.CVR <- log(sqrt(metahetero_all$sigma2_strain.CVR + metahetero_all$sigma2_center.CVR + metahetero_all$sigma2_error.CVR))
    # str(metahetero_all)
    # str(robu_all.)

    metahetero_all <- metahetero_all %>% mutate(
      var.RR = if_else(var.RR == -Inf, -7, var.RR),   # numbers chosen as limits; based on the values in the dataset
      var.VR = if_else(var.VR == -Inf, -5, var.VR),
      var.CVR = if_else(var.CVR == -Inf, -6, var.CVR)
    )

    # Combine data
    # Step1

    combinedmetahetero <- bind_rows(robu_all., metahetero_all)
    # glimpse(combinedmetahetero)

    # Steps 2&3

                  metacombohetero <- combinedmetahetero
           metacombohetero$counts <- metahetero1$par_group_size[match(metacombohetero$parameter_group, metahetero1$parameter_group)]
       metacombohetero$procedure2 <- metahetero1$procedure[match(metacombohetero$parameter_group, metahetero1$parameter_group)]
    metacombohetero$GroupingTerm2 <- metahetero1$GroupingTerm[match(metacombohetero$parameter_group, metahetero1$parameter_group)]

    # **Clean-up and rename

    metacombohetero <- metacombohetero %>% select(parameter_group, var.CVR, N.CVR, var.VR, N.VR, var.RR, N.RR, counts, procedure = procedure2, GroupingTerm = GroupingTerm2)

    # Invidual table for heterogeneity across all traits
     kable(metacombohetero,, digits = 3) %>%
      kable_styling() %>%
      scroll_box(width = "100%", height = "200px")

parameter_group | var.CVR | N.CVR | var.VR | N.VR | var.RR | N.RR | counts | procedure | GroupingTerm
--- | --- | --- | --- | --- | --- | --- | --- | --- | ---
pre-pulse inhibition | 0.004 | 5 | 0.000 | 5 | 0.000 | 5 | 5 | Acoustic Startle and Pre-pulse Inhibition (PPI) | Behaviour
B cells | 0.000 | 4 | 0.000 | 4 | 0.002 | 4 | 4 | Immunophenotyping | Immunology
cd4 nkt | -0.014 | 6 | -0.012 | 6 | 0.018 | 6 | 6 | Immunophenotyping | Immunology
cd4 t | 0.005 | 7 | 0.003 | 7 | 0.001 | 7 | 7 | Immunophenotyping | Immunology
cd8 nkt | -0.008 | 6 | 0.003 | 6 | -0.003 | 6 | 6 | Immunophenotyping | Immunology
cd8 t | 0.004 | 7 | 0.001 | 7 | 0.000 | 7 | 7 | Immunophenotyping | Immunology
cdcs | -0.002 | 2 | -0.009 | 2 | -0.002 | 2 | 2 | Immunophenotyping | Immunology
dn nkt | -0.001 | 6 | 0.001 | 6 | 0.006 | 6 | 6 | Immunophenotyping | Immunology
dn t | 0.005 | 7 | 0.000 | 7 | 0.000 | 7 | 7 | Immunophenotyping | Immunology
eosinophils | 0.000 | 3 | 0.013 | 3 | 0.005 | 3 | 3 | Hematology | Hematology
follicular b cells | -0.002 | 2 | -0.007 | 2 | 0.000 | 2 | 2 | Immunophenotyping | Immunology
luc | 0.000 | 2 | 0.021 | 2 | 0.027 | 2 | 2 | Hematology | Hematology
lymphocytes | 0.059 | 2 | 0.008 | 2 | 0.013 | 2 | 2 | Hematology | Hematology
monocytes | 0.002 | 3 | 0.006 | 3 | 0.009 | 3 | 3 | Hematology | Hematology
neutrophils | -0.010 | 3 | 0.030 | 3 | 0.018 | 3 | 3 | Hematology | Hematology
nk cells | 0.000 | 6 | 0.002 | 6 | 0.002 | 6 | 6 | Immunophenotyping | Immunology
nkt cells | 0.000 | 4 | -0.004 | 4 | -0.001 | 4 | 4 | Immunophenotyping | Immunology
percentage of live gated events | -0.007 | 2 | -0.004 | 2 | -0.001 | 2 | 2 | Immunophenotyping | Immunology
response amplitude | 0.003 | 10 | 0.004 | 10 | 0.013 | 10 | 10 | Acoustic Startle and Pre-pulse Inhibition (PPI) | Behaviour
t cells | 0.001 | 3 | 0.011 | 3 | 0.002 | 3 | 3 | Immunophenotyping | Immunology
number events | 0.001 | 2 | 0.000 | 2 | -0.001 | 2 | 2 | Immunophenotyping | Immunology
12khz-evoked abr threshold | -2.660 | 14 | -2.144 | 14 | -3.788 | 14 | 1 | Auditory Brain Stem Response | Hearing
18khz-evoked abr threshold | -2.711 | 14 | -2.193 | 14 | -3.514 | 14 | 1 | Auditory Brain Stem Response | Hearing
24khz-evoked abr threshold | -2.467 | 14 | -1.182 | 14 | -3.555 | 14 | 1 | Auditory Brain Stem Response | Hearing
30khz-evoked abr threshold | -2.419 | 14 | -2.456 | 14 | -3.401 | 14 | 1 | Auditory Brain Stem Response | Hearing
6khz-evoked abr threshold | -6.000 | 14 | -4.570 | 14 | -4.675 | 14 | 1 | Auditory Brain Stem Response | Hearing
alanine aminotransferase | -1.677 | 16 | -1.116 | 16 | -2.143 | 16 | 1 | Clinical Chemistry | Physiology
albumin | -2.221 | 16 | -2.302 | 16 | -3.575 | 16 | 1 | Clinical Chemistry | Physiology
alkaline phosphatase | -2.480 | 16 | -2.125 | 16 | -2.553 | 16 | 1 | Clinical Chemistry | Physiology
alpha-amylase | -2.316 | 9 | -1.846 | 9 | -3.002 | 9 | 1 | Clinical Chemistry | Physiology
area under glucose response curve | -2.146 | 16 | -2.088 | 16 | -2.192 | 16 | 1 | Intraperitoneal glucose tolerance test (IPGTT) | Metabolism
aspartate aminotransferase | -1.487 | 16 | -1.111 | 16 | -2.123 | 16 | 1 | Clinical Chemistry | Physiology
basophil cell count | -2.440 | 8 | -1.278 | 8 | -1.527 | 8 | 1 | Hematology | Hematology
basophil differential count | -2.659 | 8 | -1.253 | 8 | -2.318 | 8 | 1 | Hematology | Hematology
bmc/body weight | -1.871 | 17 | -2.057 | 17 | -2.580 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
body length | -2.749 | 16 | -2.780 | 16 | -4.854 | 16 | 1 | Body Composition (DEXA lean/fat) | Morphology
body temp | -6.000 | 3 | -5.000 | 3 | -6.576 | 3 | 1 | Echo | Heart
body weight | -2.210 | 18 | -2.251 | 18 | -3.638 | 18 | 1 | Body Weight | Morphology
body weight after experiment | -2.837 | 9 | -2.949 | 9 | -3.828 | 9 | 1 | Indirect Calorimetry | Metabolism
body weight before experiment | -2.597 | 9 | -2.688 | 9 | -3.677 | 9 | 1 | Indirect Calorimetry | Metabolism
bone area | -2.037 | 17 | -1.964 | 17 | -3.014 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
bone mineral content (excluding skull) | -1.799 | 17 | -1.777 | 17 | -2.565 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
bone mineral density (excluding skull) | -1.324 | 17 | -1.230 | 17 | -3.506 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
calcium | -2.447 | 16 | -2.448 | 16 | -5.041 | 16 | 1 | Clinical Chemistry | Physiology
cardiac output | -2.671 | 5 | -2.847 | 5 | -3.575 | 5 | 1 | Echo | Heart
center average speed | -2.564 | 12 | -2.981 | 12 | -2.753 | 12 | 1 | Open Field | Behaviour
center distance travelled | -2.777 | 12 | -1.859 | 12 | -1.859 | 12 | 1 | Open Field | Behaviour
center permanence time | -2.852 | 13 | -2.078 | 13 | -1.972 | 13 | 1 | Open Field | Behaviour
center resting time | -2.418 | 10 | -1.752 | 10 | -1.576 | 10 | 1 | Open Field | Behaviour
chloride | -1.617 | 10 | -1.565 | 10 | -5.111 | 10 | 1 | Clinical Chemistry | Physiology
click-evoked abr threshold | -2.317 | 11 | -1.793 | 11 | -2.889 | 11 | 1 | Auditory Brain Stem Response | Hearing
creatine kinase | -2.138 | 9 | -1.168 | 9 | -1.340 | 9 | 1 | Clinical Chemistry | Physiology
creatinine | -2.358 | 16 | -0.614 | 16 | -2.654 | 16 | 1 | Clinical Chemistry | Physiology
cv | -2.160 | 9 | -1.615 | 9 | -2.083 | 9 | 1 | Electrocardiogram (ECG) | Heart
distance travelled - total | -2.312 | 13 | -2.212 | 13 | -2.242 | 13 | 1 | Open Field | Behaviour
ejection fraction | -2.477 | 6 | -2.588 | 6 | -4.039 | 6 | 1 | Echo | Heart
end-diastolic diameter | -6.000 | 3 | -3.219 | 3 | -4.076 | 3 | 1 | Echo | Heart
end-systolic diameter | -6.000 | 3 | -5.000 | 3 | -3.875 | 3 | 1 | Echo | Heart
fasted blood glucose concentration | -1.978 | 16 | -2.041 | 16 | -2.915 | 16 | 1 | Intraperitoneal glucose tolerance test (IPGTT) | Metabolism
fat mass | -2.007 | 17 | -1.755 | 17 | -2.267 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
fat/body weight | -1.918 | 17 | -1.806 | 17 | -2.283 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
forelimb and hindlimb grip strength measurement mean | -2.365 | 16 | -2.197 | 16 | -3.135 | 16 | 1 | Grip Strength | Morphology
forelimb grip strength measurement mean | -2.601 | 16 | -2.573 | 16 | -3.234 | 16 | 1 | Grip Strength | Morphology
fractional shortening | -2.490 | 7 | -2.520 | 7 | -4.428 | 7 | 1 | Echo | Heart
free fatty acids | -2.245 | 5 | -2.204 | 5 | -4.054 | 5 | 1 | Clinical Chemistry | Physiology
fructosamine | -2.828 | 6 | -2.795 | 6 | -3.336 | 6 | 1 | Clinical Chemistry | Physiology
glucose | -2.546 | 16 | -1.963 | 16 | -2.641 | 16 | 1 | Clinical Chemistry | Physiology
hdl-cholesterol | -2.294 | 15 | -2.035 | 15 | -2.612 | 15 | 1 | Clinical Chemistry | Physiology
heart weight | -1.978 | 15 | -1.661 | 15 | -3.140 | 15 | 1 | Heart Weight | Morphology
heart weight normalised against body weight | -2.199 | 14 | -1.782 | 14 | -3.132 | 14 | 1 | Heart Weight | Morphology
hematocrit | -1.578 | 17 | -1.597 | 17 | -3.881 | 17 | 1 | Hematology | Hematology
hemoglobin | -2.449 | 17 | -2.389 | 17 | -3.875 | 17 | 1 | Hematology | Hematology
hr | -2.090 | 12 | -1.834 | 12 | -2.476 | 12 | 1 | Electrocardiogram (ECG) | Heart
hrv | -3.097 | 7 | -1.989 | 7 | -2.003 | 7 | 1 | Electrocardiogram (ECG) | Heart
initial response to glucose challenge | -2.842 | 16 | -3.311 | 16 | -2.845 | 16 | 1 | Intraperitoneal glucose tolerance test (IPGTT) | Metabolism
insulin | -1.138 | 8 | -0.657 | 8 | -1.041 | 8 | 1 | Insulin Blood Level | Metabolism
iron | -1.934 | 11 | -1.739 | 11 | -2.948 | 11 | 1 | Clinical Chemistry | Physiology
lactate dehydrogenase | -2.606 | 4 | -1.765 | 4 | -1.989 | 4 | 1 | Clinical Chemistry | Physiology
latency to center entry | -2.294 | 10 | -1.555 | 10 | -1.435 | 10 | 1 | Open Field | Behaviour
ldl-cholesterol | -1.447 | 7 | -1.147 | 7 | -0.946 | 7 | 1 | Clinical Chemistry | Physiology
lean mass | -2.219 | 17 | -2.145 | 17 | -3.418 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
lean/body weight | -1.719 | 17 | -1.798 | 17 | -3.767 | 17 | 1 | Body Composition (DEXA lean/fat) | Morphology
left anterior chamber depth | -1.863 | 2 | -1.852 | 2 | -7.000 | 2 | 1 | Eye Morphology | Eye
left corneal thickness | -6.000 | 2 | -5.000 | 2 | -7.000 | 2 | 1 | Eye Morphology | Eye
left inner nuclear layer | -6.000 | 2 | -5.000 | 2 | -7.000 | 2 | 1 | Eye Morphology | Eye
left outer nuclear layer | -6.000 | 2 | -5.000 | 2 | -7.000 | 2 | 1 | Eye Morphology | Eye
left posterior chamber depth | -6.000 | 2 | -5.000 | 2 | -4.923 | 2 | 1 | Eye Morphology | Eye
left total retinal thickness | -1.607 | 3 | -1.641 | 3 | -5.300 | 3 | 1 | Eye Morphology | Eye
locomotor activity | -2.301 | 13 | -2.882 | 13 | -2.361 | 13 | 1 | Combined SHIRPA and Dysmorphology | Behaviour
lvawd | -6.000 | 5 | -5.000 | 5 | -4.648 | 5 | 1 | Echo | Heart
lvaws | -1.627 | 4 | -1.885 | 4 | -3.427 | 4 | 1 | Echo | Heart
lvidd | -2.763 | 7 | -2.609 | 7 | -4.135 | 7 | 1 | Echo | Heart
lvids | -2.137 | 7 | -2.075 | 7 | -4.024 | 7 | 1 | Echo | Heart
lvpwd | -2.578 | 7 | -2.307 | 7 | -4.028 | 7 | 1 | Echo | Heart
lvpws | -2.718 | 6 | -2.563 | 6 | -3.982 | 6 | 1 | Echo | Heart
magnesium | -6.000 | 6 | -3.327 | 6 | -2.767 | 6 | 1 | Urinalysis | Physiology
mean cell hemoglobin concentration | -1.404 | 17 | -1.334 | 17 | -5.254 | 17 | 1 | Hematology | Hematology
mean cell volume | -1.702 | 17 | -1.737 | 17 | -5.137 | 17 | 1 | Hematology | Hematology
mean corpuscular hemoglobin | -2.334 | 17 | -2.325 | 17 | -5.730 | 17 | 1 | Hematology | Hematology
mean platelet volume | -2.727 | 13 | -2.650 | 13 | -4.354 | 13 | 1 | Hematology | Hematology
mean r amplitude | -6.000 | 8 | -2.751 | 8 | -2.519 | 8 | 1 | Electrocardiogram (ECG) | Heart
mean sr amplitude | -4.105 | 7 | -3.587 | 7 | -3.330 | 7 | 1 | Electrocardiogram (ECG) | Heart
mzb (cd21/35 high) | -2.182 | 4 | -1.822 | 4 | -7.000 | 4 | 1 | Immunophenotyping | Immunology
number of caudal vertebrae | -1.333 | 7 | -1.331 | 7 | -16.871 | 11 | 1 | X-ray | Morphology
number of center entries | -2.587 | 10 | -2.550 | 10 | -1.942 | 10 | 1 | Open Field | Behaviour
number of digits | -1.530 | 3 | -1.530 | 3 | -16.871 | 15 | 1 | X-ray | Morphology
number of lumbar vertebrae | -0.951 | 3 | -0.951 | 3 | -16.871 | 14 | 1 | X-ray | Morphology
number of pelvic vertebrae | -1.655 | 3 | -1.655 | 3 | -16.871 | 14 | 1 | X-ray | Morphology
number of rears - total | -2.017 | 8 | -1.274 | 8 | -1.892 | 8 | 1 | Open Field | Behaviour
number of ribs left | -1.234 | 3 | -1.233 | 3 | -16.871 | 15 | 1 | X-ray | Morphology
number of ribs right | -2.400 | 2 | -2.402 | 2 | -16.871 | 15 | 1 | X-ray | Morphology
number of signals | -3.327 | 9 | -3.870 | 9 | -3.407 | 11 | 1 | Electrocardiogram (ECG) | Heart
others | -3.208 | 6 | -2.931 | 6 | -5.111 | 6 | 1 | Immunophenotyping | Immunology
pdcs | -1.637 | 5 | -0.679 | 5 | -1.979 | 5 | 1 | Immunophenotyping | Immunology
percentage center time | -2.726 | 13 | -2.123 | 13 | -1.931 | 13 | 1 | Open Field | Behaviour
periphery average speed | -2.423 | 12 | -2.283 | 12 | -2.614 | 12 | 1 | Open Field | Behaviour
periphery distance travelled | -2.477 | 12 | -2.624 | 12 | -2.273 | 12 | 1 | Open Field | Behaviour
periphery permanence time | -1.889 | 13 | -2.136 | 13 | -3.343 | 13 | 1 | Open Field | Behaviour
periphery resting time | -2.418 | 10 | -2.895 | 10 | -2.544 | 10 | 1 | Open Field | Behaviour
phosphorus | -3.180 | 16 | -2.391 | 16 | -2.927 | 16 | 1 | Clinical Chemistry | Physiology
platelet count | -2.517 | 17 | -2.448 | 17 | -3.239 | 17 | 1 | Hematology | Hematology
pnn5(6>ms) | -6.000 | 6 | -1.355 | 6 | -1.160 | 6 | 1 | Electrocardiogram (ECG) | Heart
potassium | -1.788 | 10 | -1.685 | 10 | -3.393 | 10 | 1 | Clinical Chemistry | Physiology
pq | -2.365 | 7 | -2.799 | 7 | -3.489 | 7 | 1 | Electrocardiogram (ECG) | Heart
pr | -2.816 | 11 | -2.779 | 11 | -3.976 | 11 | 1 | Electrocardiogram (ECG) | Heart
qrs | -3.265 | 11 | -3.212 | 11 | -4.590 | 11 | 1 | Electrocardiogram (ECG) | Heart
qtc | -3.176 | 10 | -2.941 | 10 | -4.831 | 10 | 1 | Electrocardiogram (ECG) | Heart
qtc dispersion | -2.997 | 7 | -2.265 | 7 | -3.126 | 7 | 1 | Electrocardiogram (ECG) | Heart
red blood cell count | -2.269 | 17 | -2.315 | 17 | -3.795 | 17 | 1 | Hematology | Hematology
red blood cell distribution width | -1.512 | 13 | -1.425 | 13 | -3.810 | 13 | 1 | Hematology | Hematology
respiration rate | -2.972 | 4 | -2.405 | 4 | -3.468 | 4 | 1 | Echo | Heart
respiratory exchange ratio | -2.339 | 9 | -2.359 | 9 | -4.825 | 9 | 1 | Indirect Calorimetry | Metabolism
right anterior chamber depth | -0.465 | 2 | -0.469 | 2 | -7.000 | 2 | 1 | Eye Morphology | Eye
right corneal thickness | -2.191 | 2 | -2.409 | 2 | -4.423 | 2 | 1 | Eye Morphology | Eye
right inner nuclear layer | -1.035 | 2 | -0.935 | 2 | -3.449 | 2 | 1 | Eye Morphology | Eye
right outer nuclear layer | -6.000 | 2 | -5.000 | 2 | -7.000 | 2 | 1 | Eye Morphology | Eye
right posterior chamber depth | -6.000 | 2 | -5.000 | 2 | -4.092 | 2 | 1 | Eye Morphology | Eye
right total retinal thickness | -0.948 | 3 | -0.974 | 3 | -4.789 | 3 | 1 | Eye Morphology | Eye
rmssd | -1.249 | 7 | -0.802 | 7 | -1.936 | 7 | 1 | Electrocardiogram (ECG) | Heart
rp macrophage (cd19- cd11c-) | -1.316 | 6 | -1.285 | 6 | -2.224 | 6 | 1 | Immunophenotyping | Immunology
rr | -2.033 | 11 | -1.988 | 11 | -4.619 | 11 | 1 | Electrocardiogram (ECG) | Heart
sodium | -1.627 | 10 | -1.544 | 10 | -5.603 | 10 | 1 | Clinical Chemistry | Physiology
spleen weight | -1.306 | 10 | -1.020 | 10 | -2.743 | 10 | 1 | Immunophenotyping | Immunology
st | -2.830 | 11 | -2.548 | 11 | -4.127 | 11 | 1 | Electrocardiogram (ECG) | Heart
stroke volume | -2.179 | 5 | -2.091 | 5 | -4.875 | 5 | 1 | Echo | Heart
tibia length | -0.830 | 12 | -0.841 | 12 | -5.345 | 12 | 1 | Heart Weight | Morphology
total bilirubin | -2.143 | 16 | -1.809 | 16 | -2.661 | 16 | 1 | Clinical Chemistry | Physiology
total cholesterol | -1.302 | 16 | -1.110 | 16 | -2.930 | 16 | 1 | Clinical Chemistry | Physiology
total food intake | -1.938 | 6 | -1.755 | 6 | -3.002 | 6 | 1 | Indirect Calorimetry | Metabolism
total protein | -6.000 | 16 | -3.618 | 16 | -4.016 | 16 | 1 | Clinical Chemistry | Physiology
total water intake | -2.843 | 3 | -5.000 | 3 | -2.782 | 3 | 1 | Indirect Calorimetry | Metabolism
triglycerides | -2.116 | 16 | -1.701 | 16 | -1.947 | 16 | 1 | Clinical Chemistry | Physiology
urea (blood urea nitrogen - bun) | -1.658 | 16 | -1.442 | 16 | -2.669 | 16 | 1 | Clinical Chemistry | Physiology
uric acid | -6.000 | 3 | -1.838 | 3 | -1.085 | 3 | 1 | Clinical Chemistry | Physiology
white blood cell count | -2.181 | 16 | -1.698 | 16 | -2.411 | 16 | 1 | Hematology | Hematology
whole arena average speed | -2.253 | 13 | -2.232 | 13 | -2.518 | 13 | 1 | Open Field | Behaviour
whole arena permanence | -1.055 | 2 | -1.056 | 2 | -16.871 | 13 | 1 | Open Field | Behaviour
whole arena resting time | -2.730 | 13 | -2.866 | 13 | -2.432 | 13 | 1 | Open Field | Behaviour

    #### Meta-analysis of heterogeneity

    ## Perform meta-meta-analysis (3 for each of the 9 grouping terms: var.CVR, var.VR, var.RR)

    metacombohetero_final <- metacombohetero %>%
      group_by(GroupingTerm) %>%
      nest()

    # Final random effects meta-analyses within grouping terms, with SE of the estimate

    heterog1 <- metacombohetero_final %>%

      mutate(
        model_heteroCVR = map(data, ~ metafor::rma.uni(
          yi = .x$var.CVR, sei = sqrt(1 / 2 * (.x$N.CVR - 1)),
          control = list(optimizer = "optim",
                         optmethod = "Nelder-Mead",
                         maxit = 10000,
                         stepadj = 0.5), verbose = F)),
        model_heteroVR = map(data, ~ metafor::rma.uni(
          yi = .x$var.VR, sei = sqrt(1 / 2 * (.x$N.VR - 1)),
          control = list(optimizer = "optim",
                         optmethod = "Nelder-Mead",
                         maxit = 10000,
                         stepadj = 0.5), verbose = F)),
        model_heteroRR = map(data, ~ metafor::rma.uni(
          yi = .x$var.RR, sei = sqrt(1 / 2 * (.x$N.RR - 1)),
          control = list(optimizer = "optim",
                         optmethod = "Nelder-Mead",
                         maxit = 10000,
                         stepadj = 0.5), verbose = F)) )

    # Across all grouping terms

    metacombohetero_all_final <- metacombohetero %>%
      nest(data = everything())

    # Final random effects meta-analyses ACROSS grouping terms, with SE of the estimate

    heterog1_all <- metacombohetero_all_final %>%

      mutate( model_heteroCVR = map(data, ~ metafor::rma.uni(
          yi = .x$var.CVR, sei = sqrt(1 / 2 * (.x$N.CVR - 1)),
          control = list(optimizer = "optim",
                         optmethod = "Nelder-Mead",
                         maxit = 10000,
                         stepadj = 0.5), verbose = F)),
        model_heteroVR = map(data, ~ metafor::rma.uni(yi = .x$var.VR, sei = sqrt(1 / 2 * (.x$N.VR - 1)),
          control = list(optimizer = "optim",
                         optmethod = "Nelder-Mead",
                         maxit = 10000,
                         stepadj = 0.5), verbose = F)),
        model_heteroRR = map(data, ~ metafor::rma.uni(yi = .x$var.RR, sei = sqrt(1 / 2 * (.x$N.RR - 1)),
          control = list(optimizer = "optim",
                         optmethod = "Nelder-Mead",
                         maxit = 10000,
                         stepadj = 0.5), verbose = F)))

    # Re-structure data for each grouping term; extract heterogenenity/variance terms; delete un-used variables

    extract_heterogeneity <- function(data, type){

      tmp <- data %>%
             filter(., GroupingTerm == type) %>%
             select(., -data) %>%
             mutate(heteroCVR       = .[[2]][[1]]$b,
                    heteroCVR_lower = .[[2]][[1]]$ci.lb,
                    heteroCVR_upper = .[[2]][[1]]$ci.ub,
                    heteroCVR_se    = .[[2]][[1]]$se,
                    heteroVR        = .[[3]][[1]]$b,
                    heteroVR_lower  = .[[3]][[1]]$ci.lb,
                    heteroVR_upper  = .[[3]][[1]]$ci.ub,
                    heteroVR_se     = .[[3]][[1]]$se,
                    heteroRR        = .[[4]][[1]]$b,
                    heteroRR_lower  = .[[4]][[1]]$ci.lb,
                    heteroRR_upper  = .[[4]][[1]]$ci.ub,
                    heteroRR_se     = .[[4]][[1]]$se) %>%
            select(., GroupingTerm, heteroCVR:heteroRR_se)

      return(tmp)
    }

       Behaviour. <- extract_heterogeneity(heterog1, type = "Behaviour")
      Immunology. <- extract_heterogeneity(heterog1, type = "Immunology")
      Hematology. <- extract_heterogeneity(heterog1, type = "Hematology")
         Hearing. <- extract_heterogeneity(heterog1, type = "Hearing")
      Physiology. <- extract_heterogeneity(heterog1, type = "Physiology")
      Metabolism. <- extract_heterogeneity(heterog1, type = "Metabolism")
      Morphology. <- extract_heterogeneity(heterog1, type = "Morphology")
           Heart. <- extract_heterogeneity(heterog1, type = "Heart")
             Eye. <- extract_heterogeneity(heterog1, type = "Eye")

    #Reorder to be able to keep cell referencing
    heterog1_all <- heterog1_all %>%
                    mutate(GroupingTerm = "All") %>%
                    select(GroupingTerm, everything())

    All. <- heterog1_all %>%
            select(., -data) %>%
            mutate(heteroCVR       = .[[2]][[1]]$b,
                   heteroCVR_lower = .[[2]][[1]]$ci.lb,
                   heteroCVR_upper = .[[2]][[1]]$ci.ub,
                   heteroCVR_se    = .[[2]][[1]]$se,
                   heteroVR        = .[[3]][[1]]$b,
                   heteroVR_lower  = .[[3]][[1]]$ci.lb,
                   heteroVR_upper  = .[[3]][[1]]$ci.ub,
                   heteroVR_se     = .[[3]][[1]]$se,
                   heteroRR        = .[[4]][[1]]$b,
                   heteroRR_lower  = .[[4]][[1]]$ci.lb,
                   heteroRR_upper  = .[[4]][[1]]$ci.ub,
                   heteroRR_se     = .[[4]][[1]]$se) %>%
            select(., GroupingTerm, heteroCVR:heteroRR_se)

    heterog2 <- bind_rows(Behaviour., Morphology., Metabolism., Physiology., Immunology., Hematology., Heart., Hearing., Eye., All.)

    #### Preparation of the Heterogeneity PLOT

    #Restructure data for plotting

    heterog3 <- gather(heterog2, parameter, value, c(heteroCVR, heteroVR, heteroRR), factor_key = TRUE)

    heteroCVR.ci <- heterog3 %>%
                    filter(parameter == "heteroCVR") %>%
                    mutate(ci.low = heteroCVR_lower,
                           ci.high = heteroCVR_upper)
    heteroVR.ci <-  heterog3 %>%
                    filter(parameter == "heteroVR") %>%
                    mutate(ci.low = heteroVR_lower,
                           ci.high = heteroVR_upper)
    heteroRR.ci <-  heterog3 %>%
                    filter(parameter == "heteroRR") %>%
                    mutate(ci.low = heteroRR_lower,
                           ci.high = heteroRR_upper)

    heterog4 <- bind_rows(heteroCVR.ci, heteroVR.ci, heteroRR.ci) %>% select(GroupingTerm, parameter, value, ci.low, ci.high)

    # **Re-order grouping terms

    heterog4$GroupingTerm <- factor(heterog4$GroupingTerm,
                                    levels = c("Behaviour", "Morphology", "Metabolism",
                                               "Physiology", "Immunology", "Hematology", "Heart",
                                               "Hearing", "Eye", "All"))

    heterog4$GroupingTerm <- factor(heterog4$GroupingTerm,
                                    rev(levels(heterog4$GroupingTerm)))
    heterog4$label <- "All traits"
    # write.csv(heterog4, "heterog4.csv")

    #### Plot S1 C (Second-order meta analysis on heterogeneity)

    heterog5 <- heterog4
    heterog5$mean <- as.numeric(exp(heterog5$value))
    heterog5$ci.l <- as.numeric(exp(heterog5$ci.low))
    heterog5$ci.h <- as.numeric(exp(heterog5$ci.high))

    heterog6 <- heterog5

    HeteroS1 <-
      heterog6 %>%
      ggplot(aes(y = GroupingTerm, x = mean)) +
      geom_errorbarh(aes(
        xmin = ci.l,
        xmax = ci.h
      ),
      height = 0.1, show.legend = FALSE
      ) +
      geom_point(aes(shape = parameter),
        fill = "black",
        color = "black", size = 2.2,
        show.legend = FALSE
      ) +
      scale_x_continuous(
        limits = c(-0.1, 1.4),
        # breaks = c(0, 0.1, 0.2),
        name = "sigma^2"
      ) +
      # geom_vline(xintercept=0,
      # color='black',
      # linetype='dashed')+
      facet_grid(
        cols = vars(parameter), rows = vars(label),
        labeller = label_wrap_gen(width = 23),
        scales = "free",
        space = "free"
      ) +
      theme_bw() +
      theme(
        strip.text.y = element_text(angle = 270, size = 10, margin = margin(t = 15, r = 15, b = 15, l = 15)),
        strip.text.x = element_text(size = 12),
        strip.background = element_rect(colour = NULL, linetype = "blank", fill = "gray90"),
        text = element_text(size = 14),
        panel.spacing = unit(0.5, "lines"),
        panel.border = element_blank(),
        axis.line = element_line(),
        panel.grid.major.x = element_line(linetype = "solid", colour = "gray95"),
        panel.grid.major.y = element_line(linetype = "solid", color = "gray95"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        legend.title = element_blank(),
        axis.title.x = element_text(hjust = 0.5, size = 14),
        axis.title.y = element_blank())

# 4 Meta-analysis results

## 4.1 Preparing: Figure 4

### 4.1.1 Count data, based on First-order meta-analysis results

This includes all separate eligible traits. Re-ordering of grouping terms. Preparation of subplot 3A

    meta_clean$GroupingTerm <- factor(meta_clean$GroupingTerm,
                                      levels = c("Behaviour", "Morphology", "Metabolism", "Physiology", "Immunology", "Hematology", "Heart", "Hearing", "Eye"))
    meta_clean$GroupingTerm <- factor(meta_clean$GroupingTerm,
                                      rev(levels(meta_clean$GroupingTerm)))

    # *Preparing data for all traits

    meta.plot2.all <- meta_clean %>%
                      select(lnCVR, lnVR, lnRR, GroupingTerm) %>%
                      arrange(GroupingTerm)

    # lnVR has been removed here and in the steps below, as this is only included in the supplemental figure
          meta.plot2.all.b <- gather(meta.plot2.all, trait, value, c(lnCVR, lnRR))

    meta.plot2.all.b$trait <- factor(meta.plot2.all.b$trait, levels = c("lnCVR", "lnRR"))
          meta.plot2.all.c <- meta.plot2.all.b %>%
                              group_by_at(vars(trait, GroupingTerm)) %>%
                              summarise(
                                malebias = sum(value > 0), femalebias = sum(value <= 0), total = malebias + femalebias,
                                malepercent = malebias * 100 / total, femalepercent = femalebias * 100 / total
                              )

    meta.plot2.all.c$label <- "All traits"

    # Re-structure to create stacked bar plots

    meta.plot2.all.d <- as.data.frame(meta.plot2.all.c)
    meta.plot2.all.e <- gather(meta.plot2.all.d,
                               key = sex,
                               value = percent,
                               malepercent:femalepercent,
                               factor_key = TRUE)

    # Create new sample size variable

    meta.plot2.all.e$samplesize <- with(meta.plot2.all.e,
                                        ifelse(sex == "malepercent",
                                               malebias,
                                               femalebias))

    # Add summary row ('All') and re-arrange rows into correct order for plotting (warnings about coercing 'id' into character vector are ok)

    meta.plot2.all.f <- meta.plot2.all.e %>%
                        group_by(trait, sex) %>%
                          summarise(GroupingTerm = "All",
                                    malebias = sum(malebias),
                                    femalebias = sum(femalebias),
                                    total = malebias + femalebias,
                                    label = "All traits",
                                    samplesize = sum(samplesize)) %>%
                          mutate(percent = ifelse(sex == "femalepercent", femalebias*100/(malebias+femalebias), malebias*100/(malebias+femalebias))) %>%
                          bind_rows(meta.plot2.all.e, .) %>%
                          mutate(rownumber = row_number()) %>%
                          .[c(37, 1:9, 39, 10:18, 38, 19:27, 40, 28:36), ]
      #line references in previous code line corresponding to:
      #'lnCVR(male(All)), lnCVR(male('single grouping terms'), lnRR(male(All)), lnRR(male('single grouping terms')),
      #lnCVR(female(All)), lnCVR(female('single grouping terms'), lnRR(female(All)), lnRR(female('single grouping terms'))'

    meta.plot2.all.f$GroupingTerm <- factor(meta.plot2.all.f$GroupingTerm,
                                            levels = c("Behaviour", "Morphology", "Metabolism", "Physiology", "Immunology", "Hematology", "Heart", "Hearing", "Eye", "All"))
    meta.plot2.all.f$GroupingTerm <- factor(meta.plot2.all.f$GroupingTerm,
                                            rev(levels(meta.plot2.all.f$GroupingTerm)))

    malebias_Fig2_alltraits <-
        ggplot(meta.plot2.all.f) +
        aes(x = GroupingTerm, y = percent, fill = sex) +
        ylab("Percentage") +
        geom_col() +
        geom_hline(yintercept = 50, linetype = "dashed", color = "gray40") +
        geom_text(
          data = subset(meta.plot2.all.f, samplesize != 0), aes(label = samplesize), position = position_stack(vjust = .5),
          color = "white", size = 3.5
        ) +
        facet_grid(
          cols = vars(trait), rows = vars(label), labeller = label_wrap_gen(width = 18),
          scales = "free", space = "free"
        ) +
        scale_fill_brewer(palette = "Set2") +
        theme_bw(base_size = 18) +
        theme(
          strip.text.y = element_text(angle = 270, size = 10, margin = margin(t = 15, r = 15, b = 15, l = 15)),
          strip.text.x = element_text(size = 12),
          strip.background = element_rect(colour = NULL, linetype = "blank", fill = "gray90"),
          text = element_text(size = 14),
          panel.spacing = unit(0.5, "lines"),
          panel.border = element_blank(),
          axis.line = element_line(),
          panel.grid.major.x = element_line(linetype = "solid", colour = "gray95"),
          panel.grid.major.y = element_line(linetype = "solid", color = "gray95"),
          panel.grid.minor.y = element_blank(),
          panel.grid.minor.x = element_blank(),
          legend.position= "none",
          #axis.title.x = Percentage,
          axis.title.y = element_blank()
        ) +
        coord_flip()

### 4.1.2 Re-structure data for plotting

Data are re-structured, and grouping terms are being re-ordered

    overall3 <- gather(overall2, parameter, value, c(lnCVR, lnRR), factor_key = TRUE)

    lnCVR.ci <- overall3 %>% filter(parameter == "lnCVR") %>% mutate(ci.low = lnCVR_lower,
        ci.high = lnCVR_upper)
    lnVR.ci <- overall3 %>% filter(parameter == "lnVR") %>% mutate(ci.low = lnVR_lower,
        ci.high = lnVR_upper)
    lnRR.ci <- overall3 %>% filter(parameter == "lnRR") %>% mutate(ci.low = lnRR_lower,
        ci.high = lnRR_upper)

    overall4 <- bind_rows(lnCVR.ci, lnRR.ci) %>% select(GroupingTerm, parameter, value,
        ci.low, ci.high)

    # Re-order grouping terms

    overall4$GroupingTerm <- factor(overall4$GroupingTerm, levels = c("Behaviour", "Morphology",
        "Metabolism", "Physiology", "Immunology", "Hematology", "Heart", "Hearing", "Eye",
        "All"))
    overall4$GroupingTerm <- factor(overall4$GroupingTerm, rev(levels(overall4$GroupingTerm)))
    overall4$label <- "All traits"

    #### PLOT 4B

    # adding male / female symbols for Figure 4B additional packages for intergating
    # male / female symbols in plot 4
    library(png)
    library(grid)

    male <- readPNG(here("images", "MaleAqua.png"))
    female <- readPNG(here("images", "FemaleSalmon.png"))

    Metameta_Fig3_alltraits <- overall4 %>%
    ggplot(aes(y = GroupingTerm, x = value)) + geom_errorbarh(aes(xmin = ci.low, xmax = ci.high),
        height = 0.1, show.legend = FALSE) + geom_point(aes(shape = parameter), fill = "black",
        color = "black", size = 2.2, show.legend = FALSE) + scale_x_continuous(limits = c(-0.24,
        0.25), breaks = c(-0.2, -0.1, 0, 0.1, 0.2), name = "Effect size") + geom_vline(xintercept = 0,
        color = "black", linetype = "dashed") + facet_grid(cols = vars(parameter), rows = vars(label),
        labeller = label_wrap_gen(width = 23), scales = "free", space = "free") + theme_bw() +
        theme(strip.text.y = element_text(angle = 270, size = 10, margin = margin(t = 15,
            r = 15, b = 15, l = 15)), strip.text.x = element_text(size = 12), strip.background = element_rect(colour = NULL,
            linetype = "blank", fill = "gray90"), text = element_text(size = 14), panel.spacing = unit(0.5,
            "lines"), panel.border = element_blank(), axis.line = element_line(), panel.grid.major.x = element_line(linetype = "solid",
            colour = "gray95"), panel.grid.major.y = element_line(linetype = "solid",
            color = "gray95"), panel.grid.minor.y = element_blank(), panel.grid.minor.x = element_blank(),
            legend.title = element_blank(), axis.title.x = element_text(hjust = 0.5,
                size = 14), axis.title.y = element_blank()) + # addition of male & female symols
    annotation_custom(rasterGrob(female), xmin = -0.2, xmax = -0.1, ymin = 2.3, ymax = 4) +
        annotation_custom(rasterGrob(male), xmin = 0.1, xmax = 0.2, ymin = 2.3, ymax = 4)

    # Metameta_Fig3_alltraits

## 4.2 Main Manuscript: Figure 4

    Fig4 <- ggarrange(malebias_Fig2_alltraits, Metameta_Fig3_alltraits, nrow = 2, align = "v",
        heights = c(10, 9), labels = c("A", "B"))
    Fig4

Panel A shows the numbers of traits across functional groups that are either male-biased (blue-green) or female-biased (orange-red), as calculated in Figure 1.1D. Panel B shows effect sizes and 95% CI from separate meta-analysis for each functional group (step H in Figure1.1). Both panels represent results evaluated across all traits. Traits that are male biased are shown in blue, whereas female bias data is represented in orange.[Figure 4 in manuscript]

    ggsave("Figure4.pdf", plot = Fig4, width = 6, height = 7)

# 5 Supplemental Figures

## 5.1 Figure 5.1

This Figure is similar to Figure 4 in the main article except that lnVR has been added to the plot for comparison with lnCVR and lnRR.

    #### Combined Figure S1: overall Count data, Meta anlysis results, Heterogeneity

    FigS1 <- ggarrange(malebias_FigS1_alltraits, Metameta_FigS1_alltraits, HeteroS1,
        nrow = 3, align = "v", heights = c(1.1, 1, 1), labels = c("A", "B", "C"))
    FigS1

Numbers of either male (blue-green bars) or female (orange-red bars) biased traits (Panel A) across functional groups, this time for lnCVR (left hand side), lnVR (middle) and lnRR (right hand side). Panel B shows effect sizes from separate meta-analysis for each functional group, and Panel C contains results of heterogeneity analyses. All three panels represent results evaluated across all traits.[Figure 4 – figure supplement 1 in manuscript]

    ggsave("Figure4_1.pdf", plot = FigS1, width = 6, height = 7)

## 5.2 Figure 5.2

### 5.2.1 Sex-bias in traits sub-analysis

#### 5.2.1.1 Methods

Among the 147 traits (i.e. with ‘first-order’ meta-analytic results), we found 31 traits to be significantly male-biased and 30 female-biased for lnCVR (47 male-biased and 39 female- biased for lnRR; 74 male and 55 female-biased for lnVR). We used this statistical significance to categorize traits into two groups: male-based and female-biased traits. For each set of these sex-biased traits, we conducted second-order meta-analyses with trait as the unit of analysis in the same way as those using lnRR, lnVR and lnCVR (Figure 1.1J). As sensitivity analyses, we also employed another way of separating male- and female-biased traits. We quantified traits that were more than 10% biased towards either of the sexes, regardless of significance (N[male-biased] = 31, & N[female-biased] = 36 for lnCVR, N[male-biased] = 56, & N[female-biased] = 57 for lnVR , N[male-biased] = 42, & N[female-biased] = 33 for lnRR; see Figure 4 – figure supplement 1, Supplemental Table 3).

#### 5.2.1.2 Results

Having established that different sex-biases were present across trait functional groups, we explored specific patterns in male and female bias by categorizing traits into either male or female-biased, and then, quantifying the degree of sex bias. To do this, we focused on only those traits that were significantly biased towards one sex (as indicated by “CI not overlapping 0” in Figure Figure 4 – figure supplement 2). Selecting significantly sex-biased traits reduced the number of traits substantially (Figure 4 – figure supplement 2A; compared to all traits, Figure 4 and Figure 4 – figure supplement 1), and returned similar numbers of traits with either significant male- or female-bias for lnCVR (Figure 4 – figure supplement 2A).

Looking only at traits with significant male-bias, trait variability tended to be higher in males than in females overall (11.77% male bias, lnCVR = 0.111, CI = [0.093 to 0.129]; compared to 8.3% female-bias, lnCVR = -0.087, CI = [-0.109 to (-0.064)]; Figure 4 – figure supplement 2B). The degree of sex-bias varied across the functional trait groups, ranging from 4.6 % (Physiology, lnCVR = -0.047, CI = [-0.064 to (-0.030)] to 15.13 % (Eye morphology, lnCVR = -0.164, CI = [-0.249 to (-0.250)] in females, and from 8.6 % (Hematology, lnCVR = 0.082, CI = [0.048 to 0.117]) and 16.66 % (Heart, lnCVR = 0.154, CI = [0.085 to 0.223]) in males (Figure 4 – figure supplement 2B).

To put sex biases (sexual dimorphisms) in variability into perspective, we also quantify the degree of sexual dimorphism in significantly sex-biased traits (i.e., the extent of sexual dimorphisms). There were more traits that were significantly biased in mean trait value in males compared to females (Figure 4 – figure supplement 2A, lnRR). These traits were 11.73 % biased towards males Figure 4 – figure supplement 2B; lnRR = 0.111, CI = [0.083 – 0.138]) while 9.38 % towards females (lnRR = -0.094, CI = [-0.129 – (-0.068)]. The overall pattern for lnRR was more varied among the functional trait groups than for lnCVR, ranging from 1.2% (Hematology, lnRR = -0.013, CI = [-0.017 – (-0.007)] to 14.72 % (Immunology, lnRR = -0.159, CI = [-0.229 – (-0.089)] in female-based traits Figure 4 – figure supplement 2B, whereas in males-based traits from 1.8 % (Hearing, lnRR = 0.018, CI = [0.006 – 0.031]) to 25 % (Metabolism, lnRR = 0.223, CI = [0.110 – 0.335]) (Figure 4 – figure supplement 2B).

    library(ggpubr)
    FigS2b <- ggarrange(Metameta_FigS2_female.sig, Metameta_FigS2_male.sig, ncol = 2,
        nrow = 1, widths = c(1, 1.3), heights = c(1.2, 1.2))

    FigS2d <- ggarrange(Metameta_Fig3_female.perc, Metameta_Fig3_male.perc, ncol = 2,
        nrow = 1, widths = c(1, 1.2), heights = c(1.1, 1.1))

    # end combination Figure 5
    FigS2 <- ggarrange(malebias_FigS2_sigtraits, malebias_Fig2_over10, FigS2b, FigS2d,
        ncol = 1, nrow = 4, heights = c(2.5, 2.2, 2.1, 2), labels = c("A", " ", "B",
            " "))
    FigS2

1. Differences in numbers of affected traits, in variance (lnCVR and lnVR) and means (lnRR), where there is a significant difference between the sexes (i.e CI not overlapping zero), and where the sex bias is greater than 10% difference (regardless of significance). Panel B depicts results for the sex bias in those traits that differ between the sexes (second-order meta-analysis). Triangles represent sex bias in means (response ratio) and black circles differences in the coefficient of variation ratio (mean-adjusted variability). The orange-red bars represent trait groups with a female bias, blue-green bars male-biased traits.[Figure 4 – figure supplement 2 in manuscript]

    ggsave("Figure4_2.pdf", plot = FigS2, width = 6, height = 7)

# 6 .Supplemental Tables

## 6.1 Table 1

    # PDF pander(ST1, caption = 'Table 6.1: Summary of the available numbers of male
    # and female mice from each strain and originating institution.') HTML
    kable(ST1, caption = "Summary of the available numbers of male and female mice from each strain and originating institution.") %>%
        kable_styling() %>% scroll_box(width = "100%", height = "200px")

Summary of the available numbers of male and female mice from each strain and originating institution.
production_center | strain_name | sex | n
--- | --- | --- | ---
BCM | C57BL/6N | female | 653
BCM | C57BL/6N | male | 639
BCM | C57BL/6N;C57BL/6NTac | female | 47
BCM | C57BL/6N;C57BL/6NTac | male | 52
BCM | C57BL/6NCrl | female | 4
BCM | C57BL/6NCrl | male | 2
BCM | C57BL/6NJ | female | 6
BCM | C57BL/6NJ | male | 6
BCM | C57BL/6NTac | female | 1
BCM | C57BL/6NTac | male | 5
HMGU | C57BL/6NCrl | female | 313
HMGU | C57BL/6NCrl | male | 311
HMGU | C57BL/6NTac | female | 1045
HMGU | C57BL/6NTac | male | 1062
ICS | C57BL/6N | female | 1025
ICS | C57BL/6N | male | 1050
JAX | C57BL/6NJ | female | 2025
JAX | C57BL/6NJ | male | 2022
KMPC | C57BL/6N;C57BL/6NTac | female | 271
KMPC | C57BL/6N;C57BL/6NTac | male | 266
MARC | C57BL/6N | female | 936
MARC | C57BL/6N | male | 926
MRC Harwell | C57BL/6NTac | female | 2639
MRC Harwell | C57BL/6NTac | male | 2661
MRC Harwell | C57BL/6NTac | no_data | 3
RBRC | C57BL/6NJcl | female | 222
RBRC | C57BL/6NJcl | male | 222
RBRC | C57BL/6NTac | female | 526
RBRC | C57BL/6NTac | male | 523
TCP | C57BL/6NCrl | female | 552
TCP | C57BL/6NCrl | male | 524
TCP | C57BL6/NCrl | female | 2
TCP | C57BL6/NCrl | male | 2
UC Davis | C57BL/6N | male | 1
UC Davis | C57BL/6NCrl | female | 1155
UC Davis | C57BL/6NCrl | male | 1158
WTSI | B6Brd;B6Dnk;B6N-Tyr<c-Brd> | female | 97
WTSI | B6Brd;B6Dnk;B6N-Tyr<c-Brd> | male | 87
WTSI | C57BL/6J-Tyr<c-Brd> or C57BL/6NTac/USA | male | 3
WTSI | C57BL/6N | female | 1951
WTSI | C57BL/6N | male | 2008
WTSI | C57BL/6N;C57BL/6NTac | female | 41
WTSI | C57BL/6N;C57BL/6NTac | male | 7
WTSI | C57BL/6NCrl | male | 13
WTSI | C57BL/6NTac | female | 49
WTSI | C57BL/6NTac | male | 34

## 6.2 Table 2

    # Table detailing numbers of correlated and uncorrelated traits

    # PDF pander(cbind(meta1 %>% count(parameter_group)), caption = 'Table 6.2: Trait
    # categories (parameter_group) and the number of correlated traits within these
    # categories. Traits were meta-analysed using robumeta')

    kable(cbind(meta1 %>% count(parameter_group)), caption = "Trait categories (parameter_group) and the number of correlated traits within these categories. Traits were meta-analysed using robumeta") %>%
        kable_styling() %>% scroll_box(width = "100%", height = "200px")

Trait categories (parameter_group) and the number of correlated traits within these categories. Traits were meta-analysed using robumeta
parameter_group | n
--- | ---
12khz-evoked abr threshold | 1
18khz-evoked abr threshold | 1
24khz-evoked abr threshold | 1
30khz-evoked abr threshold | 1
6khz-evoked abr threshold | 1
alanine aminotransferase | 1
albumin | 1
alkaline phosphatase | 1
alpha-amylase | 1
area under glucose response curve | 1
aspartate aminotransferase | 1
B cells | 4
basophil cell count | 1
basophil differential count | 1
bmc/body weight | 1
body length | 1
body temp | 1
body weight | 1
body weight after experiment | 1
body weight before experiment | 1
bone area | 1
bone mineral content (excluding skull) | 1
bone mineral density (excluding skull) | 1
calcium | 1
cardiac output | 1
cd4 nkt | 6
cd4 t | 7
cd8 nkt | 6
cd8 t | 7
cdcs | 2
center average speed | 1
center distance travelled | 1
center permanence time | 1
center resting time | 1
chloride | 1
click-evoked abr threshold | 1
creatine kinase | 1
creatinine | 1
cv | 1
distance travelled - total | 1
dn nkt | 6
dn t | 7
ejection fraction | 1
end-diastolic diameter | 1
end-systolic diameter | 1
eosinophils | 3
fasted blood glucose concentration | 1
fat mass | 1
fat/body weight | 1
follicular b cells | 2
forelimb and hindlimb grip strength measurement mean | 1
forelimb grip strength measurement mean | 1
fractional shortening | 1
free fatty acids | 1
fructosamine | 1
glucose | 1
hdl-cholesterol | 1
heart weight | 1
heart weight normalised against body weight | 1
hematocrit | 1
hemoglobin | 1
hr | 1
hrv | 1
initial response to glucose challenge | 1
insulin | 1
iron | 1
lactate dehydrogenase | 1
latency to center entry | 1
ldl-cholesterol | 1
lean mass | 1
lean/body weight | 1
left anterior chamber depth | 1
left corneal thickness | 1
left inner nuclear layer | 1
left outer nuclear layer | 1
left posterior chamber depth | 1
left total retinal thickness | 1
locomotor activity | 1
luc | 2
lvawd | 1
lvaws | 1
lvidd | 1
lvids | 1
lvpwd | 1
lvpws | 1
lymphocytes | 2
magnesium | 1
mean cell hemoglobin concentration | 1
mean cell volume | 1
mean corpuscular hemoglobin | 1
mean platelet volume | 1
mean r amplitude | 1
mean sr amplitude | 1
monocytes | 3
neutrophils | 3
nk cells | 6
nkt cells | 4
number of center entries | 1
number of rears - total | 1
others | 1
pdcs | 1
percentage center time | 1
percentage of live gated events | 2
periphery average speed | 1
periphery distance travelled | 1
periphery permanence time | 1
periphery resting time | 1
phosphorus | 1
platelet count | 1
pnn5(6>ms) | 1
potassium | 1
pq | 1
pr | 1
pre-pulse inhibition | 5
qrs | 1
qtc | 1
qtc dispersion | 1
red blood cell count | 1
red blood cell distribution width | 1
respiration rate | 1
respiratory exchange ratio | 1
response amplitude | 10
right anterior chamber depth | 1
right corneal thickness | 1
right inner nuclear layer | 1
right outer nuclear layer | 1
right posterior chamber depth | 1
right total retinal thickness | 1
rmssd | 1
rp macrophage (cd19- cd11c-) | 1
rr | 1
sodium | 1
spleen weight | 1
st | 1
stroke volume | 1
t cells | 3
tibia length | 1
total bilirubin | 1
total cholesterol | 1
total food intake | 1
total protein | 1
total water intake | 1
triglycerides | 1
urea (blood urea nitrogen - bun) | 1
uric acid | 1
white blood cell count | 1
whole arena average speed | 1
whole arena resting time | 1

## 6.3 Table 3

    # PDF pander(metacombo, caption = 'Table 6.3: We use this corrected (for
    # correlated traits) results table, which contains each of the meta-analytic
    # means for all effect sizes of interest, for further analyses.  We further use
    # this table as part of the Shiny App, which is able to provide the percentage
    # differences between males and females for mean, variance and coefficient of
    # variance.')

    # HTML
    kable(metacombo, digits = 3, caption = "We use this corrected (for correlated traits) results table, which contains each of the meta-analytic means for all effect sizes of interest, for further analyses.  We further use this table as part of the Shiny App, which is able to provide the percentage differences between males and females for mean, variance and coefficient of variance.") %>%
        kable_styling() %>% scroll_box(width = "100%", height = "200px")

We use this corrected (for correlated traits) results table, which contains each of the meta-analytic means for all effect sizes of interest, for further analyses. We further use this table as part of the Shiny App, which is able to provide the percentage differences between males and females for mean, variance and coefficient of variance.
parameter_group | counts | procedure | GroupingTerm | lnCVR | lnCVR_lower | lnCVR_upper | lnCVR_se | lnVR | lnVR_lower | lnVR_upper | lnVR_se | lnRR | lnRR_lower | lnRR_upper | lnRR_se
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
pre-pulse inhibition | 5 | Acoustic Startle and Pre-pulse Inhibition (PPI) | Behaviour | 0.023 | -0.080 | 0.127 | 0.037 | 0.009 | -0.036 | 0.055 | 0.014 | -0.005 | -0.043 | 0.032 | 0.013
B cells | 4 | Immunophenotyping | Immunology | -0.094 | -0.250 | 0.062 | 0.043 | -0.100 | -0.207 | 0.008 | 0.025 | -0.003 | -0.130 | 0.125 | 0.039
cd4 nkt | 6 | Immunophenotyping | Immunology | -0.029 | -0.057 | -0.001 | 0.010 | -0.202 | -0.310 | -0.094 | 0.033 | -0.234 | -0.401 | -0.068 | 0.063
cd4 t | 7 | Immunophenotyping | Immunology | -0.151 | -0.243 | -0.059 | 0.036 | -0.170 | -0.263 | -0.077 | 0.035 | -0.003 | -0.041 | 0.035 | 0.015
cd8 nkt | 6 | Immunophenotyping | Immunology | -0.042 | -0.078 | -0.007 | 0.012 | -0.030 | -0.182 | 0.122 | 0.053 | 0.004 | -0.057 | 0.064 | 0.021
cd8 t | 7 | Immunophenotyping | Immunology | -0.122 | -0.218 | -0.027 | 0.036 | -0.158 | -0.234 | -0.082 | 0.027 | -0.042 | -0.051 | -0.032 | 0.002
cdcs | 2 | Immunophenotyping | Immunology | -0.036 | -0.359 | 0.286 | 0.025 | 0.108 | -0.057 | 0.273 | 0.013 | 0.164 | -0.170 | 0.499 | 0.026
dn nkt | 6 | Immunophenotyping | Immunology | -0.062 | -0.136 | 0.012 | 0.026 | -0.157 | -0.281 | -0.033 | 0.045 | -0.173 | -0.291 | -0.055 | 0.044
dn t | 7 | Immunophenotyping | Immunology | -0.080 | -0.184 | 0.025 | 0.042 | -0.242 | -0.343 | -0.141 | 0.041 | -0.230 | -0.252 | -0.208 | 0.007
eosinophils | 3 | Hematology | Hematology | -0.066 | -0.281 | 0.148 | 0.033 | -0.015 | -0.405 | 0.374 | 0.087 | -0.004 | -0.241 | 0.232 | 0.051
follicular b cells | 2 | Immunophenotyping | Immunology | -0.116 | -0.726 | 0.494 | 0.048 | -0.105 | -0.695 | 0.485 | 0.046 | 0.005 | -0.187 | 0.198 | 0.015
luc | 2 | Hematology | Hematology | 0.018 | -0.204 | 0.240 | 0.017 | 0.266 | -1.225 | 1.757 | 0.117 | 0.222 | -1.414 | 1.857 | 0.129
lymphocytes | 2 | Hematology | Hematology | 0.081 | -2.262 | 2.423 | 0.184 | 0.155 | -1.089 | 1.399 | 0.098 | 0.060 | -1.013 | 1.134 | 0.084
monocytes | 3 | Hematology | Hematology | -0.021 | -0.203 | 0.160 | 0.042 | 0.078 | -0.181 | 0.338 | 0.059 | 0.103 | -0.148 | 0.353 | 0.057
neutrophils | 3 | Hematology | Hematology | 0.259 | 0.013 | 0.504 | 0.056 | 0.380 | -0.206 | 0.966 | 0.132 | 0.132 | -0.267 | 0.531 | 0.092
nk cells | 6 | Immunophenotyping | Immunology | -0.041 | -0.096 | 0.013 | 0.020 | 0.016 | -0.070 | 0.102 | 0.032 | 0.047 | -0.016 | 0.111 | 0.023
nkt cells | 4 | Immunophenotyping | Immunology | 0.003 | -0.107 | 0.114 | 0.029 | -0.246 | -0.445 | -0.047 | 0.043 | -0.182 | -0.323 | -0.041 | 0.031
percentage of live gated events | 2 | Immunophenotyping | Immunology | -0.093 | -0.304 | 0.117 | 0.017 | -0.041 | -0.141 | 0.059 | 0.008 | 0.050 | 0.008 | 0.092 | 0.003
response amplitude | 10 | Acoustic Startle and Pre-pulse Inhibition (PPI) | Behaviour | 0.033 | -0.013 | 0.079 | 0.020 | 0.255 | 0.197 | 0.313 | 0.026 | 0.202 | 0.111 | 0.292 | 0.040
t cells | 3 | Immunophenotyping | Immunology | -0.134 | -0.275 | 0.007 | 0.033 | -0.124 | -0.412 | 0.164 | 0.067 | -0.001 | -0.166 | 0.165 | 0.037
12khz-evoked abr threshold | 1 | Auditory Brain Stem Response | Hearing | 0.054 | -0.006 | 0.113 | 0.030 | 0.087 | 0.007 | 0.167 | 0.041 | 0.002 | -0.021 | 0.026 | 0.012
18khz-evoked abr threshold | 1 | Auditory Brain Stem Response | Hearing | 0.024 | -0.033 | 0.081 | 0.029 | 0.025 | -0.049 | 0.099 | 0.038 | -0.020 | -0.043 | 0.003 | 0.012
24khz-evoked abr threshold | 1 | Auditory Brain Stem Response | Hearing | 0.052 | -0.015 | 0.118 | 0.034 | -0.089 | -0.332 | 0.154 | 0.124 | -0.022 | -0.044 | 0.000 | 0.011
30khz-evoked abr threshold | 1 | Auditory Brain Stem Response | Hearing | 0.017 | -0.053 | 0.088 | 0.036 | -0.034 | -0.102 | 0.033 | 0.034 | -0.050 | -0.075 | -0.025 | 0.013
6khz-evoked abr threshold | 1 | Auditory Brain Stem Response | Hearing | -0.008 | -0.042 | 0.026 | 0.017 | 0.014 | -0.019 | 0.047 | 0.017 | 0.018 | 0.006 | 0.031 | 0.006
alanine aminotransferase | 1 | Clinical Chemistry | Physiology | -0.068 | -0.190 | 0.053 | 0.062 | 0.059 | -0.132 | 0.249 | 0.097 | 0.107 | 0.032 | 0.182 | 0.038
albumin | 1 | Clinical Chemistry | Physiology | 0.113 | 0.045 | 0.181 | 0.035 | 0.056 | -0.008 | 0.120 | 0.033 | -0.057 | -0.073 | -0.040 | 0.008
alkaline phosphatase | 1 | Clinical Chemistry | Physiology | 0.104 | 0.045 | 0.164 | 0.030 | -0.311 | -0.398 | -0.224 | 0.044 | -0.422 | -0.469 | -0.374 | 0.024
alpha-amylase | 1 | Clinical Chemistry | Physiology | 0.038 | -0.042 | 0.119 | 0.041 | 0.280 | 0.162 | 0.398 | 0.060 | 0.225 | 0.179 | 0.270 | 0.023
area under glucose response curve | 1 | Intraperitoneal glucose tolerance test (IPGTT) | Metabolism | -0.153 | -0.221 | -0.085 | 0.035 | 0.275 | 0.195 | 0.355 | 0.041 | 0.436 | 0.366 | 0.506 | 0.036
aspartate aminotransferase | 1 | Clinical Chemistry | Physiology | 0.012 | -0.123 | 0.147 | 0.069 | -0.057 | -0.246 | 0.132 | 0.096 | -0.059 | -0.133 | 0.016 | 0.038
basophil cell count | 1 | Hematology | Hematology | -0.092 | -0.202 | 0.019 | 0.056 | 0.203 | -0.013 | 0.419 | 0.110 | 0.268 | 0.064 | 0.471 | 0.104
basophil differential count | 1 | Hematology | Hematology | -0.093 | -0.179 | -0.008 | 0.044 | -0.064 | -0.283 | 0.155 | 0.112 | -0.016 | -0.110 | 0.079 | 0.048
bmc/body weight | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.131 | 0.033 | 0.230 | 0.050 | -0.045 | -0.134 | 0.044 | 0.045 | -0.172 | -0.221 | -0.124 | 0.025
body length | 1 | Body Composition (DEXA lean/fat) | Morphology | -0.035 | -0.082 | 0.013 | 0.024 | -0.006 | -0.053 | 0.041 | 0.024 | 0.028 | 0.023 | 0.033 | 0.003
body temp | 1 | Echo | Heart | -0.033 | -0.107 | 0.042 | 0.038 | -0.030 | -0.104 | 0.044 | 0.038 | 0.002 | -0.001 | 0.004 | 0.001
body weight | 1 | Body Weight | Morphology | 0.025 | -0.042 | 0.091 | 0.034 | 0.234 | 0.169 | 0.298 | 0.033 | 0.210 | 0.194 | 0.225 | 0.008
body weight after experiment | 1 | Indirect Calorimetry | Metabolism | 0.085 | 0.030 | 0.141 | 0.028 | 0.285 | 0.233 | 0.337 | 0.027 | 0.203 | 0.186 | 0.220 | 0.009
body weight before experiment | 1 | Indirect Calorimetry | Metabolism | 0.105 | 0.041 | 0.169 | 0.033 | 0.304 | 0.244 | 0.364 | 0.031 | 0.201 | 0.182 | 0.220 | 0.010
bone area | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.098 | 0.027 | 0.169 | 0.036 | 0.129 | 0.053 | 0.204 | 0.038 | 0.032 | 0.000 | 0.063 | 0.016
bone mineral content (excluding skull) | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.171 | 0.063 | 0.279 | 0.055 | 0.209 | 0.102 | 0.317 | 0.055 | 0.037 | -0.013 | 0.088 | 0.026
bone mineral density (excluding skull) | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.054 | -0.088 | 0.197 | 0.073 | 0.049 | -0.109 | 0.207 | 0.081 | 0.001 | -0.019 | 0.021 | 0.010
calcium | 1 | Clinical Chemistry | Physiology | 0.010 | -0.046 | 0.066 | 0.029 | 0.014 | -0.042 | 0.070 | 0.029 | 0.004 | 0.000 | 0.007 | 0.002
cardiac output | 1 | Echo | Heart | 0.013 | -0.080 | 0.107 | 0.048 | 0.102 | 0.021 | 0.183 | 0.041 | 0.093 | 0.058 | 0.129 | 0.018
center average speed | 1 | Open Field | Behaviour | 0.017 | -0.040 | 0.074 | 0.029 | -0.059 | -0.100 | -0.017 | 0.021 | -0.072 | -0.115 | -0.030 | 0.022
center distance travelled | 1 | Open Field | Behaviour | -0.016 | -0.073 | 0.041 | 0.029 | -0.106 | -0.202 | -0.010 | 0.049 | -0.094 | -0.195 | 0.007 | 0.051
center permanence time | 1 | Open Field | Behaviour | -0.025 | -0.083 | 0.032 | 0.029 | -0.026 | -0.101 | 0.050 | 0.039 | -0.004 | -0.090 | 0.083 | 0.044
center resting time | 1 | Open Field | Behaviour | 0.024 | -0.074 | 0.123 | 0.050 | -0.023 | -0.155 | 0.109 | 0.067 | -0.063 | -0.222 | 0.095 | 0.081
chloride | 1 | Clinical Chemistry | Physiology | 0.032 | -0.127 | 0.191 | 0.081 | 0.024 | -0.144 | 0.192 | 0.086 | -0.013 | -0.018 | -0.008 | 0.003
click-evoked abr threshold | 1 | Auditory Brain Stem Response | Hearing | -0.053 | -0.153 | 0.048 | 0.051 | -0.056 | -0.183 | 0.071 | 0.065 | -0.015 | -0.058 | 0.027 | 0.022
creatine kinase | 1 | Clinical Chemistry | Physiology | 0.024 | -0.107 | 0.155 | 0.067 | -0.132 | -0.397 | 0.133 | 0.135 | -0.134 | -0.384 | 0.115 | 0.127
creatinine | 1 | Clinical Chemistry | Physiology | 0.035 | -0.023 | 0.093 | 0.030 | 0.107 | -0.220 | 0.433 | 0.167 | -0.084 | -0.132 | -0.037 | 0.024
cv | 1 | Electrocardiogram (ECG) | Heart | 0.187 | 0.072 | 0.303 | 0.059 | -0.090 | -0.248 | 0.069 | 0.081 | -0.240 | -0.341 | -0.139 | 0.051
distance travelled - total | 1 | Open Field | Behaviour | -0.019 | -0.086 | 0.048 | 0.034 | -0.127 | -0.200 | -0.055 | 0.037 | -0.112 | -0.182 | -0.043 | 0.035
ejection fraction | 1 | Echo | Heart | -0.030 | -0.135 | 0.074 | 0.053 | -0.053 | -0.148 | 0.043 | 0.049 | -0.028 | -0.049 | -0.008 | 0.011
end-diastolic diameter | 1 | Echo | Heart | 0.112 | 0.043 | 0.181 | 0.035 | 0.174 | 0.088 | 0.261 | 0.044 | 0.060 | 0.035 | 0.085 | 0.013
end-systolic diameter | 1 | Echo | Heart | -0.008 | -0.078 | 0.061 | 0.036 | 0.067 | -0.002 | 0.135 | 0.035 | 0.076 | 0.045 | 0.108 | 0.016
fasted blood glucose concentration | 1 | Intraperitoneal glucose tolerance test (IPGTT) | Metabolism | -0.018 | -0.126 | 0.090 | 0.055 | 0.070 | -0.030 | 0.171 | 0.051 | 0.087 | 0.049 | 0.124 | 0.019
fat mass | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.041 | -0.043 | 0.125 | 0.043 | 0.371 | 0.270 | 0.473 | 0.052 | 0.328 | 0.267 | 0.390 | 0.031
fat/body weight | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.078 | -0.012 | 0.167 | 0.046 | 0.202 | 0.108 | 0.296 | 0.048 | 0.124 | 0.064 | 0.183 | 0.030
forelimb and hindlimb grip strength measurement mean | 1 | Grip Strength | Morphology | 0.058 | 0.004 | 0.112 | 0.027 | 0.115 | 0.053 | 0.176 | 0.031 | 0.054 | 0.029 | 0.079 | 0.013
forelimb grip strength measurement mean | 1 | Grip Strength | Morphology | 0.027 | -0.019 | 0.072 | 0.023 | 0.100 | 0.054 | 0.145 | 0.023 | 0.070 | 0.044 | 0.096 | 0.013
fractional shortening | 1 | Echo | Heart | -0.015 | -0.116 | 0.086 | 0.052 | -0.058 | -0.156 | 0.041 | 0.050 | -0.041 | -0.057 | -0.026 | 0.008
free fatty acids | 1 | Clinical Chemistry | Physiology | 0.028 | -0.100 | 0.157 | 0.066 | 0.055 | -0.074 | 0.185 | 0.066 | 0.019 | -0.009 | 0.048 | 0.015
fructosamine | 1 | Clinical Chemistry | Physiology | -0.040 | -0.120 | 0.040 | 0.041 | -0.068 | -0.151 | 0.016 | 0.043 | -0.028 | -0.069 | 0.013 | 0.021
glucose | 1 | Clinical Chemistry | Physiology | 0.069 | 0.018 | 0.120 | 0.026 | 0.128 | 0.042 | 0.214 | 0.044 | 0.065 | 0.022 | 0.108 | 0.022
hdl-cholesterol | 1 | Clinical Chemistry | Physiology | -0.065 | -0.126 | -0.004 | 0.031 | 0.172 | 0.070 | 0.275 | 0.052 | 0.261 | 0.218 | 0.303 | 0.022
heart weight | 1 | Heart Weight | Morphology | 0.177 | 0.067 | 0.286 | 0.056 | 0.365 | 0.217 | 0.513 | 0.076 | 0.174 | 0.141 | 0.207 | 0.017
heart weight normalised against body weight | 1 | Heart Weight | Morphology | 0.079 | -0.006 | 0.165 | 0.044 | 0.036 | -0.097 | 0.168 | 0.068 | -0.050 | -0.084 | -0.016 | 0.017
hematocrit | 1 | Hematology | Hematology | 0.057 | -0.052 | 0.165 | 0.055 | 0.074 | -0.033 | 0.180 | 0.054 | 0.017 | 0.004 | 0.031 | 0.007
hemoglobin | 1 | Hematology | Hematology | 0.087 | 0.027 | 0.146 | 0.030 | 0.087 | 0.019 | 0.154 | 0.034 | 0.005 | -0.008 | 0.018 | 0.007
hr | 1 | Electrocardiogram (ECG) | Heart | -0.063 | -0.173 | 0.047 | 0.056 | -0.014 | -0.149 | 0.121 | 0.069 | 0.041 | -0.014 | 0.095 | 0.028
hrv | 1 | Electrocardiogram (ECG) | Heart | 0.172 | 0.109 | 0.235 | 0.032 | -0.081 | -0.213 | 0.050 | 0.067 | -0.250 | -0.366 | -0.135 | 0.059
initial response to glucose challenge | 1 | Intraperitoneal glucose tolerance test (IPGTT) | Metabolism | -0.097 | -0.150 | -0.043 | 0.027 | 0.043 | 0.014 | 0.072 | 0.015 | 0.118 | 0.085 | 0.151 | 0.017
insulin | 1 | Insulin Blood Level | Metabolism | -0.099 | -0.372 | 0.174 | 0.139 | 0.177 | -0.194 | 0.549 | 0.189 | 0.445 | 0.094 | 0.795 | 0.179
iron | 1 | Clinical Chemistry | Physiology | -0.097 | -0.214 | 0.019 | 0.060 | -0.253 | -0.396 | -0.111 | 0.073 | -0.153 | -0.193 | -0.113 | 0.021
lactate dehydrogenase | 1 | Clinical Chemistry | Physiology | 0.094 | -0.021 | 0.210 | 0.059 | 0.141 | -0.062 | 0.344 | 0.104 | 0.032 | -0.141 | 0.205 | 0.088
latency to center entry | 1 | Open Field | Behaviour | 0.125 | 0.033 | 0.218 | 0.047 | 0.364 | 0.206 | 0.523 | 0.081 | 0.273 | 0.074 | 0.473 | 0.102
ldl-cholesterol | 1 | Clinical Chemistry | Physiology | 0.423 | 0.155 | 0.691 | 0.137 | 0.267 | -0.096 | 0.630 | 0.185 | -0.162 | -0.601 | 0.278 | 0.224
lean mass | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.144 | 0.076 | 0.211 | 0.035 | 0.338 | 0.266 | 0.410 | 0.037 | 0.193 | 0.175 | 0.211 | 0.009
lean/body weight | 1 | Body Composition (DEXA lean/fat) | Morphology | 0.195 | 0.091 | 0.300 | 0.053 | 0.184 | 0.086 | 0.282 | 0.050 | -0.012 | -0.026 | 0.001 | 0.007
left anterior chamber depth | 1 | Eye Morphology | Eye | -0.185 | -0.431 | 0.060 | 0.125 | -0.153 | -0.401 | 0.094 | 0.126 | 0.033 | 0.028 | 0.038 | 0.002
left corneal thickness | 1 | Eye Morphology | Eye | -0.145 | -0.234 | -0.055 | 0.046 | -0.135 | -0.223 | -0.047 | 0.045 | 0.008 | -0.006 | 0.021 | 0.007
left inner nuclear layer | 1 | Eye Morphology | Eye | 0.048 | -0.036 | 0.132 | 0.043 | 0.049 | -0.035 | 0.132 | 0.043 | 0.001 | -0.010 | 0.011 | 0.005
left outer nuclear layer | 1 | Eye Morphology | Eye | -0.068 | -0.151 | 0.016 | 0.043 | -0.062 | -0.145 | 0.022 | 0.043 | 0.006 | 0.001 | 0.012 | 0.003
left posterior chamber depth | 1 | Eye Morphology | Eye | -0.263 | -0.473 | -0.053 | 0.107 | -0.269 | -0.479 | -0.058 | 0.107 | -0.003 | -0.015 | 0.009 | 0.006
left total retinal thickness | 1 | Eye Morphology | Eye | -0.198 | -0.439 | 0.044 | 0.123 | -0.193 | -0.427 | 0.040 | 0.119 | 0.003 | -0.003 | 0.009 | 0.003
locomotor activity | 1 | Combined SHIRPA and Dysmorphology | Behaviour | 0.096 | 0.022 | 0.170 | 0.038 | -0.016 | -0.058 | 0.026 | 0.021 | -0.111 | -0.176 | -0.045 | 0.033
lvawd | 1 | Echo | Heart | 0.023 | -0.025 | 0.070 | 0.024 | 0.045 | -0.001 | 0.092 | 0.024 | 0.025 | 0.011 | 0.038 | 0.007
lvaws | 1 | Echo | Heart | -0.002 | -0.252 | 0.248 | 0.128 | 0.023 | -0.178 | 0.224 | 0.103 | 0.011 | -0.031 | 0.053 | 0.021
lvidd | 1 | Echo | Heart | 0.045 | -0.024 | 0.115 | 0.035 | 0.098 | 0.021 | 0.175 | 0.039 | 0.053 | 0.038 | 0.068 | 0.008
lvids | 1 | Echo | Heart | -0.064 | -0.199 | 0.072 | 0.069 | 0.008 | -0.134 | 0.150 | 0.072 | 0.076 | 0.053 | 0.099 | 0.012
lvpwd | 1 | Echo | Heart | -0.032 | -0.126 | 0.062 | 0.048 | -0.010 | -0.127 | 0.106 | 0.060 | 0.030 | 0.013 | 0.047 | 0.009
lvpws | 1 | Echo | Heart | -0.019 | -0.101 | 0.063 | 0.042 | 0.009 | -0.082 | 0.100 | 0.047 | 0.027 | 0.006 | 0.047 | 0.010
magnesium | 1 | Urinalysis | Physiology | 0.016 | -0.023 | 0.055 | 0.020 | -0.051 | -0.117 | 0.014 | 0.033 | -0.041 | -0.114 | 0.031 | 0.037
mean cell hemoglobin concentration | 1 | Hematology | Hematology | 0.038 | -0.088 | 0.164 | 0.064 | 0.025 | -0.109 | 0.159 | 0.068 | -0.011 | -0.015 | -0.008 | 0.002
mean cell volume | 1 | Hematology | Hematology | 0.004 | -0.096 | 0.104 | 0.051 | -0.003 | -0.096 | 0.090 | 0.048 | -0.006 | -0.010 | -0.003 | 0.002
mean corpuscular hemoglobin | 1 | Hematology | Hematology | -0.003 | -0.065 | 0.060 | 0.032 | -0.019 | -0.082 | 0.044 | 0.032 | -0.017 | -0.020 | -0.014 | 0.001
mean platelet volume | 1 | Hematology | Hematology | 0.049 | -0.004 | 0.102 | 0.027 | 0.035 | -0.021 | 0.092 | 0.029 | -0.017 | -0.028 | -0.007 | 0.005
mean r amplitude | 1 | Electrocardiogram (ECG) | Heart | 0.008 | -0.028 | 0.045 | 0.019 | -0.095 | -0.163 | -0.027 | 0.035 | -0.084 | -0.150 | -0.017 | 0.034
mean sr amplitude | 1 | Electrocardiogram (ECG) | Heart | 0.028 | -0.013 | 0.070 | 0.021 | -0.088 | -0.127 | -0.048 | 0.020 | -0.113 | -0.156 | -0.070 | 0.022
number of center entries | 1 | Open Field | Behaviour | 0.015 | -0.053 | 0.084 | 0.035 | -0.036 | -0.095 | 0.023 | 0.030 | -0.059 | -0.168 | 0.050 | 0.056
number of rears - total | 1 | Open Field | Behaviour | -0.001 | -0.114 | 0.112 | 0.058 | 0.187 | -0.039 | 0.413 | 0.115 | 0.179 | 0.057 | 0.302 | 0.063
others | 1 | Immunophenotyping | Immunology | -0.168 | -0.260 | -0.077 | 0.047 | -0.152 | -0.244 | -0.059 | 0.047 | 0.020 | 0.005 | 0.034 | 0.007
pdcs | 1 | Immunophenotyping | Immunology | -0.173 | -0.400 | 0.054 | 0.116 | -0.257 | -0.719 | 0.204 | 0.235 | -0.092 | -0.252 | 0.069 | 0.082
percentage center time | 1 | Open Field | Behaviour | -0.022 | -0.086 | 0.042 | 0.033 | -0.019 | -0.091 | 0.053 | 0.037 | -0.006 | -0.097 | 0.085 | 0.046
periphery average speed | 1 | Open Field | Behaviour | -0.044 | -0.108 | 0.019 | 0.033 | -0.140 | -0.212 | -0.068 | 0.037 | -0.096 | -0.145 | -0.048 | 0.025
periphery distance travelled | 1 | Open Field | Behaviour | -0.031 | -0.092 | 0.029 | 0.031 | -0.134 | -0.187 | -0.081 | 0.027 | -0.104 | -0.171 | -0.036 | 0.035
periphery permanence time | 1 | Open Field | Behaviour | -0.037 | -0.128 | 0.054 | 0.046 | -0.029 | -0.101 | 0.042 | 0.036 | 0.008 | -0.014 | 0.029 | 0.011
periphery resting time | 1 | Open Field | Behaviour | -0.054 | -0.127 | 0.019 | 0.037 | -0.057 | -0.107 | -0.007 | 0.025 | 0.003 | -0.056 | 0.061 | 0.030
phosphorus | 1 | Clinical Chemistry | Physiology | -0.049 | -0.084 | -0.013 | 0.018 | -0.083 | -0.158 | -0.008 | 0.038 | -0.042 | -0.081 | -0.003 | 0.020
platelet count | 1 | Hematology | Hematology | 0.074 | 0.021 | 0.127 | 0.027 | 0.242 | 0.187 | 0.296 | 0.028 | 0.164 | 0.137 | 0.191 | 0.014
pnn5(6>ms) | 1 | Electrocardiogram (ECG) | Heart | 0.291 | 0.172 | 0.410 | 0.061 | -0.293 | -0.527 | -0.058 | 0.120 | -0.600 | -0.924 | -0.277 | 0.165
potassium | 1 | Clinical Chemistry | Physiology | -0.071 | -0.221 | 0.080 | 0.077 | -0.007 | -0.173 | 0.158 | 0.084 | 0.070 | 0.048 | 0.093 | 0.012
pq | 1 | Electrocardiogram (ECG) | Heart | -0.065 | -0.154 | 0.024 | 0.045 | -0.065 | -0.127 | -0.003 | 0.032 | 0.002 | -0.026 | 0.029 | 0.014
pr | 1 | Electrocardiogram (ECG) | Heart | -0.056 | -0.105 | -0.008 | 0.025 | -0.075 | -0.124 | -0.027 | 0.025 | -0.018 | -0.032 | -0.005 | 0.007
qrs | 1 | Electrocardiogram (ECG) | Heart | 0.073 | 0.035 | 0.110 | 0.019 | 0.068 | 0.030 | 0.106 | 0.019 | -0.005 | -0.015 | 0.005 | 0.005
qtc | 1 | Electrocardiogram (ECG) | Heart | 0.033 | -0.010 | 0.076 | 0.022 | 0.031 | -0.021 | 0.083 | 0.026 | -0.001 | -0.009 | 0.008 | 0.004
qtc dispersion | 1 | Electrocardiogram (ECG) | Heart | 0.003 | -0.052 | 0.059 | 0.028 | -0.005 | -0.106 | 0.097 | 0.052 | -0.008 | -0.051 | 0.036 | 0.022
red blood cell count | 1 | Hematology | Hematology | 0.077 | 0.007 | 0.147 | 0.036 | 0.100 | 0.032 | 0.168 | 0.035 | 0.023 | 0.009 | 0.037 | 0.007
red blood cell distribution width | 1 | Hematology | Hematology | 0.125 | -0.004 | 0.253 | 0.065 | 0.135 | -0.004 | 0.274 | 0.071 | 0.010 | -0.003 | 0.024 | 0.007
respiration rate | 1 | Echo | Heart | -0.138 | -0.218 | -0.059 | 0.041 | -0.070 | -0.180 | 0.039 | 0.056 | 0.061 | 0.023 | 0.099 | 0.020
respiratory exchange ratio | 1 | Indirect Calorimetry | Metabolism | -0.012 | -0.090 | 0.066 | 0.040 | -0.011 | -0.088 | 0.067 | 0.039 | 0.002 | -0.006 | 0.009 | 0.004
right anterior chamber depth | 1 | Eye Morphology | Eye | -0.449 | -1.329 | 0.431 | 0.449 | -0.416 | -1.292 | 0.460 | 0.447 | 0.032 | 0.026 | 0.037 | 0.003
right corneal thickness | 1 | Eye Morphology | Eye | -0.036 | -0.228 | 0.157 | 0.098 | -0.031 | -0.196 | 0.135 | 0.085 | -0.001 | -0.024 | 0.021 | 0.011
right inner nuclear layer | 1 | Eye Morphology | Eye | -0.255 | -0.763 | 0.254 | 0.260 | -0.279 | -0.837 | 0.280 | 0.285 | -0.018 | -0.066 | 0.031 | 0.025
right outer nuclear layer | 1 | Eye Morphology | Eye | 0.006 | -0.078 | 0.090 | 0.043 | 0.011 | -0.073 | 0.095 | 0.043 | 0.006 | 0.000 | 0.011 | 0.003
right posterior chamber depth | 1 | Eye Morphology | Eye | -0.078 | -0.291 | 0.135 | 0.109 | -0.076 | -0.289 | 0.136 | 0.109 | 0.007 | -0.018 | 0.032 | 0.013
right total retinal thickness | 1 | Eye Morphology | Eye | -0.199 | -0.646 | 0.248 | 0.228 | -0.193 | -0.629 | 0.243 | 0.222 | 0.005 | -0.005 | 0.015 | 0.005
rmssd | 1 | Electrocardiogram (ECG) | Heart | 0.180 | -0.088 | 0.448 | 0.137 | -0.016 | -0.411 | 0.379 | 0.202 | -0.118 | -0.245 | 0.009 | 0.065
rp macrophage (cd19- cd11c-) | 1 | Immunophenotyping | Immunology | -0.077 | -0.340 | 0.187 | 0.134 | -0.075 | -0.335 | 0.186 | 0.133 | -0.075 | -0.207 | 0.058 | 0.068
rr | 1 | Electrocardiogram (ECG) | Heart | -0.076 | -0.188 | 0.035 | 0.057 | -0.090 | -0.206 | 0.027 | 0.060 | -0.013 | -0.021 | -0.004 | 0.005
sodium | 1 | Clinical Chemistry | Physiology | 0.026 | -0.117 | 0.170 | 0.073 | 0.034 | -0.134 | 0.201 | 0.085 | 0.010 | 0.007 | 0.013 | 0.002
spleen weight | 1 | Immunophenotyping | Immunology | 0.187 | -0.050 | 0.425 | 0.121 | 0.113 | -0.160 | 0.387 | 0.140 | -0.154 | -0.210 | -0.098 | 0.029
st | 1 | Electrocardiogram (ECG) | Heart | 0.003 | -0.054 | 0.061 | 0.029 | -0.005 | -0.081 | 0.070 | 0.039 | -0.003 | -0.018 | 0.011 | 0.007
stroke volume | 1 | Echo | Heart | 0.059 | -0.078 | 0.197 | 0.070 | 0.157 | 0.009 | 0.306 | 0.076 | 0.094 | 0.078 | 0.110 | 0.008
tibia length | 1 | Heart Weight | Morphology | -0.148 | -0.440 | 0.145 | 0.149 | -0.137 | -0.426 | 0.151 | 0.147 | 0.010 | 0.006 | 0.013 | 0.002
total bilirubin | 1 | Clinical Chemistry | Physiology | 0.061 | -0.010 | 0.131 | 0.036 | 0.002 | -0.086 | 0.091 | 0.045 | -0.055 | -0.098 | -0.012 | 0.022
total cholesterol | 1 | Clinical Chemistry | Physiology | 0.094 | -0.075 | 0.264 | 0.086 | 0.314 | 0.113 | 0.516 | 0.103 | 0.203 | 0.175 | 0.230 | 0.014
total food intake | 1 | Indirect Calorimetry | Metabolism | -0.119 | -0.254 | 0.016 | 0.069 | -0.096 | -0.256 | 0.064 | 0.082 | 0.027 | -0.023 | 0.077 | 0.026
total protein | 1 | Clinical Chemistry | Physiology | -0.042 | -0.062 | -0.022 | 0.010 | -0.036 | -0.062 | -0.009 | 0.013 | 0.009 | -0.001 | 0.019 | 0.005
total water intake | 1 | Indirect Calorimetry | Metabolism | -0.146 | -0.237 | -0.054 | 0.047 | -0.210 | -0.268 | -0.151 | 0.030 | -0.065 | -0.137 | 0.007 | 0.037
triglycerides | 1 | Clinical Chemistry | Physiology | -0.032 | -0.123 | 0.059 | 0.047 | 0.327 | 0.209 | 0.445 | 0.060 | 0.347 | 0.259 | 0.436 | 0.045
urea (blood urea nitrogen - bun) | 1 | Clinical Chemistry | Physiology | -0.141 | -0.266 | -0.015 | 0.064 | -0.095 | -0.251 | 0.061 | 0.079 | 0.040 | 0.005 | 0.075 | 0.018
uric acid | 1 | Clinical Chemistry | Physiology | 0.037 | -0.066 | 0.139 | 0.052 | 0.363 | 0.091 | 0.634 | 0.138 | 0.447 | -0.080 | 0.975 | 0.269
white blood cell count | 1 | Hematology | Hematology | -0.091 | -0.170 | -0.011 | 0.041 | 0.117 | -0.002 | 0.236 | 0.061 | 0.198 | 0.137 | 0.259 | 0.031
whole arena average speed | 1 | Open Field | Behaviour | -0.016 | -0.086 | 0.054 | 0.036 | -0.114 | -0.184 | -0.044 | 0.036 | -0.100 | -0.152 | -0.048 | 0.027
whole arena resting time | 1 | Open Field | Behaviour | -0.053 | -0.101 | -0.005 | 0.025 | -0.059 | -0.108 | -0.011 | 0.025 | 0.005 | -0.051 | 0.061 | 0.029

## 6.4 Table 4

    # PDF pander(overall2, caption = 'Table 6.4: Summary of overall meta-analyses on
    # the functional trait group level (GroupingTerm). Results for lnCVR, lnVR and
    # lnRR and their respective upper and lower 95 percent CI's, standard error and
    # I2 values are provided.')

    # HTML
    kable(overall2, digits = 3, caption = "Summary of overall meta-analyses on the functional trait group level (GroupingTerm). Results for lnCVR, lnVR and lnRR and their respective upper and lower 95 percent CI's, standard error and I2 values are provided.") %>%
        kable_styling() %>% scroll_box(width = "100%", height = "200px")

Summary of overall meta-analyses on the functional trait group level (GroupingTerm). Results for lnCVR, lnVR and lnRR and their respective upper and lower 95 percent CI’s, standard error and I2 values are provided.
GroupingTerm | lnCVR | lnCVR_lower | lnCVR_upper | lnCVR_se | lnCVR_I2 | lnVR | lnVR_lower | lnVR_upper | lnVR_se | lnVR_I2 | lnRR | lnRR_lower | lnRR_upper | lnRR_se | lnRR_I2
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
Behaviour | -0.003504851 | -0.024 | 0.017 | 0.010 | 39.033 | -0.017834494 | -0.074 | 0.038 | 0.029 | 92.799 | -0.019920623 | -0.063 | 0.024 | 0.022 | 89.238
Morphology | 0.077445258 | 0.041 | 0.113 | 0.018 | 67.282 | 0.151417093 | 0.082 | 0.221 | 0.035 | 90.626 | 0.067816016 | 0.007 | 0.128 | 0.031 | 99.698
Metabolism | -0.043083123 | -0.113 | 0.026 | 0.035 | 84.801 | 0.091060863 | -0.034 | 0.216 | 0.064 | 96.783 | 0.142257666 | 0.036 | 0.248 | 0.054 | 99.377
Physiology | 0.012679189 | -0.014 | 0.039 | 0.014 | 70.566 | 0.035982066 | -0.028 | 0.100 | 0.033 | 90.643 | 0.016369478 | -0.044 | 0.077 | 0.031 | 99.817
Immunology | -0.068181720 | -0.098 | -0.038 | 0.015 | 40.769 | -0.111238224 | -0.162 | -0.060 | 0.026 | 58.020 | -0.057483957 | -0.107 | -0.008 | 0.025 | 96.148
Hematology | 0.021786545 | -0.017 | 0.060 | 0.020 | 60.753 | 0.080211058 | 0.032 | 0.129 | 0.025 | 64.943 | 0.038853737 | -0.002 | 0.080 | 0.021 | 99.630
Heart | 0.018383893 | -0.013 | 0.050 | 0.016 | 81.298 | -0.005081004 | -0.036 | 0.026 | 0.016 | 74.058 | -0.004893287 | -0.032 | 0.023 | 0.014 | 98.985
Hearing | 0.015730245 | -0.011 | 0.043 | 0.014 | 21.756 | 0.010685756 | -0.023 | 0.044 | 0.017 | 22.797 | -0.013236557 | -0.034 | 0.007 | 0.010 | 79.486
Eye | -0.081793214 | -0.148 | -0.016 | 0.034 | 49.810 | -0.074449703 | -0.138 | -0.011 | 0.032 | 48.309 | 0.009118622 | 0.001 | 0.017 | 0.004 | 90.528
All | 0.004655293 | -0.009 | 0.018 | 0.007 | 76.486 | 0.015663425 | -0.008 | 0.039 | 0.012 | 90.760 | 0.012433191 | -0.006 | 0.031 | 0.009 | 99.716

## 6.5 Table 5

    female <- (overall.female.plot3)
    female$Sex <- c("female")

    male <- (overall.male.plot3)
    male$Sex <- c("male")

    mf_sig <- rbind(female, male)
    mf_sig <- mf_sig[, c(1, 17, 2:16)]

    mf_sig <- mf_sig %>% arrange(desc(GroupingTerm))

    # PDF pander(mf_sig, caption = 'Table 6.5: Provides an overview of meta-analysis
    # results performed on traits that were significantly biased towards either sex.
    # This table summarizes findings for both sexes and the respective functional
    # trait groups.') HTML
    kable(mf_sig, digits = 3, caption = "Provides an overview of meta-analysis results performed on traits that were significantly biased towards either sex. This table summarizes findings for both sexes and the respective functional trait groups.") %>%
        kable_styling() %>% scroll_box(width = "100%", height = "300px")

Provides an overview of meta-analysis results performed on traits that were significantly biased towards either sex. This table summarizes findings for both sexes and the respective functional trait groups.
GroupingTerm | Sex | lnCVR | lnCVR_lower | lnCVR_upper | lnCVR_se | lnCVR_I2 | lnVR | lnVR_lower | lnVR_upper | lnVR_se | lnVR_I2 | lnRR | lnRR_lower | lnRR_upper | lnRR_se | lnRR_I2
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
Behaviour | female | -0.0531306677583289 | -0.10116634434615 | -0.00509499117050839 | 0.0245084486075866 | 0 | -0.093 | -0.120 | -0.066 | 0.014 | 42.139 | -0.094 | -0.117 | -0.072 | 0.011 | 0.000
Behaviour | male | 0.107425386004474 | 0.0498612598581906 | 0.164989512150757 | 0.0293699917959421 | 0 | 0.284 | 0.189 | 0.378 | 0.048 | 37.808 | 0.203 | 0.135 | 0.272 | 0.035 | 0.000
Morphology | female | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | -0.110 | -0.230 | 0.011 | 0.061 | 93.935
Morphology | male | 0.127090489394448 | 0.0862246543567037 | 0.167956324432192 | 0.0208502989647201 | 41.5455829967945 | 0.217 | 0.155 | 0.280 | 0.032 | 84.378 | 0.120 | 0.058 | 0.183 | 0.032 | 99.674
Metabolism | female | -0.124528442700272 | -0.166238188601915 | -0.0828186967986301 | 0.0212808736439258 | 13.3240051843312 | -0.210 | -0.268 | -0.151 | 0.030 | 0.000 | NA | NA | NA | NA | NA
Metabolism | male | 0.0939139877083824 | 0.0519967178471424 | 0.135831257569622 | 0.0213867551607469 | 0 | 0.224 | 0.100 | 0.348 | 0.063 | 95.761 | 0.223 | 0.110 | 0.335 | 0.058 | 98.647
Physiology | female | -0.0471784060711916 | -0.0640130203407664 | -0.0303437918016167 | 0.00858924674247289 | 0.610173416856019 | -0.163 | -0.296 | -0.030 | 0.068 | 92.843 | -0.117 | -0.221 | -0.014 | 0.053 | 99.301
Physiology | male | 0.0963807810373384 | 0.0630394628524223 | 0.129722099222255 | 0.0170111892095509 | 0.109385385672822 | 0.240 | 0.159 | 0.320 | 0.041 | 56.006 | 0.145 | 0.071 | 0.218 | 0.038 | 98.518
Immunology | female | -0.090049165480264 | -0.148089504316368 | -0.03200882664416 | 0.0296129619186469 | 81.548126644761 | -0.180 | -0.218 | -0.141 | 0.020 | 0.000 | -0.159 | -0.229 | -0.090 | 0.035 | 94.729
Immunology | male | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.028 | 0.001 | 0.055 | 0.014 | 44.587
Hematology | female | -0.0920412245428073 | -0.150194500108348 | -0.0338879489772661 | 0.0296705837577868 | 0 | NA | NA | NA | NA | NA | -0.013 | -0.018 | -0.007 | 0.003 | 84.990
Hematology | male | 0.0824361263690887 | 0.0482041476047534 | 0.116668105133424 | 0.0174656162227229 | 0.0891085957461772 | 0.144 | 0.046 | 0.243 | 0.050 | 86.350 | 0.116 | 0.027 | 0.205 | 0.046 | 98.445
Heart | female | -0.0911918652110014 | -0.17059643028354 | -0.0117873001384624 | 0.0405132776412588 | 66.5484415178361 | -0.084 | -0.109 | -0.059 | 0.013 | 0.528 | -0.108 | -0.179 | -0.036 | 0.036 | 98.814
Heart | male | 0.154171112862654 | 0.0851206825435622 | 0.223221543181746 | 0.0352304587552389 | 78.6992182184384 | 0.104 | 0.062 | 0.146 | 0.022 | 36.412 | 0.058 | 0.041 | 0.075 | 0.009 | 84.650
Hearing | female | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | -0.035 | -0.062 | -0.009 | 0.014 | 61.362
Hearing | male | NA | NA | NA | NA | NA | 0.087 | 0.007 | 0.167 | 0.041 | 0.000 | 0.018 | 0.006 | 0.031 | 0.006 | 0.000
Eye | female | -0.164020755634996 | -0.249857309297975 | -0.0781842019720163 | 0.0437949647748872 | 3.0711311258835 | -0.166 | -0.277 | -0.056 | 0.057 | 24.076 | NA | NA | NA | NA | NA
Eye | male | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.019 | 0.004 | 0.034 | 0.008 | 97.046
All | female | -0.0866926973518197 | -0.109246316493763 | -0.0641390782098763 | 0.011507159988573 | 66.4794235647273 | -0.130 | -0.157 | -0.103 | 0.014 | 75.919 | -0.098 | -0.129 | -0.068 | 0.016 | 99.612
All | male | 0.111341926196315 | 0.0930887040851348 | 0.129595148307494 | 0.00931303955335859 | 34.0698195492482 | 0.198 | 0.162 | 0.233 | 0.018 | 86.497 | 0.111 | 0.083 | 0.138 | 0.014 | 99.639

## 6.6 Table 6

    het <- heterog6 %>% pivot_wider(names_from = parameter, values_from = value:ci.h) %>%
        rename(lnCVR = value_heteroCVR, lnCVR_lower = ci.low_heteroCVR, lnCVR_higher = ci.high_heteroCVR,
            lnVR = value_heteroVR, lnVR_lower = ci.low_heteroVR, lnVR_higher = ci.high_heteroVR,
            lnRR = value_heteroRR, lnRR_lower = ci.low_heteroRR, lnRR_higher = ci.high_heteroRR,
            CVR = mean_heteroCVR, VR = mean_heteroVR, RR = mean_heteroRR, CVR_lower = ci.l_heteroCVR,
            VR_lower = ci.l_heteroVR, RR_lower = ci.l_heteroRR, CVR_higher = ci.h_heteroCVR,
            VR_higher = ci.h_heteroVR, RR_higher = ci.h_heteroRR) %>% select(1, 2, 5,
        8, 3, 6, 9, 4, 7, 10, 14, 17, 20, 15, 18, 21, 16, 19, 22)

    # PDF pander(het, caption = 'Table 6.6: Summarizes our findings on heterogeneity
    # due to institutions and mouse strains. These results are based on meta-analyses
    # on sigma^2 and errors for mouse strains and centers (Institutions), following
    # the identical workflow from above.')

    # HTML
    kable(het, digits = 3, caption = "Summarizes our findings on heterogeneity due to institutions and mouse strains. These results are based on meta-analyses on sigma^2 and errors for mouse strains and centers (Institutions), following the identical workflow from above.") %>%
        kable_styling() %>% scroll_box(width = "100%", height = "200px")

Summarizes our findings on heterogeneity due to institutions and mouse strains. These results are based on meta-analyses on sigma^2 and errors for mouse strains and centers (Institutions), following the identical workflow from above.
GroupingTerm | lnCVR | lnCVR_lower | lnCVR_higher | lnVR | lnVR_lower | lnVR_higher | lnRR | lnRR_lower | lnRR_higher | CVR | CVR_lower | CVR_higher | VR | VR_lower | VR_higher | RR | RR_lower | RR_higher
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
Behaviour | -1.649 | -2.464 | -0.835 | -1.545 | -2.359 | -0.730 | -2.640 | -4.135 | -1.145 | 0.192 | 0.085 | 0.434 | 0.213 | 0.094 | 0.482 | 0.071 | 0.016 | 0.318
Morphology | -1.745 | -2.424 | -1.066 | -1.725 | -2.404 | -1.046 | -7.277 | -10.004 | -4.549 | 0.175 | 0.089 | 0.344 | 0.178 | 0.090 | 0.351 | 0.001 | 0.000 | 0.011
Metabolism | -2.404 | -3.567 | -1.240 | -2.897 | -4.260 | -1.534 | -2.970 | -4.134 | -1.806 | 0.090 | 0.028 | 0.289 | 0.055 | 0.014 | 0.216 | 0.051 | 0.016 | 0.164
Physiology | -2.899 | -3.762 | -2.037 | -1.922 | -2.688 | -1.155 | -2.664 | -3.430 | -1.898 | 0.055 | 0.023 | 0.130 | 0.146 | 0.068 | 0.315 | 0.070 | 0.032 | 0.150
Immunology | -0.291 | -0.795 | 0.213 | -0.231 | -0.735 | 0.273 | -0.895 | -1.790 | 0.001 | 0.748 | 0.452 | 1.238 | 0.793 | 0.479 | 1.313 | 0.409 | 0.167 | 1.001
Hematology | -0.463 | -1.120 | 0.193 | -0.369 | -1.026 | 0.287 | -1.287 | -2.220 | -0.354 | 0.629 | 0.326 | 1.213 | 0.691 | 0.359 | 1.333 | 0.276 | 0.109 | 0.702
Heart | -3.573 | -4.272 | -2.874 | -2.980 | -3.549 | -2.411 | -3.877 | -4.443 | -3.311 | 0.028 | 0.014 | 0.056 | 0.051 | 0.029 | 0.090 | 0.021 | 0.012 | 0.036
Hearing | -3.059 | -5.049 | -1.068 | -2.361 | -4.352 | -0.370 | -3.602 | -5.592 | -1.611 | 0.047 | 0.006 | 0.344 | 0.094 | 0.013 | 0.690 | 0.027 | 0.004 | 0.200
Eye | -3.707 | -5.102 | -2.312 | -3.228 | -4.331 | -2.125 | -5.773 | -6.574 | -4.972 | 0.025 | 0.006 | 0.099 | 0.040 | 0.013 | 0.119 | 0.003 | 0.001 | 0.007
All | -2.290 | -2.648 | -1.932 | -1.936 | -2.254 | -1.619 | -3.452 | -3.944 | -2.960 | 0.101 | 0.071 | 0.145 | 0.144 | 0.105 | 0.198 | 0.032 | 0.019 | 0.052

# 7 R Session Information

    sessionInfo() %>% pander()

R version 4.0.3 (2020-10-10)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: splus2R(v.1.2-2), pander(v.0.6.3), knitr(v.1.30), here(v.0.1), png(v.0.1-7), ggpubr(v.0.4.0), robumeta(v.2.0), kableExtra(v.1.2.1), forcats(v.0.5.0), stringr(v.1.4.0), tidyr(v.1.1.2), tibble(v.3.0.4), ggplot2(v.3.3.2), tidyverse(v.1.3.0), purrr(v.0.3.4), devtools(v.2.3.2), usethis(v.1.6.3), metafor(v.2.4-0), Matrix(v.1.2-18), dplyr(v.1.0.2), readr(v.1.4.0) and pacman(v.0.5.1)

loaded via a namespace (and not attached): nlme(v.3.1-149), fs(v.1.5.0), lubridate(v.1.7.9), RColorBrewer(v.1.1-2), webshot(v.0.5.2), httr(v.1.4.2), rprojroot(v.1.3-2), tools(v.4.0.3), backports(v.1.1.10), utf8(v.1.1.4), R6(v.2.4.1), mgcv(v.1.8-33), DBI(v.1.1.0), colorspace(v.1.4-1), withr(v.2.3.0), tidyselect(v.1.1.0), prettyunits(v.1.1.1), processx(v.3.4.4), curl(v.4.3), compiler(v.4.0.3), cli(v.2.0.2), rvest(v.0.3.6), formatR(v.1.7), xml2(v.1.3.2), desc(v.1.2.0), labeling(v.0.3), scales(v.1.1.1), callr(v.3.5.0), digest(v.0.6.25), foreign(v.0.8-80), rmarkdown(v.2.4), rio(v.0.5.16), pkgconfig(v.2.0.3), htmltools(v.0.5.0), sessioninfo(v.1.1.1), highr(v.0.8), dbplyr(v.1.4.4), rlang(v.0.4.8), readxl(v.1.3.1), rstudioapi(v.0.11), farver(v.2.0.3), generics(v.0.0.2), jsonlite(v.1.7.1), zip(v.2.1.1), car(v.3.0-10), magrittr(v.1.5), Rcpp(v.1.0.5), munsell(v.0.5.0), fansi(v.0.4.1), abind(v.1.4-5), lifecycle(v.0.2.0), stringi(v.1.5.3), yaml(v.2.2.1), carData(v.3.0-4), pkgbuild(v.1.1.0), blob(v.1.2.1), crayon(v.1.3.4), lattice(v.0.20-41), splines(v.4.0.3), cowplot(v.1.1.0), haven(v.2.3.1), hms(v.0.5.3), ps(v.1.4.0), pillar(v.1.4.6), ggsignif(v.0.6.0), codetools(v.0.2-16), pkgload(v.1.1.0), reprex(v.0.3.0), glue(v.1.4.2), evaluate(v.0.14), data.table(v.1.13.0), remotes(v.2.2.0), modelr(v.0.1.8), vctrs(v.0.3.4), testthat(v.2.3.2), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), xfun(v.0.18), openxlsx(v.4.2.2), broom(v.0.7.1), rstatix(v.0.6.0), viridisLite(v.0.3.0), memoise(v.1.1.0) and ellipsis(v.0.3.1)