Pre-midterm practice assignment

ESS 103

Question 1

Download and import the Nitrate.csv file from Moodle. The dataset includes results of an agricultural researcher who measured soil nitrate levels (\(\text{mg/kg}\)) in 15 agricultural plots before (Pre_Nitrate) and after (Post_Nitrate) planting a cover crop. We are interested whether nitrates level changed significantly before and after the plantation.

1. State your hypotheses.
2. Choose the relevent test after testing all required hypotheses.
3. Run the test and discuss the results.

Solution

I use a significance level of 0.05 throughout. A small p-value provides evidence against the null hypothesis. A large p-value does not prove that the null hypothesis is true. Independence between sampling units depends on the study design and cannot be established using a normality test. I assume that different plots, samples, and lakes are independent as appropriate to each question.

Hypotheses. Let the difference be post-planting minus pre-planting nitrate. My null hypothesis is that the population mean difference equals zero. My alternative hypothesis is that the mean difference is not zero (a two-sided test).

nitrate <- read.csv("Nitrate.csv", row.names = 1)
stopifnot(nrow(nitrate) == 15, !anyNA(nitrate),
          !anyDuplicated(nitrate$Plot_ID))
nitrate$Difference <- nitrate$Post_Nitrate - nitrate$Pre_Nitrate
summary(nitrate$Difference)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -8.280  -4.610  -3.140  -2.025   0.645   8.430 

shapiro.test(nitrate$Difference)

    Shapiro-Wilk normality test

data:  nitrate$Difference
W = 0.94607, p-value = 0.4648

qqnorm(nitrate$Difference, main = "Normal Q-Q plot of nitrate differences")
qqline(nitrate$Difference)

Test choice. The measurements are paired because the same 15 plots were measured twice. For a paired t-test, normality concerns the differences, rather than the two sets of measurements separately. The Shapiro-Wilk test gives p = 0.4648, so it does not indicate a departure from normality. Equal pre- and post-planting variances are not required. I use a paired t-test, assuming that the differences from separate plots are independent.

colMeans(nitrate[c("Pre_Nitrate", "Post_Nitrate")])

 Pre_Nitrate Post_Nitrate 
    25.68133     23.65667 

t.test(nitrate$Post_Nitrate, nitrate$Pre_Nitrate, paired = TRUE)

    Paired t-test

data:  nitrate$Post_Nitrate and nitrate$Pre_Nitrate
t = -1.8384, df = 14, p-value = 0.08731
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -4.3867546  0.3374213
sample estimates:
mean difference 
      -2.024667 

Conclusion. Mean nitrate decreased from 25.68 to 23.66 mg/kg, a mean change of -2.02 mg/kg. However, the paired t-test gives t(14) = -1.84, p = 0.0873. The 95% confidence interval for the mean change is [-4.39, 0.34] mg/kg. I fail to reject the null hypothesis: these data do not provide sufficient evidence of a change in mean nitrate at the 5% level.

Question 2

Microplastic particle concentrations (particles/L) were sampled from two independent rivers: River A and River B. Import the dataset Microplastics.csv from Moodle and test whether the Microplastic particle concentrations differ significantly in the 2 sites.

1. State your hypotheses.
2. Choose the relevent test after testing all required hypotheses.
3. Run the test and discuss the results.

Solution

Hypotheses. H0: the population mean microplastic concentration is equal in River A and River B. H1: the population means differ.

micro <- read.csv("Microplastics.csv", row.names = 1)
micro$River <- factor(micro$River)
stopifnot(!anyNA(micro))
table(micro$River)

River_A River_B 
     15      15 

by(micro$Microplastics_per_L, micro$River, shapiro.test)

micro$River: River_A

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.92683, p-value = 0.2445

------------------------------------------------------------ 
micro$River: River_B

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.89249, p-value = 0.07314

fligner.test(Microplastics_per_L ~ River, data = micro)

    Fligner-Killeen test of homogeneity of variances

data:  Microplastics_per_L by River
Fligner-Killeen:med chi-squared = 6.7388, df = 1, p-value = 0.009434

par(mfrow = c(1, 2))
for (river in levels(micro$River)) {
  values <- micro$Microplastics_per_L[micro$River == river]
  qqnorm(values, main = river)
  qqline(values)
}

