thesis_code

Honors Thesis Data Summary

The impact of variable versus constant winter snow cover on maple leaf litter decomposition

Outline of Data Summary

1. Leaf Litter Mass Loss
2. Leaf Litter C:N, %C, %N
3. Respiration (Rs)

• Variable Temperature
• Variable Moisture

4. Physical Conditions

• Temperature
• Field Moisture

Directory

knitr::opts_knit$set(root.dir = "//Users//abbey//Desktop//thesis//data//REPORT//data_summary")

Libraries

library(ggplot2)
library(ggpubr)

## Loading required package: magrittr

library(plyr)

## 
## Attaching package: 'plyr'

## The following object is masked from 'package:ggpubr':
## 
##     mutate

library(car)

## Loading required package: carData

library(gridExtra)
library(knitr)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:plyr':
## 
##     here

## The following object is masked from 'package:base':
## 
##     date

library(pbkrtest)

## Loading required package: lme4

## Loading required package: Matrix

## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car

library(ggpubr)
library(rstatix)

## 
## Attaching package: 'rstatix'

## The following object is masked from 'package:plyr':
## 
##     mutate

## The following object is masked from 'package:stats':
## 
##     filter

library(formattable)
library(reshape2)
library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:formattable':
## 
##     area

## The following object is masked from 'package:rstatix':
## 
##     select

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

library(lme4)
library(wesanderson)
library(pander)
library(emmeans)

## Welcome to emmeans.
## NOTE -- Important change from versions <= 1.41:
##     Indicator predictors are now treated as 2-level factors by default.
##     To revert to old behavior, use emm_options(cov.keep = character(0))

library(lsmeans)

## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.

1. Leaf Litter Mass Loss

Preliminary Organization of Data:

#data upload, organize variables and define categorical variables 
mass_loss <- read.csv("mass_loss_compiled.csv")
#class(mass_loss$day_since_deploy)
mass_loss$day_since_deploy <- as.numeric(mass_loss$day_since_deploy)
#str(mass_loss)
mass_loss$ML_harvest=as.factor(mass_loss$ML_harvest)

#create preliminary boxplot describing percent mass loss from original dry mass 
mass_loss_plot <- ggplot(mass_loss, aes(x = factor(day_since_deploy), y = mass_loss_from_OG, fill = factor(treatment)))+
  geom_boxplot()+
  ggtitle("Mean Leaf Litter Percent Mass Remaining by Treatment") +
  scale_y_continuous(name = "Mass Remaining (%)",
                           breaks = seq(70, 145, 10),
                           limits=c(70, 145)) +
  scale_x_discrete(name = "Days in Field")+
  scale_fill_discrete(name = "Treatment", labels = c("CC", "CS"))+
  theme_classic()
mass_loss_plot

ANOVA <- mass_loss
ANOVA_ml <- ANOVA[c(1,3,5,15)]
ANOVA_ml$harvest <- ANOVA_ml$ML_harvest
ANOVA_final<- ANOVA_ml[c(2,3,4,5)]
ANOVA_final$harvest <- factor(ANOVA_final$harvest)
library(lubridate)
ANOVA_final$datefinal <- mdy(ANOVA_final$date)

#str(ANOVA_final)

table(ANOVA_final$treatment, ANOVA_final$harvest)

##     
##      0 1 2 3
##   C  5 5 5 4
##   SC 5 5 6 5

Finalized Figures:

#data statistical summary 
data_sum <- ddply(ANOVA_final, c("datefinal", "treatment"), summarise,
               N    = length(mass_loss_from_OG),
               mean = mean(mass_loss_from_OG),
               sd   = sd(mass_loss_from_OG),
               se   = sd / sqrt(N))
formattable(data_sum)

datefinal	treatment	N	mean	sd	se
2020-01-14	C	5	100.00000	0.000000	0.000000
2020-01-14	SC	5	100.00000	0.000000	0.000000
2020-02-26	C	5	102.98434	27.310703	12.213717
2020-02-26	SC	5	93.79780	15.473915	6.920145
2020-03-16	C	5	107.01334	6.765307	3.025537
2020-03-16	SC	6	88.59614	13.522174	5.520405
2020-04-07	C	4	93.16938	13.925831	6.962915
2020-04-07	SC	5	96.57595	25.003275	11.181805

data_sum$datefinal <- mdy(data_sum$datefinal)

