Step 1: Clean Data

setwd('C:\\Users\\Yang\\Desktop\\Research\\Work\\CH4 data')
GHG<-read.csv('CH4_Data.csv',head=TRUE)
GHG$F <- NA
GHG$F[which(GHG$Treatment== "Control"|GHG$Treatment == "Rainout")] <- 0
GHG$F[which(GHG$Treatment == "Fertilized"|GHG$Treatment == "Rainout+Fertilized")] <- 1
GHG$R <- NA
GHG$R[which(GHG$Treatment == "Control"|GHG$Treatment == "Fertilized")] <- 0
GHG$R[which(GHG$Treatment == "Rainout"|GHG$Treatment == "Rainout+Fertilized")] <- 1

GHG$Block <- NA
GHG$Block[which(GHG$Plot== 1|GHG$Plot == 2|GHG$Plot == 3|GHG$Plot == 4)] <- 1
GHG$Block[which(GHG$Plot== 5|GHG$Plot == 6|GHG$Plot == 7|GHG$Plot == 8)] <- 2
GHG$Block[which(GHG$Plot== 9|GHG$Plot == 10|GHG$Plot == 11|GHG$Plot == 12)] <- 3
GHG$Block[which(GHG$Plot== 13|GHG$Plot == 14|GHG$Plot == 15|GHG$Plot == 16)] <- 4

GHG[,2]<-as.factor(GHG[,2])
GHG[,3]<-as.factor(GHG[,3])
GHG[,4]<-as.factor(GHG[,4])
GHG[,7]<-as.factor(GHG[,7])
GHG[,8]<-as.factor(GHG[,8])
GHG[,9]<-as.factor(GHG[,9])
head(GHG)
##        Date Plot  Treatment Round       CO2 CH4 F R Block
## 1 11/4/2015    1    Control     1        NA  NA 0 0     1
## 2 11/4/2015    1    Control     1        NA  NA 0 0     1
## 3 11/4/2015    1    Control     1        NA  NA 0 0     1
## 4 11/4/2015    2 Fertilized     1 13.290000  NA 1 0     1
## 5 11/4/2015    2 Fertilized     1 10.950000  NA 1 0     1
## 6 11/4/2015    2 Fertilized     1  8.646667  NA 1 0     1

step 2: transform CO2 and CH4 if they are not normal

qqnorm(GHG$CO2)

GHG$CO2_log <- log(GHG$CO2)
qqnorm(GHG$CO2_log)

qqnorm(GHG$CH4)

GHG$CH4_log <- log(-GHG$CH4)
## Warning in log(-GHG$CH4): NaNs produced
qqnorm(GHG$CH4_log)

Step 3: ANOVA for CO2 and CH4

library(lme4)
## Loading required package: Matrix
CO2.aov<-aov(CO2_log~F*R*Round+Block, data=GHG)
summary(CO2.aov)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## F             1   0.47    0.47   0.922  0.339    
## R             1   0.12    0.12   0.234  0.630    
## Round         2  68.05   34.03  66.262 <2e-16 ***
## Block         3   1.78    0.59   1.153  0.331    
## F:R           1   0.17    0.17   0.335  0.564    
## F:Round       2   0.08    0.04   0.081  0.922    
## R:Round       2   0.21    0.11   0.206  0.814    
## F:R:Round     2   0.64    0.32   0.626  0.537    
## Residuals   101  51.86    0.51                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 44 observations deleted due to missingness
TukeyHSD(CO2.aov, 'Round')
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = CO2_log ~ F * R * Round + Block, data = GHG)
## 
## $Round
##           diff         lwr        upr    p adj
## 2-1 -1.8020743 -2.18098204 -1.4231666 0.000000
## 3-1 -1.3265748 -1.74640698 -0.9067425 0.000000
## 3-2  0.4754996  0.08960812  0.8613911 0.011541
CH4.aov<-aov(CH4_log~F*R*Round+Block, data=GHG)
summary(CH4.aov)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## F            1   5.09   5.089   7.050 0.00946 **
## R            1   0.12   0.123   0.170 0.68095   
## Round        2   1.23   0.616   0.853 0.42963   
## Block        3  10.09   3.365   4.662 0.00459 **
## F:R          1   3.02   3.017   4.179 0.04402 * 
## F:Round      2   0.65   0.325   0.450 0.63919   
## R:Round      2   1.54   0.772   1.069 0.34790   
## F:R:Round    2   1.29   0.646   0.895 0.41230   
## Residuals   85  61.35   0.722                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 60 observations deleted due to missingness
TukeyHSD(CH4.aov, 'F')
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = CH4_log ~ F * R * Round + Block, data = GHG)
## 
## $F
##           diff        lwr        upr     p adj
## 1-0 -0.4556639 -0.7968704 -0.1144573 0.0094612
CH4.aov2 <- aov(CH4_log~F*R+Block, data=GHG)
TukeyHSD(CH4.aov2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = CH4_log ~ F * R + Block, data = GHG)
## 
## $F
##           diff        lwr        upr    p adj
## 1-0 -0.4556639 -0.7930024 -0.1183254 0.008652
## 
## $R
##           diff        lwr       upr     p adj
## 1-0 0.07012015 -0.2639629 0.4042032 0.6777876
## 
## $Block
##           diff        lwr         upr     p adj
## 2-1  0.5284069 -0.1412233  1.19803702 0.1725217
## 3-1  0.3083208 -0.3288983  0.94553988 0.5868760
## 4-1 -0.3180162 -0.9936436  0.35761125 0.6086212
## 3-2 -0.2200861 -0.8073697  0.36719758 0.7609398
## 4-2 -0.8464230 -1.4751729 -0.21767321 0.0036643
## 4-3 -0.6263370 -1.2204497 -0.03222423 0.0347189
## 
## $`F:R`
##               diff        lwr        upr     p adj
## 1:0-0:0 -0.1152233 -0.7372804  0.5068337 0.9623445
## 0:1-0:0  0.3689031 -0.2140248  0.9518311 0.3530533
## 1:1-0:0 -0.4412554 -1.0716781  0.1891673 0.2654425
## 0:1-1:0  0.4841265 -0.1427040  1.1109570 0.1879026
## 1:1-1:0 -0.3260321 -0.9972580  0.3451938 0.5837787
## 1:1-0:1 -0.8101586 -1.4452918 -0.1750253 0.0065787

