Lab 09

Homework

A| Species Diversity

A.1| Hypothesis & Test

Null hypothesis:

\(\mu_{low}=\mu_{ctrl}=\mu{high}\)

IV: Disturbance - factor with 3 levels

DV: Species Diversity

I will use a one factor ANOVA because our independent variable is categorical (1 factor) with 3 levels and our dependent variable is continuous numerical.

A.2| Data Read

#data load
plants <- read.csv("Data/disturbanceplants.csv")
glimpse(plants)

## Rows: 90
## Columns: 5
## $ Treatment         <chr> "Control", "Control", "Control", "Control", "Control…
## $ Plant_growth_cm   <dbl> 32.48357, 29.30868, 33.23844, 37.61515, 28.82923, 28…
## $ Soil_nutrient_ppm <dbl> 151.9416, 169.3729, 135.9589, 143.4468, 142.1578, 12…
## $ Species_diversity <dbl> 3.812834, 3.071421, 2.964554, 3.741236, 3.388269, 3.…
## $ Biomass_g         <dbl> 572.0637, 428.2069, 558.1582, 500.5117, 450.9246, 52…

A.3| Assumptions Check

#fit model to obtain residuals if needed
disturbanceModel <- lm(Species_diversity~Treatment, data=plants)

#normality check
par(mfrow=c(2,2))
plot(disturbanceModel)