## Warning: All formats failed to parse. No formats found.

#str(data_sum)

#figure
ML <-ggline(ANOVA_final, x = "datefinal", y = "mass_loss_from_OG", 
       add = c("mean_se", "jitter"),
       color = "treatment", palette = "jco",  error.plot = "errorbar")+
  theme_classic()+
  labs(title = "Percent Mass Loss of Maple Leaf Litter",
              subtitle = "Plot of Mass Loss by Treatment", x = "Date", y = "Mass Remaining (%)")+
  geom_vline(xintercept = as.numeric(as.Date("2020-03-18")),linetype="solid",  color = "black", size=0.5)+
   scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

ML

Statistics (ANOVA):

#ANOVA 
res.aov2 <- aov(mass_loss_from_OG ~ treatment+harvest , data = ANOVA_final)
Anova(res.aov2, type = "III") # type III sum of squares to address unbalanced sampling

## Anova Table (Type III tests)
## 
## Response: mass_loss_from_OG
##             Sum Sq Df  F value Pr(>F)    
## (Intercept)  85277  1 337.2531 <2e-16 ***
## treatment      433  1   1.7118 0.1993    
## harvest        105  3   0.1390 0.9360    
## Residuals     8850 35                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#no significance observed

Assumptions Check:

# 1. Homogeneity of variances
plot(res.aov2, 1)

leveneTest(mass_loss_from_OG ~ treatment*harvest, data = ANOVA_final)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  7   2.723 0.02466 *
##       32                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 2. Normality
plot(res.aov2, 2)

# Extract the residuals
aov_residuals <- residuals(object = res.aov2)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals )

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.94991, p-value = 0.07526

2. Leaf Litter C:N, %C, %N

Preliminary Organization of Data:

# organizing data - gleaning columns and defining C:N 
cn <- read.csv("CN_harvest.csv")
cn <- cn[c(1,4,7,9,10,12)]
cn$C_N <- cn$perc_C/cn$perc_N
str(cn)

## 'data.frame':    40 obs. of  7 variables:
##  $ Sample_ID: Factor w/ 35 levels "5697","5698",..: 31 32 33 34 35 31 32 33 34 35 ...
##  $ perc_N   : num  1.58 1.48 1.19 1.44 1.27 ...
##  $ perc_C   : num  45.6 44.8 45.2 45.5 46 ...
##  $ harvest  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ treatment: Factor w/ 2 levels "C","SC": 2 2 2 2 2 1 1 1 1 1 ...
##  $ date     : Factor w/ 4 levels "01.14.20","02.26.20",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ C_N      : num  28.9 30.4 38 31.5 36.1 ...

cn$harvest<- as.factor(cn$harvest)
cn$datefinal <- mdy(cn$date)

## boxplot for C:N throughout harvests 
CN <- ggplot(cn, aes(x = harvest, y = C_N, color = treatment))+
  geom_boxplot()+
  theme_classic()+
  labs(title = "C:N of Maple Leaf Litter", x = "Harvest", y = "C:N")+
  scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))
CN

# boxplot for C throughout harvests 
C <- ggplot(cn, aes(x = harvest, y = perc_C, color = treatment))+
  geom_boxplot()+
  theme_classic()+
  labs(title = "%C of Maple Leaf Litter", x = "Harvest", y = "%C")+
  scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))
C

# boxplot for N throughout harvests 
N <-ggplot(cn, aes(x = harvest, y = perc_N, color = treatment))+
  geom_boxplot()+
  theme_classic()+
  labs(title = "%N of Maple Leaf Litter", x = "Harvest", y = "%N")+
  scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))
N

grid.arrange(CN, C, N, nrow = 3)

Finalized Figures:

# C:N plot 
C_N<-ggline(cn, x = "datefinal", y = "C_N", 
       add = c("mean_se", "jitter"),
       color = "treatment", palette = "jco",  error.plot = "errorbar")+
  theme_classic()+
  labs(title = "C:N of Maple Leaf Litter",
               x = "Date", y = "C:N")+
  geom_vline(xintercept = as.numeric(as.Date("2020-03-18")),linetype="solid",  color = "black", size=0.5)+
   scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

C_N

