Lecture 4: topics

• Introduction to MetaboAnalyst software
  • A useful R-based resources for metabolomics
• Elements of statistical Power Analysis

• We will use functions from packages base, utils, and stats (pre-installed and pre-loaded)
• We may also use the packages below (specifying package::function for clarity).

# Load pckgs for this R session

# --- General 
library(here)     # tools find your project's files, based on working directory
library(dplyr)    # A Grammar of Data Manipulation
library(skimr)    # Compact and Flexible Summaries of Data
library(magrittr) # A Forward-Pipe Operator for R 
library(readr)    # A Forward-Pipe Operator for R 

# Plotting & data visualization
library(ggplot2)      # Create Elegant Data Visualisations Using the Grammar of Graphics
library(ggfortify)     # Data Visualization Tools for Statistical Analysis Results
library(scatterplot3d) # 3D Scatter Plot

# --- Statistics
library(MASS)       # Support Functions and Datasets for Venables and Ripley's MASS
library(factoextra) # Extract and Visualize the Results of Multivariate Data Analyses
library(FactoMineR) # Multivariate Exploratory Data Analysis and Data Mining
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests

# --- Tidymodels (meta package)
library(rsample)    # General Resampling Infrastructure  
library(broom)      # Convert Statistical Objects into Tidy Tibbles

Name: Biopsy Data on Breast Cancer Patients
Documentation: See reference on the data downloaded and conditioned for R here https://cran.r-project.org/web/packages/MASS/MASS.pdf
Sampling details: This breast cancer database was obtained from the University of Wisconsin Hospitals, Madison from Dr. William H. Wolberg. He assessed biopsies of breast tumours for 699 patients up to 15 July 1992; each of nine attributes has been scored on a scale of 1 to 10, and the outcome is also known. The dataset contains the original Wisconsin breast cancer data with 699 observations on 11 variables.

The biopsy data contains 699 observations of 11 variables.

The dataset also contains a character variable: ID, and a factor variable: class, with two levels (“benign” and “malignant”).

# check variable types
str(biopsy)

'data.frame':   699 obs. of  11 variables:
 $ ID   : chr  "1000025" "1002945" "1015425" "1016277" ...
 $ V1   : int  5 5 3 6 4 8 1 2 2 4 ...
 $ V2   : int  1 4 1 8 1 10 1 1 1 2 ...
 $ V3   : int  1 4 1 8 1 10 1 2 1 1 ...
 $ V4   : int  1 5 1 1 3 8 1 1 1 1 ...
 $ V5   : int  2 7 2 3 2 7 2 2 2 2 ...
 $ V6   : int  1 10 2 4 1 10 10 1 1 1 ...
 $ V7   : int  3 3 3 3 3 9 3 3 1 2 ...
 $ V8   : int  1 2 1 7 1 7 1 1 1 1 ...
 $ V9   : int  1 1 1 1 1 1 1 1 5 1 ...
 $ class: Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...

There is also one incomplete variable V6

• remember the package skimr for exploring a dataframe?

# check if vars have missing values
biopsy %>% 
  # select only variables starting with "V"
  skimr::skim(starts_with("V")) %>%
  dplyr::select(skim_variable, 
                n_missing)

# A tibble: 9 × 2
  skim_variable n_missing
  <chr>             <int>
1 V1                    0
2 V2                    0
3 V3                    0
4 V4                    0
5 V5                    0
6 V6                   16
7 V7                    0
8 V8                    0
9 V9                    0

A 3D scatterplot of observations shows the first 3 principal components’ scores.

• For this one, we need the scatterplot3d() function of the scatterplot3d package;
• The color argument assigned to the Label variable;
• To add a legend, we use the legend() function and specify its coordinates via the xyz.convert() function.

# 3D scatterplot ...
plot_3d <- with(PC_scores, 
                scatterplot3d::scatterplot3d(PC_scores$PC1, 
                                             PC_scores$PC2, 
                                             PC_scores$PC3, 
                                             color = as.numeric(Label), 
                                             pch = 19, 
                                             main ="Figure 3: 3D Scatter Plot", 
                                             xlab="PC1",
                                             ylab="PC2",
                                             zlab="PC3"))