par(mfrow = c(1, 1))

Test choice. There are 15 observations per river, and the two groups are independent according to the question. Shapiro-Wilk p-values are 0.2445 and 0.0731. Neither rejects normality, although the small sample sizes mean these checks have limited power. The Fligner-Killeen variance test gives p = 0.00943, indicating unequal variances. I therefore use Welch’s independent-samples t-test, which does not require equal variances.

aggregate(Microplastics_per_L ~ River, data = micro, FUN = mean)

t.test(Microplastics_per_L ~ River, data = micro, var.equal = FALSE)

    Welch Two Sample t-test

data:  Microplastics_per_L by River
t = 2.8754, df = 19.58, p-value = 0.009488
alternative hypothesis: true difference in means between group River_A and group River_B is not equal to 0
95 percent confidence interval:
 1.376911 8.689756
sample estimates:
mean in group River_A mean in group River_B 
             9.513333              4.480000 

Conclusion. River A has a mean of 9.51 particles/L and River B has a mean of 4.48 particles/L. Welch’s test gives t = 2.88, df = 19.58, p = 0.00949. I reject H0. River A has a significantly higher mean concentration, with an estimated difference of 5.03 particles/L (95% CI: 1.38 to 8.69 particles/L). This comparison concerns the sampled rivers; it does not establish why they differ.

Question 3

Dissolved Organic Carbon (DOC) was sampled across four distinct wetland habitat types: Bog, Fen, Marsh, and Swamp. Choose and perform the correct test (after checking all assumptions) to determine whether DOC levels differ significantly across wetland types. Dataset (DOC.csv) is available on Moodle.

Solution

Hypotheses. For the rank-based test, H0 is that DOC distributions are the same across Bog, Fen, Marsh, and Swamp. H1 is that at least one distribution differs. If the distributions have similar shapes, this can be interpreted as a difference in location or medians.

doc <- read.csv("DOC.csv", row.names = 1)
doc$Wetland_Type <- factor(doc$Wetland_Type)
stopifnot(!anyNA(doc))
table(doc$Wetland_Type)

  Bog   Fen Marsh Swamp 
   10    10    10    10 

doc_model <- lm(DOC_mgL ~ Wetland_Type, data = doc)
shapiro.test(residuals(doc_model))

    Shapiro-Wilk normality test

data:  residuals(doc_model)
W = 0.84492, p-value = 6.788e-05

fligner.test(DOC_mgL ~ Wetland_Type, data = doc)

    Fligner-Killeen test of homogeneity of variances

data:  DOC_mgL by Wetland_Type
Fligner-Killeen:med chi-squared = 4.6936, df = 3, p-value = 0.1957

boxplot(DOC_mgL ~ Wetland_Type, data = doc,
        main = "DOC across wetland types", xlab = "Wetland type",
        ylab = "DOC (mg/L)", col = "lightblue")

qqnorm(residuals(doc_model), main = "Normal Q-Q plot of DOC residuals")
qqline(residuals(doc_model))

Test choice. Each habitat has 10 observations. The normality test on one-way ANOVA residuals gives p = 0.0000679, indicating non-normal residuals. The variance test gives p = 0.1957 and does not detect unequal variances. Given the non-normality and high DOC observations, I use the Kruskal-Wallis test for four independent groups. It does not require normally distributed residuals. I interpret it as a test of distributions because equal distribution shapes have not been established.

aggregate(DOC_mgL ~ Wetland_Type, data = doc, FUN = median)

kruskal.test(DOC_mgL ~ Wetland_Type, data = doc)

    Kruskal-Wallis rank sum test

data:  DOC_mgL by Wetland_Type
Kruskal-Wallis chi-squared = 16.639, df = 3, p-value = 0.0008385

pairwise.wilcox.test(doc$DOC_mgL, doc$Wetland_Type,
                     p.adjust.method = "holm", exact = FALSE)

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  doc$DOC_mgL and doc$Wetland_Type 

      Bog    Fen    Marsh 
Fen   0.6776 -      -     
Marsh 0.0026 0.0364 -     
Swamp 0.0453 0.2270 0.2270

P value adjustment method: holm 

