Introduction to R

Integrative mini-project

Dr. Yu Cheng Hsu

Learning goal

  • Design a simple analysis pipeline from an experiment description
  • Apply R to clean data and run (with guidance) a paired comparison
  • Interpret what the test is asking—not memorize every formula

Background: Darwin’s Experiment

  • Source: The Effects of Cross and Self Fertilization in the Vegetable Kingdom by Charles Darwin in 1876.
  • Hypothesis: Cross-pollination could make their offspring stronger (in terms of height).
  • Experiment design
    • 16 pairs of corn,
      • one by self-pollination
      • one by cross-pollination.
    • Each pair of observations was grown in the same pots,
    • 4 pots in this study.

Illustration of experiment design

Illustration of experiment design

The Data Science Workflow

Before diving into the code, let’s look at the roadmap. Every robust data science project follows a structured life cycle.

As a professional scientist to

  • Translate the experiment into an actionable analysis plan
    • Outline the potential plan for
      • Processing
      • Analysis
      • Visualization
  • Convey the scientific findings

A full picture of analytics

Griffin et al. (2018) proposed best practices in the life sciences. In this module, so far, we have covered a simple traditional approach to research data.

Data lifecycle, figure from Griffin et al. (2018)

Data lifecycle, figure from Griffin et al. (2018)
Life Cycle Stage Example
Collecting Read the data
Integrating -
Processing Manipulating tables
Analyzing Reporting statistics
Publishing Generating report and figures

Project Roadmap

Some potential steps we need to perform:

  1. Data Import & Inspection: Load the experiment data and ensure data integrity.
  2. Data Manipulation (Processing): Calculate the height difference, capturing the paired experimental design.
  3. Exploratory Data Analysis (EDA): Visualize the paired nature of the plants and the distribution of their height differences.
  4. Statistical test: Formulate our hypothesis and conduct statistical inference
  5. Conclusion & Reporting: Interpret the statistical output (p-value, mean difference).

Data Import and Inspection

Before any analysis, the absolute first step is to inspect the data to ensure its integrity.

Something you need check before you do anything:

  • Missing values
  • Typos and human errors
  • Data structures
# Load necessary libraries and the dataset
library(tidyverse)
zea_mays =read_csv("https://raw.githubusercontent.com/ychsu987/ychsu987.github.io/refs/heads/main/assets/data/real-example/zea_mays_clean.csv")
head(zea_mays)
# A tibble: 6 × 4
  Pots   Pair Cross_fertilized_plant Self_fertilized_plant
  <chr> <dbl>                  <dbl>                 <dbl>
1 Pot1      1                   23.5                  17.4
2 Pot1      2                   12                    20.4
3 Pot1      3                   21                    20  
4 Pot2      4                   22                    20  
5 Pot2      5                   19.1                  18.4
6 Pot2      6                   21.5                  18.6

Inspecting Data

