BI 303 Lab 5 Data Analysis I Guide
written by Max Helmberger
2024-02-11
Welcome to lab!
Your goals for today are to:
- Conduct t-tests and ANOVA to detect significant differences between groups
- Determine effect sizes and understand their importance
- Modify figures to include statistical information
Loading Packages and Readying Data
First, remember to load the necessary packages.
library(ggplot2)
library(lterdatasampler)We’ll be using the and_vertebrates dataset. Take a quick
a look at it! View(), str() and ?
are your friends.
Exercise 1: What does the
sectionvariable refer to and how many levels does it have? Comment your answers.
t-tests
A t-test is one of the simplest statistical tests you can use.
They’re used to determine whether the means of two groups are
significantly different from one another. The code below will run a
t-test comparing the length_1_mm variable between
“sections” of the creek.
t.test(length_1_mm~section,data=and_vertebrates)##
## Welch Two Sample t-test
##
## data: length_1_mm by section
## t = 18.43, df = 32189, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group CC and group OG is not equal to 0
## 95 percent confidence interval:
## 6.098594 7.550174
## sample estimates:
## mean in group CC mean in group OG
## 77.16313 70.33874Exercise 2: What’s wrong with this analysis? Take a look at the
and_vertebratesdataset again if you have to. Comment your answer.
That’s right! The length_1_mm variable contains
measurements for both trout and salamanders, which you don’t want to be
mixing. It doesn’t make sense to compare them both with the same
analysis. Now would be a good time to split up the dataset. You could
achieve this by subsetting every time you reference a length variable,
but it’d be easier to make two new versions of the
and_vertebrates dataset.
and_trout = and_vertebrates[which(and_vertebrates$species=="Cutthroat trout"),]
and_salamander = and_vertebrates[which(and_vertebrates$species=="Coastal giant salamander"),]Now we can compare trout and salamander lengths between different
sections of the creek separately (we’ll just do trout here so as not to
repeat ourselves too much). It’s usually a good idea to start with a
quick visualization of your data. For this, you can use base R’s
boxplot() function to make a simple boxplot with a little
less typing than it would take with ggplot2.
boxplot(length_1_mm~section,data=and_trout)
Looks like lengths are about the same, maybe a little higher in the
clear-cut forest section. Well, let’s go ahead and run the test. See the
code below. The ~, in R-speak, basically means “as a
function of”, so the code below “translates” to “Using a t-test, test
for a difference in length as a function of section using the and_trout
dataset”!
t.test(length_1_mm~section,data=and_trout)##
## Welch Two Sample t-test
##
## data: length_1_mm by section
## t = 7.8578, df = 20116, p-value = 4.105e-15
## alternative hypothesis: true difference in means between group CC and group OG is not equal to 0
## 95 percent confidence interval:
## 2.914639 4.851980
## sample estimates:
## mean in group CC mean in group OG
## 85.30988 81.42657R is telling you a few different things here, not just the p-value!
- t = 7.8578 - The t-statistic.
- df = 20116 - The degrees of freedom, equal to the total number of values in each group minus 2 (at least that’s the definition for a two-sample t-test).
- 95 percent confidence interval = 2.91 to 4.85 - The test estimates a confidence interval for the difference between the two means (mean length in both forest sections). In other words, there is a 95% chance that the difference in mean trout length between the two forest sections is between 2.91 and 4.85 mm. If this confidence interval included zero (that is, if it went from a negative number to a positive number), then you couldn’t assume there was a difference between those means.
- mean in group CC = 85.31 - The mean snout-fork length for trout in the CC forest section.
- mean in group OG = 81.43 - The mean snout-fork length for trout in the OG forest section.
Here in lab, just know that you should report the t-statistic and degrees of freedom in your writeups along with your p-values.
R displays the p-value here as 4.105e-15.
e-15 is equivalent to “times 10^-15”. The p-value, written
out, would thus be 0.000000000000004105. Highly significant! So
significant, in fact, that when you report such a p-value, you don’t
have to specify the exact number, just that it’s below a certain very
small threshold (usually 0.001 or 0.0001). See below:
Cutthroat trout snout-fork length was significantly higher in the clear-cut section of the creek comapared to the old-growth forest section (two-sample t-test, t = 7.858, df = 20116, p < 0.001).
Despite the difference being statistically significant, it doesn’t look that large, does it? I mean, you made the boxplot and saw for yourself. This is why you also want to calculate effect size.
Cohen’s d is the most common effect size measurement for the difference between two groups. It divides the difference in two groups’ means by the square root of the mean of their variances. In general:
- a Cohen’s d of 0.2 or smaller = small effect size
- a Cohen’s d of around 0.5 = medium effect size
- a Cohen’s d of 0.8 or larger = large effect size
You can calculate Cohen’s d manually in R, but the lsr
package contains a simple function to do it automatically. The
cohensD function calculates the effect size of the
difference between any two variables.
install.packages("lsr")library(lsr)
cohensD(subset(and_trout,section=="CC")$length_1_mm,subset(and_trout,section=="OG")$length_1_mm)## [1] 0.1099336A Cohen’s d value of .11 is very small. The difference is there, but it would be only telling half the story if you only reported the significant difference and didn’t indicate the size of the difference.
A better way to report the results of this t-test would thus be something like this:
Cutthroat trout in the clear-cut forest section had a mean snout-fork length of 85.31 mm, significantly higher than the mean snout-fork length of 81.43 mm for trout in the old growth forest section (two-sample t-test, t = 7.858, df = 20116, p < 0.001), but with a small effect size (Cohen’s d = 0.11).
Exercise 3: Calculate Cohen’s d for the difference in trout weight between pools and isolated pools (hint: these are two levels of the
unittypevariable).
t-test assumptions
For a t-test, the only assumption you really need to worry about testing in R is the assumption of normality. Both variables you’re comparing need to be normally distributed, or else a t-test’s result can’t be trusted. One common way to test this assumption in R is the Shapiro-Wilk test.
shapiro.test(subset(and_trout,year=="1987")$weight_g)##
## Shapiro-Wilk normality test
##
## data: subset(and_trout, year == "1987")$weight_g
## W = 0.8159, p-value < 2.2e-16The Shapiro-Wilk test tests the null hypothesis that the data are normally distributed. If the p-value obtained is less than 0.05, the null hypothesis is rejected, indicating that the data are not normally distributed.
At large sample sizes (as in, above 50), the Shapiro-Wilk test
becomes oversensitive and unreliable. Because of this, the
shapiro.test() function in R doesn’t even work for
variables with more than 5000 values. The Mack Creek dataset we’re using
here is huge, with over 30000 replicates (which is why the line of code
above only used one year’s data), so the quantitative methods aren’t
super helpful. That said, for your experiments here in lab, with much
smaller sample sizes, they would be.
So what do you do if you can’t use a quantitative test? That’s where
visual methods come in, which are somewhat subjective, but work for any
size of dataset. You’ll generate what’s called a
quantile-quantile plot or qqplot for
short. Not to be confused with a ggplot! A qqplot plots a theoretical
distribution (in this case, normal) against your data. The
qqline() function adds a 1-to-1 line to the plot.
qqnorm(subset(and_trout,section=="CC")$length_1_mm)
qqline(subset(and_trout,section=="CC")$length_1_mm)
If your data are normally distributed, then all or most of the points should fall on or very close to the line. This can be subjective! A couple of outliers close to the minimum or maximum of your data are okay, but if the points veer hard off the line or don’t align with it, the data are not normally distributed.
The qqplot above is clear, though; our trout lengths are not normally distributed, at least not within the clear-cut forests.
Exercise 4: Make a qqplot for the trout lengths from the old-growth forest section. Are those data normally distributed? Comment your answer and explain why.
Another way to visually assess normality is with a density plot. The density plot of normally-distributed data should resemble a simple bell curve.
ggplot(and_trout,aes(x=length_1_mm,fill=section)) +
geom_density(alpha=0.4)## Warning: Removed 5 rows containing non-finite values (`stat_density()`).
In this case, it’s very clear that the data aren’t normally distributed in either section. They’re bimodal.
Were you to test the assumptions of a t-test and find them unmet, you
would want to use the nonparametric equivalent of a t-test instead, a
Wilcoxon rank sum test (wilcox.test())! Luckily, the code
is very similar.
wilcox.test(length_1_mm~section,data=and_trout)##
## Wilcoxon rank sum test with continuity correction
##
## data: length_1_mm by section
## W = 55178843, p-value = 7.697e-16
## alternative hypothesis: true location shift is not equal to 0The test returns two numbers, a “W” value, which the Wilcoxon test statistic (see here for a more detailed explanation if you’re interested), and the p-value. You would report both in your Results section.
ANOVA
ANOVA (ANalysis Of VAriance) is another common statistical test used to assess differences between 3 or more groups, either 3 or more levels of a single factor or 3 or more treatments found in a factorial experiment.
Let’s say you wanted to test whether trout weight differs between creek section and between channel units (different microhabitats within the creek). First, let’s visualize.
ggplot(data=and_trout,mapping=aes(x=unittype,y=weight_g)) +
geom_boxplot() +
facet_wrap(~section)## Warning: Removed 7839 rows containing non-finite values (`stat_boxplot()`).
It definitely looks like weights vary between channel units! Though not that not all channel units are represented within each section. That’s okay. Now, let’s test!
First, you have to create an ANOVA model.
troutweight.aov=aov(weight_g~section+unittype+section*unittype,data=and_trout)This line of code uses the aov() function to create an
ANOVA model (called troutweight.aov in R), a model of how
weight_g differs “as a function of” section,
as a function of unittype, and as a function of
section*unittype. This latter term is the
interaction of section and
unitytpe, which allows the test to detect whether the
effect of one variable (say, unittype) on trout weight
depends on the other variable. More on that later.
ANOVA interpretation and multiple comparisons
Once you’ve made your ANOVA model, just use the
summary() function to check its results.
summary(troutweight.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## section 1 2780 2780 31.19 2.39e-08 ***
## unittype 6 60399 10066 112.94 < 2e-16 ***
## section:unittype 3 3706 1235 13.86 5.09e-09 ***
## Residuals 11973 1067204 89
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 8449 observations deleted due to missingnessThere’s a lot to see here! You’ll notice you have rows for each one of your model terms (section, unittype, and their interaction).
- Df: Degrees of freedom. Note that there are separate values for each model term and the model residuals.
- F value: A measure of the variation between group means relative to the variation within each group. Larger is better; it means your groups are more different from each other.
- Pr(>F): The p-value for that model term.
Indicates that particular variable as having a significant effect. The
section:unittypeterm is the interaction (see above).
Were you to report these results, you would write something like this, providing the F value, degrees of freedom, residual degrees of freedom, and p-value for each term like so:
Trout weight significantly differed by creek section (two-way ANOVA, F1,11973 = 31.19, p < 0.0001), channel unit (two-way ANOVA, F6,11973 = 112.94, p < 0.0001), and their interaction (two-way ANOVA, F3,11973 = 13.86, p < 0.0001).
Knowing that unittype has a significant effect on trout
weight (statistically speaking, not necessarily biologically)
doesn’t tell you which channel units significantly differed from each
other. That’s where multiple comparisons come in. The
TukeyHSD() function runs the Tukey’s Honestly Significant
Difference test on your ANOVA model.
TukeyHSD(troutweight.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weight_g ~ section + unittype + section * unittype, data = and_trout)
##
## $section
## diff lwr upr p adj
## OG-CC -0.9656505 -1.304595 -0.6267063 0
##
## $unittype
## diff lwr upr p adj
## I-C 2.2714112 -3.543093147 8.0859156 0.9117485
## IP-C -6.1458797 -9.083048794 -3.2087107 0.0000000
## P-C 3.8614589 3.249057933 4.4738599 0.0000000
## R-C -0.7438777 -2.320967077 0.8332117 0.8069744
## S-C -4.1636182 -14.691401556 6.3641652 0.9068012
## SC-C -3.1490786 -3.929934897 -2.3682223 0.0000000
## IP-I -8.4172910 -14.914678339 -1.9199036 0.0025636
## P-I 1.5900477 -4.237826331 7.4179217 0.9845924
## R-I -3.0152889 -9.021643178 2.9910653 0.7567086
## S-I -6.4350294 -18.452654108 5.5825953 0.6959176
## SC-I -5.4204898 -11.268465287 0.4274857 0.0901996
## P-IP 10.0073387 7.043790709 12.9708866 0.0000000
## R-IP 5.4020021 2.101276012 8.7027281 0.0000290
## S-IP 1.9822616 -8.937528012 12.9020512 0.9983265
## SC-IP 2.9968012 -0.006083795 5.9996861 0.0508658
## R-P -4.6053366 -6.231025775 -2.9796474 0.0000000
## S-P -8.0250771 -18.560250428 2.5100963 0.2708557
## SC-P -7.0105375 -7.885403731 -6.1356713 0.0000000
## S-R -3.4197405 -14.054685401 7.2152044 0.9645850
## SC-R -2.4052009 -4.101540167 -0.7088616 0.0005810
## SC-S 1.0145396 -9.531766838 11.5608460 0.9999578
##
## $`section:unittype`
## diff lwr upr p adj
## OG:C-CC:C -0.8927032 -1.65052896 -0.1348774 0.0060194
## CC:I-CC:C NA NA NA NA
## OG:I-CC:C 1.3383320 -5.28437426 7.9610383 0.9999920
## CC:IP-CC:C NA NA NA NA
## OG:IP-CC:C -7.0789590 -10.43710887 -3.7208090 0.0000000
## CC:P-CC:C 4.9652389 3.97894490 5.9515328 0.0000000
## OG:P-CC:C 1.9423888 0.98542064 2.8993570 0.0000000
## CC:R-CC:C -0.4083959 -3.66346509 2.8466733 1.0000000
## OG:R-CC:C -1.8052640 -3.95836625 0.3478383 0.2181271
## CC:S-CC:C -4.1310469 -16.11128444 7.8491907 0.9970168
## OG:S-CC:C NA NA NA NA
## CC:SC-CC:C -3.6657197 -4.78791181 -2.5435276 0.0000000
## OG:SC-CC:C -3.0835897 -4.52576153 -1.6414179 0.0000000
## CC:I-OG:C NA NA NA NA
## OG:I-OG:C 2.2310352 -4.39630841 8.8583788 0.9976551
## CC:IP-OG:C NA NA NA NA
## OG:IP-OG:C -6.1862558 -9.55354190 -2.8189696 0.0000001
## CC:P-OG:C 5.8579420 4.84097557 6.8749085 0.0000000
## OG:P-OG:C 2.8350920 1.84654098 3.8236431 0.0000000
## CC:R-OG:C 0.4843073 -2.78018662 3.7488012 0.9999998
## OG:R-OG:C -0.9125608 -3.07988505 1.2547635 0.9807114
## CC:S-OG:C -3.2383437 -15.22114540 8.7444581 0.9997720
## OG:S-OG:C NA NA NA NA
## CC:SC-OG:C -2.7730165 -3.92225981 -1.6237732 0.0000000
## OG:SC-OG:C -2.1908865 -3.65420624 -0.7275668 0.0000458
## OG:I-CC:I NA NA NA NA
## CC:IP-CC:I NA NA NA NA
## OG:IP-CC:I NA NA NA NA
## CC:P-CC:I NA NA NA NA
## OG:P-CC:I NA NA NA NA
## CC:R-CC:I NA NA NA NA
## OG:R-CC:I NA NA NA NA
## CC:S-CC:I NA NA NA NA
## OG:S-CC:I NA NA NA NA
## CC:SC-CC:I NA NA NA NA
## OG:SC-CC:I NA NA NA NA
## CC:IP-OG:I NA NA NA NA
## OG:IP-OG:I -8.4172910 -15.80813358 -1.0264484 0.0099778
## CC:P-OG:I 3.6269068 -3.03043193 10.2842456 0.8640004
## OG:P-OG:I 0.6040568 -6.04900050 7.2571142 1.0000000
## CC:R-OG:I -1.7467279 -9.09130817 5.5978523 0.9999461
## OG:R-OG:I -3.1435960 -10.07058905 3.7833971 0.9638768
## CC:S-OG:I -5.4693789 -19.13954625 8.2007885 0.9877243
## OG:S-OG:I NA NA NA NA
## CC:SC-OG:I -5.0040517 -11.68287639 1.6747729 0.3986653
## OG:SC-OG:I -4.4219217 -11.16189527 2.3180518 0.6291371
## OG:IP-CC:IP NA NA NA NA
## CC:P-CC:IP NA NA NA NA
## OG:P-CC:IP NA NA NA NA
## CC:R-CC:IP NA NA NA NA
## OG:R-CC:IP NA NA NA NA
## CC:S-CC:IP NA NA NA NA
## OG:S-CC:IP NA NA NA NA
## CC:SC-CC:IP NA NA NA NA
## OG:SC-CC:IP NA NA NA NA
## CC:P-OG:IP 12.0441978 8.61825388 15.4701417 0.0000000
## OG:P-OG:IP 9.0213478 5.60373106 12.4389645 0.0000000
## CC:R-OG:IP 6.6705630 2.04889761 11.2922285 0.0001140
## OG:R-OG:IP 5.2736950 1.34939053 9.1979995 0.0005600
## CC:S-OG:IP 2.9479121 -9.47345693 15.3692811 0.9999474
## OG:S-OG:IP NA NA NA NA
## CC:SC-OG:IP 3.4132392 -0.05427156 6.8807500 0.0588275
## OG:SC-OG:IP 3.9953692 0.41149206 7.5792464 0.0134703
## OG:P-CC:P -3.0228500 -4.19575624 -1.8499438 0.0000000
## CC:R-CC:P -5.3736348 -8.69860037 -2.0486691 0.0000054
## OG:R-CC:P -6.7705028 -9.02788412 -4.5131215 0.0000000
## CC:S-CC:P -9.0962857 -21.09570292 2.9031315 0.3782829
## OG:S-CC:P NA NA NA NA
## CC:SC-CC:P -8.6309586 -9.94215496 -7.3197622 0.0000000
## OG:SC-CC:P -8.0488286 -9.64249743 -6.4551597 0.0000000
## CC:R-OG:P -2.3507848 -5.66716964 0.9656001 0.4967277
## OG:R-OG:P -3.7476528 -5.99237612 -1.5029295 0.0000020
## CC:S-OG:P -6.0734357 -18.07047808 5.9236067 0.9171122
## OG:S-OG:P NA NA NA NA
## CC:SC-OG:P -5.6081086 -6.89739068 -4.3188265 0.0000000
## OG:SC-OG:P -5.0259786 -6.60166661 -3.4502905 0.0000000
## OG:R-CC:R -1.3968681 -5.23333388 2.4395978 0.9949139
## CC:S-CC:R -3.7226510 -16.11654917 8.6712473 0.9992818
## OG:S-CC:R NA NA NA NA
## CC:SC-CC:R -3.2573238 -6.62510285 0.1104552 0.0699683
## OG:SC-CC:R -2.6751938 -6.16266856 0.8122809 0.3579123
## CC:S-OG:R -2.3257829 -14.47687697 9.8253112 0.9999958
## OG:S-OG:R NA NA NA NA
## CC:SC-OG:R -1.8604558 -4.18043631 0.4595248 0.2846161
## OG:SC-OG:R -1.2783258 -3.76887875 1.2122272 0.9086714
## OG:S-CC:S NA NA NA NA
## CC:SC-CC:S 0.4653271 -11.54602384 12.4766781 1.0000000
## OG:SC-CC:S 1.0474571 -10.99800246 13.0929167 1.0000000
## CC:SC-OG:S NA NA NA NA
## OG:SC-OG:S NA NA NA NA
## OG:SC-CC:SC 0.5821300 -1.09903619 2.2632962 0.9968902It gives you confidence intervals and p-values for every single pairwise difference in your dataset (e.g. riffle vs. cascade, isolated pool vs. rapid, etc.). You don’t need to report all the pairwise p-values, but knowing which ones are different from which other ones is important for when you make figures or discuss results in more detail.
ANOVA effect size
The effect size metric we’ll use for our ANOVAs is called
Cohen’s f. The function to calculate it is in the
package sjstats, so go ahead and install that.
install.packages("sjstats")library(sjstats)
effectsize::cohens_f(troutweight.aov)## # Effect Size for ANOVA (Type I)
##
## Parameter | Cohen's f (partial) | 95% CI
## ----------------------------------------------------
## section | 0.05 | [0.04, Inf]
## unittype | 0.24 | [0.22, Inf]
## section:unittype | 0.06 | [0.04, Inf]
##
## - One-sided CIs: upper bound fixed at [Inf].Cohen’s f is a bit different from Cohen’s d. The value you get for each independent variable is an estimate of what percentage of variation in the dependent variable (in this case, trout weight) is explained by variation in that independent variable.
- a Cohen’s f of 0.1 or smaller = small effect size
- a Cohen’s f of around 0.25 = medium effect size
- a Cohen’s f of 0.4 or larger = large effect size
So, looking at the output above, we see that section and
the section:unittype interaction term had small effects on
trout weight, whereas unittype has a medium effect. That
tracks, doesn’t it? There didn’t seem to be a big difference between
sections (even though it was statistically significant), but the
difference between channel units was more noticeable.
Trout weight significantly differed by creek section (two-way ANOVA, F1,11973 = 31.19, p < 0.0001), channel unit (two-way ANOVA, F6,11973 = 112.94, p < 0.0001), and their interaction (two-way ANOVA, F3,11973 = 13.86, p < 0.0001), though channel unit had the largest effect on trout weight (Cohen’s f = 0.24).
Exercise 5: Using ANOVA, test for differences in salamander snout-tail length between different channel units. Report your results.
ANOVA assumptions and nonparametric equivalent
Testing the assumptions of an ANOVA requires a different procedure
than a t-test does. The assumption of normality applies instead to the
residuals of the ANOVA model, the differences between
each value in the dataset and the mean of that value’s group. This means
you have to fit the model in order to test whether it fulfills the
assumptions in the first place. To test whether the residuals are
normally distributed, use shapiro.test() or make a qqplot
(or histogram, use the hist() function.
qqnorm(resid(troutweight.aov))
qqline(resid(troutweight.aov))
Super not normal! ANOVA also assumes equal variance between groups,
which is assessed via the bartlett.test() function.
bartlett.test(weight_g~interaction(section,unittype),data=and_trout)##
## Bartlett test of homogeneity of variances
##
## data: weight_g by interaction(section, unittype)
## Bartlett's K-squared = 717.29, df = 10, p-value < 2.2e-16Note that since we’re conducting a two-way ANOVA here (as we have two
factors, section and unittype), we need to put
the two factors inside the interaction() function. If we
were just doing a one-way ANOVA, we wouldn’t need the
interaction().
Like the Shapiro-Wilk test, a p-value below 0.05 means the assumption (of equal variance) is violated. That’s the case here.
Nonparametric equivalent - Kruskal-Wallis test
The nonparametric equivalent of an ANOVA is known as a Kruskal-Wallis
test. It’s simple to run, though you have to use the
interaction() function if you want to make it mimic a
two-way ANOVA.
kruskal.test(weight_g~section,data=and_trout)##
## Kruskal-Wallis rank sum test
##
## data: weight_g by section
## Kruskal-Wallis chi-squared = 36.527, df = 1, p-value = 1.506e-09kruskal.test(weight_g~unittype,data=and_trout)##
## Kruskal-Wallis rank sum test
##
## data: weight_g by unittype
## Kruskal-Wallis chi-squared = 849.15, df = 6, p-value < 2.2e-16kruskal.test(weight_g~interaction(section,unittype),data=and_trout)##
## Kruskal-Wallis rank sum test
##
## data: weight_g by interaction(section, unittype)
## Kruskal-Wallis chi-squared = 917.42, df = 10, p-value < 2.2e-16You’d want to report the chi-squared value, df, and p-value.
Just like the Tukey’s HSD test lets you determine which ANOVA groups are significantly different from each other, the Kruskal post-hoc Nemenyi test lets you do so after conducted a Kruskal-Wallis test.
install.packages("PMCMRplus")library(PMCMRplus)kwAllPairsNemenyiTest(weight_g~section,data=and_trout)
kwAllPairsNemenyiTest(weight_g~unittype,data=and_trout)
kwAllPairsNemenyiTest(weight_g~interaction(section,unittype),data=and_trout)Uh oh! You’ll notice the first two lines of code didn’t work. This is
because the kwAllPairsNemenyiTest() function specifically
requires factors, which if you look at the section and
unittype variables using class(), you’ll see
for yourself. (The third line worked because of the
interaction() function.)
You can make a new variable in the and_trout dataset
that’s a factor version of one of your variables of interest, or you can
put as.factor() directly into the test function - this
avoids having to alter the original dataset. See examples of both
approaches below.
and_trout$section.f=as.factor(and_trout$section)
kwAllPairsNemenyiTest(weight_g~section.f,data=and_trout)## Warning in kwAllPairsNemenyiTest.default(c(1.75, 1.95, 5.6, 2.15, 6.9, 5.9, :
## Ties are present, p-values are not corrected.##
## Pairwise comparisons using Tukey-Kramer-Nemenyi all-pairs test with Tukey-Dist approximation## data: weight_g by section.f## CC
## OG 1.5e-09##
## P value adjustment method: single-step## alternative hypothesis: two.sidedkwAllPairsNemenyiTest(weight_g~as.factor(unittype),data=and_trout)## Warning in kwAllPairsNemenyiTest.default(c(1.75, 1.95, 5.6, 2.15, 6.9, 5.9, :
## Ties are present, p-values are not corrected.##
## Pairwise comparisons using Tukey-Kramer-Nemenyi all-pairs test with Tukey-Dist approximation## data: weight_g by as.factor(unittype)## C I IP P R S
## I 0.5241 - - - - -
## IP 9.3e-15 3.2e-10 - - - -
## P < 2e-16 1.0000 < 2e-16 - - -
## R 1.0000 0.5134 6.2e-14 9.1e-11 - -
## S 0.7742 0.3185 0.6572 0.1681 0.8055 -
## SC < 2e-16 0.0021 3.3e-11 < 2e-16 2.7e-10 0.9999##
## P value adjustment method: single-step## alternative hypothesis: two.sidedFrom all these yests, you’ll have gotten warnings saying that “ties are present” - this is okay to ignore here. The way these nonparametric tests work is by assigning “ranks” to all the values in your dataset, and so if any values are the same, you have “ties” and the test can’t calculate the exact p-value, but will still calculate it with enough accuracy to trust ones that are appreciably below 0.05.
Representing Statistical Results in Figures
ggplot(data=and_salamander,mapping=aes(x=unittype, y=length_1_mm,fill=unittype)) +
geom_boxplot() +
theme_bw() +
theme(legend.position = "none") +
labs(x="Channel Unit Type",y="Salamander Snout-Vent Length (mm)")## Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Consider this plot. There are a few different ways the results of
statistical tests can be represented on figures. We’ll focus on the most
common, the use of “Tukey letters” to denote significantly different
groups on a boxplot or barplot. First, you need to install the
ggtukey package.
install.packages("ggpubr")
install.packages("devtools")
devtools::install_github("https://github.com/ethanbass/ggtukey/")library(ggtukey)
ggplot(data=and_salamander,mapping=aes(x=unittype, y=length_1_mm,fill=unittype)) +
geom_boxplot() +
theme_bw() +
theme(legend.position = "none") +
labs(x="Channel Unit Type",y="Salamander Snout-Vent Length (mm)") +
geom_tukey()## Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
A few things to note here.
If two plots share a letter, it means they don’t significantly differ from each other.
In case you’re wondering, the reason that S (steps, small waterfalls) is not significantly different from any other group is because its sample size is tiny (just 2 salamanders found in steps).
Typically, you want the letters placed fully above each plot. This can be done with the code below. A less elegant but much more efficient way to make this type of plot is to simply save the R plot as a PDF or image and add the letters manually in a PDF editor or graphics program.
Regardless, in many cases, adding letters requires one small tweak to the code, adjusting the scale of the y-axis to ensure that there’s enough room. This can be done like so, using the
ylim()function:
ggplot(data=and_salamander,mapping=aes(x=unittype, y=length_1_mm,fill=unittype)) +
geom_boxplot() +
theme_bw() +
theme(legend.position = "none") +
labs(x="Channel Unit Type",y="Salamander Snout-Vent Length (mm)") +
geom_tukey(where="whisker",vjust=-14,hjust=0.5) +
ylim(0,210)## Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Note that the geom_tukey() now contains some arguments
to specify the placement of the labels. Remember that you can always use
? to view the help file for a function and see what the
available arguments are.
?geom_tukeyExercise 6: Make a boxplot comparing trout snout-fork length between channel units. Give it proper axis labels and add Tukey letters in a sensible location.