.
├── Descriptive Statistics
| └── Describes an individual dataset by e.g.
| ├── Central Tendency (Mean, Median, Mode)
| ├── Measures of Variability (Range, Variance, Standard Deviation)
| └── Frequency distribution
|
|
└── Inferential Statistics
Session 4 :: Inferential Statistics
.
├── Descriptive Statistics
| └── Describes an individual dataset by e.g.
| ├── Central Tendency (Mean, Median, Mode)
| ├── Measures of Variability (Range, Variance, Standard Deviation)
| └── Frequency distribution
|
|
└── Inferential Statistics
└── Uses a small sample to make predictions about an entire population
├── Hypothesis testing ( --> p-value)
└── Confidence intervals for inferred values of numerical parameters
Categorical Tests
The Chi-Square[d] Test – Idea
Given
a set of individuals (anything: mice, days of the year, soil samples …)
one or more categorical trait(s), measured on each individual
(symptom severity, ice day or not, whether species X was found in a soil sample …)
a pre-specified hypothesis about the frequencies of the categories
The test rejects our hypothesis if this chi-square distributed quantity it “too large”:
Chi-Square Test of Goodness of Fit
Chi-Square Test of Goodness of Fit
Chi-Square Test of Goodness of Fit
Chi-Square Test of Goodness of Fit
Chi-Square Test of Independence [or: of Association]
Chi-Square Test of Independence – Expected Counts
Chi-Square Test of Independence – Expected Counts
Chi-Square Test of Independence – Expected Counts
Chi-Square Test of Independence – Expected Counts
Chi-Square Test of Independence – Expected Counts
Important
For a chi-square test to be valid, ideally all expected counts, but minimally 80% of the expected counts, should be >= 5 !
Example Dataset
The data matrix used in the example is in the repository in the data/ folder. Copy it to your project directory and read it into R like so:
mouse_data <-readRDS("mouse_data.rds")# Have look:mouse_data
genotype severity
1 white none
2 white strong
3 white strong
4 white strong
5 white medium
6 white none
7 white strong
8 white none
9 white medium
10 white strong
11 white medium
12 white weak
13 white none
14 white medium
15 white none
16 white weak
17 white weak
18 white none
19 white weak
20 white medium
21 white none
22 white medium
23 white none
24 white none
25 white none
26 white none
27 white medium
28 white weak
29 white none
30 white none
31 white medium
32 white medium
33 white none
34 white weak
35 white medium
36 white medium
37 white weak
38 white weak
39 white strong
40 white strong
41 white strong
42 white medium
43 white weak
44 white medium
45 white weak
46 white strong
47 white medium
48 white weak
49 white strong
50 white medium
51 black strong
52 black medium
53 black medium
54 black strong
55 black strong
56 black medium
57 black medium
58 black strong
59 black medium
60 black strong
61 black medium
62 black medium
63 black medium
64 black weak
65 black strong
66 black none
67 black strong
68 black strong
69 black medium
70 black weak
71 black weak
72 black weak
73 black none
74 black weak
75 black strong
76 black strong
77 black strong
78 black medium
79 black strong
80 black strong
81 black strong
82 black weak
83 black strong
84 black strong
85 black strong
86 black strong
87 black medium
88 black weak
89 black weak
90 black weak
91 black strong
92 black strong
93 black weak
94 black medium
95 black strong
96 black medium
97 black medium
98 black strong
99 black medium
100 black strong
# test the white genotype's frequencies # against the uniform distributionchisq.test(x =table(mouse_data)[1,])
Chi-squared test for given probabilities
data: table(mouse_data)[1, ]
X-squared = 1.36, df = 3, p-value = 0.7149
Chi-square Tests in R – the chisq.test Function
# test the black genotype's frequencies # against the uniform distributionchisq.test(x =table(mouse_data)[2,])
Chi-squared test for given probabilities
data: table(mouse_data)[2, ]
X-squared = 18.64, df = 3, p-value = 0.0003245
Chi-square Tests in R – the chisq.test Function
# Are the black genotype's frequencies # compatible with the white's probabilities?tbl <-table(mouse_data)chisq.test(x = tbl[2,], p = tbl[1,]/sum(tbl[1,]))
Chi-squared test for given probabilities
data: tbl[2, ]
X-squared = 27.277, df = 3, p-value = 5.151e-06
An alternative For Sparse Tables – Fisher’s Exact Test
For 2x2 contingency tables, it is possible to exactly calculate the total probability P of all tables with
the same row and column sums as the observed table
the same or a stronger non-uniform count distribution as the observed table
This p-value calculation is not affected by small cell counts}
For tables larger than 2x2, an extension exists which can still compute an exact p-value
Function fisher.test implements both the 2x2 and the extended case
?fisher.test
:::
An alternative For Sparse Tables – Fisher’s Exact Test
fisher.test(x=mouse_data[,1], y=mouse_data[,2])
Fisher's Exact Test for Count Data
data: mouse_data[, 1] and mouse_data[, 2]
p-value = 0.002082
alternative hypothesis: two.sided
# As with chisq.test, one can get the same result# by specifying the pre-computed contingency table:fisher.test(x=table(mouse_data))
Fisher's Exact Test for Count Data
data: table(mouse_data)
p-value = 0.002082
alternative hypothesis: two.sided
Quantitative Tests (Numerical Data)
Entry point: Our mouse categorical data and the Normal Distribution
Quantitative Tests (Numerical Data)
A similar principle applies when numerical data are compared to a reference distribution:
Quantitative Tests (Numerical Data)
From here on we will use Fisher's Iris dataset as a (numerical) example.
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50
Median :5.800 Median :3.000 Median :4.350 Median :1.300 virginica :50
Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
Quantitative Tests (Numerical Data)
We will restrict ourselves here to the classical tests:
t-test
Anova
linear model
These tests are conceptually closely related andbuild upon each other hierarchically.
Test Functions and Their Assumptions
Do the Assumptions Hold for the Iris Flower Traits?
Test Normality of Sepal.Length using shapiro.test():
A Quick Tour of Concepts: Independent Samples T-Test
It tests whether two population means are identical.
A Quick Tour of Concepts: Paired Samples T-Test
It tests whether the mean difference between pairs of measurements is zero or not.
Equivalent to an unpaired t-test on the differences between the conditions
Note that the two conditions can be anything, not necessarily pre/post!
Running T-Tests in R – the t.test() function
The t.test() function is defined in the stats package like chisq.test(). It, too, is multi-functional, in the sense that it can run all variants of t-tests described above.
See its help page:
?t.test
Running T-Tests in R – the t.test() function
Running a One-Sample T-Test in R
# Is Sepal.Width > 2 in all species jointly?t.test(x = iris$Sepal.Width,mu =2)
One Sample t-test
data: iris$Sepal.Width
t = 29.71, df = 149, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 2
95 percent confidence interval:
2.987010 3.127656
sample estimates:
mean of x
3.057333
Running a One-Sample T-Test in R
s <-"virginica"t.test(x = iris$Sepal.Width[ iris$Species == s ],mu =2) # Can you express it in dplyr style?
One Sample t-test
data: iris$Sepal.Width[iris$Species == s]
t = 21.356, df = 49, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 2
95 percent confidence interval:
2.882347 3.065653
sample estimates:
mean of x
2.974
Two Sample t-test
data: iris$Sepal.Width[iris$Species == s1] and iris$Sepal.Width[iris$Species == s2]
t = 9.455, df = 98, p-value = 1.845e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.519895 0.796105
sample estimates:
mean of x mean of y
3.428 2.770
Paired t-test
data: iris$Sepal.Width[iris$Species == s] and iris$Sepal.Length[iris$Species == s]
t = -42.323, df = 49, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-1.652926 -1.503074
sample estimates:
mean difference
-1.578
Combine Testing and Pretty Printing
## Define a wrapper function around ggboxplotplot ## (slightly simpler: It computes all possible pairs of group_var's levels and passes them ## (as the list of comparisons to draw) on to stat_compare_means()require(dplyr)require(ggplot2)require(ggpubr)ggbp <-function(data, group_var, test_var, color, palette,var.equal=FALSE,paired =FALSE) {# Construct a list, containing all# possible pairs of "group_var"'s levels: cmb <-utils::combn(data |> dplyr::pull(group_var) |>levels(),2,# insist to return a *list*!simplify=FALSE ) ggboxplot(data=data, x = group_var, y = test_var,color = group_var, palette = palette,order =levels(data |>pull(group_var)),ylab = test_var, xlab = group_var) +# pass this list to stat_compare: stat_compare_means(method ="t.test",paired = paired,method.args =list(var.equal = var.equal),comparisons = cmb ) }
Anova tests the hypothesis: All of >= 2 group means are identical.
For two groups, this is the same as a standard 2-sample t-test.
A Quick Tour of Concepts: ANOVA
A Quick Tour of Concepts: ANOVA
A Quick Tour of Concepts: ANOVA
Important
A significant F value rejects the Null "hypothesis that all means are equal"
One needs to do individual pairwise post-hoc t-tests to find the actual differences!
Multiple post-hoc tests need multiple testing correction
ANOVA – A Beast in Two Deguises
Anova as described above answers questions of the kind: “I have data from N sources: Is there any evidence that not all are actually
samples from the same underlying population?“
This type of anova test is implemented by the R functions stats::aov and stats::oneway.test.
The second kind is implemented by functions stats::anova and car::Anova.
It asks: “Given a linear model like Sepal.Width ~ Sepal.Length + Species,
do I actually need both Sepal.Length and Species to describe Sepal.Width,
or can I drop a predictor?”
Note that in both cases the goal is a parsimonious description of the data, which nevertheless identifies and keeps true, data-inherent structure.
?aov
oneway.test(): An ANOVA Test Which Can Deal With Unequal Variances
# run_formula_test is defined in R/functions.Rlapply(colnames(iris)[1:4], ## run oneway.test() for each trait \(x)run_formula_test(iris, "Species",x,test_f=oneway.test) )|>bind_rows() |>## concatenate the results by rowselect(statistic,p.value,group_var,test_var) ## filter outputs to keep
# A tibble: 4 × 4
statistic p.value group_var test_var
<dbl> <chr> <chr> <chr>
1 139. 1.5e-28 Species Sepal.Length
2 45.0 1.4e-14 Species Sepal.Width
3 1828. 2.7e-66 Species Petal.Length
4 1277. 4.1e-64 Species Petal.Width
The Null Hypothesis of Equality of all Species Means is strongly rejected for all 4 traits. However the Petal traits are clearly more species specific than the Sepal trait.
oneway.test(): An ANOVA Test Which Can Deal With Unequal Variances
Important
Note: The oneway.test() function expects a formula lhs ~ rhs, where lhs is a column of normally distributed values, and rhs is a column of categorical group labels. The groups so specified need not have equal variance.
oneway.test(Petal.Width ~ Species, data=iris)
One-way analysis of means (not assuming equal variances)
data: Petal.Width and Species
F = 1276.9, num df = 2.000, denom df = 84.951, p-value < 2.2e-16
?oneway.test
A Quick Tour of Concepts: Linear Regression
Multiple Linear Regression
?lm
T-test, Anova, and Linear Model
Comparing Groups: Sepal.Width ~ Species
Sepal.Width meets all assumptions, so we can run a t.test() with var.equal=TRUE :
Important
At least in stats::t.test v. 4.4.3, var.equal=FALSE by default!
If you insist to assume equal variance, you must set it to TRUE. Better check ?t.test!
Species estimate statistic p.value parameter conf.low conf.high method alternative
1 setosa 3.428 63.94592 7.167044e-49 49 3.320271 3.535729 One Sample t-test two.sided
2 versicolor 2.770 62.41862 2.310251e-48 49 2.680820 2.859180 One Sample t-test two.sided
3 virginica 2.974 65.20798 2.780789e-49 49 2.882347 3.065653 One Sample t-test two.sided
Sepal.Width is different from zero with very high probability in all species, but the exact value of the means seem to differ: at least setosa seems to have a clearly higher mean.
T-test, Anova, and Linear Model
Comparing Groups: Sepal.Width ~ Species
Does Anova confirm that the means are not equal?
We can use the standard function aov, because normality and equal variance hold for this trait.
aov_res <-aov(Sepal.Width ~ Species, data=iris) # keep it for later use aov_res |>summary() # output the summary only
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 11.35 5.672 49.16 <2e-16 ***
Residuals 147 16.96 0.115
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
T-test, Anova, and Linear Model
Comparing Groups: Sepal.Width ~ Species
The standard function for post-hoc tests is stats::TukeyHSD:
it is tailored to avoid inflation of significance due to multiple testing.
It assumes normality and equal variance (holds for Sepal.Width). If these assumptions are violated, rstatix::games_howell_test can be used as an alternative.
T-test, Anova, and Linear Model
Comparing Groups: Sepal.Width ~ Species
You feed TukeyHSD() with the Anova output object:
TukeyHSD(aov_res)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sepal.Width ~ Species, data = iris)
$Species
diff lwr upr p adj
versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802
We see that that all pairs of means are significantly different. Notably, setosa has a clearly larger Sepal.Width than both virginica and versicolor.
T-test, Anova, and Linear Models
Beyond Bulk Comparison: Model the Data Points Individually
(See also this and that resource on the relation of the t-test, Anova, and the linear model)
A More Complex Model – Include Sepal.Length (Numeric!)
mod <-lm(Sepal.Width ~ Sepal.Length + Species, data=iris)
mod_sum <-summary(mod) # This produces an output in the style of the previous slide
mod_sum
Call:lm(formula = Sepal.Width ~ Sepal.Length + Species, data = iris)Residuals: Min 1Q Median 3Q Max -0.95096 -0.16522 0.00171 0.18416 0.72918 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.67650 0.23536 7.123 4.46e-11 ***Sepal.Length 0.34988 0.04630 7.557 4.19e-12 ***Speciesversicolor -0.98339 0.07207 -13.644 < 2e-16 ***Speciesvirginica -1.00751 0.09331 -10.798 < 2e-16 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.289 on 146 degrees of freedomMultiple R-squared: 0.5693, Adjusted R-squared: 0.5604 F-statistic: 64.32 on 3 and 146 DF, p-value: < 2.2e-16
T-test, Anova, and Linear Models
Nearly everything which appears in an lm() output can be extracted:
Get the Table of Model Coefficient Estimates and T-Test Results:
# The t-test compares the individual estimates # against zeromod_sum$coefficients
Nearly everything which appears in an lm() output can be extracted:
Get the F Statistic (Compares the Entire Model to the Null Model):
mod_sum$fstatistic # lazy style -- just copy it
value numdf dendf
64.31529 3.00000 146.00000
# Or do it yourself: The total significance of the model is the significance gained # when comparing to the NULL model (without any predictors):anova(lm(Sepal.Width ~1, data = iris), mod)
Analysis of Variance Table
Model 1: Sepal.Width ~ 1
Model 2: Sepal.Width ~ Sepal.Length + Species
Res.Df RSS Df Sum of Sq F Pr(>F)
1 149 28.307
2 146 12.193 3 16.114 64.315 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
T-test, Anova, and Linear Models
Nearly everything which appears in an lm() output can be extracted:
lm(Sepal.Width ~ Sepal.Length, data = iris) |>summary()
No significant relation ?
Refining a Model: The Role of Sepal.Length
lm(Sepal.Width ~ Sepal.Length + Species, data = iris) |>summary()
There is a (strong!) pattern, driven by the Species Factor: When taking the Species-dependent slopes into account, Sepal.Width clearly depends on Sepal.Length within each individual Species.
Refining a Model: The Role of Sepal.Length – and Species!
ggplot2::ggplot( iris,aes(y=Sepal.Width,x=Sepal.Length,color = Species ))+geom_point()+geom_smooth(method="lm")
NOTE: Written like this, it seems to fit the line of each species individually!
Refining a Model: The Role of Sepal.Length – and Species!
Last Topic Just Scratched: Anova for Comparing Models
The base::anova() and car::Anova() functions can do more than checking
whether group means are equal:
they can help to streamline an entire model,
by testing the variance explained by the individual predictors of the model (e.g. “Sepal.Length” and “Species” above).
Predictors with a non-significant contribution can be dropped
Last Topic Just Scratched: Anova for Comparing Models
anova() is meant for balanced data (same number of observations for each group).
When passed
a single result object of a modelling function, it will
start with the Null model,
then add the predictors in the order in which they appear in the model
report the significance gained or lost at each step
# "Sum Sq" (Sum of Squares) is the amount of variance# explained by adding the predictor to the model, # given all other predictors already added. anova(lm(Sepal.Width ~ Sepal.Length + Species, data=iris))
Analysis of Variance Table
Response: Sepal.Width
Df Sum Sq Mean Sq F value Pr(>F)
Sepal.Length 1 0.3913 0.3913 4.6851 0.03205 *
Species 2 15.7225 7.8613 94.1304 < 2e-16 ***
Residuals 146 12.1931 0.0835
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Last Topic Just Scratched: Anova for Comparing Models
anova() is meant for balanced data (same number of observations for each group).
When passed
a single result object of a modelling function, it will
start with the Null model,
then add the predictors in the order in which they appear in the model
report the significance gained or lost at each step
Analysis of Variance Table
Response: Sepal.Width
Df Sum Sq Mean Sq F value Pr(>F)
Sepal.Length 1 0.3913 0.3913 4.6851 0.03205 *
Species 2 15.7225 7.8613 94.1304 < 2e-16 ***
Residuals 146 12.1931 0.0835
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Last Topic Just Scratched: Anova for Comparing Models
anova() is meant for balanced data (same number of observations for each group).
When passed
two or more separate objects, it will
compare them sequentially: 2 to 1, 3 to 2, …
This makes most sense if the models are nested: models later in the list contains all predictors of earlier models, plus some more.
# "Sum of Sq" (Sum of Squares) is the amount of variance # explained by adding the predictor to the model, # given all other predictors already added. anova(lm(Sepal.Width ~1, # means: "Intercept only" data=iris), lm(Sepal.Width ~ Species, data=iris),lm(Sepal.Width ~ Species + Sepal.Length , data=iris))
Analysis of Variance Table
Model 1: Sepal.Width ~ 1
Model 2: Sepal.Width ~ Species
Model 3: Sepal.Width ~ Species + Sepal.Length
Res.Df RSS Df Sum of Sq F Pr(>F)
1 149 28.307
2 147 16.962 2 11.3449 67.922 < 2.2e-16 ***
3 146 12.193 1 4.7689 57.102 4.187e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Last Topic Just Scratched: Anova for Comparing Models
anova() is meant for balanced data (same number of observations for each group).
When passed
two or more separate objects, it will
compare them sequentially: 2 to 1, 3 to 2, …
This makes most sense if the models are nested: models later in the list contains all predictors of earlier models, plus some more.
This gives you essentially the same result:
anova(lm(Sepal.Width ~ Species + Sepal.Length , data=iris))
Analysis of Variance Table
Response: Sepal.Width
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 11.3449 5.6725 67.922 < 2.2e-16 ***
Sepal.Length 1 4.7689 4.7689 57.102 4.187e-12 ***
Residuals 146 12.1931 0.0835
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
?anova.lm
Last Topic Just Scratched: Anova for Comparing Models
car::Anova() has less restrictions than stats::anova():
It can deal with unbalanced data
It computes “Type-II” or “Type-III” analysis-of-variance tables
Type-II [default]: Test each term after all others, but “omit a term’s higher-order relatives”
Type-III: Test each term after all others [the help page cautions against this]
These strategies assure that each test asks the question: Will this added term detect additional structure, not already detected by the other terms?
Anova Table (Type II tests)
Response: Sepal.Width
Sum Sq Df F value Pr(>F)
Sepal.Length 4.7689 1 57.102 4.187e-12 ***
Species 15.7225 2 94.130 < 2.2e-16 ***
Residuals 12.1931 146
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
If tested against a model which already contains Species, it becomes evident that Sepal.Length helps “things fall in place”: Its relative contribution is now much higher.
