Set up
1) Load packages
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(janitor)
library(grid)
library(ggthemes)
library(tidyverse)
library(extrafont)
#font_import()
loadfonts(device = "win")
windowsFonts(Times = windowsFont("TT Times New Roman"))
library(lsmeans)
library(nlme)
library(car)
library(lme4)
library(multcomp)
library(lmerTest)
library(multcompView)
library(semPlot)
library(lavaan)
2) Import data
#import data
infiltrat <- read_excel("3-29-21 part 8 nrcs soil health.xlsx", sheet=3)
3) Data Wrangling
#change column names
infiltrat1 <- infiltrat %>%
clean_names()
#str(soilhealth)
#View(soilhealth)
infiltrat2 <- infiltrat1
infiltrat2 <- infiltrat2[c(1,5,9,13,17,21,25,29,33,37),]
infiltrat2 <- infiltrat2 %>%
dplyr::filter(location!="Ottawa")
infiltrat1 <- as.data.frame(infiltrat1)
# keeps aggregate data
inf <- infiltrat1 %>%
dplyr::select(location, treatment, reps, precip, single, cornell)
inf1 <- inf %>%
filter(location!="Ottawa", treatment!="IR")
inf1$precip <- as.factor(inf1$precip)
#vartable(inf1)
inf1$location_f =factor(inf1$location, levels=c('Tribune', 'Hays', 'Manhattan'))
4) Theme James
theme_James <- function(base_size=14, base_family="TT Times New Roman") {
(theme_foundation(base_size=base_size, base_family=base_family)+
theme_bw()+
theme(panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
axis.title = element_text(color="black",size=rel(1.2)),
axis.text = element_text(color="black", size = 12),
legend.key = element_rect(colour = NA),
legend.spacing = unit(0, "cm"),
legend.text = element_text(size=12),
legend.title = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(color="Black",size = rel(1.5),face = "bold",hjust = 0.5),
strip.text = element_text(color="Black",size = rel(1),face="bold")
))
}
Infiltration
Infiltration function
#function to calculate the mean and standard error
mean_se_fx_inf <- function(select_variable, aggregate_size, aggregate_size_conversion){
subset_data <- inf1[, select_variable]
subset_data_1 <- gather(subset_data, aggregate_size , value, - treatment, - location)
subset_data_1$aggregate_size <- factor(subset_data_1$aggregate_size,
levels = aggregate_size)
subset_data_1$aggregate_size <- plyr::mapvalues(subset_data_1$aggregate_size,
from=aggregate_size, to=aggregate_size_conversion)
se <- function(x, na.rm=TRUE) {
if (na.rm) x <- na.omit(x)
sqrt(var(x)/length(x))
}
subset_data_2 <- subset_data_1 %>% group_by(location, treatment, aggregate_size) %>%
summarize(mean_data = mean(value, na.rm = TRUE), standard_error = se(value))
#Reposition location labels
subset_data_2$location <- factor(subset_data_2$location, levels = c("Tribune", "Hays", "Manhattan"))
return(subset_data_2)
}
Infiltration
method_inf <- mean_se_fx_inf(c("treatment", "location", "single", "cornell"),
c("single", "cornell"),
c("Single-Ring", "Cornell"))
method_inf$location_f =factor(method_inf$location, levels=c('Tribune', 'Hays', 'Manhattan'))
colsnp <- c( "AG" = "black", "EA" = "grey40", "NP" = "grey90")
ggplot(data=method_inf, aes(x=aggregate_size, y=mean_data, fill = treatment)) +
geom_bar(position=position_dodge(), stat="identity", colour = "black") +
geom_errorbar(aes(ymin=mean_data-standard_error, ymax=mean_data+standard_error),
width=.2, # Width of the error bars
position=position_dodge(.9)) +
scale_y_continuous(limits=c(0,35)) +
facet_wrap(facets=vars(location_f), strip.position="bottom") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position="top",
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
strip.background=element_rect(size=0.5, colour = "black"),
axis.text.x=element_text(colour="black", size=10)) +
ylab(expression(Infiltration~(cm ~h^{-1}))) +
xlab("")+
facet_wrap(facets=vars(location_f), strip.position="bottom") +
scale_fill_manual(values = colsnp) +
theme_James()
ggsave("infil.png", height=6, width=10)
Infiltration statistics
s<-"Single Ring"
s
## [1] "Single Ring"
single <- lmer(single ~ treatment*precip + (1|reps), data=inf1, na.action=na.omit)
anova(single, type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 951.71 475.85 2 23.179 46.197 8.191e-09 ***
## precip 464.18 232.09 2 23.178 22.532 3.673e-06 ***
## treatment:precip 476.75 119.19 4 23.172 11.571 2.612e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c<- "Cornell Sprinkle Infiltrometer"
c
## [1] "Cornell Sprinkle Infiltrometer"
cornell <- lmer(cornell ~ treatment*precip + (1|reps), data=inf1, na.action=na.omit)
anova(cornell, type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 36.15 18.077 2 25 4.4263 0.02261 *
## precip 17.86 8.928 2 25 2.1862 0.13336
## treatment:precip 319.70 79.924 4 25 19.5697 2.086e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Extracting model means and pairwise comparisons
#The lsmeans part in my code is post-hoc analysis. "specs = c("Tillage"), by = c("fertilizer")" means comparing the lsmeans of target properties in treatment "tillage" by different fertilizer application (i.e. compare bG activity in NT and CT when manure was applied).
s
## [1] "Single Ring"
#means
single_means <- lsmeans(single, specs = "treatment", by = "precip")
#single_means
#PWC
single_pwc <- cld(single_means, adjust = "none", Letters = letters, reversed = T)
#transforming it in a data frame to use on ggplot
single_pwc <- as.data.frame(single_pwc)
single_pwc
## treatment precip lsmean SE df lower.CL upper.CL .group
## 1 NP 472 28.062914 1.986313 24.71643 23.9696447 32.156183 a
## 2 EA 472 6.401788 1.720460 23.01103 2.8428386 9.960737 b
## 3 AG 472 4.356709 1.720460 23.01103 0.7977601 7.915658 b
## 5 NP 579 8.509000 1.720460 23.01103 4.9500509 12.067949 a
## 4 AG 579 2.667000 1.720460 23.01103 -0.8919491 6.225949 b
## 6 EA 579 2.540000 1.720460 23.01103 -1.0189491 6.098949 b
## 7 NP 850 8.892644 1.720460 23.01103 5.3336945 12.451593 a
## 8 EA 850 6.032500 1.720460 23.01103 2.4735509 9.591449 ab
## 9 AG 850 1.843352 1.720460 23.01103 -1.7155968 5.402301 b
#Determining the real SE
real_se_single <- inf1 %>%
na.omit() %>%
group_by(precip, treatment) %>%
summarise(
n=n(),
mean=mean(single),
sd=sd(single)
) %>%
mutate( se=sd/sqrt(n))
c
## [1] "Cornell Sprinkle Infiltrometer"
#means
cornell_means <- lsmeans(cornell, specs = "treatment", by = "precip")
#cornell_means
#PWC
cornell_pwc <- cld(cornell_means, adjust = "none", Letters = letters, reversed = T)
cornell_pwc
## precip = 472:
## treatment lsmean SE df lower.CL upper.CL .group
## NP 22.2 1.01 25 20.1 24.2 a
## AG 18.5 1.01 25 16.4 20.5 b
## EA 18.0 1.19 25 15.6 20.4 b
##
## precip = 579:
## treatment lsmean SE df lower.CL upper.CL .group
## AG 21.7 1.01 25 19.7 23.8 a
## NP 20.5 1.19 25 18.0 22.9 a
## EA 14.0 1.01 25 11.9 16.1 b
##
## precip = 850:
## treatment lsmean SE df lower.CL upper.CL .group
## EA 22.5 1.01 25 20.5 24.6 a
## NP 17.8 1.01 25 15.7 19.8 b
## AG 13.0 1.01 25 10.9 15.1 c
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#transforming it in a data frame to use on ggplot
cornell_pwc <- as.data.frame(cornell_pwc)
#cornell_pwc
#Determining the real SE
real_se_cornell <- inf1 %>%
na.omit() %>%
group_by(precip, treatment) %>%
summarise(
n=n(),
mean=mean(cornell),
sd=sd(cornell)
) %>%
mutate( se=sd/sqrt(n))
#totinf <- merge(single_pwc, cornell_pwc, by=c("treatment", "precip"))
#totinf
single_pwc$smean <- single_pwc$lsmean
mean_inf<- merge(inf1, single_pwc, by=c("treatment", "precip"))
cornell_pwc$cmean <- cornell_pwc$lsmean
mean_inf2<- merge(mean_inf, cornell_pwc, by=c("treatment", "precip"))
#mean_inf2
Single Ring Plotting
#add locations
df <- data.frame (location = c("Tribune", "Hays", "Manhattan"),
precip = c("472", "579", "850"))
df<- merge(df, real_se_single, by=c("precip"))
single_pwc_1 <- merge(single_pwc, df, by=c("precip", "treatment"))
single_pwc_1 <- as.data.frame(single_pwc_1)
single_pwc_1$location_f =factor(single_pwc_1$location, levels=c('Tribune', 'Hays', 'Manhattan'))
colinf <- c( "AG" = "#ff6500", "EA" = "deepskyblue1", "NP" = "green4")
ggplot(data=single_pwc_1, aes(x=treatment, y=lsmean, fill = treatment)) +
geom_bar(position=position_dodge(), stat="identity", colour = "black") +
geom_errorbar(aes(ymin=lsmean-se, ymax=lsmean+se),
width=.2, # Width of the error bars
position=position_dodge(.9)) +
scale_y_continuous(limits=c(0,35)) +
facet_wrap(facets=vars(location_f), strip.position="bottom") +
xlab("")+
facet_wrap(facets=vars(location_f), strip.position="bottom") +
scale_fill_manual(values = colinf) +
theme_James() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
strip.background=element_rect(size=0.5, colour = "black"),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
legend.position= c(0.92,0.8)) +
ylab(expression(~Single-Ring ~Infiltration~(cm ~h^{-1})))+
geom_label(aes(label=trimws(.group), y = lsmean+se+3),
label.padding = unit(.3,"lines"), show.legend=NA, label.size = NA, fill=NA, font="bold") +
guides(fill = guide_legend(override.aes = aes(label="")))
# ggsave("infil.png", height=6, width=10)
Cornell Sprinkle Infiltrometer Plotting
#add locations
df <- data.frame (location = c("Tribune", "Hays", "Manhattan"),
precip = c("472", "579", "850"))
df<- merge(df, real_se_cornell, by=c("precip"))
cornell_pwc_1 <- merge(cornell_pwc, df, by=c("precip", "treatment"))
cornell_pwc_1 <- as.data.frame(cornell_pwc_1)
cornell_pwc_1$location_f =factor(cornell_pwc_1$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ggplot(data=cornell_pwc_1, aes(x=treatment, y=lsmean, fill = treatment)) +
geom_bar(position=position_dodge(), stat="identity", colour = "black") +
geom_errorbar(aes(ymin=lsmean-se, ymax=lsmean+se),
width=.2, # Width of the error bars
position=position_dodge(.9)) +
scale_y_continuous(limits=c(0,35)) +
facet_wrap(facets=vars(location_f), strip.position="bottom") +
xlab("")+
facet_wrap(facets=vars(location_f), strip.position="bottom") +
scale_fill_manual(values = colinf) +
theme_James() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
strip.background=element_rect(size=0.5, colour = "black"),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
legend.position= c(0.92,0.87)) +
ylab(expression(~Cornell ~Infiltration~(cm ~h^{-1})))+
geom_label(aes(label=trimws(.group), y = lsmean+se+2),
label.padding = unit(.3,"lines"), show.legend=NA, label.size = NA, fill=NA, font="bold" ) +
guides(fill = guide_legend(override.aes = aes(label="")))
# ggsave("infil.png", height=6, width=10)
Combined graph
#mean_inf
single_pwc1<- single_pwc %>%
dplyr::select( treatment, precip, .group) %>%
dplyr::mutate(aggregate_size= "Single-Ring")
#add locations
df <- data.frame (location = c("Tribune", "Hays", "Manhattan"),
precip = c("472", "579", "850"))
singlef <- merge(df, single_pwc1, by=c("precip"))
cornell_pwc1 <- cornell_pwc %>%
dplyr::select( treatment, precip, .group) %>%
dplyr::mutate(aggregate_size= "Cornell")
#add locations
df <- data.frame (location = c("Tribune", "Hays", "Manhattan"),
precip = c("472", "579", "850"))
cornellf <- merge(df, cornell_pwc1, by=c("precip"))
inf <- rbind(singlef, cornellf)
infilf <- merge(method_inf, inf, by= c("location", "treatment", "aggregate_size"))
Final graph
ggplot(data=infilf, aes(x=aggregate_size, y=mean_data, fill = treatment)) +
geom_bar(position=position_dodge(), stat="identity", colour = "black") +
geom_errorbar(aes(ymin=mean_data-standard_error, ymax=mean_data+standard_error),
width=.2, # Width of the error bars
position=position_dodge(.9)) +
scale_y_continuous(limits=c(0,35)) +
facet_wrap(facets=vars(location_f), strip.position="bottom") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.position="top",
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
strip.background=element_rect(size=0.5, colour = "black"),
axis.text.x=element_text(colour="black", size=10)) +
ylab(expression(Infiltration~(cm ~h^{-1}))) +
xlab("")+
facet_wrap(facets=vars(location_f), strip.position="bottom") +
scale_fill_manual(values = colsnp) +
theme_James() +
geom_text(data = infilf,
aes(x = aggregate_size, group=treatment, y = mean_data+standard_error +1,
label = format(.group, nsmall = 0, digits=1, scientific = FALSE)),
color="black", position=position_dodge(.9), hjust="center", inherit.aes = TRUE)
Right post hoc test
Single Ring Infiltrometer
s<-"Single Ring"
s
## [1] "Single Ring"
single <- lmer(single ~ treatment*precip + (1|reps), data=inf1, na.action=na.omit)
anova(single, type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 951.71 475.85 2 23.179 46.197 8.191e-09 ***
## precip 464.18 232.09 2 23.178 22.532 3.673e-06 ***
## treatment:precip 476.75 119.19 4 23.172 11.571 2.612e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#means
single_means_int <- lsmeans(single, ~treatment*precip, adjust="tukey")
#single_means
#PWC
single_pwc_int <- cld(single_means_int, adjust = "none", Letters = letters, reversed = T)
single_pwc_int
## treatment precip lsmean SE df lower.CL upper.CL .group
## NP 472 28.06 1.99 24.7 23.970 32.16 a
## NP 850 8.89 1.72 23.0 5.334 12.45 b
## NP 579 8.51 1.72 23.0 4.950 12.07 b
## EA 472 6.40 1.72 23.0 2.843 9.96 bc
## EA 850 6.03 1.72 23.0 2.474 9.59 bc
## AG 472 4.36 1.72 23.0 0.798 7.92 bc
## AG 579 2.67 1.72 23.0 -0.892 6.23 c
## EA 579 2.54 1.72 23.0 -1.019 6.10 c
## AG 850 1.84 1.72 23.0 -1.716 5.40 c
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#transforming it in a data frame to use on ggplot
single_pwc_int <- as.data.frame(single_pwc)
#single_pwc
#Determining the real SE
real_se_single <- inf1 %>%
na.omit() %>%
group_by(precip, treatment) %>%
summarise(
n=n(),
mean=mean(single),
sd=sd(single)
) %>%
mutate( se=sd/sqrt(n))
#add locations
df <- data.frame (location = c("Tribune", "Hays", "Manhattan"),
precip = c("472", "579", "850"))
df<- merge(df, real_se_single, by=c("precip"))
single_pwc_2 <- merge(single_pwc_int, df, by=c("precip", "treatment"))
single_pwc_2 <- as.data.frame(single_pwc_2)
single_pwc_2$location_f =factor(single_pwc_2$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ggplot(data=single_pwc_2, aes(x=treatment, y=lsmean, fill = treatment)) +
geom_bar(position=position_dodge(), stat="identity", colour = "black") +
geom_errorbar(aes(ymin=lsmean-se, ymax=lsmean+se),
width=.2, # Width of the error bars
position=position_dodge(.9)) +
scale_y_continuous(limits=c(0,35)) +
facet_wrap(facets=vars(location_f), strip.position="bottom") +
xlab("")+
facet_wrap(facets=vars(location_f), strip.position="bottom") +
scale_fill_manual(values = colinf) +
theme_James() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
strip.background=element_rect(size=0.5, colour = "black"),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
legend.position= c(0.92,0.8)) +
ylab(expression(~Single-Ring ~Infiltration~(cm ~h^{-1})))+
geom_label(aes(label=trimws(.group), y = lsmean+se+2),
label.padding = unit(.3,"lines"), show.legend=NA,label.size = NA, fill=NA, font="bold" ) +
guides(fill = guide_legend(override.aes = aes(label="")))
Cornell Sprinkle Infiltrometer
c<- "Cornell Sprinkle Infiltrometer"
c
## [1] "Cornell Sprinkle Infiltrometer"
cornell <- lmer(cornell ~ treatment*precip + (1|reps), data=inf1, na.action=na.omit)
anova(cornell, type=3)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 36.15 18.077 2 25 4.4263 0.02261 *
## precip 17.86 8.928 2 25 2.1862 0.13336
## treatment:precip 319.70 79.924 4 25 19.5697 2.086e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#means
cornell_means_int <- lsmeans(cornell, ~treatment*precip, adjust="tukey")
#cornell_means
#PWC
cornell_pwc_int <- cld(cornell_means_int, adjust = "none", Letters = letters, reversed = T)
cornell_pwc_int
## treatment precip lsmean SE df lower.CL upper.CL .group
## EA 850 22.5 1.01 25 20.5 24.6 a
## NP 472 22.2 1.01 25 20.1 24.2 a
## AG 579 21.7 1.01 25 19.7 23.8 a
## NP 579 20.5 1.19 25 18.0 22.9 ab
## AG 472 18.5 1.01 25 16.4 20.5 b
## EA 472 18.0 1.19 25 15.6 20.4 b
## NP 850 17.8 1.01 25 15.7 19.8 b
## EA 579 14.0 1.01 25 11.9 16.1 c
## AG 850 13.0 1.01 25 10.9 15.1 c
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#transforming it in a data frame to use on ggplot
cornell_pwc_int <- as.data.frame(cornell_pwc_int)
#Determining the real SE
real_se_cornell <- inf1 %>%
na.omit() %>%
group_by(precip, treatment) %>%
summarise(
n=n(),
mean=mean(cornell),
sd=sd(cornell)
) %>%
mutate( se=sd/sqrt(n))
#add locations
df <- data.frame (location = c("Tribune", "Hays", "Manhattan"),
precip = c("472", "579", "850"))
df<- merge(df, real_se_cornell, by=c("precip"))
cornell_pwc_2 <- merge(cornell_pwc_int, df, by=c("precip", "treatment"))
cornell_pwc_2 <- as.data.frame(cornell_pwc_2)
cornell_pwc_2$location_f =factor(cornell_pwc_2$location, levels=c('Tribune', 'Hays', 'Manhattan'))
ggplot(data=cornell_pwc_2, aes(x=treatment, y=lsmean, fill = treatment)) +
geom_bar(position=position_dodge(), stat="identity", colour = "black") +
geom_errorbar(aes(ymin=lsmean-se, ymax=lsmean+se),
width=.2, # Width of the error bars
position=position_dodge(.9)) +
scale_y_continuous(limits=c(0,35)) +
facet_wrap(facets=vars(location_f), strip.position="bottom") +
xlab("")+
facet_wrap(facets=vars(location_f), strip.position="bottom") +
scale_fill_manual(values = colinf) +
theme_James() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=0.5),
strip.background=element_rect(size=0.5, colour = "black"),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
legend.position= c(0.92,0.87)) +
ylab(expression(~Cornell ~Infiltration~(cm ~h^{-1})))+
geom_label(aes(label=trimws(.group), y = lsmean+se+2),
label.padding = unit(.3,"lines"), show.legend=NA, label.size = NA, fill=NA, font="bold") +
guides(fill = guide_legend(override.aes = aes(label="")))