#%C plot 
C_ <- ggline(cn, x = "datefinal", y = "perc_C", 
       add = c("mean_se", "jitter"),
       color = "treatment", palette = "jco",  error.plot = "errorbar")+
  theme_classic()+
  labs(title = "%C of Maple Leaf Litter",
               x = "Date", y = "%C")+
  geom_vline(xintercept = as.numeric(as.Date("2020-03-18")),linetype="solid",  color = "black", size=0.5)+
   scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

C_

#%N plot
N_<-ggline(cn, x = "datefinal", y = "perc_N", 
       add = c("mean_se", "jitter"),
       color = "treatment", palette = "jco",  error.plot = "errorbar")+
  theme_classic()+
  labs(title = "%N of Maple Leaf Litter",
               x = "Date", y = "%N")+
  geom_vline(xintercept = as.numeric(as.Date("2020-03-18")),linetype="solid",  color = "black", size=0.5)+
   scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

N_

ggarrange(C_N,N_,C_, ML)

Statistics (ANOVA):

#preliminary organization
ANOVA_CN <- cn
ANOVA_CN <- ANOVA_CN[c(2,3,4,5,7)]
ANOVA_CN$harvest <- factor(ANOVA_CN$harvest)
#str(ANOVA_CN)
table(ANOVA_CN$treatment, ANOVA_CN$harvest)

##     
##      0 1 2 3
##   C  5 5 5 4
##   SC 5 5 6 5

#ANOVA for C:N
res.aov_CN <- aov(C_N ~ treatment+harvest , data = ANOVA_CN)
Anova(res.aov_CN, type = "III") 

## Anova Table (Type III tests)
## 
## Response: C_N
##             Sum Sq Df  F value  Pr(>F)    
## (Intercept) 9778.1  1 245.8125 < 2e-16 ***
## treatment    161.0  1   4.0462 0.05202 .  
## harvest      314.8  3   2.6379 0.06482 .  
## Residuals   1392.3 35                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# type III sum of squares to address unbalanced sampling
# summary(res.aov_CN)
# treatment has significant effect on C:N

#ANOVA for %C
res.aov_C <- aov(perc_C ~ treatment+harvest , data = ANOVA_CN)
Anova(res.aov_C, type = "III") 

## Anova Table (Type III tests)
## 
## Response: perc_C
##              Sum Sq Df    F value    Pr(>F)    
## (Intercept) 16816.7  1 16122.9554 < 2.2e-16 ***
## treatment       8.4  1     8.0211  0.007617 ** 
## harvest        62.4  3    19.9335 1.039e-07 ***
## Residuals      36.5 35                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# type III sum of squares to address unbalanced sampling
# summary(res.aov_C)
# both treatment and harvest are significant for %C
TukeyHSD(res.aov_C)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = perc_C ~ treatment + harvest, data = ANOVA_CN)
## 
## $treatment
##           diff       lwr        upr     p adj
## SC-C -1.023188 -1.679651 -0.3667256 0.0032115
## 
## $harvest
##           diff       lwr        upr     p adj
## 1-0 -0.5881892 -1.819954  0.6435754 0.5766123
## 2-0 -1.1752680 -2.378712  0.0281764 0.0576811
## 3-0 -3.4153813 -4.680899 -2.1498634 0.0000001
## 2-1 -0.5870788 -1.790523  0.6163656 0.5594240
## 3-1 -2.8271921 -4.092710 -1.5616743 0.0000042
## 3-2 -2.2401133 -3.478083 -1.0021432 0.0001307

#ANOVA for %N
res.aov_N <- aov(perc_N ~ treatment+harvest , data = ANOVA_CN)
Anova(res.aov_N, type = "III") 

## Anova Table (Type III tests)
## 
## Response: perc_N
##              Sum Sq Df  F value  Pr(>F)    
## (Intercept) 14.3870  1 365.6558 < 2e-16 ***
## treatment    0.1012  1   2.5720 0.11776    
## harvest      0.3198  3   2.7097 0.05989 .  
## Residuals    1.3771 35                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# type III sum of squares to address unbalanced sampling
# summary(res.aov_N)
# no significance observed 

Assumptions Check:

# C:N
# 1. Homogeneity of variances
plot(res.aov_CN, 1)

leveneTest(C_N ~ treatment*harvest, data = ANOVA_CN)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  0.8689 0.5412
##       32

# 2. Normality
plot(res.aov_CN, 2)

