1 Take Home Test 2

1.1 Ecological Background

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.


1.2 Descriptions of Experiments

1.2.1 Experiment 1: (dataframe = experiment1)

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.


1.2.2 Experiment 2: (dataframe = experiment2)

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.


1.2.3 Experiment 3: (dataframe = experiment3)

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.


1.3 Data Setup

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")

1.4 Generate data

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

1.5 Sample Answer Question 1: 2-sample t-test

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

1.5.1 Question 1 2-sample t-test

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
)

1.6 Sample Answer Question 2: 2-sample PAIRED t-test

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.

1.6.1 Question 2: paired t-test

1.6.1.1 Correct analysis: paired t-test

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

1.6.2 Incorrect analysis

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

1.6.2.1 Output from correct analysis

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]]

1.6.2.2 Output from INcorrect analysis

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]]



1.7 Sample Answer Question 3: ANOVA on log-transformed data

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

1.7.1 ANOVA

#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]