Earthworms are very beneficial to agricultural soils but are highly detrimental to forest soils. If you find a worm in the dirt in the woods, it’s almost definitely an invasive species. There are several invasive species that can occur in soils, from small red ones to big night crawlers. All of these invasive worm species degrade soils by increasing decomposition rates, eating fungi, and reducing organic matter content. Changes in forest soil ecology can impact plant physiology such as photosynthesis by impacting the ability of plants to access water. To determine the impact of earthworms on forest plants we carried out 3 experimental studies. The 1st two experiments were conducted with potted plants in a greenhouse. The 2nd experiment was also conducted with potted plants, but the plants were placed in the forest so that they experience more natural light conditions.

For this study we grew 40 Trillium in pots in a greenhouse. 20 of the plants were grown without worms, and 20 had 5 earthworms added to their pot. One week after the experiment began, we measured photosynthesis on all of the plants.

For this study we grew 40 Trillium in pots in a greenhouse. We 1st measured photosynthetic rates on all of the plants. We then introduced 5 earthworms to each pot. After 1 week, we re-measured photosynthesis of each plant.

For our 3rd study, we investigated the impact of 2 different species of worms. We used 60 potted Trillium that were placed in the forest. 20 pots were controls and had no worms added. The other 40 plots had one of two different invasive worm species added. 20 pots had “species 1” added, and 20 pots had “species 2” added. After 1 week, photosynthesis was measured on all pots. Preliminary analyses indicate that these data should be log transformed.

**Preliminaries** Clear workspace and load functions to make the data

```
rm(list = ls())
setwd("~/1_R/NLBs_tutorials_Rpubs_etc")
source("make_my_data.R")
```

`make.my.day.ta(first.initial = "N",second.initial = "B")`

```
## [1] "Data for experiment 1 is in the object named experiment1"
## [1] "The mean of the response variable for experiment 1 is 52.2"
## [1] "Type summary(experiment1) to see it"
## [1] "It should look like this:"
## plant treatment photosynth
## Min. : 1.00 control:20 Min. :46.67
## 1st Qu.:10.75 worms :20 1st Qu.:50.91
## Median :20.50 Median :52.25
## Mean :20.50 Mean :52.22
## 3rd Qu.:30.25 3rd Qu.:53.80
## Max. :40.00 Max. :57.55
## [1] ""
## [1] "Data for experiment 2 is in the object experiment2"
## [1] "Data for experiment 3 is in the object experiment3"
```

Was data made?

`ls()`

```
## [1] "experiment1" "experiment2" "experiment3"
## [4] "make.my.day.ta" "simulate.anova.l" "simulate.t.test.a"
## [7] "simulate.t.test.b"
```

`dim(experiment1)`

`## [1] 40 3`

`dim(experiment2)`

`## [1] 40 3`

`dim(experiment3)`

`## [1] 60 3`

The questions say “20 of the plants were grown without worms, and 20 had 5 earthworms added to their pot. One week after the experiment began, we measured photosynthesis on all of the plants.”

This is a 2-sample t-test situation. One set of plants were a control, one set of plants had the experimental worm treatment. ALl were measured at the same time. One set of measurements was made on each plants; plants were not measured twice.

Look at the data

`summary(experiment1)`

```
## plant treatment photosynth
## Min. : 1.00 control:20 Min. :46.67
## 1st Qu.:10.75 worms :20 1st Qu.:50.91
## Median :20.50 Median :52.25
## Mean :20.50 Mean :52.22
## 3rd Qu.:30.25 3rd Qu.:53.80
## Max. :40.00 Max. :57.55
```

```
exp1.output <- t.test(photosynth ~ treatment, data = experiment1)
exp1.output
```

```
##
## Welch Two Sample t-test
##
## data: photosynth by treatment
## t = -5.8385, df = 37.633, p-value = 9.802e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.668953 -2.264230
## sample estimates:
## mean in group control mean in group worms
## 50.49017 53.95676
```

Extract output from t.test object

```
Q1.type.of.test <- exp1.output$method
Q1.test.statistic.type <- "t"
Q1.test.statistic.value <- exp1.output$statistic
Q1.p.value <- exp1.output$p.value
Q1.conclusion <- ifelse(exp1.output$p.value < 0.05,
"reject Ho", "Accept Ho")
Q1.effect.size.value <- exp1.output$estimate[[1]] -
exp1.output$estimate[[2]]
Q1.CI.upper.value <- exp1.output$conf.int[[1]]
Q1.CI.lower.value <- exp1.output$conf.int[[2]]
```

Make dataframe of output

```
Q1.df <- data.frame(type.of.test = Q1.type.of.test
,test.statistic.type = Q1.test.statistic.type
,test.statistic.value = Q1.test.statistic.value
,p.value = Q1.p.value
,conclusion = Q1.conclusion
,effect.size.value = Q1.effect.size.value
,CI.upper.value = Q1.CI.upper.value
,CI.lower.value = Q1.CI.lower.value
)
```

This is a paired-test situation. The question says “We 1st measured photosynthetic rates on all of the plants. We then introduced 5 earthworms to each pot. After 1 week, we re-measured photosynthesis of each plant.” All plants were measured w/o worms, then all plants we re-measurd with worms. There are therefore two seperate measurement on the exact sample plant.