# Extract the residuals
aov_residuals_CN <- residuals(object = res.aov_CN)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals_CN)

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals_CN
## W = 0.94767, p-value = 0.06307

#%C
# 1. Homogeneity of variances
plot(res.aov_C, 1)

leveneTest(perc_C ~ treatment*harvest, data = ANOVA_CN)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  0.8341 0.5672
##       32

# 2. Normality
plot(res.aov_C, 2)

# Extract the residuals
aov_residuals_C <- residuals(object = res.aov_C)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals_C)

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals_C
## W = 0.95761, p-value = 0.1387

#%N
# 1. Homogeneity of variances
plot(res.aov_N, 1)

leveneTest(perc_N ~ treatment*harvest, data = ANOVA_CN)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  7  0.3779 0.9084
##       32

# 2. Normality
plot(res.aov_N, 2)

# Extract the residuals
aov_residuals_N <- residuals(object = res.aov_N)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals_N)

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals_N
## W = 0.96586, p-value = 0.264

Post Hoc Analysis:

# wut 

3. Respiration (Rs)

Variable Temperature

Preliminary Organization of Data:

# loading in data sets from 2 experimental trials 
sc.t1<- read.csv("07.03_SC_deltaT.csv")
sc.t1$timept <- 1
sc.t1$trt <- "CS"
sc.t1$watercont <- 76.8
sc.t1$bag <- 1

sc.t2<- read.csv("10.03_SC_deltaT.csv")
sc.t2$timept <- 2
sc.t2$trt <- "CS"
sc.t2$watercont <- 77.7
sc.t2$bag <- 2

cc.t1 <- read.csv("07.03_CC_deltaT.csv")
cc.t1$timept <- 1
cc.t1$trt <- "CC"
cc.t1$watercont <- 67.1
cc.t1$bag <- 3

cc.t2 <- read.csv("10.03_CC_deltaT.csv")
cc.t2$timept <- 2
cc.t2$trt <- "CC"
cc.t2$watercont <- 69.6
cc.t2$bag <- 4

# bind all together
all <- rbind(sc.t1, sc.t2, cc.t1, cc.t2)
all$time <- hms(all$HHMMSS)

## Warning in .parse_hms(..., order = "HMS", quiet = quiet): Some strings failed to
## parse, or all strings are NAs

# normalize for water content differences, divide Rs by wc% 
all$Rs_Final <- all$Rs/all$watercont
all$bag <- as.factor(all$bag)
all$trt<-as.factor(all$trt) 

Finalized Figures:

# preliminary figure for all trials, CS v CC with variable temperature
prelim_rs <- ggplot(all, aes(x = Tsoil_C, y = Rs_Final, color = factor(trt)))+
  geom_point()+
  xlab("Leaf Litter Temperature (°C)")+
  ylab(bquote('Respiration ('*mu~'g' ~CO[2]~ g^-1~s^-1~'%WC'^-1*')'))+
  labs(title="Leaf Litter Respiration by Treatment Over Variable Temperature")+  
  geom_smooth(method=lm)+
  theme_classic()+
  scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))
prelim_rs

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

# preliminary graph of Rs by treatment type 
 prelim_rs_boxo <- ggplot(all, aes(x = trt, y = Rs_Final, color = trt))+
  geom_boxplot()+
  xlab("Treatment")+
  ylab(bquote('Respiration ('*mu~'g' ~CO[2]~ g^-1~s^-1~'%WC'^-1*')'))+
  labs(title="Boxplot of Leaf Litter Respiration by Treatment Over Variable Temperature")+
  scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))+
  theme_classic()
prelim_rs_boxo

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

Statistics (ANCOVA)

# without random effects being considered 
mod2_new <- lm(Rs_Final~Tsoil_C*trt, data=all )
Anova(mod2_new, type="III")

## Anova Table (Type III tests)
## 
## Response: Rs_Final
##                 Sum Sq   Df   F value Pr(>F)    
## (Intercept) 0.00007096    1  900.9448 <2e-16 ***
## Tsoil_C     0.00054879    1 6967.3752 <2e-16 ***
## trt         0.00000005    1    0.6805 0.4095    
## Tsoil_C:trt 0.00001066    1  135.3819 <2e-16 ***
## Residuals   0.00012043 1529                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.lm(mod2_new)