# ... + legend
legend(plot_3d$xyz.convert(0.5, 0.7, 0.5), 
       pch = 19, 
       yjust=-0.6,
       xjust=-0.9,
       legend = levels(PC_scores$Label), 
       col = seq_along(levels(PC_scores$Label)))

Biplots have two key elements: scores (the 2 axes) and loadings (the vectors). As in the scores plot, each point represents an observation projected in the space of principal components where:

• Biopsies of the same class are located closer to each other, which indicates that they have similar scores referred to the 2 main principal components;
• The loading vectors show strength and direction of association of original variables with new PC variables.

As expected from PCA, the single PC1 accounts for variance in almost all original variables, while V9 has the major projection along PC2.

• Power analysis helps with the key question How many observations/subjects do I need for my experiment? (= \(n\))

  • Too small of a sample size can under detect the effect of interest in your experiment
  • Too large of a sample size may lead to unnecessary wasting of resources
  • We strive to have just the sufficient number of observations needed to have a good chance of detecting the effect researched. (Even more so in a very time-consuming or expensive experiment.)

• When should we do power analysis?

  • (Ideally), before the experiment: a priori power analysis allows to determine the necessary sample size \(n\) of a test, given a desired \(\alpha\) level, a desired power level (\(1- \beta\)), and the size of the effect to be detected (a measure of difference between \(H_o\) and \(H_1\))
  • In reality, sometimes you can only do post-hoc power analysis after the experiment, so the sample size \(n\) is already given.
    • In this case, given \(n\), \(\alpha\), and a specified effect size, the analysis will return the power (\(1- \beta\)) of the test, or \(\beta\) (i.e. the probability of Type II error = incorrectly retaining \(H_o\)).

• A specified effect size (i.e. the minimum deviation from \(H_o\) that you hope to detect for a meaningful result)
  • The larger the effect size, the easier it is to detect an effect and require fewer obs
• • As \(standard deviation\) gets bigger, it is harder to detect a significant difference, so you’ll need a bigger sample size.

