Chapter 5 - Introducing Statistics in R
Omar Masoud
Chapter 5: Introducing Statistics in R
5.1 Introduction
In this chapter, we will practice four statistical analyses:
- Chi-square contingency table analysis
- Two-sample t-test
- Simple linear regression
- One-way ANOVA
Setup
5.2 Chi-square Contingency Table Analysis
5.2.1 Importing and Exploring the Data
# Load required packages
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(ggfortify)
# Import Chapter 5 datasets
lady <- read.csv("ladybirds_morph_colour.csv")
ozone <- read.csv("ozone.csv")
plant_gr <- read.csv("plant.growth.rate.csv")
daphnia <- read.csv("Daphniagrowth.csv")
# Check datasets
str(lady)## 'data.frame': 20 obs. of 4 variables:
## $ Habitat : chr "Rural" "Rural" "Rural" "Rural" ...
## $ Site : chr "R1" "R2" "R3" "R4" ...
## $ morph_colour: chr "black" "black" "black" "black" ...
## $ number : int 10 3 4 7 6 15 18 9 12 16 ...str(ozone)## 'data.frame': 20 obs. of 3 variables:
## $ Ozone : num 61.7 64 72.4 56.8 52.4 44.8 70.4 67.6 68.8 53.7 ...
## $ Garden.location: chr "West" "West" "West" "West" ...
## $ Garden.ID : chr "G1" "G2" "G3" "G4" ...str(plant_gr)## 'data.frame': 50 obs. of 2 variables:
## $ soil.moisture.content: num 0.47 0.541 1.698 0.826 0.857 ...
## $ plant.growth.rate : num 21.3 27 39 30.2 37.1 ...str(daphnia)## 'data.frame': 40 obs. of 3 variables:
## $ parasite : chr "control" "control" "control" "control" ...
## $ rep : int 1 2 3 4 5 6 7 8 9 10 ...
## $ growth.rate: num 1.07 1.27 1.32 1.08 1.2 ...library(dplyr)
totals <- lady %>%
group_by(Habitat, morph_colour) %>%
summarise(total.number = sum(number), .groups = "drop")
totals## # A tibble: 4 × 3
## Habitat morph_colour total.number
## <chr> <chr> <int>
## 1 Industrial black 115
## 2 Industrial red 85
## 3 Rural black 30
## 4 Rural red 705.2.3 Visualizing the Data
library(ggplot2)
ggplot(totals, aes(x = Habitat,
y = total.number,
fill = morph_colour)) +
geom_bar(stat = "identity", position = "dodge")
ggplot(totals, aes(x = Habitat,
y = total.number,
fill = morph_colour)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c(black = "black",
red = "red"))
5.2.5 Interpreting the Graph
Null Hypothesis (H0): There is no association between ladybird colour and habitat.
Alternative Hypothesis (H1): There is an association between ladybird colour and habitat.
Interpretation:
The graph shows that black ladybirds are relatively more common in industrial habitats, while red ladybirds are relatively more common in rural habitats.
This suggests a possible association between ladybird colour and habitat.
We expect to reject the null hypothesis, but a Chi-square test is needed to evaluate the statistical evidence. ### 5.2.6 Chi-square Test
lady.mat <- xtabs(number ~ Habitat + morph_colour,
data = lady)
lady.mat## morph_colour
## Habitat black red
## Industrial 115 85
## Rural 30 70lady.chi <- chisq.test(lady.mat)
lady.chi##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: lady.mat
## X-squared = 19.103, df = 1, p-value = 1.239e-05Explore the Chi-square test results
names(lady.chi)## [1] "statistic" "parameter" "p.value" "method" "data.name" "observed"
## [7] "expected" "residuals" "stdres"Observed frequencies
lady.chi$observed## morph_colour
## Habitat black red
## Industrial 115 85
## Rural 30 70Expected frequencies
lady.chi$expected## morph_colour
## Habitat black red
## Industrial 96.66667 103.33333
## Rural 48.33333 51.666675.2.7 Interpretation and Conclusion
Statistical Results:
A Chi-square contingency table analysis was performed to examine the association between ladybird colour morph and habitat type.
The results showed a statistically significant association between ladybird colour and habitat (χ² = 19.103, df = 1, p < 0.001).
Black ladybirds were relatively more common in industrial habitats (57.5%), whereas red ladybirds were relatively more common in rural habitats (70%).
Conclusion:
We reject the null hypothesis (H0) and conclude that ladybird colour morph is significantly associated with habitat type. ## 5.3 Two-sample t-test
5.3.1 Visualizing Ozone Data
library(ggplot2)
ggplot(ozone, aes(x = Ozone)) +
geom_histogram(binwidth = 10) +
facet_wrap(~ Garden.location, ncol = 1) +
theme_bw()
### 5.3.3 Performing the Two-sample t-test
t.test(Ozone ~ Garden.location, data = ozone)##
## Welch Two Sample t-test
##
## data: Ozone by Garden.location
## t = 4.2363, df = 17.656, p-value = 0.0005159
## alternative hypothesis: true difference in means between group East and group West is not equal to 0
## 95 percent confidence interval:
## 8.094171 24.065829
## sample estimates:
## mean in group East mean in group West
## 77.34 61.265.5 Simple Linear Regression
5.5.1 Plotting the Data
library(ggplot2)
ggplot(plant_gr,
aes(x = soil.moisture.content,
y = plant.growth.rate)) +
geom_point() +
ylab("Plant Growth Rate (mm/week)") +
theme_bw()
5.5.3 Fitting the Linear Regression Model
5.5.4 Checking Model Assumptions
# Fit the linear regression model
model_pgr <- lm(plant.growth.rate ~ soil.moisture.content,
data = plant_gr)
# Diagnostic plots
library(ggfortify)
autoplot(model_pgr, smooth.colour = NA)## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 50 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 50 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 50 rows containing missing values or values outside the scale range
## (`geom_line()`).
### 5.5.5 Linear Regression Results
# ANOVA table
anova(model_pgr)## Analysis of Variance Table
##
## Response: plant.growth.rate
## Df Sum Sq Mean Sq F value Pr(>F)
## soil.moisture.content 1 2521.15 2521.15 156.08 < 2.2e-16 ***
## Residuals 48 775.35 16.15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Regression model summary
summary(model_pgr)##
## Call:
## lm(formula = plant.growth.rate ~ soil.moisture.content, data = plant_gr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.9089 -3.0747 0.2261 2.6567 8.9406
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.348 1.283 15.08 <2e-16 ***
## soil.moisture.content 12.750 1.021 12.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.019 on 48 degrees of freedom
## Multiple R-squared: 0.7648, Adjusted R-squared: 0.7599
## F-statistic: 156.1 on 1 and 48 DF, p-value: < 2.2e-165.5.6 Plotting the Regression Line
ggplot(plant_gr, aes(x = soil.moisture.content,
y = plant.growth.rate)) +
geom_point() +
geom_smooth(method = "lm") +
ylab("Plant Growth Rate (mm/week)") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
### 5.5.7 Linear Regression Summary
A simple linear regression was performed to examine the relationship between soil moisture content and plant growth rate.
Results:
- Intercept = 19.348
- Slope = 12.750
- F(1, 48) = 156.08
- p < 0.001
- R² = 0.7648
- Adjusted R² = 0.7599
5.6 One-way ANOVA
5.6.1 Importing and Plotting Daphnia Data
library(ggplot2)
ggplot(daphnia, aes(x = parasite, y = growth.rate)) +
geom_boxplot() +
theme_bw() +
coord_flip()
### 5.6.2 Constructing the ANOVA Model
5.6.3 Checking Model Assumptions
# Fit the one-way ANOVA model
model_grow <- lm(growth.rate ~ parasite, data = daphnia)
# Check model assumptions
library(ggfortify)
autoplot(model_grow, smooth.colour = NA)## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).
### 5.6.4 One-way ANOVA Results
# ANOVA table
anova(model_grow)## Analysis of Variance Table
##
## Response: growth.rate
## Df Sum Sq Mean Sq F value Pr(>F)
## parasite 3 3.1379 1.04597 32.325 2.571e-10 ***
## Residuals 36 1.1649 0.03236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model summary and treatment contrasts
summary(model_grow)##
## Call:
## lm(formula = growth.rate ~ parasite, data = daphnia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41930 -0.09696 0.01408 0.12267 0.31790
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.21391 0.05688 21.340 < 2e-16 ***
## parasiteMetschnikowia bicuspidata -0.41275 0.08045 -5.131 1.01e-05 ***
## parasitePansporella perplexa -0.13755 0.08045 -1.710 0.0959 .
## parasitePasteuria ramosa -0.73171 0.08045 -9.096 7.34e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1799 on 36 degrees of freedom
## Multiple R-squared: 0.7293, Adjusted R-squared: 0.7067
## F-statistic: 32.33 on 3 and 36 DF, p-value: 2.571e-105.6.5 Treatment Contrasts and Group Means
library(dplyr)
# Calculate mean growth rate for each treatment
sumDat <- daphnia %>%
group_by(parasite) %>%
summarise(meanGR = mean(growth.rate))
# Display group means
sumDat## # A tibble: 4 × 2
## parasite meanGR
## <chr> <dbl>
## 1 Metschnikowia bicuspidata 0.801
## 2 Pansporella perplexa 1.08
## 3 Pasteuria ramosa 0.482
## 4 control 1.215.6.6 Final ANOVA Plot
ggplot(daphnia, aes(x = parasite, y = growth.rate)) +
geom_point(position = position_jitter(width = 0.1),
alpha = 0.6) +
geom_point(data = sumDat,
aes(x = parasite, y = meanGR),
inherit.aes = FALSE,
shape = 23, size = 4, fill = "red") +
coord_flip() +
theme_bw() +
labs(x = "Parasite",
y = "Growth Rate")
## 5.7 Chapter Summary
Four statistical methods were applied using R:
- Chi-square contingency table analysis.
- Welch two-sample t-test.
- Simple linear regression.
- One-way ANOVA.
Data visualization, model fitting, diagnostic plots, and statistical summaries were performed using R, ggplot2, dplyr, and ggfortify.