## 
## Call:
## lm(formula = Rs_Final ~ Tsoil_C * trt, data = all)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0013627 -0.0001769 -0.0000231  0.0001507  0.0013300 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.277e-04  2.758e-05  30.016   <2e-16 ***
## Tsoil_C        1.452e-04  1.739e-06  83.471   <2e-16 ***
## trtCS         -2.875e-05  3.485e-05  -0.825     0.41    
## Tsoil_C:trtCS -2.647e-05  2.275e-06 -11.635   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0002807 on 1529 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.9101, Adjusted R-squared:  0.9099 
## F-statistic:  5161 on 3 and 1529 DF,  p-value: < 2.2e-16

# summary(mod2_new)
# shows significant interaction between temperature of leaf litter and treatment type 

# to include random effect of bag
# https://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
# http://home.cc.umanitoba.ca/~krussll/stats/mixed-effects.html
# https://stats.stackexchange.com/questions/211642/mixed-effect-model-ancova-with-lmer-in-r
# this is usign the lmer function - not sure if it crosses over to the aov function?
# aov gets tricky with random efects, better to move to using lme4 funciton

#mod2 <- lmer(Rs_Final~Tsoil_C*trt + (1|bag), data=all )
#Anova(mod2)
# analysis of deviance table
#summary(mod2)
#ranef(mod2) # random effects of each bag 

Variable Moisture

Preliminary Organization of Data:

deltaM <- read.csv("deltaM_compiled.csv")
# change wc% to categorical variable 
deltaM$wc_percent=as.factor(deltaM$wc_percent)
# str(deltaM)
# create subset particular to each treatment (SC and CC)
deltaM_SC <- subset(deltaM, treatment == "SC")
deltaM_SC_troubleshoot <- read.csv("deltaM_SC_troubleshoot.csv")
#str(deltaM_SC)
#str(deltaM_SC_troubleshoot)
deltaM_SC_troubleshoot$wc_percent=as.factor(deltaM_SC_troubleshoot$wc_percent)
deltaM_CC <- subset(deltaM, treatment == "CC")
deltaM_SC <- subset(deltaM, treatment == "SC")

#scaling for differences in temperatures (divide Rs by Tsoil_C to create Rs_final)
deltaM$Rs_final <- (deltaM$Rs / deltaM$Tsoil_C)
#str(deltaM)
deltaM$WC <- as.factor(deltaM$WC)

Finalized Figures: (still need to figure out continuous x-axis)

#str(deltaM)
deltaM$WC <- as.numeric(deltaM$WC)

deltaM_summary_plot <- ggplot(deltaM, aes(x = WC, y = Rs_final, fill = treatment, group = WC))+
  geom_boxplot()+
  xlab("Leaf Litter Water Content (%)")+
  ylab(bquote('Respiration ('*mu~'g' ~CO[2]~ g^-1~s^-1~'°C'^-1*')'))+
  labs(title="Boxplot of Leaf Litter Respiration by Treatment Over Variable Moisture")+
  theme_classic()+
  scale_fill_discrete(name = "Treatment", labels = c("CC", "CS"))
  
deltaM_summary_plot

# linear plot 
wc<- read.csv("deltaM_compiled.csv")
# normalize for differences in tempertaure of leaf litter 
wc$Rs_Final <- wc$Rs/wc$Tsoil_C
wc$trt <- wc$treatment

# graph Rs_Final by wc litter for cc/sc
prelim_rs_wc <- ggplot(wc, aes(x = wc_percent, y = Rs_Final, color = trt))+
  geom_point()+
  labs(title="Leaf Litter Respiration by Treatment Over Variable Moisture")+  
  geom_smooth(method=lm)+
  xlab("Leaf Litter Water Content (%)")+
  ylab(bquote('Respiration ('*mu~'g' ~CO[2]~ g^-1~s^-1~'°C'^-1*')'))+
  theme_classic()+
  scale_color_discrete(name = "Treatment", labels = c("CC", "CS"))
  
prelim_rs_wc

Statistics (ANCOVA):

# ANCOVA with Rs normalized for leaf litter temperature differences 
# no random effects being considered - unecessary
#mod1_wc <- aov(Rs_Final~wc_percent*trt, data=wc)
#Anova(mod1_wc, type="III") 
#summary(mod1_wc)
# significant effect of water content, not significant for treatment
# no significant interaction - therefore get rid of interaction term in model 