Conclusion. The Kruskal-Wallis test gives chi-squared = 16.64, df = 3, p = 0.000839. I reject H0: DOC differs across wetland types. Pairwise Wilcoxon tests with Holm correction identify differences between Bog and Marsh (adjusted p = 0.0026), Fen and Marsh (0.0364), and Bog and Swamp (0.0453). The sample DOC values tend to be higher in Bog than in Marsh and Swamp, and higher in Fen than in Marsh. The remaining pairs are not significant after correction.

Question 4

An ecology experiment measured native grass biomass under two temperature conditions (Ambient vs. Warmed) and two nitrogen levels (Control vs. Enriched). You can find the results in biomass.csv on Moodle.

1. Create a ggplot2 scatter plot showing biomass across the 2 temperature states, and add Nitrogen levels in different colors. Make sure you add the respective graph and axes titles.

2. Perform the correct test (based on your hypotheses and assumption testing) to check for impact of temperature, nitrogen, and their interaction on plant biomass. If needed, perform a post-hoc test to determine exactly which combinations produce significant results.

Solution

Scatter plot

biomass <- read.csv("biomass.csv", row.names = 1)
biomass$Temperature <- factor(biomass$Temperature,
                              levels = c("Ambient", "Warmed"))
biomass$Nitrogen <- factor(biomass$Nitrogen,
                           levels = c("Control", "Enriched"))
stopifnot(!anyNA(biomass))
ggplot(biomass, aes(x = Temperature, y = Plant_Biomass, color = Nitrogen)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15,
                                            dodge.width = 0.55, seed = 103),
             size = 2.6, alpha = 0.85) +
  labs(title = "Native grass biomass under temperature and nitrogen treatments",
       x = "Temperature condition", y = "Plant biomass (reported units)",
       color = "Nitrogen treatment") +
  theme_minimal(base_size = 12)

The file does not specify a biomass unit, so I do not assign an unreported unit to the axis.

Hypotheses and assumptions

The temperature null hypothesis is that marginal mean biomass is equal under ambient and warmed conditions. The nitrogen null hypothesis is that marginal mean biomass is equal under control and enriched conditions. The interaction null hypothesis is that the temperature effect is the same at both nitrogen levels. Each alternative states that the corresponding effect is present.

with(biomass, table(Temperature, Nitrogen))

           Nitrogen
Temperature Control Enriched
    Ambient      10       10
    Warmed       10       10

biomass_model <- aov(Plant_Biomass ~ Temperature * Nitrogen, data = biomass)
shapiro.test(residuals(biomass_model))

    Shapiro-Wilk normality test

data:  residuals(biomass_model)
W = 0.94163, p-value = 0.03925

biomass$Group <- interaction(biomass$Temperature, biomass$Nitrogen)
fligner.test(Plant_Biomass ~ Group, data = biomass)

    Fligner-Killeen test of homogeneity of variances

data:  Plant_Biomass by Group
Fligner-Killeen:med chi-squared = 1.6302, df = 3, p-value = 0.6526

par(mfrow = c(1, 2))
plot(biomass_model, which = 1)
plot(biomass_model, which = 2)

par(mfrow = c(1, 1))

Test choice. I use a two-way ANOVA because there are two categorical factors and a continuous response, and I need to test their interaction. The design is balanced with 10 observations in each combination. Equal variances are plausible (Fligner-Killeen p = 0.6526). However, residual normality is rejected at the 5% level (Shapiro-Wilk p = 0.0393). Thus, the normality assumption is imperfect. I retain the original-scale ANOVA because a balanced design with similar variances is reasonably robust to moderate non-normality, but its p-values should be interpreted with that limitation. I assume each row represents a separate experimental unit.

summary(biomass_model)

                     Df Sum Sq Mean Sq F value  Pr(>F)    
Temperature           1  855.2   855.2  269.01 < 2e-16 ***
Nitrogen              1 1170.4  1170.4  368.18 < 2e-16 ***
Temperature:Nitrogen  1   37.2    37.2   11.69 0.00158 ** 
Residuals            36  114.4     3.2                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aggregate(Plant_Biomass ~ Temperature + Nitrogen, data = biomass, FUN = mean)

TukeyHSD(aov(Plant_Biomass ~ Group, data = biomass))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Plant_Biomass ~ Group, data = biomass)

