Models for Soil Nutrients
Jonathan Spiess
January 31, 2021
SoilTukey21 <-read.csv("D:/R/data/SRMSoils 2021Tukey.csv", head=TRUE, stringsAsFactors = FALSE)
SoilTukey21$Variable <- factor(SoilTukey21$Variable)
SoilTukey21$Contrast <- factor(SoilTukey21$Contrast)
SoilTukey21$Years <- factor(SoilTukey21$Years)
print(levels(SoilTukey21$Variable))## [1] "Ammonium" "Calcium" "Decomposition" "Magnesium"
## [5] "Nitrate" "Phosphorus" "Potassium" "Soil Moisture"
## [9] "Total Abundance" "Total C:N" "Total Carbon" "Total Nitrogen"SoilTukey21$Variable = factor(SoilTukey21$Variable,
levels(SoilTukey21$Variable)[c(2,7,4,6,11,12,10,9,3,8,5,1)])
ggplot(SoilTukey21, aes(x=Variable, fill= Contrast, y=Estimate),color=Contrast)+
geom_hline(yintercept=0, color="black")+
geom_linerange(aes(ymin=Lower,ymax=Upper, color=Contrast), size=1.5, position=position_dodge(0.5))+
geom_point(size=3, pch=21, position=position_dodge(0.5), color="gray38", stroke=1)+
scale_colour_viridis(discrete = TRUE)+
scale_fill_viridis(discrete = TRUE, name="")+
xlab ("")+ylab ("Standardized Tukey's Estimate with 95% CI")+
# ggtitle(label = " Recent Burn Low No Difference Recent Burn Higher")+
guides(color=FALSE)+
coord_flip()+
theme_bw()+
theme(legend.position = c(0.89, 0.575),
axis.title=element_text(size=8, face="bold"),
axis.text.y =element_text(size=8, face="bold"),
axis.text.x = element_text(size=8, face="bold"),
legend.text=element_text(size=5),
legend.title = element_text(size=5),
plot.title = element_text(size=7, face="bold"),
#plot.margin = margin(6,2,2,2),
panel.grid.minor = element_blank())
Bnorm <- fitdist(SoilM$Total_Biomass, "norm")
plot(Bnorm) #not terrible, qq is a little wonky at the ends
Blognorm <- fitdist(log(SoilM$Total_Biomass), "norm")
plot(Blognorm) #use log distribution
Bgamma <- fitdist(SoilM$Total_Biomass, "gamma")
plot(Bgamma) #does not look too different than log normal
#untransformed
ggplot(SoilM, aes(x=Location, y=Total_Biomass, color=TSF))+
geom_boxplot(size=1)+
labs(y="Total Abundance", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1, scales = "free")
#with log
ggplot(SoilM, aes(x=Location, y=log(Total_Biomass), color=TSF))+
geom_boxplot(size=1)+
labs(y="log(Total Abundance)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
Total abundance stats
TAnull <-lmer(log(Total_Biomass) ~ 1 + (1|Location/Year), data=SoilM, REML = FALSE) #singular fit issues## boundary (singular) fit: see ?isSingularTAT <- lmer(log(Total_Biomass) ~ TSF + (1|Location/Year), data=SoilM, REML = FALSE) #singular fit## boundary (singular) fit: see ?isSingularTAE <- lmer(log(Total_Biomass) ~ ESD + (1|Location/Year), data=SoilM, REML = FALSE)## boundary (singular) fit: see ?isSingularTATE <- lmer(log(Total_Biomass) ~ TSF + ESD + (1|Location/Year), data=SoilM, REML = FALSE)
anova(TAnull, TAT, TAE, TATE)## Data: SoilM
## Models:
## TAnull: log(Total_Biomass) ~ 1 + (1 | Location/Year)
## TAT: log(Total_Biomass) ~ TSF + (1 | Location/Year)
## TAE: log(Total_Biomass) ~ ESD + (1 | Location/Year)
## TATE: log(Total_Biomass) ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## TAnull 4 105.577 118.607 -48.788 97.577
## TAT 8 100.213 126.273 -42.106 84.213 13.364 4 0.009626 **
## TAE 8 70.079 96.139 -27.040 54.079 30.133 0 < 2.2e-16 ***
## TATE 12 67.777 106.867 -21.888 43.777 10.303 4 0.035623 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TA.tscale <- lmer(scale(log(Total_Biomass)) ~ TSF + ESD + (1|Location/Year), data=SoilM, REML = FALSE)
Mult_TA.TSF <- glht(TA.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_TA.TSF) #NYB lower than 1-2yr##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Total_Biomass)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilM, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.21440 0.13545 1.583 0.5002
## 2yr - 3yr - 1yr < == 0 -0.08088 0.13535 -0.598 0.9746
## 3yr - 4yr - 1yr < == 0 -0.02133 0.17406 -0.123 0.9999
## Unburned - 1yr < == 0 -0.31391 0.17252 -1.820 0.3544
## 2yr - 3yr - 1yr - 2yr == 0 -0.29528 0.13497 -2.188 0.1789
## 3yr - 4yr - 1yr - 2yr == 0 -0.23572 0.17607 -1.339 0.6596
## Unburned - 1yr - 2yr == 0 -0.52831 0.17338 -3.047 0.0187 *
## 3yr - 4yr - 2yr - 3yr == 0 0.05955 0.17408 0.342 0.9969
## Unburned - 2yr - 3yr == 0 -0.23303 0.17356 -1.343 0.6572
## Unburned - 3yr - 4yr == 0 -0.29258 0.21842 -1.340 0.6591
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TA.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Total_Biomass)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilM, REML = FALSE)
##
## Quantile = 2.7137
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.21440 -0.15317 0.58196
## 2yr - 3yr - 1yr < == 0 -0.08088 -0.44817 0.28642
## 3yr - 4yr - 1yr < == 0 -0.02133 -0.49367 0.45102
## Unburned - 1yr < == 0 -0.31391 -0.78209 0.15427
## 2yr - 3yr - 1yr - 2yr == 0 -0.29528 -0.66154 0.07099
## 3yr - 4yr - 1yr - 2yr == 0 -0.23572 -0.71352 0.24208
## Unburned - 1yr - 2yr == 0 -0.52831 -0.99882 -0.05780
## 3yr - 4yr - 2yr - 3yr == 0 0.05955 -0.41284 0.53194
## Unburned - 2yr - 3yr == 0 -0.23303 -0.70403 0.23797
## Unburned - 3yr - 4yr == 0 -0.29258 -0.88533 0.30016Mult_TA.ESD <- glht(TA.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_TA.ESD) #sandy lower than all except clayey##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Total_Biomass)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilM, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.1552 0.1984 0.782 0.9278
## Saline Lowland - Clayey == 0 0.3455 0.1547 2.233 0.1492
## Sandy - Clayey == 0 -0.3609 0.1410 -2.560 0.0677 .
## Thin Claypan - Clayey == 0 1.1060 0.3791 2.918 0.0248 *
## Saline Lowland - Loamy == 0 0.1903 0.1999 0.952 0.8624
## Sandy - Loamy == 0 -0.5161 0.1778 -2.902 0.0264 *
## Thin Claypan - Loamy == 0 0.9508 0.3954 2.405 0.0997 .
## Sandy - Saline Lowland == 0 -0.7064 0.1239 -5.703 <0.001 ***
## Thin Claypan - Saline Lowland == 0 0.7605 0.3580 2.125 0.1880
## Thin Claypan - Sandy == 0 1.4669 0.3616 4.056 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TA.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Total_Biomass)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilM, REML = FALSE)
##
## Quantile = 2.6751
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.15517 -0.37550 0.68583
## Saline Lowland - Clayey == 0 0.34549 -0.06848 0.75945
## Sandy - Clayey == 0 -0.36090 -0.73796 0.01616
## Thin Claypan - Clayey == 0 1.10601 0.09189 2.12013
## Saline Lowland - Loamy == 0 0.19032 -0.34446 0.72510
## Sandy - Loamy == 0 -0.51607 -0.99171 -0.04042
## Thin Claypan - Loamy == 0 0.95085 -0.10684 2.00854
## Sandy - Saline Lowland == 0 -0.70639 -1.03774 -0.37503
## Thin Claypan - Saline Lowland == 0 0.76053 -0.19708 1.71813
## Thin Claypan - Sandy == 0 1.46691 0.49949 2.43434Variables with single sampling event
P
#P Overall
Pnorm <- fitdist(SoilJuly181920$P_ppm, "norm")
plot(Pnorm) #qq, empirical, and pp not good
Plognorm <- fitdist(log(SoilJuly181920$P_ppm), "norm")
plot(Plognorm) #slightly better
Pgamma <- fitdist(SoilJuly181920$P_ppm, "gamma")
plot(Pgamma) #qq wonky on right end
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=P_ppm, color=TSF))+
geom_boxplot(size=1)+
labs(y="P_ppm", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
#with log
ggplot(SoilJuly181920, aes(x=Location, y=log(P_ppm), color=TSF))+
geom_boxplot(size=1)+
labs(y="log(P_ppm)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
Pnull <-lmer(log(P_ppm) ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
PT <- lmer(log(P_ppm) ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
PE <- lmer(log(P_ppm) ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
PTE <- lmer(log(P_ppm) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
anova(Pnull, PT, PE, PTE)## Data: SoilJuly181920
## Models:
## Pnull: log(P_ppm) ~ 1 + (1 | Location/Year)
## PT: log(P_ppm) ~ TSF + (1 | Location/Year)
## PE: log(P_ppm) ~ ESD + (1 | Location/Year)
## PTE: log(P_ppm) ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Pnull 4 455.36 470.01 -223.68 447.36
## PT 8 457.04 486.35 -220.52 441.04 6.3142 4 0.1769
## PE 8 370.37 399.67 -177.18 354.37 86.6789 0 <2e-16 ***
## PTE 12 373.56 417.51 -174.78 349.56 4.8058 4 0.3078
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1P.tscale <- lmer(scale(log(P_ppm)) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
Mult_P.TSF <- glht(P.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_P.TSF) #no difference!##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(P_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.12394 0.12050 1.029 0.837
## 2yr - 3yr - 1yr < == 0 0.03890 0.13824 0.281 0.999
## 3yr - 4yr - 1yr < == 0 -0.19467 0.18281 -1.065 0.819
## Unburned - 1yr < == 0 -0.08501 0.12618 -0.674 0.961
## 2yr - 3yr - 1yr - 2yr == 0 -0.08505 0.13802 -0.616 0.972
## 3yr - 4yr - 1yr - 2yr == 0 -0.31861 0.18214 -1.749 0.395
## Unburned - 1yr - 2yr == 0 -0.20895 0.12624 -1.655 0.453
## 3yr - 4yr - 2yr - 3yr == 0 -0.23356 0.18855 -1.239 0.722
## Unburned - 2yr - 3yr == 0 -0.12391 0.15028 -0.825 0.920
## Unburned - 3yr - 4yr == 0 0.10965 0.19644 0.558 0.980
## (Adjusted p values reported -- single-step method)confint(Mult_P.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(P_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.7107
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.12394 -0.20269 0.45057
## 2yr - 3yr - 1yr < == 0 0.03890 -0.33582 0.41361
## 3yr - 4yr - 1yr < == 0 -0.19467 -0.69021 0.30088
## Unburned - 1yr < == 0 -0.08501 -0.42703 0.25701
## 2yr - 3yr - 1yr - 2yr == 0 -0.08505 -0.45917 0.28908
## 3yr - 4yr - 1yr - 2yr == 0 -0.31861 -0.81233 0.17512
## Unburned - 1yr - 2yr == 0 -0.20895 -0.55114 0.13324
## 3yr - 4yr - 2yr - 3yr == 0 -0.23356 -0.74465 0.27753
## Unburned - 2yr - 3yr == 0 -0.12391 -0.53127 0.28346
## Unburned - 3yr - 4yr == 0 0.10965 -0.42284 0.64215Mult_P.ESD <- glht(P.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_P.ESD) #sandy lowest, salo lower than thin clay but not lower than clayey##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(P_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 -0.18971 0.17686 -1.073 0.8023
## SaLo - Clayey == 0 -0.06339 0.13696 -0.463 0.9892
## Sandy - Clayey == 0 -0.92992 0.12636 -7.359 <0.001 ***
## ThCl - Clayey == 0 0.86409 0.33768 2.559 0.0681 .
## SaLo - Loamy == 0 0.12631 0.17887 0.706 0.9494
## Sandy - Loamy == 0 -0.74021 0.15956 -4.639 <0.001 ***
## ThCl - Loamy == 0 1.05379 0.35340 2.982 0.0205 *
## Sandy - SaLo == 0 -0.86653 0.11085 -7.817 <0.001 ***
## ThCl - SaLo == 0 0.92748 0.31985 2.900 0.0264 *
## ThCl - Sandy == 0 1.79401 0.32233 5.566 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_P.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(P_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6757
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 -0.18971 -0.66292 0.28351
## SaLo - Clayey == 0 -0.06339 -0.42985 0.30307
## Sandy - Clayey == 0 -0.92992 -1.26802 -0.59182
## ThCl - Clayey == 0 0.86409 -0.03945 1.76763
## SaLo - Loamy == 0 0.12631 -0.35229 0.60492
## Sandy - Loamy == 0 -0.74021 -1.16716 -0.31327
## ThCl - Loamy == 0 1.05379 0.10822 1.99937
## Sandy - SaLo == 0 -0.86653 -1.16313 -0.56993
## ThCl - SaLo == 0 0.92748 0.07165 1.78331
## ThCl - Sandy == 0 1.79401 0.93156 2.65645K
#K Overall
Knorm <- fitdist(SoilJuly181920$K_ppm, "norm")
plot(Knorm) #qq wonky
Klognorm <- fitdist(log(SoilJuly181920$K_ppm), "norm")
plot(Klognorm) #seems good
Kgamma <- fitdist(SoilJuly181920$K_ppm, "gamma")
plot(Kgamma) #qq wonky on right end
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=K_ppm, color=TSF))+
geom_boxplot(size=1)+
labs(y="K_ppm", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
#with log
ggplot(SoilJuly181920, aes(x=Location, y=log(K_ppm), color=TSF))+
geom_boxplot(size=1)+
labs(y="log(K_ppm)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
Knull <-lmer(log(K_ppm) ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #singular## boundary (singular) fit: see ?isSingularKT <- lmer(log(K_ppm) ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
KE <- lmer(log(K_ppm) ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #singular## boundary (singular) fit: see ?isSingularKTE <- lmer(log(K_ppm) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #singular## boundary (singular) fit: see ?isSingularanova(Knull, KT, KE, KTE)## Data: SoilJuly181920
## Models:
## Knull: log(K_ppm) ~ 1 + (1 | Location/Year)
## KT: log(K_ppm) ~ TSF + (1 | Location/Year)
## KE: log(K_ppm) ~ ESD + (1 | Location/Year)
## KTE: log(K_ppm) ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Knull 4 206.57 221.22 -99.287 198.57
## KT 8 202.34 231.65 -93.171 186.34 12.232 4 0.015704 *
## KE 8 147.69 176.99 -65.842 131.69 54.656 0 < 2.2e-16 ***
## KTE 12 140.73 184.69 -58.366 116.73 14.952 4 0.004801 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1K.tscale <- lmer(scale(log(K_ppm)) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)## boundary (singular) fit: see ?isSingularMult_K.TSF <- glht(K.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_K.TSF) #RB, 1-2, and 3-4 are higher than unburned##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(K_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.002468 0.133650 0.018 1.00000
## 2yr - 3yr - 1yr < == 0 -0.092252 0.152797 -0.604 0.97362
## 3yr - 4yr - 1yr < == 0 0.122725 0.200961 0.611 0.97250
## Unburned - 1yr < == 0 -0.480132 0.139076 -3.452 0.00471 **
## 2yr - 3yr - 1yr - 2yr == 0 -0.094720 0.152575 -0.621 0.97080
## 3yr - 4yr - 1yr - 2yr == 0 0.120258 0.200260 0.601 0.97414
## Unburned - 1yr - 2yr == 0 -0.482600 0.139126 -3.469 0.00464 **
## 3yr - 4yr - 2yr - 3yr == 0 0.214977 0.208113 1.033 0.83553
## Unburned - 2yr - 3yr == 0 -0.387880 0.164365 -2.360 0.12216
## Unburned - 3yr - 4yr == 0 -0.602857 0.213909 -2.818 0.03723 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_K.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(K_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.7124
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.002468 -0.360043 0.364979
## 2yr - 3yr - 1yr < == 0 -0.092252 -0.506695 0.322192
## 3yr - 4yr - 1yr < == 0 0.122725 -0.422358 0.667808
## Unburned - 1yr < == 0 -0.480132 -0.857360 -0.102904
## 2yr - 3yr - 1yr - 2yr == 0 -0.094720 -0.508562 0.319123
## 3yr - 4yr - 1yr - 2yr == 0 0.120258 -0.422925 0.663440
## Unburned - 1yr - 2yr == 0 -0.482600 -0.859963 -0.105236
## 3yr - 4yr - 2yr - 3yr == 0 0.214977 -0.349504 0.779458
## Unburned - 2yr - 3yr == 0 -0.387880 -0.833701 0.057941
## Unburned - 3yr - 4yr == 0 -0.602857 -1.183060 -0.022655Mult_K.ESD <- glht(K.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_K.ESD) #sandy lower than clayey, loamy, and salo##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(K_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.12937 0.19239 0.672 0.958
## SaLo - Clayey == 0 0.14269 0.14835 0.962 0.859
## Sandy - Clayey == 0 -0.77311 0.13836 -5.588 <0.001 ***
## ThCl - Clayey == 0 -0.63529 0.36660 -1.733 0.385
## SaLo - Loamy == 0 0.01332 0.18966 0.070 1.000
## Sandy - Loamy == 0 -0.90248 0.17327 -5.209 <0.001 ***
## ThCl - Loamy == 0 -0.76466 0.38200 -2.002 0.241
## Sandy - SaLo == 0 -0.91580 0.12054 -7.598 <0.001 ***
## ThCl - SaLo == 0 -0.77798 0.35071 -2.218 0.155
## ThCl - Sandy == 0 0.13782 0.35181 0.392 0.994
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_K.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(K_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6756
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.12937 -0.38540 0.64414
## SaLo - Clayey == 0 0.14269 -0.25425 0.53964
## Sandy - Clayey == 0 -0.77311 -1.14331 -0.40291
## ThCl - Clayey == 0 -0.63529 -1.61617 0.34559
## SaLo - Loamy == 0 0.01332 -0.49415 0.52080
## Sandy - Loamy == 0 -0.90248 -1.36608 -0.43887
## ThCl - Loamy == 0 -0.76466 -1.78674 0.25743
## Sandy - SaLo == 0 -0.91580 -1.23831 -0.59329
## ThCl - SaLo == 0 -0.77798 -1.71636 0.16040
## ThCl - Sandy == 0 0.13782 -0.80350 1.07915Ca
#Ca Overall
Canorm <- fitdist(SoilJuly181920$Ca_ppm, "norm")
plot(Canorm) #all pretty bad!
Calognorm <- fitdist(log(SoilJuly181920$Ca_ppm), "norm")
plot(Calognorm) #qq and pp are off more than other log's, but it's what we have 
#Cagamma <- fitdist(SoilJuly181920$Ca_ppm, "gamma") #errors out
#plot(Cagamma) #qq wonky on right end
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=Ca_ppm, color=TSF))+
geom_boxplot(size=1)+
labs(y="Ca_ppm", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
#with log
ggplot(SoilJuly181920, aes(x=Location, y=log(Ca_ppm), color=TSF))+
geom_boxplot(size=1)+
labs(y="log(Ca_ppm)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
ggplot(SoilJuly181920, aes(x=Location, y=log(Ca_ppm), color=ESD))+
geom_boxplot(size=1)+
labs(y="log(Ca_ppm)", x="Location")+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
Canull <-lmer(log(Ca_ppm) ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
CaT <- lmer(log(Ca_ppm) ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
CaE <- lmer(log(Ca_ppm) ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
CaTE <- lmer(log(Ca_ppm) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
anova(Canull, CaT, CaE, CaTE)## Data: SoilJuly181920
## Models:
## Canull: log(Ca_ppm) ~ 1 + (1 | Location/Year)
## CaT: log(Ca_ppm) ~ TSF + (1 | Location/Year)
## CaE: log(Ca_ppm) ~ ESD + (1 | Location/Year)
## CaTE: log(Ca_ppm) ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Canull 4 432.76 447.41 -212.38 424.76
## CaT 8 438.48 467.78 -211.24 422.48 2.2761 4 0.6851
## CaE 8 417.61 446.91 -200.80 401.61 20.8714 0 <2e-16 ***
## CaTE 12 421.98 465.94 -198.99 397.98 3.6262 4 0.4590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ca.tscale <- lmer(scale(log(Ca_ppm)) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
Mult_Ca.TSF <- glht(Ca.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_Ca.TSF) #No TSF difference##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Ca_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.01238 0.13255 0.093 1.000
## 2yr - 3yr - 1yr < == 0 0.04106 0.15094 0.272 0.999
## 3yr - 4yr - 1yr < == 0 -0.05178 0.19738 -0.262 0.999
## Unburned - 1yr < == 0 -0.20908 0.13694 -1.527 0.537
## 2yr - 3yr - 1yr - 2yr == 0 0.02868 0.15072 0.190 1.000
## 3yr - 4yr - 1yr - 2yr == 0 -0.06416 0.19666 -0.326 0.997
## Unburned - 1yr - 2yr == 0 -0.22146 0.13698 -1.617 0.478
## 3yr - 4yr - 2yr - 3yr == 0 -0.09284 0.20527 -0.452 0.991
## Unburned - 2yr - 3yr == 0 -0.25014 0.16035 -1.560 0.515
## Unburned - 3yr - 4yr == 0 -0.15730 0.20763 -0.758 0.941
## (Adjusted p values reported -- single-step method)confint(Mult_Ca.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Ca_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.7137
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.01238 -0.34733 0.37210
## 2yr - 3yr - 1yr < == 0 0.04106 -0.36855 0.45067
## 3yr - 4yr - 1yr < == 0 -0.05178 -0.58742 0.48386
## Unburned - 1yr < == 0 -0.20908 -0.58069 0.16254
## 2yr - 3yr - 1yr - 2yr == 0 0.02868 -0.38033 0.43768
## 3yr - 4yr - 1yr - 2yr == 0 -0.06416 -0.59785 0.46953
## Unburned - 1yr - 2yr == 0 -0.22146 -0.59317 0.15025
## 3yr - 4yr - 2yr - 3yr == 0 -0.09284 -0.64989 0.46421
## Unburned - 2yr - 3yr == 0 -0.25014 -0.68528 0.18501
## Unburned - 3yr - 4yr == 0 -0.15730 -0.72075 0.40615Mult_Ca.ESD <- glht(Ca.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_Ca.ESD) #sandy is higher than clayey and thin claypan? Loamy is highest##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Ca_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.5916 0.1948 3.037 0.01716 *
## SaLo - Clayey == 0 0.3907 0.1509 2.590 0.06280 .
## Sandy - Clayey == 0 0.4438 0.1391 3.191 0.01040 *
## ThCl - Clayey == 0 -0.8073 0.3720 -2.170 0.17079
## SaLo - Loamy == 0 -0.2008 0.1973 -1.018 0.83098
## Sandy - Loamy == 0 -0.1477 0.1758 -0.840 0.90813
## ThCl - Loamy == 0 -1.3988 0.3895 -3.592 0.00250 **
## Sandy - SaLo == 0 0.0531 0.1221 0.435 0.99149
## ThCl - SaLo == 0 -1.1980 0.3520 -3.403 0.00506 **
## ThCl - Sandy == 0 -1.2511 0.3549 -3.525 0.00326 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_Ca.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Ca_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6751
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.59155 0.07041 1.11269
## SaLo - Clayey == 0 0.39072 -0.01284 0.79428
## Sandy - Clayey == 0 0.44382 0.07170 0.81594
## ThCl - Clayey == 0 -0.80730 -1.80242 0.18782
## SaLo - Loamy == 0 -0.20083 -0.72874 0.32708
## Sandy - Loamy == 0 -0.14773 -0.61797 0.32252
## ThCl - Loamy == 0 -1.39884 -2.44066 -0.35703
## Sandy - SaLo == 0 0.05310 -0.27343 0.37963
## ThCl - SaLo == 0 -1.19802 -2.13975 -0.25628
## ThCl - Sandy == 0 -1.25112 -2.20058 -0.30165Mg
#Mg Overall
Mgnorm <- fitdist(SoilJuly181920$Mg_ppm, "norm")
plot(Mgnorm) #pp is wonky
Mglognorm <- fitdist(log(SoilJuly181920$Mg_ppm), "norm")
plot(Mglognorm) #better 
Mggamma <- fitdist(SoilJuly181920$Mg_ppm, "gamma") #errors out
plot(Mggamma) #better than normal, not better than log
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=Mg_ppm, color=TSF))+
geom_boxplot(size=1)+
labs(y="Mg_ppm (scales free)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1, scales = "free")
ggplot(SoilJuly181920, aes(x=Location, y=Mg_ppm, color=ESD))+
geom_boxplot(size=1)+
labs(y="Mg_ppm (scales free)", x="Location")+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1, scales = "free")
#with log
ggplot(SoilJuly181920, aes(x=Location, y=log(Mg_ppm), color=TSF))+
geom_boxplot(size=1)+
labs(y="log(Mg_ppm)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
ggplot(SoilJuly181920, aes(x=Location, y=log(Mg_ppm), color=ESD))+
geom_boxplot(size=1)+
labs(y="log(Mg_ppm) (scales free)", x="Location")+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1, scales = "free")
Mgnull <-lmer(log(Mg_ppm) ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
MgT <- lmer(log(Mg_ppm) ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
MgE <- lmer(log(Mg_ppm) ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
MgTE <- lmer(log(Mg_ppm) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
anova(Mgnull, MgT, MgE, MgTE)## Data: SoilJuly181920
## Models:
## Mgnull: log(Mg_ppm) ~ 1 + (1 | Location/Year)
## MgT: log(Mg_ppm) ~ TSF + (1 | Location/Year)
## MgE: log(Mg_ppm) ~ ESD + (1 | Location/Year)
## MgTE: log(Mg_ppm) ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Mgnull 4 242.41 257.06 -117.21 234.41
## MgT 8 245.24 274.55 -114.62 229.24 5.1688 4 0.2704
## MgE 8 228.63 257.94 -106.32 212.63 16.6122 0 <2e-16 ***
## MgTE 12 232.64 276.60 -104.32 208.64 3.9907 4 0.4073
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mg.tscale <- lmer(scale(log(Mg_ppm)) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
Mult_Mg.TSF <- glht(Mg.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_Mg.TSF) #No TSF difference##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Mg_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.06333 0.14587 0.434 0.992
## 2yr - 3yr - 1yr < == 0 0.04985 0.16692 0.299 0.998
## 3yr - 4yr - 1yr < == 0 -0.10525 0.21987 -0.479 0.989
## Unburned - 1yr < == 0 -0.21089 0.15204 -1.387 0.628
## 2yr - 3yr - 1yr - 2yr == 0 -0.01348 0.16666 -0.081 1.000
## 3yr - 4yr - 1yr - 2yr == 0 -0.16858 0.21908 -0.769 0.937
## Unburned - 1yr - 2yr == 0 -0.27423 0.15210 -1.803 0.364
## 3yr - 4yr - 2yr - 3yr == 0 -0.15510 0.22744 -0.682 0.959
## Unburned - 2yr - 3yr == 0 -0.26075 0.18005 -1.448 0.588
## Unburned - 3yr - 4yr == 0 -0.10565 0.23461 -0.450 0.991
## (Adjusted p values reported -- single-step method)confint(Mult_Mg.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Mg_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.7133
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.06333 -0.33244 0.45911
## 2yr - 3yr - 1yr < == 0 0.04985 -0.40304 0.50274
## 3yr - 4yr - 1yr < == 0 -0.10525 -0.70182 0.49132
## Unburned - 1yr < == 0 -0.21089 -0.62341 0.20162
## 2yr - 3yr - 1yr - 2yr == 0 -0.01348 -0.46568 0.43872
## 3yr - 4yr - 1yr - 2yr == 0 -0.16858 -0.76300 0.42584
## Unburned - 1yr - 2yr == 0 -0.27423 -0.68691 0.13845
## 3yr - 4yr - 2yr - 3yr == 0 -0.15510 -0.77220 0.46201
## Unburned - 2yr - 3yr == 0 -0.26075 -0.74928 0.22779
## Unburned - 3yr - 4yr == 0 -0.10565 -0.74221 0.53091Mult_Mg.ESD <- glht(Mg.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_Mg.ESD) #sandy is lower than thin clay and salo##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Mg_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 -0.03693 0.21267 -0.174 0.9998
## SaLo - Clayey == 0 0.22015 0.16443 1.339 0.6403
## Sandy - Clayey == 0 -0.32282 0.15230 -2.120 0.1903
## ThCl - Clayey == 0 0.76271 0.40580 1.880 0.3011
## SaLo - Loamy == 0 0.25708 0.21320 1.206 0.7253
## Sandy - Loamy == 0 -0.28589 0.19174 -1.491 0.5392
## ThCl - Loamy == 0 0.79964 0.42400 1.886 0.2982
## Sandy - SaLo == 0 -0.54297 0.13326 -4.075 <0.001 ***
## ThCl - SaLo == 0 0.54256 0.38567 1.407 0.5954
## ThCl - Sandy == 0 1.08553 0.38804 2.797 0.0356 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_Mg.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(log(Mg_ppm)) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6751
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 -0.03693 -0.60587 0.53200
## SaLo - Clayey == 0 0.22015 -0.21973 0.66003
## Sandy - Clayey == 0 -0.32282 -0.73024 0.08460
## ThCl - Clayey == 0 0.76271 -0.32285 1.84827
## SaLo - Loamy == 0 0.25708 -0.31327 0.82743
## Sandy - Loamy == 0 -0.28589 -0.79882 0.22704
## ThCl - Loamy == 0 0.79964 -0.33463 1.93391
## Sandy - SaLo == 0 -0.54297 -0.89945 -0.18648
## ThCl - SaLo == 0 0.54256 -0.48917 1.57430
## ThCl - Sandy == 0 1.08553 0.04747 2.12359Total Nitrogen
#TN Overall
TNnorm <- fitdist(SoilJuly181920$TotalN_percent, "norm")
plot(TNnorm) #qq could be better
TNlognorm <- fitdist(log(SoilJuly181920$TotalN_percent), "norm")
plot(TNlognorm) #better, but lots of negative values? 
TNgamma <- fitdist(SoilJuly181920$TotalN_percent, "gamma") #errors out
plot(TNgamma) #about the same as normal
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=TotalN_percent, color=TSF))+
geom_boxplot(size=1)+
labs(y="TotalN_percent", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
ggplot(SoilJuly181920, aes(x=Location, y=TotalN_percent, color=ESD))+
geom_boxplot(size=1)+
labs(y="TotalN_percent", x="Location")+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
TNnull <-lmer(TotalN_percent ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
TNT <- lmer(TotalN_percent ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
TNE <- lmer(TotalN_percent ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
TNTE <- lmer(TotalN_percent ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
anova(TNnull, TNT, TNE, TNTE)## Data: SoilJuly181920
## Models:
## TNnull: TotalN_percent ~ 1 + (1 | Location/Year)
## TNT: TotalN_percent ~ TSF + (1 | Location/Year)
## TNE: TotalN_percent ~ ESD + (1 | Location/Year)
## TNTE: TotalN_percent ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## TNnull 4 -859.54 -844.89 433.77 -867.54
## TNT 8 -858.41 -829.11 437.21 -874.41 6.8724 4 0.1428
## TNE 8 -920.65 -891.35 468.32 -936.65 62.2386 0 <2e-16 ***
## TNTE 12 -919.64 -875.68 471.82 -943.64 6.9901 4 0.1364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TN.tscale <- lmer(scale(TotalN_percent) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
Mult_TN.TSF <- glht(TN.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_TN.TSF) #RB higher than NYB##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalN_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 -0.03529 0.12916 -0.273 0.9987
## 2yr - 3yr - 1yr < == 0 -0.07438 0.14801 -0.503 0.9866
## 3yr - 4yr - 1yr < == 0 -0.05738 0.19540 -0.294 0.9983
## Unburned - 1yr < == 0 -0.33564 0.13497 -2.487 0.0898 .
## 2yr - 3yr - 1yr - 2yr == 0 -0.03909 0.14778 -0.265 0.9989
## 3yr - 4yr - 1yr - 2yr == 0 -0.02209 0.19468 -0.113 1.0000
## Unburned - 1yr - 2yr == 0 -0.30035 0.13503 -2.224 0.1650
## 3yr - 4yr - 2yr - 3yr == 0 0.01700 0.20179 0.084 1.0000
## Unburned - 2yr - 3yr == 0 -0.26126 0.16035 -1.629 0.4696
## Unburned - 3yr - 4yr == 0 -0.27826 0.20931 -1.329 0.6651
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TN.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalN_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.7108
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 -0.03529 -0.38542 0.31484
## 2yr - 3yr - 1yr < == 0 -0.07438 -0.47560 0.32684
## 3yr - 4yr - 1yr < == 0 -0.05738 -0.58708 0.47231
## Unburned - 1yr < == 0 -0.33564 -0.70151 0.03023
## 2yr - 3yr - 1yr - 2yr == 0 -0.03909 -0.43969 0.36151
## 3yr - 4yr - 1yr - 2yr == 0 -0.02209 -0.54984 0.50565
## Unburned - 1yr - 2yr == 0 -0.30035 -0.66639 0.06569
## 3yr - 4yr - 2yr - 3yr == 0 0.01700 -0.53002 0.56401
## Unburned - 2yr - 3yr == 0 -0.26126 -0.69593 0.17341
## Unburned - 3yr - 4yr == 0 -0.27826 -0.84565 0.28913Mult_TN.ESD <- glht(TN.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_TN.ESD) #sandy is lower than clayey, salo, and loamy; salo is higher than clayey##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalN_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.08182 0.18959 0.432 0.9918
## SaLo - Clayey == 0 0.38875 0.14681 2.648 0.0538 .
## Sandy - Clayey == 0 -0.61022 0.13545 -4.505 <0.001 ***
## ThCl - Clayey == 0 0.04003 0.36199 0.111 1.0000
## SaLo - Loamy == 0 0.30693 0.19176 1.601 0.4666
## Sandy - Loamy == 0 -0.69204 0.17105 -4.046 <0.001 ***
## ThCl - Loamy == 0 -0.04179 0.37884 -0.110 1.0000
## Sandy - SaLo == 0 -0.99896 0.11882 -8.407 <0.001 ***
## ThCl - SaLo == 0 -0.34872 0.34285 -1.017 0.8313
## ThCl - Sandy == 0 0.65024 0.34552 1.882 0.2997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TN.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalN_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6747
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.081818 -0.425269 0.588905
## SaLo - Clayey == 0 0.388746 -0.003932 0.781423
## Sandy - Clayey == 0 -0.610219 -0.972501 -0.247936
## ThCl - Clayey == 0 0.040026 -0.928189 1.008240
## SaLo - Loamy == 0 0.306928 -0.205973 0.819829
## Sandy - Loamy == 0 -0.692037 -1.149541 -0.234532
## ThCl - Loamy == 0 -0.041792 -1.055079 0.971495
## Sandy - SaLo == 0 -0.998964 -1.316777 -0.681151
## ThCl - SaLo == 0 -0.348720 -1.265744 0.568304
## ThCl - Sandy == 0 0.650244 -0.273909 1.574397Total Carbon
#TC Overall
TCnorm <- fitdist(SoilJuly181920$TotalC_percent, "norm")
plot(TCnorm) #qq could be better
TClognorm <- fitdist(log(SoilJuly181920$TotalC_percent), "norm")
plot(TClognorm) #better, but some of negative values? 
TCgamma <- fitdist(SoilJuly181920$TotalC_percent, "gamma") #errors out
plot(TCgamma) #this is better
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=TotalC_percent, color=TSF))+
geom_boxplot(size=1)+
labs(y="TotalC_percent", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
ggplot(SoilJuly181920, aes(x=Location, y=TotalC_percent, color=ESD))+
geom_boxplot(size=1)+
labs(y="TotalC_percent", x="Location")+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
TCnull <-lmer(TotalC_percent ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
TCT <- lmer(TotalC_percent ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
TCE <- lmer(TotalC_percent ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
TCTE <- lmer(TotalC_percent ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
anova(TCnull, TCT, TCE, TCTE)## Data: SoilJuly181920
## Models:
## TCnull: TotalC_percent ~ 1 + (1 | Location/Year)
## TCT: TotalC_percent ~ TSF + (1 | Location/Year)
## TCE: TotalC_percent ~ ESD + (1 | Location/Year)
## TCTE: TotalC_percent ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## TCnull 4 555.43 570.08 -273.71 547.43
## TCT 8 554.26 583.56 -269.13 538.26 9.1676 4 0.05704 .
## TCE 8 491.74 521.04 -237.87 475.74 62.5231 0 < 2e-16 ***
## TCTE 12 489.61 533.57 -232.81 465.61 10.1276 4 0.03833 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TC.tscale <- lmer(scale(TotalC_percent) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
Mult_TC.TSF <- glht(TC.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_TC.TSF) #1-2 higher than NYB##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalC_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.06101 0.13302 0.459 0.9906
## 2yr - 3yr - 1yr < == 0 -0.12858 0.15043 -0.855 0.9109
## 3yr - 4yr - 1yr < == 0 -0.03924 0.19467 -0.202 0.9996
## Unburned - 1yr < == 0 -0.34344 0.13576 -2.530 0.0815 .
## 2yr - 3yr - 1yr - 2yr == 0 -0.18959 0.15022 -1.262 0.7087
## 3yr - 4yr - 1yr - 2yr == 0 -0.10025 0.19398 -0.517 0.9853
## Unburned - 1yr - 2yr == 0 -0.40445 0.13577 -2.979 0.0232 *
## 3yr - 4yr - 2yr - 3yr == 0 0.08934 0.20391 0.438 0.9921
## Unburned - 2yr - 3yr == 0 -0.21486 0.15635 -1.374 0.6376
## Unburned - 3yr - 4yr == 0 -0.30420 0.20057 -1.517 0.5443
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TC.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalC_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.7146
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.06101 -0.30008 0.42210
## 2yr - 3yr - 1yr < == 0 -0.12858 -0.53693 0.27977
## 3yr - 4yr - 1yr < == 0 -0.03924 -0.56770 0.48921
## Unburned - 1yr < == 0 -0.34344 -0.71197 0.02508
## 2yr - 3yr - 1yr - 2yr == 0 -0.18959 -0.59737 0.21819
## 3yr - 4yr - 1yr - 2yr == 0 -0.10025 -0.62682 0.42631
## Unburned - 1yr - 2yr == 0 -0.40445 -0.77301 -0.03590
## 3yr - 4yr - 2yr - 3yr == 0 0.08934 -0.46420 0.64287
## Unburned - 2yr - 3yr == 0 -0.21486 -0.63929 0.20956
## Unburned - 3yr - 4yr == 0 -0.30420 -0.84867 0.24028Mult_TC.ESD <- glht(TC.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_TC.ESD) #sandy is lower than clayey, salo, and loamy; salo is higher than clayey##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalC_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.33490 0.19522 1.716 0.39434
## SaLo - Clayey == 0 0.58938 0.15108 3.901 < 0.001 ***
## Sandy - Clayey == 0 -0.45994 0.13945 -3.298 0.00729 **
## ThCl - Clayey == 0 0.28571 0.37270 0.767 0.93268
## SaLo - Loamy == 0 0.25448 0.19735 1.289 0.67210
## Sandy - Loamy == 0 -0.79484 0.17612 -4.513 < 0.001 ***
## ThCl - Loamy == 0 -0.04919 0.39006 -0.126 0.99993
## Sandy - SaLo == 0 -1.04932 0.12227 -8.582 < 0.001 ***
## ThCl - SaLo == 0 -0.30367 0.35290 -0.860 0.90072
## ThCl - Sandy == 0 0.74565 0.35572 2.096 0.19949
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TC.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TotalC_percent) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6728
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.33490 -0.18688 0.85668
## SaLo - Clayey == 0 0.58938 0.18558 0.99318
## Sandy - Clayey == 0 -0.45994 -0.83266 -0.08722
## ThCl - Clayey == 0 0.28571 -0.71045 1.28187
## SaLo - Loamy == 0 0.25448 -0.27301 0.78197
## Sandy - Loamy == 0 -0.79484 -1.26558 -0.32410
## ThCl - Loamy == 0 -0.04919 -1.09177 0.99338
## Sandy - SaLo == 0 -1.04932 -1.37612 -0.72251
## ThCl - SaLo == 0 -0.30367 -1.24692 0.63959
## ThCl - Sandy == 0 0.74565 -0.20513 1.69642Total C:N
#TCNR Overall
TCNRnorm <- fitdist(SoilJuly181920$TCNR, "norm")
plot(TCNRnorm) #qq and pp are off
TCNRlognorm <- fitdist(log(SoilJuly181920$TCNR), "norm")
plot(TCNRlognorm) #not really any better
TCNRgamma <- fitdist(SoilJuly181920$TCNR, "gamma")
plot(TCNRgamma) #none are great, use normal?
#untransformed
ggplot(SoilJuly181920, aes(x=Location, y=TCNR, color=TSF))+
geom_boxplot(size=1)+
labs(y="TCNR (scales free)", x="Location")+
scale_color_brewer(palette = "Set1")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1, scales = "free")
ggplot(SoilJuly181920, aes(x=Location, y=TCNR, color=ESD))+
geom_boxplot(size=1)+
labs(y="TCNR", x="Location")+
scale_color_brewer(palette = "Dark2")+
theme_bw()+
theme(legend.position = "bottom",
axis.title=element_text(size=12, face="bold"),
axis.text=element_text(size=12),
axis.text.x=element_text(angle=33, hjust=1),
legend.text=element_text(size=8),
legend.title = element_text(size=8),
panel.grid.major.x = element_blank())+
facet_wrap( ~ Year,ncol = 1)
TCNRnull <-lmer(TCNR ~ 1 + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #failed to converge?## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00273451 (tol = 0.002, component 1)TCNRT <- lmer(TCNR ~ TSF + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #failed to converge## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00332473 (tol = 0.002, component 1)TCNRE <- lmer(TCNR ~ ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #converged
TCNRTE <- lmer(TCNR ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE) #converged
anova(TCNRnull, TCNRT, TCNRE, TCNRTE) #none are really good## Data: SoilJuly181920
## Models:
## TCNRnull: TCNR ~ 1 + (1 | Location/Year)
## TCNRT: TCNR ~ TSF + (1 | Location/Year)
## TCNRE: TCNR ~ ESD + (1 | Location/Year)
## TCNRTE: TCNR ~ TSF + ESD + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## TCNRnull 4 1109.8 1124.5 -550.91 1101.8
## TCNRT 8 1111.2 1140.5 -547.61 1095.2 6.6050 4 0.1583
## TCNRE 8 1103.6 1132.9 -543.80 1087.6 7.6138 0 <2e-16 ***
## TCNRTE 12 1106.4 1150.3 -541.18 1082.4 5.2481 4 0.2628
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(TCNRTE)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: TCNR ~ TSF + ESD + (1 | Location/Year)
## Data: SoilJuly181920
##
## AIC BIC logLik deviance df.resid
## 1106.4 1150.3 -541.2 1082.4 276
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0956 -0.4277 -0.0565 0.2591 7.1247
##
## Random effects:
## Groups Name Variance Std.Dev.
## Year:Location (Intercept) 0.02591 0.1610
## Location (Intercept) 0.10646 0.3263
## Residual 2.43284 1.5598
## Number of obs: 288, groups: Year:Location, 18; Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 9.14068 0.30959 29.525
## TSF1yr - 2yr 0.27314 0.26025 1.050
## TSF2yr - 3yr -0.27395 0.29233 -0.937
## TSF3yr - 4yr -0.48320 0.37442 -1.291
## TSFUnburned -0.05413 0.26258 -0.206
## ESDLoamy 0.81454 0.37458 2.175
## ESDSaLo 0.78682 0.28861 2.726
## ESDSandy 0.14726 0.26931 0.547
## ESDThCl 0.83139 0.71362 1.165
##
## Correlation of Fixed Effects:
## (Intr) TSF1-2 TSF2-3 TSF3-4 TSFUnb ESDLmy ESDSaL ESDSnd
## TSF1yr-2yr -0.430
## TSF2yr-3yr -0.386 0.447
## TSF3yr-4yr -0.336 0.352 0.319
## TSFUnburned -0.428 0.496 0.430 0.338
## ESDLoamy -0.425 0.004 -0.006 0.040 0.032
## ESDSaLo -0.594 0.034 0.042 0.098 0.016 0.404
## ESDSandy -0.597 -0.003 0.005 0.010 -0.006 0.491 0.649
## ESDThCl -0.277 0.006 -0.022 0.068 0.071 0.181 0.308 0.294TCNR.tscale <- lmer(scale(TCNR) ~ TSF + ESD + (1|Location/Year), data=SoilJuly181920, REML = FALSE)
Mult_TCNR.TSF <- glht(TCNR.tscale, linfct=mcp(TSF = "Tukey"))
summary(Mult_TCNR.TSF) #No difference##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TCNR) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.16525 0.15745 1.050 0.828
## 2yr - 3yr - 1yr < == 0 -0.16574 0.17686 -0.937 0.880
## 3yr - 4yr - 1yr < == 0 -0.29233 0.22652 -1.291 0.691
## Unburned - 1yr < == 0 -0.03275 0.15886 -0.206 1.000
## 2yr - 3yr - 1yr - 2yr == 0 -0.33099 0.17664 -1.874 0.325
## 3yr - 4yr - 1yr - 2yr == 0 -0.45758 0.22576 -2.027 0.247
## Unburned - 1yr - 2yr == 0 -0.19800 0.15884 -1.247 0.719
## 3yr - 4yr - 2yr - 3yr == 0 -0.12659 0.23882 -0.530 0.984
## Unburned - 2yr - 3yr == 0 0.13299 0.17986 0.739 0.946
## Unburned - 3yr - 4yr == 0 0.25958 0.22853 1.136 0.783
## (Adjusted p values reported -- single-step method)confint(Mult_TCNR.TSF)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TCNR) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.717
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.16525 -0.26253 0.59303
## 2yr - 3yr - 1yr < == 0 -0.16574 -0.64626 0.31478
## 3yr - 4yr - 1yr < == 0 -0.29233 -0.90778 0.32311
## Unburned - 1yr < == 0 -0.03275 -0.46436 0.39886
## 2yr - 3yr - 1yr - 2yr == 0 -0.33099 -0.81090 0.14893
## 3yr - 4yr - 1yr - 2yr == 0 -0.45758 -1.07096 0.15580
## Unburned - 1yr - 2yr == 0 -0.19800 -0.62955 0.23356
## 3yr - 4yr - 2yr - 3yr == 0 -0.12659 -0.77546 0.52228
## Unburned - 2yr - 3yr == 0 0.13299 -0.35569 0.62167
## Unburned - 3yr - 4yr == 0 0.25958 -0.36133 0.88049Mult_TCNR.ESD <- glht(TCNR.tscale, linfct=mcp(ESD = "Tukey"))
summary(Mult_TCNR.ESD) #Salo higher than clayey, sandy lower than salo##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TCNR) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.49279 0.22662 2.175 0.1699
## SaLo - Clayey == 0 0.47602 0.17461 2.726 0.0434 *
## Sandy - Clayey == 0 0.08909 0.16293 0.547 0.9800
## ThCl - Clayey == 0 0.50299 0.43173 1.165 0.7507
## SaLo - Loamy == 0 -0.01677 0.22328 -0.075 1.0000
## Sandy - Loamy == 0 -0.40370 0.20408 -1.978 0.2522
## ThCl - Loamy == 0 0.01020 0.44990 0.023 1.0000
## Sandy - SaLo == 0 -0.38693 0.14186 -2.728 0.0434 *
## ThCl - SaLo == 0 0.02697 0.41288 0.065 1.0000
## ThCl - Sandy == 0 0.41390 0.41429 0.999 0.8411
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_TCNR.ESD)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = scale(TCNR) ~ TSF + ESD + (1 | Location/Year),
## data = SoilJuly181920, REML = FALSE)
##
## Quantile = 2.6768
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.492790 -0.113811 1.099390
## SaLo - Clayey == 0 0.476022 0.008632 0.943412
## Sandy - Clayey == 0 0.089093 -0.347030 0.525215
## ThCl - Clayey == 0 0.502988 -0.652670 1.658646
## SaLo - Loamy == 0 -0.016768 -0.614442 0.580907
## Sandy - Loamy == 0 -0.403697 -0.949978 0.142584
## ThCl - Loamy == 0 0.010198 -1.194090 1.214487
## Sandy - SaLo == 0 -0.386929 -0.766659 -0.007200
## ThCl - SaLo == 0 0.026966 -1.078237 1.132169
## ThCl - Sandy == 0 0.413895 -0.695060 1.522850Moisture
#Moisture
Mnorm <- fitdist(Soil1920$Moisture, "norm")
plot(Mnorm) #ok
Mlognorm <- fitdist(log(Soil1920$Moisture), "norm")
plot(Mlognorm) #best
Mgamma <- fitdist(Soil1920$Moisture, "gamma")
plot(Mgamma) 
#Full nesting
Mnull1 <-lmer(log(Moisture) ~ 1 + (1|Location/Month/Year), data=Soil1920, REML = FALSE)
MT1 <-lmer(log(Moisture) ~ TSF + (1|Location/Month/Year), data=Soil1920, REML = FALSE)## boundary (singular) fit: see ?isSingularME1 <-lmer(log(Moisture) ~ ESD + (1|Location/Month/Year), data=Soil1920, REML = FALSE)
MTE1 <-lmer(log(Moisture) ~ TSF + ESD + (1|Location/Month/Year), data=Soil1920, REML = FALSE)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00426226 (tol = 0.002, component 1)anova(Mnull1, MT1, ME1, MTE1)## Data: Soil1920
## Models:
## Mnull1: log(Moisture) ~ 1 + (1 | Location/Month/Year)
## MT1: log(Moisture) ~ TSF + (1 | Location/Month/Year)
## ME1: log(Moisture) ~ ESD + (1 | Location/Month/Year)
## MTE1: log(Moisture) ~ TSF + ESD + (1 | Location/Month/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Mnull1 5 341.16 364.38 -165.582 331.16
## MT1 9 316.71 358.50 -149.355 298.71 32.454 4 1.545e-06 ***
## ME1 9 212.40 254.19 -97.198 194.40 104.314 0 < 2.2e-16 ***
## MTE1 13 187.65 248.02 -80.827 161.65 32.741 4 1.350e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mult_MT1 <- glht(MTE1, linfct=mcp(TSF = "Tukey"))
summary(Mult_MT1) #Everybody higher than NYB! ##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil1920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 0.064146 0.024817 2.585 0.07061 .
## 2yr - 3yr - 1yr < == 0 0.008855 0.024806 0.357 0.99640
## 3yr - 4yr - 1yr < == 0 0.061268 0.031893 1.921 0.29863
## Unburned - 1yr < == 0 -0.114635 0.031660 -3.621 0.00268 **
## 2yr - 3yr - 1yr - 2yr == 0 -0.055291 0.024643 -2.244 0.15877
## 3yr - 4yr - 1yr - 2yr == 0 -0.002878 0.032149 -0.090 0.99998
## Unburned - 1yr - 2yr == 0 -0.178781 0.031793 -5.623 < 0.001 ***
## 3yr - 4yr - 2yr - 3yr == 0 0.052413 0.031893 1.643 0.46118
## Unburned - 2yr - 3yr == 0 -0.123489 0.031821 -3.881 < 0.001 ***
## Unburned - 3yr - 4yr == 0 -0.175903 0.040103 -4.386 < 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_MT1)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil1920, REML = FALSE)
##
## Quantile = 2.7135
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 0.064146 -0.003194 0.131486
## 2yr - 3yr - 1yr < == 0 0.008855 -0.058456 0.076166
## 3yr - 4yr - 1yr < == 0 0.061268 -0.025273 0.147809
## Unburned - 1yr < == 0 -0.114635 -0.200545 -0.028724
## 2yr - 3yr - 1yr - 2yr == 0 -0.055291 -0.122161 0.011579
## 3yr - 4yr - 1yr - 2yr == 0 -0.002878 -0.090115 0.084358
## Unburned - 1yr - 2yr == 0 -0.178781 -0.265051 -0.092510
## 3yr - 4yr - 2yr - 3yr == 0 0.052413 -0.034129 0.138956
## Unburned - 2yr - 3yr == 0 -0.123489 -0.209835 -0.037144
## Unburned - 3yr - 4yr == 0 -0.175903 -0.284722 -0.067083Mult_ME1 <- glht(MTE1, linfct=mcp(ESD = "Tukey"))
summary(Mult_ME1) #Sandy is the lowest; Salo Thin Claypan > Clayey Loamy##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil1920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.02318 0.03640 0.637 0.96485
## SaLo - Clayey == 0 0.12904 0.02833 4.555 < 0.001 ***
## Sandy - Clayey == 0 -0.13082 0.02588 -5.056 < 0.001 ***
## ThCl - Clayey == 0 0.25007 0.06973 3.586 0.00281 **
## SaLo - Loamy == 0 0.10586 0.03692 2.868 0.02877 *
## Sandy - Loamy == 0 -0.15401 0.03275 -4.703 < 0.001 ***
## ThCl - Loamy == 0 0.22689 0.07281 3.116 0.01348 *
## Sandy - SaLo == 0 -0.25987 0.02277 -11.412 < 0.001 ***
## ThCl - SaLo == 0 0.12102 0.06572 1.842 0.32102
## ThCl - Sandy == 0 0.38089 0.06645 5.732 < 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_ME1)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil1920, REML = FALSE)
##
## Quantile = 2.6739
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.023182 -0.074139 0.120503
## SaLo - Clayey == 0 0.129044 0.053294 0.204793
## Sandy - Clayey == 0 -0.130824 -0.200014 -0.061633
## ThCl - Clayey == 0 0.250068 0.063615 0.436521
## SaLo - Loamy == 0 0.105862 0.007148 0.204576
## Sandy - Loamy == 0 -0.154005 -0.241563 -0.066448
## ThCl - Loamy == 0 0.226887 0.032211 0.421562
## Sandy - SaLo == 0 -0.259868 -0.320757 -0.198978
## ThCl - SaLo == 0 0.121024 -0.054699 0.296748
## ThCl - Sandy == 0 0.380892 0.203223 0.558561#For testing diiferences between months
Mnull <-lmer(log(Moisture) ~ 1 + (1|Location/Year), data=Soil1920, REML = FALSE)
MMTE <-lmer(log(Moisture) ~ Month + TSF + ESD + (1|Location/Year), data=Soil1920, REML = FALSE)
Mult_MM <- glht(MMTE, linfct=mcp(Month = "Tukey"))
summary(Mult_MM) #September > July > June > August##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ Month + TSF + ESD + (1 | Location/Year),
## data = Soil1920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## July - June == 0 0.26885 0.03219 8.351 <1e-04 ***
## August - June == 0 -0.17190 0.03219 -5.340 <1e-04 ***
## September - June == 0 0.40820 0.03219 12.680 <1e-04 ***
## August - July == 0 -0.44075 0.03219 -13.692 <1e-04 ***
## September - July == 0 0.13935 0.03219 4.329 <1e-04 ***
## September - August == 0 0.58010 0.03219 18.020 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_MM)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ Month + TSF + ESD + (1 | Location/Year),
## data = Soil1920, REML = FALSE)
##
## Quantile = 2.5699
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## July - June == 0 0.26885 0.18612 0.35158
## August - June == 0 -0.17190 -0.25463 -0.08917
## September - June == 0 0.40820 0.32547 0.49093
## August - July == 0 -0.44075 -0.52348 -0.35802
## September - July == 0 0.13935 0.05662 0.22208
## September - August == 0 0.58010 0.49737 0.66283Soil1920 <-
Soil1920 %>%
unite("MonthTSF", c("Month", "TSF"), remove = F ) %>%
as_tibble()
Soil1920 <-
Soil1920 %>%
unite("MonthESD", c("Month", "ESD"), remove = F ) %>%
as_tibble()This is a mess!
MTMint <- lmer(log(Moisture) ~ MonthTSF + (1|Location/Year), data=Soil1920, REML = FALSE)
Mult_MTMint <- glht(MTMint, linfct=mcp(MonthTSF = "Tukey"))
summary(Mult_MTMint) #a mess!## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(Moisture) ~ MonthTSF + (1 | Location/Year),
## data = Soil1920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value
## August_1yr < - August_1yr - 2yr == 0 -0.1491571 0.0646880 -2.306
## August_2yr - 3yr - August_1yr - 2yr == 0 -0.1068955 0.0646880 -1.652
## August_3yr - 4yr - August_1yr - 2yr == 0 0.1027790 0.0801984 1.282
## August_Unburned - August_1yr - 2yr == 0 -0.5139975 0.0801984 -6.409
## July_1yr - 2yr - August_1yr - 2yr == 0 0.3963020 0.0646880 6.126
## July_1yr < - August_1yr - 2yr == 0 0.2861533 0.0646880 4.424
## July_2yr - 3yr - August_1yr - 2yr == 0 0.3528516 0.0646880 5.455
## July_3yr - 4yr - August_1yr - 2yr == 0 0.0550594 0.0801984 0.687
## July_Unburned - August_1yr - 2yr == 0 0.4770062 0.0801984 5.948
## June_1yr - 2yr - August_1yr - 2yr == 0 0.0935199 0.0646880 1.446
## June_1yr < - August_1yr - 2yr == 0 0.1041125 0.0650360 1.601
## June_2yr - 3yr - August_1yr - 2yr == 0 0.0127489 0.0646880 0.197
## June_3yr - 4yr - August_1yr - 2yr == 0 0.1906136 0.0791613 2.408
## June_Unburned - August_1yr - 2yr == 0 -0.1660102 0.0801984 -2.070
## September_1yr - 2yr - August_1yr - 2yr == 0 0.5311686 0.0646880 8.211
## September_1yr < - August_1yr - 2yr == 0 0.4302591 0.0646880 6.651
## September_2yr - 3yr - August_1yr - 2yr == 0 0.4763019 0.0646880 7.363
## September_3yr - 4yr - August_1yr - 2yr == 0 0.3890918 0.0801984 4.852
## September_Unburned - August_1yr - 2yr == 0 0.4529313 0.0801984 5.648
## August_2yr - 3yr - August_1yr < == 0 0.0422616 0.0646880 0.653
## August_3yr - 4yr - August_1yr < == 0 0.2519361 0.0801984 3.141
## August_Unburned - August_1yr < == 0 -0.3648404 0.0801984 -4.549
## July_1yr - 2yr - August_1yr < == 0 0.5454591 0.0646880 8.432
## July_1yr < - August_1yr < == 0 0.4353104 0.0646880 6.729
## July_2yr - 3yr - August_1yr < == 0 0.5020087 0.0646880 7.760
## July_3yr - 4yr - August_1yr < == 0 0.2042165 0.0801984 2.546
## July_Unburned - August_1yr < == 0 0.6261633 0.0801984 7.808
## June_1yr - 2yr - August_1yr < == 0 0.2426770 0.0646880 3.751
## June_1yr < - August_1yr < == 0 0.2532696 0.0650360 3.894
## June_2yr - 3yr - August_1yr < == 0 0.1619060 0.0646880 2.503
## June_3yr - 4yr - August_1yr < == 0 0.3397707 0.0791613 4.292
## June_Unburned - August_1yr < == 0 -0.0168531 0.0801984 -0.210
## September_1yr - 2yr - August_1yr < == 0 0.6803257 0.0646880 10.517
## September_1yr < - August_1yr < == 0 0.5794162 0.0646880 8.957
## September_2yr - 3yr - August_1yr < == 0 0.6254590 0.0646880 9.669
## September_3yr - 4yr - August_1yr < == 0 0.5382489 0.0801984 6.711
## September_Unburned - August_1yr < == 0 0.6020884 0.0801984 7.507
## August_3yr - 4yr - August_2yr - 3yr == 0 0.2096746 0.0801984 2.614
## August_Unburned - August_2yr - 3yr == 0 -0.4071019 0.0801984 -5.076
## July_1yr - 2yr - August_2yr - 3yr == 0 0.5031976 0.0646880 7.779
## July_1yr < - August_2yr - 3yr == 0 0.3930489 0.0646880 6.076
## July_2yr - 3yr - August_2yr - 3yr == 0 0.4597472 0.0646880 7.107
## July_3yr - 4yr - August_2yr - 3yr == 0 0.1619549 0.0801984 2.019
## July_Unburned - August_2yr - 3yr == 0 0.5839018 0.0801984 7.281
## June_1yr - 2yr - August_2yr - 3yr == 0 0.2004155 0.0646880 3.098
## June_1yr < - August_2yr - 3yr == 0 0.2110080 0.0650360 3.244
## June_2yr - 3yr - August_2yr - 3yr == 0 0.1196445 0.0646880 1.850
## June_3yr - 4yr - August_2yr - 3yr == 0 0.2975091 0.0791613 3.758
## June_Unburned - August_2yr - 3yr == 0 -0.0591146 0.0801984 -0.737
## September_1yr - 2yr - August_2yr - 3yr == 0 0.6380641 0.0646880 9.864
## September_1yr < - August_2yr - 3yr == 0 0.5371546 0.0646880 8.304
## September_2yr - 3yr - August_2yr - 3yr == 0 0.5831974 0.0646880 9.016
## September_3yr - 4yr - August_2yr - 3yr == 0 0.4959874 0.0801984 6.185
## September_Unburned - August_2yr - 3yr == 0 0.5598268 0.0801984 6.981
## August_Unburned - August_3yr - 4yr == 0 -0.6167765 0.0948103 -6.505
## July_1yr - 2yr - August_3yr - 4yr == 0 0.2935230 0.0801984 3.660
## July_1yr < - August_3yr - 4yr == 0 0.1833743 0.0801984 2.287
## July_2yr - 3yr - August_3yr - 4yr == 0 0.2500726 0.0801984 3.118
## July_3yr - 4yr - August_3yr - 4yr == 0 -0.0477196 0.0914827 -0.522
## July_Unburned - August_3yr - 4yr == 0 0.3742272 0.0948103 3.947
## June_1yr - 2yr - August_3yr - 4yr == 0 -0.0092591 0.0801984 -0.115
## June_1yr < - August_3yr - 4yr == 0 0.0013334 0.0805208 0.017
## June_2yr - 3yr - August_3yr - 4yr == 0 -0.0900301 0.0801984 -1.123
## June_3yr - 4yr - August_3yr - 4yr == 0 0.0878346 0.0905740 0.970
## June_Unburned - August_3yr - 4yr == 0 -0.2687892 0.0948103 -2.835
## September_1yr - 2yr - August_3yr - 4yr == 0 0.4283896 0.0801984 5.342
## September_1yr < - August_3yr - 4yr == 0 0.3274800 0.0801984 4.083
## September_2yr - 3yr - August_3yr - 4yr == 0 0.3735228 0.0801984 4.657
## September_3yr - 4yr - August_3yr - 4yr == 0 0.2863128 0.0914827 3.130
## September_Unburned - August_3yr - 4yr == 0 0.3501523 0.0948103 3.693
## July_1yr - 2yr - August_Unburned == 0 0.9102995 0.0801984 11.351
## July_1yr < - August_Unburned == 0 0.8001508 0.0801984 9.977
## July_2yr - 3yr - August_Unburned == 0 0.8668491 0.0801984 10.809
## July_3yr - 4yr - August_Unburned == 0 0.5690569 0.0948103 6.002
## July_Unburned - August_Unburned == 0 0.9910037 0.0914827 10.833
## June_1yr - 2yr - August_Unburned == 0 0.6075174 0.0801984 7.575
## June_1yr < - August_Unburned == 0 0.6181099 0.0804379 7.684
## June_2yr - 3yr - August_Unburned == 0 0.5267464 0.0801984 6.568
## June_3yr - 4yr - August_Unburned == 0 0.7046111 0.0939355 7.501
## June_Unburned - August_Unburned == 0 0.3479873 0.0914827 3.804
## September_1yr - 2yr - August_Unburned == 0 1.0451661 0.0801984 13.032
## September_1yr < - August_Unburned == 0 0.9442565 0.0801984 11.774
## September_2yr - 3yr - August_Unburned == 0 0.9902993 0.0801984 12.348
## September_3yr - 4yr - August_Unburned == 0 0.9030893 0.0948103 9.525
## September_Unburned - August_Unburned == 0 0.9669288 0.0914827 10.570
## July_1yr < - July_1yr - 2yr == 0 -0.1101487 0.0646880 -1.703
## July_2yr - 3yr - July_1yr - 2yr == 0 -0.0434504 0.0646880 -0.672
## July_3yr - 4yr - July_1yr - 2yr == 0 -0.3412426 0.0801984 -4.255
## July_Unburned - July_1yr - 2yr == 0 0.0807042 0.0801984 1.006
## June_1yr - 2yr - July_1yr - 2yr == 0 -0.3027821 0.0646880 -4.681
## June_1yr < - July_1yr - 2yr == 0 -0.2921895 0.0650360 -4.493
## June_2yr - 3yr - July_1yr - 2yr == 0 -0.3835531 0.0646880 -5.929
## June_3yr - 4yr - July_1yr - 2yr == 0 -0.2056884 0.0791613 -2.598
## June_Unburned - July_1yr - 2yr == 0 -0.5623122 0.0801984 -7.012
## September_1yr - 2yr - July_1yr - 2yr == 0 0.1348666 0.0646880 2.085
## September_1yr < - July_1yr - 2yr == 0 0.0339570 0.0646880 0.525
## September_2yr - 3yr - July_1yr - 2yr == 0 0.0799999 0.0646880 1.237
## September_3yr - 4yr - July_1yr - 2yr == 0 -0.0072102 0.0801984 -0.090
## September_Unburned - July_1yr - 2yr == 0 0.0566293 0.0801984 0.706
## July_2yr - 3yr - July_1yr < == 0 0.0666983 0.0646880 1.031
## July_3yr - 4yr - July_1yr < == 0 -0.2310939 0.0801984 -2.882
## July_Unburned - July_1yr < == 0 0.1908529 0.0801984 2.380
## June_1yr - 2yr - July_1yr < == 0 -0.1926334 0.0646880 -2.978
## June_1yr < - July_1yr < == 0 -0.1820409 0.0650360 -2.799
## June_2yr - 3yr - July_1yr < == 0 -0.2734044 0.0646880 -4.227
## June_3yr - 4yr - July_1yr < == 0 -0.0955397 0.0791613 -1.207
## June_Unburned - July_1yr < == 0 -0.4521635 0.0801984 -5.638
## September_1yr - 2yr - July_1yr < == 0 0.2450153 0.0646880 3.788
## September_1yr < - July_1yr < == 0 0.1441057 0.0646880 2.228
## September_2yr - 3yr - July_1yr < == 0 0.1901485 0.0646880 2.939
## September_3yr - 4yr - July_1yr < == 0 0.1029385 0.0801984 1.284
## September_Unburned - July_1yr < == 0 0.1667780 0.0801984 2.080
## July_3yr - 4yr - July_2yr - 3yr == 0 -0.2977922 0.0801984 -3.713
## July_Unburned - July_2yr - 3yr == 0 0.1241546 0.0801984 1.548
## June_1yr - 2yr - July_2yr - 3yr == 0 -0.2593317 0.0646880 -4.009
## June_1yr < - July_2yr - 3yr == 0 -0.2487392 0.0650360 -3.825
## June_2yr - 3yr - July_2yr - 3yr == 0 -0.3401027 0.0646880 -5.258
## June_3yr - 4yr - July_2yr - 3yr == 0 -0.1622380 0.0791613 -2.049
## June_Unburned - July_2yr - 3yr == 0 -0.5188618 0.0801984 -6.470
## September_1yr - 2yr - July_2yr - 3yr == 0 0.1783170 0.0646880 2.757
## September_1yr < - July_2yr - 3yr == 0 0.0774074 0.0646880 1.197
## September_2yr - 3yr - July_2yr - 3yr == 0 0.1234503 0.0646880 1.908
## September_3yr - 4yr - July_2yr - 3yr == 0 0.0362402 0.0801984 0.452
## September_Unburned - July_2yr - 3yr == 0 0.1000797 0.0801984 1.248
## July_Unburned - July_3yr - 4yr == 0 0.4219468 0.0948103 4.450
## June_1yr - 2yr - July_3yr - 4yr == 0 0.0384605 0.0801984 0.480
## June_1yr < - July_3yr - 4yr == 0 0.0490531 0.0805208 0.609
## June_2yr - 3yr - July_3yr - 4yr == 0 -0.0423105 0.0801984 -0.528
## June_3yr - 4yr - July_3yr - 4yr == 0 0.1355542 0.0905740 1.497
## June_Unburned - July_3yr - 4yr == 0 -0.2210696 0.0948103 -2.332
## September_1yr - 2yr - July_3yr - 4yr == 0 0.4761092 0.0801984 5.937
## September_1yr < - July_3yr - 4yr == 0 0.3751997 0.0801984 4.678
## September_2yr - 3yr - July_3yr - 4yr == 0 0.4212425 0.0801984 5.253
## September_3yr - 4yr - July_3yr - 4yr == 0 0.3340324 0.0914827 3.651
## September_Unburned - July_3yr - 4yr == 0 0.3978719 0.0948103 4.197
## June_1yr - 2yr - July_Unburned == 0 -0.3834863 0.0801984 -4.782
## June_1yr < - July_Unburned == 0 -0.3728938 0.0804379 -4.636
## June_2yr - 3yr - July_Unburned == 0 -0.4642573 0.0801984 -5.789
## June_3yr - 4yr - July_Unburned == 0 -0.2863926 0.0939355 -3.049
## June_Unburned - July_Unburned == 0 -0.6430164 0.0914827 -7.029
## September_1yr - 2yr - July_Unburned == 0 0.0541624 0.0801984 0.675
## September_1yr < - July_Unburned == 0 -0.0467472 0.0801984 -0.583
## September_2yr - 3yr - July_Unburned == 0 -0.0007044 0.0801984 -0.009
## September_3yr - 4yr - July_Unburned == 0 -0.0879144 0.0948103 -0.927
## September_Unburned - July_Unburned == 0 -0.0240749 0.0914827 -0.263
## June_1yr < - June_1yr - 2yr == 0 0.0105926 0.0650360 0.163
## June_2yr - 3yr - June_1yr - 2yr == 0 -0.0807710 0.0646880 -1.249
## June_3yr - 4yr - June_1yr - 2yr == 0 0.0970937 0.0791613 1.227
## June_Unburned - June_1yr - 2yr == 0 -0.2595301 0.0801984 -3.236
## September_1yr - 2yr - June_1yr - 2yr == 0 0.4376487 0.0646880 6.766
## September_1yr < - June_1yr - 2yr == 0 0.3367391 0.0646880 5.206
## September_2yr - 3yr - June_1yr - 2yr == 0 0.3827820 0.0646880 5.917
## September_3yr - 4yr - June_1yr - 2yr == 0 0.2955719 0.0801984 3.686
## September_Unburned - June_1yr - 2yr == 0 0.3594114 0.0801984 4.482
## June_2yr - 3yr - June_1yr < == 0 -0.0913636 0.0650360 -1.405
## June_3yr - 4yr - June_1yr < == 0 0.0865011 0.0795011 1.088
## June_Unburned - June_1yr < == 0 -0.2701226 0.0804379 -3.358
## September_1yr - 2yr - June_1yr < == 0 0.4270561 0.0650360 6.566
## September_1yr < - June_1yr < == 0 0.3261466 0.0650360 5.015
## September_2yr - 3yr - June_1yr < == 0 0.3721894 0.0650360 5.723
## September_3yr - 4yr - June_1yr < == 0 0.2849794 0.0805208 3.539
## September_Unburned - June_1yr < == 0 0.3488188 0.0804379 4.336
## June_3yr - 4yr - June_2yr - 3yr == 0 0.1778647 0.0791613 2.247
## June_Unburned - June_2yr - 3yr == 0 -0.1787591 0.0801984 -2.229
## September_1yr - 2yr - June_2yr - 3yr == 0 0.5184197 0.0646880 8.014
## September_1yr < - June_2yr - 3yr == 0 0.4175102 0.0646880 6.454
## September_2yr - 3yr - June_2yr - 3yr == 0 0.4635530 0.0646880 7.166
## September_3yr - 4yr - June_2yr - 3yr == 0 0.3763429 0.0801984 4.693
## September_Unburned - June_2yr - 3yr == 0 0.4401824 0.0801984 5.489
## June_Unburned - June_3yr - 4yr == 0 -0.3566238 0.0939355 -3.796
## September_1yr - 2yr - June_3yr - 4yr == 0 0.3405550 0.0791613 4.302
## September_1yr < - June_3yr - 4yr == 0 0.2396455 0.0791613 3.027
## September_2yr - 3yr - June_3yr - 4yr == 0 0.2856883 0.0791613 3.609
## September_3yr - 4yr - June_3yr - 4yr == 0 0.1984782 0.0905740 2.191
## September_Unburned - June_3yr - 4yr == 0 0.2623177 0.0939355 2.793
## September_1yr - 2yr - June_Unburned == 0 0.6971788 0.0801984 8.693
## September_1yr < - June_Unburned == 0 0.5962692 0.0801984 7.435
## September_2yr - 3yr - June_Unburned == 0 0.6423120 0.0801984 8.009
## September_3yr - 4yr - June_Unburned == 0 0.5551020 0.0948103 5.855
## September_Unburned - June_Unburned == 0 0.6189415 0.0914827 6.766
## September_1yr < - September_1yr - 2yr == 0 -0.1009095 0.0646880 -1.560
## September_2yr - 3yr - September_1yr - 2yr == 0 -0.0548667 0.0646880 -0.848
## September_3yr - 4yr - September_1yr - 2yr == 0 -0.1420768 0.0801984 -1.772
## September_Unburned - September_1yr - 2yr == 0 -0.0782373 0.0801984 -0.976
## September_2yr - 3yr - September_1yr < == 0 0.0460428 0.0646880 0.712
## September_3yr - 4yr - September_1yr < == 0 -0.0411672 0.0801984 -0.513
## September_Unburned - September_1yr < == 0 0.0226722 0.0801984 0.283
## September_3yr - 4yr - September_2yr - 3yr == 0 -0.0872100 0.0801984 -1.087
## September_Unburned - September_2yr - 3yr == 0 -0.0233706 0.0801984 -0.291
## September_Unburned - September_3yr - 4yr == 0 0.0638395 0.0948103 0.673
## Pr(>|z|)
## August_1yr < - August_1yr - 2yr == 0 0.7218
## August_2yr - 3yr - August_1yr - 2yr == 0 0.9837
## August_3yr - 4yr - August_1yr - 2yr == 0 0.9993
## August_Unburned - August_1yr - 2yr == 0 <0.01 ***
## July_1yr - 2yr - August_1yr - 2yr == 0 <0.01 ***
## July_1yr < - August_1yr - 2yr == 0 <0.01 **
## July_2yr - 3yr - August_1yr - 2yr == 0 <0.01 ***
## July_3yr - 4yr - August_1yr - 2yr == 0 1.0000
## July_Unburned - August_1yr - 2yr == 0 <0.01 ***
## June_1yr - 2yr - August_1yr - 2yr == 0 0.9966
## June_1yr < - August_1yr - 2yr == 0 0.9884
## June_2yr - 3yr - August_1yr - 2yr == 0 1.0000
## June_3yr - 4yr - August_1yr - 2yr == 0 0.6462
## June_Unburned - August_1yr - 2yr == 0 0.8656
## September_1yr - 2yr - August_1yr - 2yr == 0 <0.01 ***
## September_1yr < - August_1yr - 2yr == 0 <0.01 ***
## September_2yr - 3yr - August_1yr - 2yr == 0 <0.01 ***
## September_3yr - 4yr - August_1yr - 2yr == 0 <0.01 ***
## September_Unburned - August_1yr - 2yr == 0 <0.01 ***
## August_2yr - 3yr - August_1yr < == 0 1.0000
## August_3yr - 4yr - August_1yr < == 0 0.1558
## August_Unburned - August_1yr < == 0 <0.01 ***
## July_1yr - 2yr - August_1yr < == 0 <0.01 ***
## July_1yr < - August_1yr < == 0 <0.01 ***
## July_2yr - 3yr - August_1yr < == 0 <0.01 ***
## July_3yr - 4yr - August_1yr < == 0 0.5371
## July_Unburned - August_1yr < == 0 <0.01 ***
## June_1yr - 2yr - August_1yr < == 0 0.0236 *
## June_1yr < - August_1yr < == 0 0.0142 *
## June_2yr - 3yr - August_1yr < == 0 0.5719
## June_3yr - 4yr - August_1yr < == 0 <0.01 **
## June_Unburned - August_1yr < == 0 1.0000
## September_1yr - 2yr - August_1yr < == 0 <0.01 ***
## September_1yr < - August_1yr < == 0 <0.01 ***
## September_2yr - 3yr - August_1yr < == 0 <0.01 ***
## September_3yr - 4yr - August_1yr < == 0 <0.01 ***
## September_Unburned - August_1yr < == 0 <0.01 ***
## August_3yr - 4yr - August_2yr - 3yr == 0 0.4823
## August_Unburned - August_2yr - 3yr == 0 <0.01 ***
## July_1yr - 2yr - August_2yr - 3yr == 0 <0.01 ***
## July_1yr < - August_2yr - 3yr == 0 <0.01 ***
## July_2yr - 3yr - August_2yr - 3yr == 0 <0.01 ***
## July_3yr - 4yr - August_2yr - 3yr == 0 0.8896
## July_Unburned - August_2yr - 3yr == 0 <0.01 ***
## June_1yr - 2yr - August_2yr - 3yr == 0 0.1761
## June_1yr < - August_2yr - 3yr == 0 0.1175
## June_2yr - 3yr - August_2yr - 3yr == 0 0.9486
## June_3yr - 4yr - August_2yr - 3yr == 0 0.0225 *
## June_Unburned - August_2yr - 3yr == 0 1.0000
## September_1yr - 2yr - August_2yr - 3yr == 0 <0.01 ***
## September_1yr < - August_2yr - 3yr == 0 <0.01 ***
## September_2yr - 3yr - August_2yr - 3yr == 0 <0.01 ***
## September_3yr - 4yr - August_2yr - 3yr == 0 <0.01 ***
## September_Unburned - August_2yr - 3yr == 0 <0.01 ***
## August_Unburned - August_3yr - 4yr == 0 <0.01 ***
## July_1yr - 2yr - August_3yr - 4yr == 0 0.0316 *
## July_1yr < - August_3yr - 4yr == 0 0.7347
## July_2yr - 3yr - August_3yr - 4yr == 0 0.1675
## July_3yr - 4yr - August_3yr - 4yr == 0 1.0000
## July_Unburned - August_3yr - 4yr == 0 0.0108 *
## June_1yr - 2yr - August_3yr - 4yr == 0 1.0000
## June_1yr < - August_3yr - 4yr == 0 1.0000
## June_2yr - 3yr - August_3yr - 4yr == 0 0.9999
## June_3yr - 4yr - August_3yr - 4yr == 0 1.0000
## June_Unburned - August_3yr - 4yr == 0 0.3217
## September_1yr - 2yr - August_3yr - 4yr == 0 <0.01 ***
## September_1yr < - August_3yr - 4yr == 0 <0.01 **
## September_2yr - 3yr - August_3yr - 4yr == 0 <0.01 ***
## September_3yr - 4yr - August_3yr - 4yr == 0 0.1616
## September_Unburned - August_3yr - 4yr == 0 0.0290 *
## July_1yr - 2yr - August_Unburned == 0 <0.01 ***
## July_1yr < - August_Unburned == 0 <0.01 ***
## July_2yr - 3yr - August_Unburned == 0 <0.01 ***
## July_3yr - 4yr - August_Unburned == 0 <0.01 ***
## July_Unburned - August_Unburned == 0 <0.01 ***
## June_1yr - 2yr - August_Unburned == 0 <0.01 ***
## June_1yr < - August_Unburned == 0 <0.01 ***
## June_2yr - 3yr - August_Unburned == 0 <0.01 ***
## June_3yr - 4yr - August_Unburned == 0 <0.01 ***
## June_Unburned - August_Unburned == 0 0.0200 *
## September_1yr - 2yr - August_Unburned == 0 <0.01 ***
## September_1yr < - August_Unburned == 0 <0.01 ***
## September_2yr - 3yr - August_Unburned == 0 <0.01 ***
## September_3yr - 4yr - August_Unburned == 0 <0.01 ***
## September_Unburned - August_Unburned == 0 <0.01 ***
## July_1yr < - July_1yr - 2yr == 0 0.9776
## July_2yr - 3yr - July_1yr - 2yr == 0 1.0000
## July_3yr - 4yr - July_1yr - 2yr == 0 <0.01 **
## July_Unburned - July_1yr - 2yr == 0 1.0000
## June_1yr - 2yr - July_1yr - 2yr == 0 <0.01 ***
## June_1yr < - July_1yr - 2yr == 0 <0.01 ***
## June_2yr - 3yr - July_1yr - 2yr == 0 <0.01 ***
## June_3yr - 4yr - July_1yr - 2yr == 0 0.4955
## June_Unburned - July_1yr - 2yr == 0 <0.01 ***
## September_1yr - 2yr - July_1yr - 2yr == 0 0.8591
## September_1yr < - July_1yr - 2yr == 0 1.0000
## September_2yr - 3yr - July_1yr - 2yr == 0 0.9996
## September_3yr - 4yr - July_1yr - 2yr == 0 1.0000
## September_Unburned - July_1yr - 2yr == 0 1.0000
## July_2yr - 3yr - July_1yr < == 0 1.0000
## July_3yr - 4yr - July_1yr < == 0 0.2920
## July_Unburned - July_1yr < == 0 0.6679
## June_1yr - 2yr - July_1yr < == 0 0.2346
## June_1yr < - July_1yr < == 0 0.3461
## June_2yr - 3yr - July_1yr < == 0 <0.01 **
## June_3yr - 4yr - July_1yr < == 0 0.9997
## June_Unburned - July_1yr < == 0 <0.01 ***
## September_1yr - 2yr - July_1yr < == 0 0.0204 *
## September_1yr < - July_1yr < == 0 0.7762
## September_2yr - 3yr - July_1yr < == 0 0.2566
## September_3yr - 4yr - July_1yr < == 0 0.9993
## September_Unburned - July_1yr < == 0 0.8611
## July_3yr - 4yr - July_2yr - 3yr == 0 0.0276 *
## July_Unburned - July_2yr - 3yr == 0 0.9923
## June_1yr - 2yr - July_2yr - 3yr == 0 <0.01 **
## June_1yr < - July_2yr - 3yr == 0 0.0183 *
## June_2yr - 3yr - July_2yr - 3yr == 0 <0.01 ***
## June_3yr - 4yr - July_2yr - 3yr == 0 0.8750
## June_Unburned - July_2yr - 3yr == 0 <0.01 ***
## September_1yr - 2yr - July_2yr - 3yr == 0 0.3740
## September_1yr < - July_2yr - 3yr == 0 0.9997
## September_2yr - 3yr - July_2yr - 3yr == 0 0.9315
## September_3yr - 4yr - July_2yr - 3yr == 0 1.0000
## September_Unburned - July_2yr - 3yr == 0 0.9995
## July_Unburned - July_3yr - 4yr == 0 <0.01 **
## June_1yr - 2yr - July_3yr - 4yr == 0 1.0000
## June_1yr < - July_3yr - 4yr == 0 1.0000
## June_2yr - 3yr - July_3yr - 4yr == 0 1.0000
## June_3yr - 4yr - July_3yr - 4yr == 0 0.9948
## June_Unburned - July_3yr - 4yr == 0 0.7025
## September_1yr - 2yr - July_3yr - 4yr == 0 <0.01 ***
## September_1yr < - July_3yr - 4yr == 0 <0.01 ***
## September_2yr - 3yr - July_3yr - 4yr == 0 <0.01 ***
## September_3yr - 4yr - July_3yr - 4yr == 0 0.0338 *
## September_Unburned - July_3yr - 4yr == 0 <0.01 **
## June_1yr - 2yr - July_Unburned == 0 <0.01 ***
## June_1yr < - July_Unburned == 0 <0.01 ***
## June_2yr - 3yr - July_Unburned == 0 <0.01 ***
## June_3yr - 4yr - July_Unburned == 0 0.1978
## June_Unburned - July_Unburned == 0 <0.01 ***
## September_1yr - 2yr - July_Unburned == 0 1.0000
## September_1yr < - July_Unburned == 0 1.0000
## September_2yr - 3yr - July_Unburned == 0 1.0000
## September_3yr - 4yr - July_Unburned == 0 1.0000
## September_Unburned - July_Unburned == 0 1.0000
## June_1yr < - June_1yr - 2yr == 0 1.0000
## June_2yr - 3yr - June_1yr - 2yr == 0 0.9995
## June_3yr - 4yr - June_1yr - 2yr == 0 0.9996
## June_Unburned - June_1yr - 2yr == 0 0.1224
## September_1yr - 2yr - June_1yr - 2yr == 0 <0.01 ***
## September_1yr < - June_1yr - 2yr == 0 <0.01 ***
## September_2yr - 3yr - June_1yr - 2yr == 0 <0.01 ***
## September_3yr - 4yr - June_1yr - 2yr == 0 0.0295 *
## September_Unburned - June_1yr - 2yr == 0 <0.01 **
## June_2yr - 3yr - June_1yr < == 0 0.9977
## June_3yr - 4yr - June_1yr < == 0 0.9999
## June_Unburned - June_1yr < == 0 0.0838 .
## September_1yr - 2yr - June_1yr < == 0 <0.01 ***
## September_1yr < - June_1yr < == 0 <0.01 ***
## September_2yr - 3yr - June_1yr < == 0 <0.01 ***
## September_3yr - 4yr - June_1yr < == 0 0.0482 *
## September_Unburned - June_1yr < == 0 <0.01 **
## June_3yr - 4yr - June_2yr - 3yr == 0 0.7610
## June_Unburned - June_2yr - 3yr == 0 0.7734
## September_1yr - 2yr - June_2yr - 3yr == 0 <0.01 ***
## September_1yr < - June_2yr - 3yr == 0 <0.01 ***
## September_2yr - 3yr - June_2yr - 3yr == 0 <0.01 ***
## September_3yr - 4yr - June_2yr - 3yr == 0 <0.01 ***
## September_Unburned - June_2yr - 3yr == 0 <0.01 ***
## June_Unburned - June_3yr - 4yr == 0 0.0202 *
## September_1yr - 2yr - June_3yr - 4yr == 0 <0.01 **
## September_1yr < - June_3yr - 4yr == 0 0.2092
## September_2yr - 3yr - June_3yr - 4yr == 0 0.0381 *
## September_3yr - 4yr - June_3yr - 4yr == 0 0.7967
## September_Unburned - June_3yr - 4yr == 0 0.3505
## September_1yr - 2yr - June_Unburned == 0 <0.01 ***
## September_1yr < - June_Unburned == 0 <0.01 ***
## September_2yr - 3yr - June_Unburned == 0 <0.01 ***
## September_3yr - 4yr - June_Unburned == 0 <0.01 ***
## September_Unburned - June_Unburned == 0 <0.01 ***
## September_1yr < - September_1yr - 2yr == 0 0.9915
## September_2yr - 3yr - September_1yr - 2yr == 0 1.0000
## September_3yr - 4yr - September_1yr - 2yr == 0 0.9656
## September_Unburned - September_1yr - 2yr == 0 1.0000
## September_2yr - 3yr - September_1yr < == 0 1.0000
## September_3yr - 4yr - September_1yr < == 0 1.0000
## September_Unburned - September_1yr < == 0 1.0000
## September_3yr - 4yr - September_2yr - 3yr == 0 0.9999
## September_Unburned - September_2yr - 3yr == 0 1.0000
## September_Unburned - September_3yr - 4yr == 0 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(MTMint)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 0.07158387 0.245583932
## .sig02 0.00000000 0.249924267
## .sigma 0.30157583 0.333578914
## (Intercept) 2.10573708 2.401303183
## MonthTSFAugust_1yr < -0.27610457 -0.022209649
## MonthTSFAugust_2yr - 3yr -0.23384301 0.020051918
## MonthTSFAugust_3yr - 4yr -0.05506473 0.260519575
## MonthTSFAugust_Unburned -0.67173804 -0.356153730
## MonthTSFJuly_1yr - 2yr 0.26935455 0.523249477
## MonthTSFJuly_1yr < 0.15920587 0.413100798
## MonthTSFJuly_2yr - 3yr 0.22590416 0.479799081
## MonthTSFJuly_3yr - 4yr -0.10278437 0.212799945
## MonthTSFJuly_Unburned 0.31926566 0.634849973
## MonthTSFJune_1yr - 2yr -0.03342755 0.220467373
## MonthTSFJune_1yr < -0.02351644 0.231745339
## MonthTSFJune_2yr - 3yr -0.11419856 0.139696370
## MonthTSFJune_3yr - 4yr 0.03477529 0.346338166
## MonthTSFJune_Unburned -0.32375075 -0.008166436
## MonthTSFSeptember_1yr - 2yr 0.40422113 0.658116052
## MonthTSFSeptember_1yr < 0.30331160 0.557206523
## MonthTSFSeptember_2yr - 3yr 0.34935441 0.603249332
## MonthTSFSeptember_3yr - 4yr 0.23124807 0.546832373
## MonthTSFSeptember_Unburned 0.29519072 0.610775035Nitrate
#Nitrate
Nnorm <- fitdist(Soil181920$NO3_ppm, "norm")
plot(Nnorm) #Not good!
Nlognorm <- fitdist(log(Soil181920$NO3_ppm), "norm")
plot(Nlognorm) #best, but not perfect
Ngamma <- fitdist(Soil181920$NO3_ppm, "gamma")
plot(Ngamma) #not good
#Full nesting
Nnull1 <-lmer(log(NO3_ppm) ~ 1 + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## boundary (singular) fit: see ?isSingularNT1 <-lmer(log(NO3_ppm) ~ TSF + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## boundary (singular) fit: see ?isSingularNE1 <-lmer(log(NO3_ppm) ~ ESD + (1|Location/Month/Year), data=Soil181920, REML = FALSE)
NTE1 <-lmer(log(NO3_ppm) ~ TSF + ESD + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## boundary (singular) fit: see ?isSingularanova(Nnull1, NT1, NE1, NTE1) #not much improvement over the null## Data: Soil181920
## Models:
## Nnull1: log(NO3_ppm) ~ 1 + (1 | Location/Month/Year)
## NT1: log(NO3_ppm) ~ TSF + (1 | Location/Month/Year)
## NE1: log(NO3_ppm) ~ ESD + (1 | Location/Month/Year)
## NTE1: log(NO3_ppm) ~ TSF + ESD + (1 | Location/Month/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Nnull1 5 2295.4 2320.6 -1142.7 2285.4
## NT1 9 2286.3 2331.7 -1134.2 2268.3 17.074 4 0.001870 **
## NE1 9 2290.7 2336.2 -1136.3 2272.7 0.000 0 1.000000
## NTE1 13 2284.1 2349.7 -1129.0 2258.1 14.650 4 0.005486 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mult_NT1 <- glht(NTE1, linfct=mcp(TSF = "Tukey"))
summary(Mult_NT1) #RB higher than 3-4 and 2-3##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NO3_ppm) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 -0.05535 0.05004 -1.106 0.7974
## 2yr - 3yr - 1yr < == 0 -0.17492 0.05767 -3.033 0.0194 *
## 3yr - 4yr - 1yr < == 0 -0.21708 0.07658 -2.835 0.0353 *
## Unburned - 1yr < == 0 -0.10953 0.05280 -2.075 0.2238
## 2yr - 3yr - 1yr - 2yr == 0 -0.11957 0.05751 -2.079 0.2218
## 3yr - 4yr - 1yr - 2yr == 0 -0.16173 0.07620 -2.122 0.2034
## Unburned - 1yr - 2yr == 0 -0.05418 0.05281 -1.026 0.8383
## 3yr - 4yr - 2yr - 3yr == 0 -0.04216 0.07848 -0.537 0.9828
## Unburned - 2yr - 3yr == 0 0.06539 0.06348 1.030 0.8363
## Unburned - 3yr - 4yr == 0 0.10755 0.08324 1.292 0.6881
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_NT1)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NO3_ppm) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Quantile = 2.712
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 -0.055348 -0.191061 0.080365
## 2yr - 3yr - 1yr < == 0 -0.174919 -0.331312 -0.018525
## 3yr - 4yr - 1yr < == 0 -0.217080 -0.424765 -0.009395
## Unburned - 1yr < == 0 -0.109530 -0.252712 0.033652
## 2yr - 3yr - 1yr - 2yr == 0 -0.119570 -0.275544 0.036404
## 3yr - 4yr - 1yr - 2yr == 0 -0.161732 -0.368392 0.044928
## Unburned - 1yr - 2yr == 0 -0.054182 -0.197413 0.089050
## 3yr - 4yr - 2yr - 3yr == 0 -0.042161 -0.254992 0.170669
## Unburned - 2yr - 3yr == 0 0.065389 -0.106763 0.237540
## Unburned - 3yr - 4yr == 0 0.107550 -0.118201 0.333301Mult_NE1 <- glht(NTE1, linfct=mcp(ESD = "Tukey"))
summary(Mult_NE1) #Sandy is lower than Salo; ##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NO3_ppm) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.009215 0.073432 0.125 0.9999
## SaLo - Clayey == 0 0.076718 0.056937 1.347 0.6343
## Sandy - Clayey == 0 -0.053391 0.052435 -1.018 0.8308
## ThCl - Clayey == 0 0.208172 0.140148 1.485 0.5424
## SaLo - Loamy == 0 0.067503 0.074226 0.909 0.8811
## Sandy - Loamy == 0 -0.062606 0.066182 -0.946 0.8653
## ThCl - Loamy == 0 0.198957 0.146623 1.357 0.6280
## Sandy - SaLo == 0 -0.130108 0.045970 -2.830 0.0321 *
## ThCl - SaLo == 0 0.131454 0.132722 0.990 0.8445
## ThCl - Sandy == 0 0.261563 0.133734 1.956 0.2620
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_NE1)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NO3_ppm) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Quantile = 2.6752
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.009215 -0.187227 0.205656
## SaLo - Clayey == 0 0.076718 -0.075597 0.229032
## Sandy - Clayey == 0 -0.053391 -0.193663 0.086881
## ThCl - Clayey == 0 0.208172 -0.166746 0.583090
## SaLo - Loamy == 0 0.067503 -0.131063 0.266068
## Sandy - Loamy == 0 -0.062606 -0.239652 0.114441
## ThCl - Loamy == 0 0.198957 -0.193283 0.591196
## Sandy - SaLo == 0 -0.130108 -0.253086 -0.007131
## ThCl - SaLo == 0 0.131454 -0.223599 0.486507
## ThCl - Sandy == 0 0.261563 -0.096197 0.619322#For testing diiferences between months
Nnull <-lmer(log(NO3_ppm) ~ 1 + (1|Location/Year), data=Soil181920, REML = FALSE)
NM <-lmer(log(NO3_ppm) ~ Month + (1|Location/Year), data=Soil181920, REML = FALSE)
NT <-lmer(log(NO3_ppm) ~ TSF + (1|Location/Year), data=Soil181920, REML = FALSE)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00204976 (tol = 0.002, component 1)NE <-lmer(log(NO3_ppm) ~ ESD + (1|Location/Year), data=Soil181920, REML = FALSE)
NMTE <-lmer(log(NO3_ppm) ~ Month + ESD + TSF + (1|Location/Year), data=Soil181920, REML = FALSE)
anova(Nnull, NM, NT, NE, NMTE)## Data: Soil181920
## Models:
## Nnull: log(NO3_ppm) ~ 1 + (1 | Location/Year)
## NM: log(NO3_ppm) ~ Month + (1 | Location/Year)
## NT: log(NO3_ppm) ~ TSF + (1 | Location/Year)
## NE: log(NO3_ppm) ~ ESD + (1 | Location/Year)
## NMTE: log(NO3_ppm) ~ Month + ESD + TSF + (1 | Location/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Nnull 4 2584.1 2604.3 -1288.0 2576.1
## NM 7 2461.7 2497.0 -1223.8 2447.7 128.38 3 <2e-16 ***
## NT 8 2580.6 2621.0 -1282.3 2564.6 0.00 1 1
## NE 8 2583.4 2623.8 -1283.7 2567.4 0.00 0 1
## NMTE 15 2456.7 2532.5 -1213.4 2426.7 140.70 7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mult_NM <- glht(NMTE, linfct=mcp(Month = "Tukey"))
summary(Mult_NM) #September > July > June > August##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NO3_ppm) ~ Month + ESD + TSF + (1 | Location/Year),
## data = Soil181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## July - June == 0 -0.14948 0.05685 -2.630 0.04246 *
## August - June == 0 -0.20509 0.05685 -3.608 0.00169 **
## September - June == 0 0.40033 0.05685 7.042 < 0.001 ***
## August - July == 0 -0.05561 0.05684 -0.978 0.76183
## September - July == 0 0.54981 0.05684 9.672 < 0.001 ***
## September - August == 0 0.60542 0.05684 10.650 < 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(NM)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 0.1359622 0.32772729
## .sig02 0.0000000 0.25595299
## .sigma 0.6617481 0.71853261
## (Intercept) 1.0804451 1.37733052
## MonthJuly -0.2611914 -0.03587699
## MonthAugust -0.3168039 -0.09148952
## MonthSeptember 0.2886149 0.51392933#Ammonium
Anorm <- fitdist(Soil181920$NH4_ppm, "norm")
plot(Anorm) #Not good!
Alognorm <- fitdist(log(Soil181920$NH4_ppm), "norm")
plot(Alognorm) #best, but not perfect
Agamma <- fitdist(Soil181920$NH4_ppm, "gamma")
plot(Agamma) #not good
#Full nesting
Anull1 <-lmer(log(NH4_ppm) ~ 1 + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## boundary (singular) fit: see ?isSingularAT1 <-lmer(log(NH4_ppm) ~ TSF + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## boundary (singular) fit: see ?isSingularAE1 <-lmer(log(NH4_ppm) ~ ESD + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## boundary (singular) fit: see ?isSingularATE1 <-lmer(log(NH4_ppm) ~ TSF + ESD + (1|Location/Month/Year), data=Soil181920, REML = FALSE)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00250854 (tol = 0.002, component 1)anova(Anull1, AT1, AE1, ATE1) #ESD does not do much here## Data: Soil181920
## Models:
## Anull1: log(NH4_ppm) ~ 1 + (1 | Location/Month/Year)
## AT1: log(NH4_ppm) ~ TSF + (1 | Location/Month/Year)
## AE1: log(NH4_ppm) ~ ESD + (1 | Location/Month/Year)
## ATE1: log(NH4_ppm) ~ TSF + ESD + (1 | Location/Month/Year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Anull1 5 2577.8 2603.0 -1283.9 2567.8
## AT1 9 2574.8 2620.3 -1278.4 2556.8 10.959 4 0.02703 *
## AE1 9 2576.1 2621.6 -1279.1 2558.1 0.000 0 1.00000
## ATE1 13 2572.3 2638.0 -1273.2 2546.3 11.785 4 0.01902 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mult_AT1 <- glht(AT1, linfct=mcp(TSF = "Tukey"))
summary(Mult_AT1) #RB higher than 3-4 ## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NH4_ppm) ~ TSF + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1yr - 2yr - 1yr < == 0 -0.13405 0.05436 -2.466 0.0941 .
## 2yr - 3yr - 1yr < == 0 -0.08826 0.06278 -1.406 0.6144
## 3yr - 4yr - 1yr < == 0 -0.22983 0.08326 -2.761 0.0437 *
## Unburned - 1yr < == 0 -0.03005 0.05762 -0.522 0.9845
## 2yr - 3yr - 1yr - 2yr == 0 0.04579 0.06271 0.730 0.9475
## 3yr - 4yr - 1yr - 2yr == 0 -0.09578 0.08317 -1.152 0.7718
## Unburned - 1yr - 2yr == 0 0.10400 0.05760 1.805 0.3607
## 3yr - 4yr - 2yr - 3yr == 0 -0.14157 0.08522 -1.661 0.4480
## Unburned - 2yr - 3yr == 0 0.05821 0.06960 0.836 0.9161
## Unburned - 3yr - 4yr == 0 0.19978 0.09175 2.177 0.1811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_AT1)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NH4_ppm) ~ TSF + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Quantile = 2.7105
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 1yr - 2yr - 1yr < == 0 -0.134048 -0.281378 0.013281
## 2yr - 3yr - 1yr < == 0 -0.088262 -0.258412 0.081888
## 3yr - 4yr - 1yr < == 0 -0.229828 -0.455490 -0.004166
## Unburned - 1yr < == 0 -0.030053 -0.186217 0.126112
## 2yr - 3yr - 1yr - 2yr == 0 0.045787 -0.124195 0.215768
## 3yr - 4yr - 1yr - 2yr == 0 -0.095780 -0.321196 0.129637
## Unburned - 1yr - 2yr == 0 0.103996 -0.052139 0.260130
## 3yr - 4yr - 2yr - 3yr == 0 -0.141566 -0.372541 0.089408
## Unburned - 2yr - 3yr == 0 0.058209 -0.130427 0.246846
## Unburned - 3yr - 4yr == 0 0.199776 -0.048914 0.448465Mult_AE1 <- glht(ATE1, linfct=mcp(ESD = "Tukey"))
summary(Mult_AE1) #Sandy is lower than Salo; ##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NH4_ppm) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Loamy - Clayey == 0 0.09771 0.07980 1.225 0.7131
## SaLo - Clayey == 0 -0.04102 0.06195 -0.662 0.9596
## Sandy - Clayey == 0 0.00516 0.05690 0.091 1.0000
## ThCl - Clayey == 0 0.37309 0.15237 2.449 0.0898 .
## SaLo - Loamy == 0 -0.13874 0.08112 -1.710 0.3973
## Sandy - Loamy == 0 -0.09255 0.07196 -1.286 0.6739
## ThCl - Loamy == 0 0.27538 0.15957 1.726 0.3880
## Sandy - SaLo == 0 0.04618 0.04998 0.924 0.8747
## ThCl - SaLo == 0 0.41411 0.14400 2.876 0.0281 *
## ThCl - Sandy == 0 0.36793 0.14524 2.533 0.0726 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)confint(Mult_AE1)##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lmer(formula = log(NH4_ppm) ~ TSF + ESD + (1 | Location/Month/Year),
## data = Soil181920, REML = FALSE)
##
## Quantile = 2.6739
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## Loamy - Clayey == 0 0.09771 -0.11566 0.31108
## SaLo - Clayey == 0 -0.04102 -0.20667 0.12462
## Sandy - Clayey == 0 0.00516 -0.14698 0.15729
## ThCl - Clayey == 0 0.37309 -0.03432 0.78050
## SaLo - Loamy == 0 -0.13874 -0.35565 0.07817
## Sandy - Loamy == 0 -0.09255 -0.28496 0.09986
## ThCl - Loamy == 0 0.27538 -0.15130 0.70206
## Sandy - SaLo == 0 0.04618 -0.08745 0.17982
## ThCl - SaLo == 0 0.41411 0.02908 0.79914
## ThCl - Sandy == 0 0.36793 -0.02041 0.75627