mod2_wc <- aov(Rs_Final~ wc_percent + trt, data=wc)
Anova(mod2_wc, type="III") 

## Anova Table (Type III tests)
## 
## Response: Rs_Final
##                 Sum Sq Df F value    Pr(>F)    
## (Intercept) 3.0100e-07  1  0.2188    0.6427    
## wc_percent  2.8253e-05  1 20.5066 5.991e-05 ***
## trt         1.5810e-06  1  1.1479    0.2909    
## Residuals   5.0977e-05 37                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#summary.lm(mod2_wc)
# wc still has a signifcant effect, treatment does not

# post hoc 

#postHocs<-glht(mod2_wc, linfct = mcp(trt = "Tukey"))
#summary(postHocs)
#confint(postHocs) 

# post hoc notes 
# https://stats.stackexchange.com/questions/51780/how-to-perform-an-ancova-in-r
#posth=glht(model.1, linfct=mcp(factorvariable="Tukey"))  ##gives the post-hoc Tukey analysis
#summary(posth) ##shows the output in a nice format.

4. Physical Conditions

Temperature

Preliminary Organization of Data:

feb_CC <-read.csv("feb_cc_compiled.csv", stringsAsFactors = FALSE)
# merging date and time columns 
feb_CC$datetime <- paste(feb_CC$date, feb_CC$time)
feb_CC$datetime2 <- mdy_hms(feb_CC$datetime)
#plot(feb_CC$datetime2,feb_CC$temp_C)

feb_SC <-read.csv("feb_sc_compiled.csv", stringsAsFactors = FALSE)
# merging date and time columns 
feb_SC$datetime <- paste(feb_SC$date, feb_SC$time)
feb_SC$datetime2 <- mdy_hms(feb_SC$datetime)
#plot(feb_SC$datetime2,feb_SC$temp_C)

jun_CC <-read.csv("jun_cc_compiled.csv", stringsAsFactors = FALSE)
# merging date and time columns 
jun_CC$datetime <- paste(jun_CC$date, jun_CC$time)
jun_CC$datetime2 <- mdy_hms(jun_CC$datetime)
#plot(jun_CC$datetime2,jun_CC$temp_C)

jun_SC <-read.csv("jun_SC_compiled.csv", stringsAsFactors = FALSE)
# merging date and time columns 
jun_SC$datetime <- paste(jun_SC$date, jun_SC$time)
jun_SC$datetime2 <- mdy_hms(jun_SC$datetime)
#plot(jun_SC$datetime2,jun_SC$temp_C)

mcgowan_amb <-read.csv("mcgowan_amb.csv", stringsAsFactors = FALSE)
# merging date and time columns 
mcgowan_amb$datetime <- paste(mcgowan_amb$date, mcgowan_amb$time)
mcgowan_amb$datetime2 <- mdy_hms(mcgowan_amb$datetime)
#plot(mcgowan_amb$datetime2,mcgowan_amb$temp_C)

gamefarm_amb <-read.csv("game_farm_road.csv", stringsAsFactors = FALSE)
# merging date and time columns 
gamefarm_amb$datetime <- paste(gamefarm_amb$date, gamefarm_amb$time)
gamefarm_amb$datetime2 <- mdy_hms(gamefarm_amb$datetime)

## Warning: 1737 failed to parse.

# plot(gamefarm_amb$datetime2,gamefarm_amb$temp_C)

#merging CC datafarmes 
feb_CC$temp_C_feb<-feb_CC$temp_C
jun_CC$temp_C_jun<-jun_CC$temp_C
CC=cbind(feb_CC,jun_CC)
CC <- CC[, !duplicated(colnames(CC))]
CC_final<- CC[,c(6,7,8)]

#merging CS datafarmes 
feb_SC$temp_C_feb_SC<-feb_SC$temp_C
jun_SC$temp_C_jun_SC<-jun_SC$temp_C
SC<-merge(feb_SC, jun_SC, by="datetime2", all = T)
SC_final<- SC[,c(1,7,13)]

both_sites<-merge(SC_final, CC_final, by="datetime2", all = T)

#merging MCG dataframes
mcgowan_amb$temp_C_mcg <- mcgowan_amb$temp_C
mcgowan_amb$precip_mcg <- mcgowan_amb$precip
mcgowan_amb_final<- mcgowan_amb[,c(7,8,9)]