# Check the data structure and look for missing values
## Apparently you can use zea_mays %>% anyNA() and so on
str(zea_mays)
spc_tbl_ [17 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ Pots                  : chr [1:17] "Pot1" "Pot1" "Pot1" "Pot2" ...
 $ Pair                  : num [1:17] 1 2 3 4 5 6 7 8 9 10 ...
 $ Cross_fertilized_plant: num [1:17] 23.5 12 21 22 19.1 ...
 $ Self_fertilized_plant : num [1:17] 17.4 20.4 20 20 18.4 ...
 - attr(*, "spec")=
  .. cols(
  ..   Pots = col_character(),
  ..   Pair = col_double(),
  ..   Cross_fertilized_plant = col_double(),
  ..   Self_fertilized_plant = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
anyNA(zea_mays) # Returns FALSE if there are no missing values
[1] TRUE
  • There is NA data being identified
zea_mays |>
  filter( is.na(Self_fertilized_plant) | is.na(Cross_fertilized_plant))
# A tibble: 1 × 4
  Pots   Pair Cross_fertilized_plant Self_fertilized_plant
  <chr> <dbl>                  <dbl>                 <dbl>
1 Pot6     NA                     NA                    NA

Data wrangling plan

After inspecting the data here are some data wrangling (manipulation) works need to be done

  1. Handling missing data
  2. Make a long table so that each observation (Zea mays) must have its own row.
    1. Currently there are two observations in a row

And you will see the reasons later on

Removing Missing Data

Let’s remove the broken pairs before calculating our differences.

Good practice of coding

# 1. Drop any incomplete pairs
zea_mays <- zea_mays |> 
  drop_na()
# 2. Verify the clean-up was successful
anyNA(zea_mays)
[1] FALSE

Not recommended way

#| echo: true
#| eval: false
# Not recommended approach
zea_mays <- zea_mays[-1,]
# 2. Verify the clean-up was successful
anyNA(zea_mays)
  • Try to avoid hard code on specific columns, rows by indexing

Reasons

  • Not robust to the change of data
  • Hard to understand the intention of the code

Data Preparation for Analysis

For our visual and statistical analysis, we need to calculate the true biological signal we care about: the height difference within each pair. We will also pivot the data for plotting.

# Calculate the paired difference
zea_mays <- zea_mays |>
  mutate(diff.height = Cross_fertilized_plant - Self_fertilized_plant)

Analysis and Visualization

Quick review

A hypothesis test asks: if there were truly no difference (or no effect), how surprising is the pattern we see in the data?

  • We state a skeptical null (\(H_0\)) and an alternative (\(H_a\)).
  • A small \(p\)-value means “this extreme a result would be unusual under \(H_0\)” — not “the hypothesis is proven.”

Why pairing matters here

  • Each cross- and self-pollinated plant shared a pot (same soil, water, light).
  • Comparing within pairs removes that shared environment; treating all plants as independent ignores it.

Keep that story in mind for the next two slides (wrong vs right test).

Scatter Plot with Pair Information

This part we would like create a plot like this. Before that there are some issues in the our data.

  • Every row must represent a single observation

Implication: Not compatibile to ggplot

In data wrangling we define

  • Wide table: a table with each observation at each column
  • Long table: a table with each observation at each row

Reshape the zea_mays

We use pivot_longer() to reshape the zea_mays

# Reshape for visualization
zea_mays_long <- zea_mays |> 
  select(Pots,Pair,Self_fertilized_plant,Cross_fertilized_plant) |>
  pivot_longer(cols = c(Cross_fertilized_plant, Self_fertilized_plant), 
               names_to = "Type", values_to = "Height")

head(zea_mays_long,3)
# A tibble: 3 × 4
  Pots   Pair Type                   Height
  <chr> <dbl> <chr>                   <dbl>
1 Pot1      1 Cross_fertilized_plant   23.5
2 Pot1      1 Self_fertilized_plant    17.4
3 Pot1      2 Cross_fertilized_plant   12
  • cols: You tell R which columns are observationss
  • names_to = "Type": The old column headers (i.e. Self_fertilized_plant and Cross_fertilized_plant) become categorical values in a new column called Type.
  • values_to = "Height": The actual measured numbers inside those cells are stacked into a new column called Height.

Visualization

# Visualization of the data
ggplot(zea_mays_long, 
  aes(x = Type, y = Height, group = Pair)) +
  geom_point(aes(color = Type), size = 3) +
  geom_line() +
  scale_x_discrete(name="Pollination Type")+
  scale_y_continuous(name="Height (inches)")

Wrong Approach: Independent t-test

We have 32 plants in total. Let’s see what happens if we adopt Darwin’s way’ (i.e. wrong statistical test). If we forget that the plants were grown in pairs (sharing the same pot), we might mistakenly assume we have two completely independent groups of 16 plants.

The Wrong Approach: Independent t-test

Hypothesis:

\[ \begin{aligned} H_0: & \mu_{\text{cross}} \leq \mu_{\text{self}} \\ H_a: & \mu_{\text{cross}} > \mu_{\text{self}} \end{aligned} \]

# INCORRECT: Running an unpaired 
# (independent) t-test
t.test(zea_mays$Cross_fertilized_plant, 
       zea_mays$Self_fertilized_plant, 
       alternative = "greater",
       paired = FALSE)

    Welch Two Sample t-test

data:  zea_mays$Cross_fertilized_plant and zea_mays$Self_fertilized_plant
t = 1.406, df = 27.048, p-value = 0.08555
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 -0.5046025        Inf
sample estimates:
mean of x mean of y 
 19.11406  16.72656

Right Approach: Paired t-test

The problem with Darwin’s fallacy is that

  • Zea mays in the same pot shared
    • Same nutrient
    • Same weather condition

Violating the independent sample assumption of t-test. The corrected way is proposed by Ronald Fisher

Hypothesis:

\[ \begin{aligned} H_0: & \mu_{\text{diff}} \leq 0 \\ H_a: & \mu_{\text{diff}} > 0 \end{aligned} \]

# CORRECT: Running a paired t-test
t.test(zea_mays$Cross_fertilized_plant, 
       zea_mays$Self_fertilized_plant, 
       alternative = "greater",
       paired = TRUE)

    Paired t-test

data:  zea_mays$Cross_fertilized_plant and zea_mays$Self_fertilized_plant
t = 2.0549, df = 15, p-value = 0.02887
alternative hypothesis: true mean difference is greater than 0
95 percent confidence interval:
 0.3507378       Inf
sample estimates:
mean difference 
         2.3875

Discussion on the results

  • Using wrong statistical test could draw opposite interpretation
    • Darwin: \(p=0.086\) (Fail to reject \(H_0\))
    • Fisher: \(p=0.029\) (Reject \(H_0\))
  • Cross fertilization could make the plant bigger
  • Limitation of the experiments
    • Control of genotypes
    • Sample size

Final remarks on the mini-project

For analysis in data

  • Always knowing the context of the data
  • Understand the research question
  • Plan the analysis and data wrangling

For coding

  • Coding help you record all of your process
  • Let your calculation transparent and
  • easy to validate
  • Clear documentation of your code

Reading list

Bibliography notes

Griffin, Philippa C, Jyoti Khadake, Kate S LeMay, Suzanna E Lewis, Sandra Orchard, Andrew Pask, Bernard Pope, et al. 2018. “Best Practice Data Life Cycle Approaches for the Life Sciences.” F1000Research 6: 1618.
Homer, Winslow. 1865. The Veteran in a New Field. Oil on canvas. The Metropolitan Museum of Art.
Ross, Sheldon M. 2017. Introductory Statistics. Academic Press.