Step 4: Calculate Mean and Se for CO2 and CH4 for each treatment and thier interaction

aggregate(CO2 ~ F, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   F     CO2.M    CO2.SE
## 1 0 3.8505552 0.6280660
## 2 1 3.5102796 0.5978107
aggregate(CO2 ~ R, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   R     CO2.M    CO2.SE
## 1 0 3.5249277 0.5857741
## 2 1 3.8533024 0.6438373
aggregate(CO2 ~ Round, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   Round      CO2.M     CO2.SE
## 1     1 8.86702858 1.01399497
## 2     2 1.26620000 0.08881436
## 3     3 1.96312500 0.16534021
aggregate(CO2 ~ F*R*Round, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##    F R Round       CO2.M      CO2.SE
## 1  0 0     1  7.55044341  2.07014459
## 2  1 0     1 10.13339821  2.17255581
## 3  0 1     1  9.42331389  1.79339441
## 4  1 1     1  7.82141205  2.50199783
## 5  0 0     2  1.16437500  0.08870348
## 6  1 0     2  1.66000000  0.28641356
## 7  0 1     2  1.09375000  0.14093993
## 8  1 1     2  1.17800000  0.15545387
## 9  0 0     3  1.90250000  0.21138785
## 10 1 0     3  1.79875000  0.30329404
## 11 0 1     3  1.86875000  0.33181017
## 12 1 1     3  2.28250000  0.46611521
aggregate(CH4 ~ F, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   F      CH4.M     CH4.SE
## 1 0 -1.7831084  0.2272607
## 2 1 -0.9147837  0.1994959
aggregate(CH4 ~ R, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   R      CH4.M     CH4.SE
## 1 0 -1.2272688  0.1991213
## 2 1 -1.5041828  0.2430105
aggregate(CH4 ~ Round, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   Round      CH4.M     CH4.SE
## 1     1 -1.4038549  0.1541584
## 2     2 -1.6497391  0.3201052
## 3     3 -0.9412202  0.2487723
aggregate(CH4 ~ F*R, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##   F R      CH4.M     CH4.SE
## 1 0 0 -1.3422970  0.2216635
## 2 1 0 -1.0916998  0.3493738
## 3 0 1 -2.2847214  0.3997153
## 4 1 1 -0.7496621  0.2078931
aggregate(CH4 ~ F*R*Round, GHG, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))
##    F R Round        CH4.M       CH4.SE
## 1  0 0     1 -1.149705636  0.318940281
## 2  1 0     1 -1.531588456  0.331468421
## 3  0 1     1 -1.884686923  0.300440432
## 4  1 1     1 -1.049438576  0.243538340
## 5  0 0     2 -1.257681699  0.388975593
## 6  1 0     2 -1.684969116  0.669592022
## 7  0 1     2 -3.218642317  0.944596450
## 8  1 1     2 -0.672291207  0.389573180
## 9  0 0     3 -1.675913857  0.396965402
## 10 1 0     3  0.007376999  0.622309973
## 11 0 1     3 -1.543297134  0.355614989
## 12 1 1     3 -0.553046739  0.409417044

Conclusion 1) For CO2, only Round has effect on CO2 flux, compared to round 1 (8.86), both round 2 (1.27) and round 3 (1.96) were significantly less than round 1. Also, round 3 was significantly higher than round 2; however, it is unusuall that round 1’s value (Nov, 2015) was much higher than round 2 (Jan, 2016) and round 3 (May, 2016), need to pay attention to the accuracy of measurements

2) For CH4, Fertilizaiton increased CH4 flux from -1.78 to -0.91; Rainout had no effect; But Fertilization had interaction with Rainout that only with Rainout treatments, Fertilizaiton increased CH4 flux from -2.28 to -0.75