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