#merging GFR dataframes
gamefarm_amb$temp_C_GFR <- gamefarm_amb$temp_C
gamefarm_amb_final <- gamefarm_amb[,c(11,12)]
allsites <- merge(both_sites, mcgowan_amb_final, by = "datetime2", all = T)
final_all_sites <- merge (allsites,gamefarm_amb_final, by = "datetime2", all = T)

#correct for beginning time record 
complete <- subset(final_all_sites, datetime2 > "2020-01-14 09:00:00")
complete$mean_CC = (complete$temp_C_feb + complete$temp_C_jun) / 2
complete$mean_SC <- rowMeans(complete[c('temp_C_feb_SC', 'temp_C_jun_SC')], na.rm=TRUE)

Finalized Figures:

reduce <- complete[c(-2,-3,-4,-5,-7)]
test <- melt(reduce, id.vars = c("datetime2"))

ggplot(test, aes(x = datetime2, y = value, group = variable, color = variable))+
  geom_line()+
  facet_wrap( ~ variable, ncol=2, labeller = labeller(variable = 
    c("temp_C_mcg" = "McGowan Weather Station",
      "temp_C_GFR" = "Game Farm Weather Station",
      "mean_CC" = "Current Conditions (CC)", 
      "mean_SC" = "Constant Snowcover (CS)")))+
  labs(title="Hourly Temperature Record of Permanant Weather Stations and Field Loggers",
       x ="Date", y = "Temperature (°C)")+ 
   theme_classic()+
  theme(legend.position = "none")+
  geom_vline(aes(xintercept = as.integer(as.POSIXct("2020-03-18 12:00:00"))), linetype = "dashed", col = "black", size = 0.5)+
  theme(plot.title = element_text(size=11))

## Warning: Removed 267 rows containing missing values (geom_path).

Temperature Analysis (via max, min, and mean): Looking at temperature range for all stations / loggers

#McGowan Woods range 
range(reduce$temp_C_mcg, na.rm = T)

## [1] -22.3  19.5

#Game Farm Road range 
range(reduce$temp_C_GFR, na.rm = T)

## [1] -23.52  20.44

#Current Conditions Loger 
range(reduce$mean_CC, na.rm = T)

## [1] -8.4545 28.7870

#Snow Cover Logger
range(reduce$mean_SC, na.rm = T)

## [1] -5.125 26.781

Field Moisture

Preliminary Organization of Data:

field_wc <- read.csv("water_content_field.csv")
#str(field_wc)
field_wc$day_since_deploy <- as.numeric(field_wc$day_since_deploy)
field_wc$ML_harvest=as.factor(field_wc$ML_harvest)
field_wc$datefinal <- mdy(field_wc$date)

#create preliminary boxplot looking at spread of water content by treatment
wc_plot <- ggplot(field_wc, aes(x = factor(day_since_deploy), y = perc_wc, fill = factor(treatment)))+
  geom_boxplot()+
  ggtitle("Leaf Litter Percent Water Content by Treatment") +
  scale_y_continuous(name = "Water Content (%)") +
  scale_x_discrete(name = "Days in Field")+
  scale_fill_discrete(name = "Treatment", labels = c("CC", "CS"))+
  theme_classic()

wc_plot

# subset data by each harvest in order to isolate comparison of water content between treatments 
wc_1 <- subset(field_wc, ML_harvest==1)
wc_2 <- subset(field_wc, ML_harvest==2)
wc_3 <- subset(field_wc, ML_harvest==3)

Finalized Figures:

#harvest 1
h1 <- ggplot(wc_1, aes(x = ML_harvest, y = perc_wc, fill = factor(treatment)))+
  geom_boxplot()+
  scale_y_continuous(name = "Water Content (%)") +
  scale_x_discrete(name = "Harvest Number")+
  scale_fill_discrete(name = "Treatment", labels = c("CC", "CS"))+
  theme_classic()
h1

#harvest 2
h2 <- ggplot(wc_2, aes(x = ML_harvest, y = perc_wc, fill = factor(treatment)))+
  geom_boxplot()+
  scale_y_continuous(name = "Water Content (%)") +
  scale_x_discrete(name = "Harvest Number")+
  scale_fill_discrete(name = "Treatment", labels = c("CC", "CS"))+
  theme_classic()
h2