$Group
                                   diff        lwr       upr     p adj
Warmed.Control-Ambient.Control    7.320  5.1725466  9.467453 0.0000000
Ambient.Enriched-Ambient.Control  8.891  6.7435466 11.038453 0.0000000
Warmed.Enriched-Ambient.Control  20.066 17.9185466 22.213453 0.0000000
Ambient.Enriched-Warmed.Control   1.571 -0.5764534  3.718453 0.2180468
Warmed.Enriched-Warmed.Control   12.746 10.5985466 14.893453 0.0000000
Warmed.Enriched-Ambient.Enriched 11.175  9.0275466 13.322453 0.0000000

Conclusion. Temperature is significant, F(1,36) = 269.01, p < 0.001, and nitrogen is significant, F(1,36) = 368.18, p < 0.001. Their interaction is also significant, F(1,36) = 11.69, p = 0.00158. Because the interaction is significant, I interpret the treatment combinations rather than relying only on the main effects.

Mean biomass is 24.43 for Ambient-Control, 31.75 for Warmed-Control, 33.32 for Ambient-Enriched, and 44.49 for Warmed-Enriched. Warming increases mean biomass by 7.32 units under control nitrogen and by 11.18 units under enriched nitrogen. This larger increase under enrichment explains the positive interaction on the original biomass scale.

Tukey’s comparisons show that all pairs differ significantly except Warmed-Control versus Ambient-Enriched (adjusted p = 0.218). Warmed-Enriched produces the highest mean biomass. These ANOVA and Tukey conclusions share the residual-normality limitation described above.

Question 5

Soil moisture content (%) was monitored across 12 forest plots over four consecutive observation periods: Week 0, Week 2, Week 4, and Week 6. Data is found in soil_mosture.csv on Moodle.

1. State your hypotheses
2. Choose the correct test bearing in mind each plot was tested multiple times. Check the assumptions.
3. Interpret your findings.

Solution

Hypotheses. H0: the population mean soil moisture is equal at Week 0, Week 2, Week 4, and Week 6. H1: at least one time-point mean differs.

moisture <- read.csv("soil_moisture.csv", row.names = 1)
moisture$Subject <- factor(moisture$Subject)
moisture$Time <- factor(moisture$Time,
                        levels = c("Week0", "Week2", "Week4", "Week6"))
stopifnot(!anyNA(moisture), all(table(moisture$Subject, moisture$Time) == 1))
# Rows are plots and columns are observation times, preserving pairing.
moisture_wide <- as.matrix(xtabs(Moisture_Pct ~ Subject + Time, data = moisture))
colMeans(moisture_wide)

   Week0    Week2    Week4    Week6 
48.85167 44.26250 39.25000 34.59750 

Test choice and assumptions. The same 12 plots were measured four times, so an ordinary independent-groups ANOVA would be inappropriate. I use a one-way repeated-measures ANOVA. This requires independent plots, approximately normal within-plot contrasts, and sphericity (equal variances of differences between time points). There are no missing or duplicated plot-time measurements.

time_pairs <- combn(ncol(moisture_wide), 2, simplify = FALSE)
normality <- do.call(rbind, lapply(time_pairs, function(j) {
  differences <- moisture_wide[, j[1]] - moisture_wide[, j[2]]
  data.frame(Comparison = paste(colnames(moisture_wide)[j], collapse = " - "),
             Shapiro_p = shapiro.test(differences)$p.value)
}))
normality

par(mfrow = c(2, 3))
for (j in time_pairs) {
  differences <- moisture_wide[, j[1]] - moisture_wide[, j[2]]
  qqnorm(differences, main = paste(colnames(moisture_wide)[j], collapse = " - "))
  qqline(differences)
}

par(mfrow = c(1, 1))
# Multivariate representation of the repeated measurements.
moisture_multivariate <- lm(moisture_wide ~ 1)
mauchly.test(moisture_multivariate, X = ~ 1)

    Mauchly's test of sphericity
    Contrasts orthogonal to
    ~1


data:  SSD matrix from lm(formula = moisture_wide ~ 1)
W = 0.79074, p-value = 0.8097