Do the **CORRECT** Analysis of paired t-test

```
exp2.output.CORRECT <- t.test(photosynth ~ treatment,
paired = TRUE, #set for paired t-test
data = experiment2)
exp2.output.CORRECT
```

```
##
## Paired t-test
##
## data: photosynth by treatment
## t = -9.9276, df = 19, p-value = 5.916e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.790967 -3.122561
## sample estimates:
## mean of the differences
## -3.956764
```

Do the **INcorect** Analysis of paired t-test

```
exp2.output.INcorrect <- t.test(photosynth ~ treatment,
data = experiment2)
exp2.output.INcorrect
```

```
##
## Welch Two Sample t-test
##
## data: photosynth by treatment
## t = -5.672, df = 36.482, p-value = 1.827e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.370905 -2.542623
## sample estimates:
## mean in group control mean in group worms
## 50.49017 54.44694
```

Extract output from paired t.test object for the CORRECT analysis

```
Q2corr.type.of.test <- exp2.output.CORRECT$method
Q2corr.test.statistic.type <- "t"
Q2corr.test.statistic.value <- exp2.output.CORRECT$statistic
Q2corr.p.value <- exp2.output.CORRECT$p.value
Q2corr.conclusion <- ifelse(exp2.output.CORRECT$p.value < 0.05,
"reject Ho", "Accept Ho")
Q2corr.effect.size.value <- exp2.output.CORRECT$estimate[[1]]
Q2corr.CI.upper.value <- exp2.output.CORRECT$conf.int[[1]]
Q2corr.CI.lower.value <- exp2.output.CORRECT$conf.int[[2]]
```

Extract output from paired t.test object for the INcorrect analysis

```
# Q2INcorr.type.of.test <- exp2.output.INcorrect$method
# Q2INcorr.test.statistic.type <- "t"
# Q2INcorr.test.statistic.value <- exp2.output.INcorrect$statistic
# Q2INcorr.p.value <- exp2.output.INcorrect$p.value
# Q2INcorr.conclusion <- ifelse(exp2.output.INcorrect$p.value < 0.05,
# "reject Ho", "Accept Ho")
# Q2INcorr.effect.size.value <- exp2.output.INcorrect$estimate[[2]]
# exp2.output.INcorrect$estimate[[1]]-
#
# Q2INcorr.CI.upper.value <- exp2.output.INcorrect$conf.int[[1]]
# Q2INcorr.CI.lower.value <- exp2.output.INcorrect$conf.int[[2]]
```

`summary(experiment3)`

```
## plant treatment photosynth
## Min. : 1.00 control :20 Min. : 30.36
## 1st Qu.:15.75 worm.spp.1:20 1st Qu.: 57.29
## Median :30.50 worm.spp.2:20 Median : 72.36
## Mean :30.50 Mean : 78.86
## 3rd Qu.:45.25 3rd Qu.: 94.89
## Max. :60.00 Max. :201.80
```

```
#Null model
Q3.m.null <- lm( log(photosynth) ~ 1, data = experiment3)
#alt model
Q3.m.alt <- lm( log(photosynth) ~ treatment, data = experiment3)
```

ANOVA Omnibus test

```
Q3.omnibus <- anova(Q3.m.null,
Q3.m.alt)
```

pairwise t-test

```
pairwise.t.default <- pairwise.t.test(x = experiment3$photosynth,
g = experiment3$treatment)#default = hom
pairwise.t.none <- pairwise.t.test(x = experiment3$photosynth,
g = experiment3$treatment,
p.adjust.method = "none")
pairwise.t.bonf <- pairwise.t.test(x = experiment3$photosynth,
g = experiment3$treatment,
p.adjust.method = "bonferroni")
min(pairwise.t.default$p.value, na.rm = T)
```

`## [1] 3.382145e-07`

`min(pairwise.t.none$p.value, na.rm = T)`

`## [1] 1.127382e-07`

`min(pairwise.t.bonf$p.value, na.rm = T)`

`## [1] 3.382145e-07`

TukeyHSD

```
Q3.m.alt.aov <- aov(log(photosynth) ~ treatment, data = experiment3)
Q3.Tukey <- TukeyHSD(Q3.m.alt.aov)
```

```
Q3.F <- Q3.omnibus$F[[2]]
Q3.P <- Q3.omnibus$"Pr(>F)"[[2]]
Q3.conclusion <- ifelse(Q3.P < 0.05, "Reject Ho", "Fail to reject Ho")
Q3.pairwise.p.uncorr <- min(pairwise.t.none$p.value, na.rm = T)
Q3.pairwise.p.bonffe <- min(pairwise.t.bonf$p.value, na.rm = T)
Q3.pairwise.p.tukey <- min(Q3.Tukey$treatment[,4])
i.max <- which.max(max(Q3.Tukey$treatment[,1]))
Q3.Tukey.ES <- Q3.Tukey$treatment[i.max, 1] #largest effect size
Q3.Tukey.CI <- Q3.Tukey$treatment[i.max, 2:3] #CI for largest ES
Q3.Tukey.CI.lwr <- Q3.Tukey.CI[1]
Q3.Tukey.CI.upr <- Q3.Tukey.CI[2]
```