• \(\alpha\) is the significance level of the test (i.e. the probability of incorrectly rejecting the null hypothesis (a false positive).
  • Understanding if the test is one-tailed (difference has a direction) or two-tailed
• \(\beta\) is the probability of accepting the null hypothesis, even though it is false (a false negative), when the real difference is equal to the minimum effect size.
  • \(1- \beta\) is the power of a test is the probability of correctly rejecting the null hypothesis (getting a significant result) when the real difference is equal to the minimum effect size.
    • a power of 80% (equivalent to a beta of 20%) is probably the most common, while some people use 50% or 90%

So (since \(\alpha\) and \(1-\beta\) are normally set) the key piece of information we need is the effect size, which is essentially a function of the difference between the means of the null and alternative hypotheses over the variation (standard deviation) in the data.

The tricky part is that effect size is related to biological/practical significance rather than statistical significance

How should you estimate a meaningful Effect Size?

• Use preliminary information in the form of pilot study
• Use background information in the form of similar studies in the literature
• (With no prior information), make an estimated guess on the effect size expected (see guidelines next)

Most R functions for sample size only allow you to enter effect size as input

GOAL: Imagine this is a pilot study, in which we tested fish is (on average) different form 20 cm in length.

The guanapo_data dataset contains information on fish lengths from the Guanapo river pilot

# Load data on river fish length 
fishlength_data <- readr::read_csv(here::here("practice", "data_input", "04_datasets", 
                                              "fishlength.csv"),
                              show_col_types = FALSE)

# Select a portion of the data (i.e. out "pilot" experiment) 
guanapo_data <- fishlength_data %>% 
  dplyr::filter(river == "Guanapo")

# Pilot experiment data 
names(guanapo_data)

[1] "id"     "river"  "length"

mean_H1 <-  mean(guanapo_data$length) # 18.29655
mean_H1

[1] 18.29655

sd_sample <- sd(guanapo_data$length)  # 2.584636
sd_sample

[1] 2.584636

Let’s compute the one sample t-test with stats::t.test against a hypothetical average fish length (\(mean\_H_o = 20\) )

# Hypothetical fish population length mean (H0)
mean_H0 <- 20
# one-sample mean t-test 
t_stat <- stats::t.test(x = guanapo_data$length,
                        mu = mean_H0,
                        alternative = "two.sided")
# one-sample t-test results
t_stat

    One Sample t-test

data:  guanapo_data$length
t = -3.5492, df = 28, p-value = 0.001387
alternative hypothesis: true mean is not equal to 20
95 percent confidence interval:
 17.31341 19.27969
sample estimates:
mean of x 
 18.29655 

• There appear to be a statistically significant result here: the mean length of the fish appears to be different from 20 cm.

QUESTION: In a new study of the same fish, what sample size n would you need to get a comparable result?

• We input Cohen’s d (after calculating it manually) following: \(effect\ size\ \approx \frac{{Mean}_{H_1}\ -{\ Mean}_{H_0}}{Std\ Dev}\)

• We use pwr::pwr.t.test to calculate the minimum sample size n required:

# Cohen's d formula 
eff_size <- (mean_H1 - mean_H0) / sd_sample # -0.6590669

# power analysis to actually calculate the minimum sample size required:
pwr::pwr.t.test(d = eff_size, 
                sig.level = 0.05, 
                power = 0.8,
                type = "one.sample")

     One-sample t test power calculation 

              n = 20.07483
              d = 0.6590669
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

We would need n = 21 (rounding up) observations for an experiment (e.g. in different river) to detect an effect size as the pilot study at a 5% significance level and 80% power.

Let’s look at the entire fishlength_data with the lengths of fish from 2 separate rivers.

# Explore complete data 
fishlength_data %>% 
  dplyr::group_by (river) %>% 
  dplyr::summarise (N = n(), 
                    mean_len = mean(length),
                    sd_len = sd(length)) 

# A tibble: 2 × 4
  river       N mean_len sd_len
  <chr>   <int>    <dbl>  <dbl>
1 Aripo      39     20.3   1.78
2 Guanapo    29     18.3   2.58

Visualize quickly the 2 samples (rivers) with a boxplot

# visualize the data
fishlength_data %>% 
  ggplot(aes(x = river, y = length)) +
  geom_boxplot()

The cortisol_data dataset contains information about cortisol levels measured on 20 participants in the morning and evening

# Load data 
cortisol_data <- read.csv(file = here::here("practice", "data_input", "04_datasets", 
                                        "cortisol.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

# Explore data 
names(cortisol_data)

[1] "patient_id" "time"       "cortisol"  

cortisol_data %>% 
  dplyr::group_by (time) %>% 
  dplyr::summarise (
    N = n(), 
    mean_cort = mean(cortisol),
    sd_cort = sd(cortisol)) 

# A tibble: 2 × 4
  time        N mean_cort sd_cort
  <chr>   <int>     <dbl>   <dbl>
1 evening    20      197.    87.5
2 morning    20      313.    73.8

Notice the difference in the paired sample means is quite large

GOAL: Flipping the question, if we know the given n (20 patients observed twice): How big should the effect size be to be detected at power of 0.8 and significance level 0.05?

• We use pwr::pwr.t.test, with the argument specification type = "paired", but this time to estimate the effect size

# power analysis to actually calculate the effect size at the desired conditions:
pwr::pwr.t.test(n = 20, 
                #d =  eff_size, 
                sig.level = 0.05, 
                power = 0.8,
                type = "paired")

     Paired t test power calculation 

              n = 20
              d = 0.6604413
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number of *pairs*

The functions returns the effect size (Cohen’s metric): d = 0.6604413. So, with this experimental design we would be able to detect a medium-large effect size.

The mussels_data dataset contains information about the length of the anterior adductor muscle scar in the mussel Mytilus trossulus across five locations around the world!

# Load data 
mussels_data <- read.csv(file = here::here("practice", "data_input", "04_datasets", 
                                        "mussels.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

# Explore data 
names(mussels_data)

[1] "length"   "location"

stats <- mussels_data %>% 
  dplyr::group_by (location) %>% 
  dplyr::summarise (
    N = n(), 
    mean_len = mean(length),
    sd_len = sd(length)) 

stats

# A tibble: 5 × 4
  location       N mean_len  sd_len
  <chr>      <int>    <dbl>   <dbl>
1 Magadan        8   0.0780 0.0129 
2 Newport        8   0.0748 0.00860
3 Petersburg     7   0.103  0.0162 
4 Tillamook     10   0.0802 0.0120 
5 Tvarminne      6   0.0957 0.0130 

There appears to be a noticeable difference in lenght at average measurements at least between some of the locations

Assuming we verified the required assumptions, let’s run the ANOVA test to confirm the visual intuition

• With the stats::aov followed by the command summary

# Summary of test outputs: 
summary_ANOVA <- summary(stats::aov(length ~ location,
                   data = mussels_data))

# From which we extract all the output elements 
# F value 
summary_ANOVA[[1]][["F value"]] # 7.121019

[1] 7.121019       NA

# p value 
summary_ANOVA[[1]][["Pr(>F)"]]  # 0.0002812242

[1] 0.0002812242           NA

# df of numerator and denominator
summary_ANOVA[[1]][["Df"]]      # 4, 34 

[1]  4 34

# Sum of Square BETWEEN groups
SSB <- summary_ANOVA[[1]][["Sum Sq"]][1]  # 0.004519674
# Sum of Square WITHIN groups
SSW <- summary_ANOVA[[1]][["Sum Sq"]][2]  # 0.005394906

• A one-way ANOVA test confirms that the mean lengths of muscle scar differed significantly between locations ( F = 7.121, with df = [4, 34], and p = 0.000281).

In ANOVA it may be tricky to decide what kind of effect size we are looking for:

• if we care about an overall significance test, the sample size needed is a function of the standard deviation of the group means
• if we’re interested in the comparisons of means, there are other ways of expressing the effect size (e.g. a difference between the smallest and largest means)

Here let’s consider an overall test in which we could reasonably collect the same n. of observations in each group

n_loc <- nrow(stats)

means_by_loc <- c(0.0780, 0.0748, 0.103, 0.0802, 0.0957)
overall_mean <-  mean(means_by_loc)
sd_by_loc <- c(0.0129, 0.00860, 0.0162, 0.0120, 0.0130)
overall_sd <-  mean(sd_by_loc)

# Effect Size f formula
Cohen_f = sqrt( sum( (1/n_loc) * (means_by_loc - overall_mean)^2) ) /overall_sd
Cohen_f # EXTREMELY BIG 

[1] 0.877622

# Power analysis with given f 
pwr::pwr.anova.test(k = n_loc,
                    n = NULL,
                    f = Cohen_f,
                    sig.level = 0.05,
                    power = 0.80)

     Balanced one-way analysis of variance power calculation 

              k = 5
              n = 4.166759
              f = 0.877622
      sig.level = 0.05
          power = 0.8

NOTE: n is number in each group

The n output ( = 5 observations per group) -as opposed to >6 per group- would be sufficient if we wanted to confidently detect the difference observed in the previous study

The ideas covered before apply also to linear models, although here:

• we use pwr.f2.test() to do the power calculation
• the effect sizes (\(f^2\)) is based on \(R^2\)

\[ f^2=\ \frac{R^2}{1-\ R^2}\]

# define the linear model
lm_mussels <- lm(length ~ location, 
                 data = mussels_data)

# summarise the model
summary(lm_mussels)

From the linear model we get that the \(R^2\) value is 0.4559 and we can use this to calculate Cohen’s \(f^2\) value using the formula

# Extract R squared
R_2 <- summary(lm_mussels)$r.squared
# compute f squared
f_2 <- R_2 / (1 - R_2)
f_2

[1] 0.837767

Our model has 5 parameters (because we have 5 groups) and so the numerator degrees of freedom \(u\) will be 4 (5−1=4).

Hence, we carry out the power analysis with the function pwr.f2.test:

# power analysis for overall linear model 
pwr::pwr.f2.test(u = 4, v = NULL, 
                 f2 = 0.8378974,
                 sig.level = 0.05 , power = 0.8)

     Multiple regression power calculation 

              u = 4
              v = 14.62182
             f2 = 0.8378974
      sig.level = 0.05
          power = 0.8

Recall that, in the F statistic evaluating the model,

• u the df for the numerator: \(df_{between} =k−1 = 5-1 = 4\)
• v the df for the denominator: \(df_{within} = n-k = ?\)
  • so \(n = v+5\)

pwr::pwr.f2.test(u = 4, f2 = 0.8378974,
            sig.level = 0.05 , power = 0.8)

     Multiple regression power calculation 

              u = 4
              v = 14.62182
             f2 = 0.8378974
      sig.level = 0.05
          power = 0.8

This tells us that the denominator degrees of freedom v should be 15 (14.62 rounded up), and this means that we would only need 20 observations n = v+5 in total across all 5 groups to detect this effect size

(Simplifying a little)

Let’s re-load a dataset from Lab # 3 (the NHANES dataset) for a quick demonstration of data splitting in an ML predictive modeling scenario

• We can try predicting BMI from age (in years), PhysActive, and gender, using linear regression model (which is a Supervised ML algorithm)

• (we already saved this dataset)

# (we already saved this dataset in our project folders)

# Use `here` in specifying all the subfolders AFTER the working directory 
nhanes <- read.csv(file = here::here("practice", "data_input", "03_datasets",
                                      "nhanes.samp.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

In this case the regression models serves for predicting numeric, continuous quantities

# fitting  linear regression model specification
lin_mod <- lm(BMI ~ Age + Gender + PhysActive, data = nhanes_train)

summary(lin_mod)

Call:
lm(formula = BMI ~ Age + Gender + PhysActive, data = nhanes_train)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.685  -4.674  -1.419   4.257  38.016 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   30.14217    1.30426  23.110  < 2e-16 ***
Age            0.01429    0.02198   0.650  0.51596    
Gendermale    -0.72960    0.72176  -1.011  0.31275    
PhysActiveYes -2.26539    0.73620  -3.077  0.00225 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.903 on 367 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.03416,   Adjusted R-squared:  0.02626 
F-statistic: 4.327 on 3 and 367 DF,  p-value: 0.005155

Using the above model, we can predict the BMI for different individuals (those left in the testing data)

• with the function predict, where we specify the argument newdata = nhanes_test)
• adding the prediction interval (the 95% CI), which gives uncertainty around a single value of the prediction

# Obtain predictions from the results of a model fitting function
pred_bmi <- stats::predict(lin_mod, 
               newdata = nhanes_test,
               interval = "confidence" )
head(pred_bmi)

       fit      lwr      upr
1 28.59148 27.33499 29.84797
2 27.70464 26.45051 28.95878
3 30.88546 29.72888 32.04203
4 28.01911 26.64955 29.38867
5 29.78421 28.04027 31.52815
6 27.60459 26.24230 28.96688

The ultimate goal of holding data back from the model training process was to evaluate its predictive performance on new data.

A common measure used is the RMSE (Root Mean Square Error) = a measure of the distance between observed values and predicted values in the testing dataset

# Computing the Root Mean Square Error
RMSE_test <- sqrt(mean((nhanes_test$BMI - predict(lin_mod, nhanes_test))^2, na.rm = T))
RMSE_test # 6.170518

[1] 6.170518

The RMSE (= 6.170518) tells us, (roughly speaking) by how much, on average, the new observed BMI values differ from those predicted by our model