Mauchly’s test gives W = 0.7907, p = 0.8097, so there is no evidence against sphericity. Five of six difference-normality tests have p > 0.05; Week 0 minus Week 2 is borderline at p = 0.0492. These small-sample checks do not prove multivariate normality. I report repeated-measures ANOVA and check whether a nonparametric Friedman test supports the overall conclusion.

moisture_model <- aov(Moisture_Pct ~ Time + Error(Subject / Time), data = moisture)
summary(moisture_model)

Error: Subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11    552   50.18               

Error: Subject:Time
          Df Sum Sq Mean Sq F value Pr(>F)    
Time       3 1369.9   456.6   185.5 <2e-16 ***
Residuals 33   81.2     2.5                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Also displays Greenhouse-Geisser and Huynh-Feldt corrections.
anova(moisture_multivariate, X = ~ 1, test = "Spherical")

friedman.test(Moisture_Pct ~ Time | Subject, data = moisture)

    Friedman rank sum test

data:  Moisture_Pct and Time and Subject
Friedman chi-squared = 36, df = 3, p-value = 7.488e-08

Conclusion. Mean moisture decreases from 48.85% to 44.26%, 39.25%, and 34.60% over the four observations. The repeated-measures ANOVA gives F(3,33) = 185.49, p < 0.001, so I reject H0. The Friedman sensitivity check also finds an overall time effect (chi-squared = 36.00, df = 3, p = 7.49e-8). Moisture differs significantly across time, with a clear decrease in the observed means. An overall test alone does not establish that every pair of weeks differs.

Question 6

A lake clarity index (scale 1 to 10) was measured across 10 monitoring lakes during three distinct seasons: Spring, Summer, and Autumn. (The dataset is called water_clarity.csv).

1. State your hypotheses.
2. Report whether seasonal variations in water clarity are statistically significant.

Solution

Hypotheses. For the Friedman test, H0 is that the seasons have the same clarity rank distribution within lakes (no systematic seasonal effect). H1 is that at least one season differs.

clarity <- read.csv("water_clarity.csv", row.names = 1)
clarity$Lake_ID <- factor(clarity$Lake_ID)
clarity$Season <- factor(clarity$Season, levels = c("Spring", "Summer", "Autumn"))
stopifnot(!anyNA(clarity), all(table(clarity$Lake_ID, clarity$Season) == 1))
by(clarity$Clarity_Index, clarity$Season, shapiro.test)

clarity$Season: Spring

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.93989, p-value = 0.5518

------------------------------------------------------------ 
clarity$Season: Summer

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.69944, p-value = 0.0008672

------------------------------------------------------------ 
clarity$Season: Autumn

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.90737, p-value = 0.2634

range(clarity$Clarity_Index)

[1] 0.2 6.1

sum(clarity$Clarity_Index < 1 | clarity$Clarity_Index > 10)

[1] 7

boxplot(Clarity_Index ~ Season, data = clarity,
        main = "Lake clarity across seasons", xlab = "Season",
        ylab = "Clarity index", col = "lightblue")

Test choice. These are repeated measurements on the same 10 lakes. I use the Friedman test, with lake as the blocking variable, because it accommodates paired measurements over three seasons and does not require normality or sphericity. The index is bounded and Summer values are strongly non-normal (Shapiro-Wilk p = 0.000867). I assume different lakes are independent, and each lake has one measurement in every season.

Data note. The question describes a 1-10 scale, but the supplied file contains values below 1 (the observed range is 0.2-6.1). I keep the supplied values for this analysis; the scale description should be checked against the original data source.

aggregate(Clarity_Index ~ Season, data = clarity, FUN = median)

friedman.test(Clarity_Index ~ Season | Lake_ID, data = clarity)

    Friedman rank sum test

data:  Clarity_Index and Season and Lake_ID
Friedman chi-squared = 11.4, df = 2, p-value = 0.003346

Conclusion. The Friedman test gives chi-squared = 11.40, df = 2, p = 0.00335. I reject H0 and conclude that seasonal variation in clarity is statistically significant. The sample medians are 2.70 in Spring, 0.65 in Summer, and 4.75 in Autumn. These summaries suggest lower clarity in Summer and higher clarity in Autumn, but the overall test does not determine which individual season pairs differ.