#harvest 3
h3 <- ggplot(wc_3, aes(x = ML_harvest, y = perc_wc, fill = factor(treatment)))+
  geom_boxplot()+
  scale_y_continuous(name = "Water Content (%)") +
  scale_x_discrete(name = "Harvest Number")+
  scale_fill_discrete(name = "Treatment", labels = c("CC", "CS"))+
  theme_classic()

h3

wc_compiled <- ggarrange(h1,h2,h3)
wc_compiled

annotate_figure(wc_compiled, top = text_grob("Field Conditions: Leaf Litter Water Content by Harvest and Treatment", size =14))

Statistical Tests:

# see if statistically significant differences between treatments across harvests
aov_wc <- aov(perc_wc ~ treatment*ML_harvest , data = field_wc)
Anova(aov_wc, type = "III")

## Anova Table (Type III tests)
## 
## Response: perc_wc
##                      Sum Sq Df   F value  Pr(>F)    
## (Intercept)           34035  1 3461.8225 < 2e-16 ***
## treatment                31  1    3.1661 0.08785 .  
## ML_harvest            17195  2  874.4846 < 2e-16 ***
## treatment:ML_harvest   9673  2  491.9142 < 2e-16 ***
## Residuals               236 24                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov_wc)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = perc_wc ~ treatment * ML_harvest, data = field_wc)
## 
## $treatment
##          diff      lwr      upr p adj
## SC-C 26.24234 23.87405 28.61064     0
## 
## $ML_harvest
##          diff       lwr       upr p adj
## 2-1 -34.36270 -37.78401 -30.94139     0
## 3-1 -76.87437 -80.47215 -73.27659     0
## 3-2 -42.51167 -46.03114 -38.99221     0
## 
## $`treatment:ML_harvest`
##                  diff        lwr        upr     p adj
## SC:1-C:1    3.5285914  -2.602963   9.660146 0.4966479
## C:2-C:1   -73.0161740 -79.147729 -66.884619 0.0000000
## SC:2-C:1    3.2699370  -2.600580   9.140454 0.5312115
## C:3-C:1   -73.2957297 -79.799225 -66.792234 0.0000000
## SC:3-C:1  -73.9373165 -80.068871 -67.805762 0.0000000
## C:2-SC:1  -76.5447654 -82.676320 -70.413211 0.0000000
## SC:2-SC:1  -0.2586544  -6.129171   5.611862 0.9999925
## C:3-SC:1  -76.8243210 -83.327817 -70.320825 0.0000000
## SC:3-SC:1 -77.4659079 -83.597463 -71.334353 0.0000000
## SC:2-C:2   76.2861110  70.415594  82.156628 0.0000000
## C:3-C:2    -0.2795557  -6.783051   6.223940 0.9999933
## SC:3-C:2   -0.9211425  -7.052697   5.210412 0.9969556
## C:3-SC:2  -76.5656666 -82.823658 -70.307675 0.0000000
## SC:3-SC:2 -77.2072535 -83.077770 -71.336737 0.0000000
## SC:3-C:3   -0.6415868  -7.145083   5.861909 0.9995970

# assumptions - all good
# 1. Homogeneity of variances
plot(aov_wc, 1)

leveneTest(perc_wc ~ treatment*ML_harvest, data = field_wc)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  1.6349  0.189
##       24

# 2. Normality
plot(aov_wc, 2)

# Extract the residuals
aov_residuals_wc <- residuals(object = aov_wc)
# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals_wc)

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals_wc
## W = 0.97796, p-value = 0.769

# summary statistics 
require("dplyr")

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union

## The following object is masked from 'package:gridExtra':
## 
##     combine

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

wc_stats<-group_by(field_wc, treatment, ML_harvest) %>%
  summarise(
    count = n(),
    mean = mean(perc_wc, na.rm = TRUE),
    sd = sd(perc_wc, na.rm = TRUE)
  )
formattable(wc_stats, col.names = c("Treatment", "Harvest Number", "Count", "Mean", "SD"), align = c("l","c", "c", "c", "c"))

Treatment	Harvest Number	Count	Mean	SD
C	1	5	82.504607	5.334972
C	2	5	9.488433	1.587738
C	3	4	9.208878	1.775822
SC	1	5	86.033199	3.953438
SC	2	6	85.774544	2.424147
SC	3	5	8.567291	1.632781