shapiro.test(disturbanceModel$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  disturbanceModel$residuals
## W = 0.9807, p-value = 0.2021

bartlett.test(Species_diversity~Treatment, data=plants)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Species_diversity by Treatment
## Bartlett's K-squared = 5.5025, df = 2, p-value = 0.06385

What do the assumption checks tell me?

The assumption checks showed no significant evidence of nonnormal residuals (Shapiro–Wilk, p = 0.202) or unequal variances (Bartlett, p = 0.064).

A.4| LM & ANOVA

#model summary
summary(disturbanceModel) #significance exists between high and control AND low and control ...

## 
## Call:
## lm(formula = Species_diversity ~ Treatment, data = plants)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.16848 -0.45755  0.02019  0.38194  1.89357 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 3.5328     0.1212  29.146  < 2e-16 ***
## TreatmentHigh_disturbance  -1.1332     0.1714  -6.611 2.95e-09 ***
## TreatmentLow_disturbance   -0.5040     0.1714  -2.940   0.0042 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6639 on 87 degrees of freedom
## Multiple R-squared:  0.3353, Adjusted R-squared:   0.32 
## F-statistic: 21.94 on 2 and 87 DF,  p-value: 1.929e-08

#R^2: 32% of variation explained

Conclusion:

Disturbance significantly affected species diversity (one-factor ANOVA: F₂,₈₇ = 21.94, p < 0.001), explaining 33.5% of the variation. Therefore, I rejected the null hypothesis that all treatment means were equal. The model summary showed that diversity was significantly lower under high disturbance than under control conditions (difference = −1.13, p < 0.001) and low disturbance (difference = −0.50, adjusted p < 0.05).

A.5| Linear Equation

Write the linear equation for species diversity based on disturbance:

\(y=\beta_0+\beta_{1Low}+\beta_{2High}\rightarrow y=3.53-0.504(Low)-1.13(High)\), where under low disturbance we plug in a 1, else 0 OR under high disturbance we plug in a 1, else 0.

A.6| Predicted Diversity for Low Disturbance

#type the calculation by hand using the model coefficients
y <- 3.53-0.504-0
y

## [1] 3.026

Predicted species diversity:

Species diversity is predicted to be 3.026 under low disturbance.

B| Pairwise comparisons

B.1| Model Read

#prepare model for means comparisons
m1 <- aov(Species_diversity~Treatment, data=plants)

B.2| Tukey Test

#Tukey test
TukeyHSD(m1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Species_diversity ~ Treatment, data = plants)
## 
## $Treatment
##                                        diff        lwr         upr     p adj
## High_disturbance-Control         -1.1331985 -1.5419448 -0.72445219 0.0000000
## Low_disturbance-Control          -0.5039744 -0.9127207 -0.09522814 0.0115951
## Low_disturbance-High_disturbance  0.6292241  0.2204778  1.03797034 0.0012015

Which groups are different?

All groups are different, including: Low-Control, High-Control, and Low-High

B.3| SNK Test

snk <- SNK.test(disturbanceModel, "Treatment")

snk$groups #woweeeee! this is lovely. 

##                  Species_diversity groups
## Control                   3.532796      a
## Low_disturbance           3.028822      b
## High_disturbance          2.399598      c

Which groups are different?

All groups are still different.

B.4| Figure 1

#add significance letter
letters <- data.frame(Treatment = c("Control", "High_disturbance", "Low_disturbance"),label = c("A", "C", "B"),height = c(4.2, 4.2, 4.7))

#boxplot it
ggplot(plants, aes(Treatment, Species_diversity, color=Treatment))+
  geom_boxplot()+
  theme_classic()+
  theme(legend.position = "none")+
  geom_text(data = letters, aes(Treatment, height, label=label))+
  labs(x="Treatment", y="Species Diversity")+
  scale_color_manual(values=c("#715cb2","#a6cd64","#f9738b"))+
  scale_x_discrete(labels = c("Control"="Control", "High_disturbance"="High Disturbance", "Low_disturbance"="Low Disturbance"))

Figure 1. Species Diversity Across Disturbance Treatments. (2, 87) = 21.94, p < .001. Treatment explained 33.5% of the variation (R² = .335; adjusted R² = .320). The fitted model was ŷ = 3.5328 − 1.1332(High) − 0.5040(Low), where treatment indicators equal 1 when applicable and 0 otherwise; Control is the reference group. Diversity was highest in Control, intermediate under low disturbance, and lowest under high disturbance. Different letters indicate significant pairwise differences according to Tukey’s HSD (all adjusted p < .05).

C| Data Simulation

C.1| Research Question & Variables

Research Q: Are there treatment differences in soil organic carbon across detrital input and removal plots?

IV: Treatments: No Litter, Control, Double Litter, No Roots | 1-Factor, 4-levels.

DV: Soil Organic Carbon (\(\frac{mg C}{g Soil}\))

C.2| Create Similarity Between Groups

I need to create data that has similar means between those two groups to force that similarity.

C.3| Simulate Data

Control <- rnorm(18, 20, 3)
DL <- rnorm(18, 18, 2.8) #double litter that will be similar in means to control
NL <- rnorm(18, 8, 2.9)
NR <- rnorm(18, 2, 2.5) #we dont have this treatment at our site, but i have read that microbes are associated in higher numbers near roots.

#add all to DF 
soil <- data.frame(Control, DL, NL, NR)

#change dataframe so that there are two columns, one with treatment and other with SOC amount
soil <- stack(data.frame(Control, DL, NL, NR))
#change column names
names(soil) <- c("SOC", "Treatment")

str(soil) #perfect!

## 'data.frame':    72 obs. of  2 variables:
##  $ SOC      : num  16 18.4 17.5 15.8 19.2 ...
##  $ Treatment: Factor w/ 4 levels "Control","DL",..: 1 1 1 1 1 1 1 1 1 1 ...

C.4| Assumptions Check

#fit model to obtain residuals if needed
carbonModel <- lm(SOC~Treatment, data=soil)

#now use model to check assumptions (which we know are all normal with similar variance because...we set them to be that way)
par(mfrow=c(2,2)) #set it up so the 4 plots are in a 2x2 grid to make it easier to see
plot(carbonModel)

#double check just in case
shapiro.test(carbonModel$residuals) #normal!

## 
##  Shapiro-Wilk normality test
## 
## data:  carbonModel$residuals
## W = 0.97957, p-value = 0.2916

bartlett.test(SOC~Treatment, data=soil)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  SOC by Treatment
## Bartlett's K-squared = 3.313, df = 3, p-value = 0.3458

The assumption checks tell me that I can proceed with using the linear model function/1-factor ANOVA because we assume that samples were independently sampled, the SOC response variable residuals follow a normal distribution, and the data is approximately equal in variance around each group mean.

C.5| LM

#model summary
summary(carbonModel) #control is the reference group beaccuse of alphabet, which is perfect. 

## 
## Call:
## lm(formula = SOC ~ Treatment, data = soil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7135 -2.3208  0.2359  2.0461  5.2544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.2571     0.6474  28.201   <2e-16 ***
## TreatmentDL  -0.1034     0.9155  -0.113     0.91    
## TreatmentNL -11.0147     0.9155 -12.031   <2e-16 ***
## TreatmentNR -16.5632     0.9155 -18.091   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.747 on 68 degrees of freedom
## Multiple R-squared:  0.8775, Adjusted R-squared:  0.8721 
## F-statistic: 162.3 on 3 and 68 DF,  p-value: < 2.2e-16

#DL is not statistically different from control with respect to mean, and that is great because it is what was set
#I want to know more, is NR different from NL? In order to know that I have to run SNK OR TUKEY

C.6| SNK Test

soilmod2 <- aov(SOC ~ Treatment, data = soil)
snk <- SNK.test(soilmod2, "Treatment")
snk$groups

##               SOC groups
## Control 18.257119      a
## DL      18.153749      a
## NL       7.242443      b
## NR       1.693937      c

Conclusion

Ahhh YEAH! There you can see what letters you can put on your boxplots (which is extremely helpful to what I am working on with my research. This is incredibly useful. Here you can see that Control and DL share similar means and NL and NR all have different means from control and one another.