code for comparing RAP estimates to field cover estimates of plant functional group cover
load in packages
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(vegan)## Warning: package 'vegan' was built under R version 4.0.5## Loading required package: permute## Warning: package 'permute' was built under R version 4.0.5## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.0.5## This is vegan 2.5-7library(RColorBrewer) # for graph color palettes
library(PairedData)## Warning: package 'PairedData' was built under R version 4.0.4## Loading required package: MASS## Warning: package 'MASS' was built under R version 4.0.5##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: gld## Warning: package 'gld' was built under R version 4.0.5## Loading required package: mvtnorm## Warning: package 'mvtnorm' was built under R version 4.0.5##
## Attaching package: 'PairedData'## The following object is masked from 'package:base':
##
## summarylibrary(ggpubr)## Warning: package 'ggpubr' was built under R version 4.0.5library(htmltools) # for RMD## Warning: package 'htmltools' was built under R version 4.0.5#install.packages('htmltools')import files
setwd("~/Research/Summer2021")
plant_traits <-read.csv("PlantTraits.csv", header = T)
#species list with plant habit
#code is USDA plant species code
#habits are PGR, AGR, AFO, PFO, and PSH
#FXGroup is combination of life form, natiity, longevity, size and regen strategy
#Count data in wide format
foliar_cover <- read.csv("FoliarCover_PCover.csv", header = T)
#PlotID and transect by species matrix
#within matrix: count of hits for each transect
plotmeta <- read.csv("Plot_meta.csv", header = T)
#treatment information for all the plotsadd in plot metadata for each column
head(plotmeta,2)## PlotID Num Trt SH_GR pasture old_point old_name
## 1 101T 101 T HSHG LRC NA HSHG_LRC_T
## 2 102T 102 T HSHG KC 31 HSHG_T_7head(foliar_cover, 1)## PLOTID TRANSECT ACMI2 ACNE9 ACTH7 AGGL AGGR ALAC4 AMMEM2 AMRE2 AMSIN ANLU2
## 1 101C 2 0 0 11.66667 0 0 0 0 0 0 0
## ARAR8 ARLU ARTRV ASTRA BASA3 BRJA BRTE CADO2 CAEU CASTI2 CHDO CHVI8
## 1 0 0 63.33333 0 0 0 1.666667 0 0 0 0 0
## CIUN COLI2 COPA3 CORA5 CRAC2 CRAT CROC DF DG DS ELEL5 EQLA ERAP
## 1 1.666667 0 0 0 0 0 0 0 0 1.666667 0 0 0
## ERCA8 ERHE2 ERNA10 EROV ERPU2 ERUM ERWR FEID GARA2 HEUN HYCA4 IPAGA3 IRIS
## 1 0 0 0 0 0 0 0 0 0 0 0 3.333333 0
## IVAX KOMA L LASE LECI4 LERE7 LILO LIRU4 LOCO4 LODI LONU2 LOTR LUAR6 LUCA
## 1 0 0 0 0 0 0 0 0 0 0 0 0 100 0
## LUPIN LUSE4 MADIA MAGR3 MEAL6 MEBU MINU N NEST5 PEHU PENST PHHE2
## 1 0 0 3.333333 0 0 0 0 18.33333 0 0 0 0
## PHLO2 POBU POSE PSSP6 PUTR2 SEIN2 TECA2 TRDU UNK100 VIPUV2 WL ZIPA2
## 1 0 0 13.33333 0 0 0 0 0 0 0 0 0
## ZIPA2.1
## 1 0try1 = merge(foliar_cover, plotmeta, by.x = "PLOTID", by.y="PlotID")
head(try1,1)## PLOTID TRANSECT ACMI2 ACNE9 ACTH7 AGGL AGGR ALAC4 AMMEM2 AMRE2 AMSIN ANLU2
## 1 101C 2 0 0 11.66667 0 0 0 0 0 0 0
## ARAR8 ARLU ARTRV ASTRA BASA3 BRJA BRTE CADO2 CAEU CASTI2 CHDO CHVI8
## 1 0 0 63.33333 0 0 0 1.666667 0 0 0 0 0
## CIUN COLI2 COPA3 CORA5 CRAC2 CRAT CROC DF DG DS ELEL5 EQLA ERAP
## 1 1.666667 0 0 0 0 0 0 0 0 1.666667 0 0 0
## ERCA8 ERHE2 ERNA10 EROV ERPU2 ERUM ERWR FEID GARA2 HEUN HYCA4 IPAGA3 IRIS
## 1 0 0 0 0 0 0 0 0 0 0 0 3.333333 0
## IVAX KOMA L LASE LECI4 LERE7 LILO LIRU4 LOCO4 LODI LONU2 LOTR LUAR6 LUCA
## 1 0 0 0 0 0 0 0 0 0 0 0 0 100 0
## LUPIN LUSE4 MADIA MAGR3 MEAL6 MEBU MINU N NEST5 PEHU PENST PHHE2
## 1 0 0 3.333333 0 0 0 0 18.33333 0 0 0 0
## PHLO2 POBU POSE PSSP6 PUTR2 SEIN2 TECA2 TRDU UNK100 VIPUV2 WL ZIPA2
## 1 0 0 13.33333 0 0 0 0 0 0 0 0 0
## ZIPA2.1 Num Trt SH_GR pasture old_point old_name
## 1 0 101 C HSHG LRC NA HSHG_LRC_Cwhich(colnames(try1)=="old_name") #8## [1] 95#reorder columns so that metadata is first, species are second
foliar_cover <- try1[, c(1:2, 90:95, 3:89)]create a species only matrix for ALL SPECIES! Also create an output dataset to plop all of our output nalues into.
head(foliar_cover, 1)## PLOTID TRANSECT Num Trt SH_GR pasture old_point old_name ACMI2 ACNE9
## 1 101C 2 101 C HSHG LRC NA HSHG_LRC_C 0 0
## ACTH7 AGGL AGGR ALAC4 AMMEM2 AMRE2 AMSIN ANLU2 ARAR8 ARLU ARTRV ASTRA
## 1 11.66667 0 0 0 0 0 0 0 0 0 63.33333 0
## BASA3 BRJA BRTE CADO2 CAEU CASTI2 CHDO CHVI8 CIUN COLI2 COPA3 CORA5
## 1 0 0 1.666667 0 0 0 0 0 1.666667 0 0 0
## CRAC2 CRAT CROC DF DG DS ELEL5 EQLA ERAP ERCA8 ERHE2 ERNA10 EROV ERPU2
## 1 0 0 0 0 0 1.666667 0 0 0 0 0 0 0 0
## ERUM ERWR FEID GARA2 HEUN HYCA4 IPAGA3 IRIS IVAX KOMA L LASE LECI4 LERE7
## 1 0 0 0 0 0 0 3.333333 0 0 0 0 0 0 0
## LILO LIRU4 LOCO4 LODI LONU2 LOTR LUAR6 LUCA LUPIN LUSE4 MADIA MAGR3 MEAL6
## 1 0 0 0 0 0 0 100 0 0 0 3.333333 0 0
## MEBU MINU N NEST5 PEHU PENST PHHE2 PHLO2 POBU POSE PSSP6 PUTR2
## 1 0 0 18.33333 0 0 0 0 0 0 13.33333 0 0
## SEIN2 TECA2 TRDU UNK100 VIPUV2 WL ZIPA2 ZIPA2.1
## 1 0 0 0 0 0 0 0 0meta_cols = c("PLOTID", "TRANSECT", "Num", "Trt", "SH_GR", "pasture", "old_name") #these are the plot metadata columns we would like to retain through the calculations
which(colnames(foliar_cover)=="ACMI2") #8## [1] 9which(colnames(foliar_cover)=="ZIPA2") #93## [1] 94total.abundance.matrix <- foliar_cover[, 9:94]
#total.abundance.matrix <- lapply(total.abundance.matrix, as.numeric)Also create an output dataset to plop all of our output nalues into.
i_wide <- foliar_cover[,meta_cols]Generate estimates of plant functional group cover
Annual forbs and grasses
codes_AF <- plant_traits %>%
filter(FXGroup == "FO.NAT.A"|
FXGroup == "FO.INV.A"|
FXGroup == "GR.INV.A")
AF_sp_list <- codes_AF$code # a string containing al the species
A.abundance.matrix <- dplyr::select(foliar_cover, any_of(AF_sp_list))
i_wide$AFG <- rowSums(A.abundance.matrix) #this is summing all of the cover values for each row, remember that each row is a plotPerennial forbs and grasses
codes_PFG <- plant_traits %>%
filter(FXGroup == "FO.NAT.P"|
FXGroup == "FO.Unk.P"|
FXGroup == "GR.NAT.P.LG."|
FXGroup == "GR.INV.P.SM." |
FXGroup == "GR.NAT.P.SM.")
PF_sp_list <- codes_PFG$code # a string containing al the PF species
PF.abundance.matrix <- dplyr::select(foliar_cover, any_of(PF_sp_list))
i_wide$PFG <- rowSums(PF.abundance.matrix) #this is summing all of the cover values for each row, remember that each row is a plotShrubs
codes_SH <- plant_traits %>%
filter(FXGroup == "SH.NAT.P.SE"|
FXGroup == "SH.NAT.P.SP" )
SH_sp_list <- codes_SH$code # a string containing al the PF species
SH.abundance.matrix <- dplyr::select(foliar_cover, any_of(SH_sp_list))
i_wide$SHR <- rowSums(SH.abundance.matrix) #this is summing all of the cover values for each row, remember that each row is a plotCreate an output dataset with means by plot
#first we need to aggregate field cover by plot
field_cover <- i_wide %>%
group_by(PLOTID)%>%
summarize(AFG =mean(AFG),
PFG =mean(PFG),
SHR =mean(SHR))
field_cover$Method = "FIELD"need to go from wide to long format
columns_rich = c("AFG", "PFG", "SHR")
field_long = pivot_longer(field_cover, columns_rich,
names_to = "FXGroup",
values_to = "cover")## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(columns_rich)` instead of `columns_rich` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.head(field_long)## # A tibble: 6 x 4
## PLOTID Method FXGroup cover
## <chr> <chr> <chr> <dbl>
## 1 101C FIELD AFG 4.44
## 2 101C FIELD PFG 119.
## 3 101C FIELD SHR 50.6
## 4 101T FIELD AFG 1.67
## 5 101T FIELD PFG 124.
## 6 101T FIELD SHR 43.3Bring in RAP data
RAP <- read.csv("RAP_Extract_30m.csv", header = T)
head(RAP,2)## system.index old_name new_name
## 1 2021_00000000000000000034 HSHG_LRC_C 101C
## 2 2021_00000000000000000035 HSHG_LRC_T 101T
## GlobalID ImageID timestamp year AFG
## 1 40fee60c-03e1-41cc-b52d-80efc314a226 2021 1.61e+12 2021 13.398359
## 2 25c31297-ec53-4788-b6c0-7459704f1173 2021 1.61e+12 2021 8.420237
## BGR LTR PFG SHR TRE
## 1 7.816773 18.57703 35.62261 22.38377 0
## 2 11.072015 17.19234 28.77758 30.99271 0RAP$Method = "RAP"
#need to get new plot name in the RAP data
RAP_withMeta = merge(RAP, plotmeta, by.x = "old_name", by.y="old_name")
head(RAP_withMeta)## old_name system.index new_name
## 1 HSHG_C_1 2021_00000000000000000002 105T
## 2 HSHG_C_5 2021_00000000000000000022 103C
## 3 HSHG_C_6 2021_00000000000000000023 104C
## 4 HSHG_C_8 2021_00000000000000000025 102C
## 5 HSHG_LRC_C 2021_00000000000000000034 101C
## 6 HSHG_LRC_T 2021_00000000000000000035 101T
## GlobalID ImageID timestamp year AFG
## 1 f1de3649-a702-4d76-bf3e-1dd8b1cc8156 2021 1.61e+12 2021 10.484047
## 2 e2202b2a-425d-4254-be3b-6ae17d0569ff 2021 1.61e+12 2021 10.819836
## 3 c49bf404-f9e1-4027-a9ea-829328976a35 2021 1.61e+12 2021 9.230419
## 4 e7bce277-baf4-44c0-983b-a00689165519 2021 1.61e+12 2021 16.653881
## 5 40fee60c-03e1-41cc-b52d-80efc314a226 2021 1.61e+12 2021 13.398359
## 6 25c31297-ec53-4788-b6c0-7459704f1173 2021 1.61e+12 2021 8.420237
## BGR LTR PFG SHR TRE Method PlotID Num Trt SH_GR pasture
## 1 6.909754 11.73291 33.38560 34.31905 0 RAP 105T 105 T HSHG HG
## 2 14.453139 19.41310 23.21656 28.28571 0 RAP 103C 103 C HSHG KC
## 3 11.810565 18.80237 28.59563 27.59199 0 RAP 104C 104 C HSHG KC
## 4 7.000000 16.23744 31.52785 24.29041 0 RAP 102C 102 C HSHG KC
## 5 7.816773 18.57703 35.62261 22.38377 0 RAP 101C 101 C HSHG LRC
## 6 11.072015 17.19234 28.77758 30.99271 0 RAP 101T 101 T HSHG LRC
## old_point
## 1 3
## 2 35
## 3 36
## 4 38
## 5 NA
## 6 NAcolumns_rich = c("AFG", "PFG", "SHR")
RAP_long = pivot_longer(RAP_withMeta, columns_rich,
names_to = "FXGroup",
values_to = "cover")
RAP_long <- RAP_long %>% dplyr::select(Method, PlotID, FXGroup, cover)
head(RAP_long,3)## # A tibble: 3 x 4
## Method PlotID FXGroup cover
## <chr> <chr> <chr> <dbl>
## 1 RAP 105T AFG 10.5
## 2 RAP 105T PFG 33.4
## 3 RAP 105T SHR 34.3RAP_long$PLOTID <- RAP_long$PlotID
head(field_long)## # A tibble: 6 x 4
## PLOTID Method FXGroup cover
## <chr> <chr> <chr> <dbl>
## 1 101C FIELD AFG 4.44
## 2 101C FIELD PFG 119.
## 3 101C FIELD SHR 50.6
## 4 101T FIELD AFG 1.67
## 5 101T FIELD PFG 124.
## 6 101T FIELD SHR 43.3some data cleaning and merge the field estimates and RAP
field_long$FieldCover <-
field_long$cover
RAP_long$RAPCover <-
RAP_long$cover
comp = dplyr::full_join(field_long, RAP_long, by = c("PLOTID", "FXGroup"))
head(comp,3)## # A tibble: 3 x 9
## PLOTID Method.x FXGroup cover.x FieldCover Method.y PlotID cover.y RAPCover
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <dbl> <dbl>
## 1 101C FIELD AFG 4.44 4.44 RAP 101C 13.4 13.4
## 2 101C FIELD PFG 119. 119. RAP 101C 35.6 35.6
## 3 101C FIELD SHR 50.6 50.6 RAP 101C 22.4 22.4comp <- na.omit(comp)
cover <- comp %>% dplyr::select(PLOTID, FXGroup, FieldCover, RAPCover)
cover_long <- pivot_longer(cover, c("FieldCover", "RAPCover"),
names_to = "method",
values_to = "cover")
head(cover_long)## # A tibble: 6 x 4
## PLOTID FXGroup method cover
## <chr> <chr> <chr> <dbl>
## 1 101C AFG FieldCover 4.44
## 2 101C AFG RAPCover 13.4
## 3 101C PFG FieldCover 119.
## 4 101C PFG RAPCover 35.6
## 5 101C SHR FieldCover 50.6
## 6 101C SHR RAPCover 22.4cover_long <- na.omit(cover_long)data viz
View all functional groups on one plto
p <- ggplot(cover_long, aes(x=FXGroup, y=cover, fill = method))+
geom_boxplot()+
scale_fill_brewer(palette="Dark2")+
xlab("Plant functional group") + ylab("Cover") +
labs(fill = "Method")+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.text = element_text(size=11, colour = "black"),
axis.line = element_line(colour = "black"), axis.title = element_text(size=16, colour = "black"),
legend.title = element_text(size = 16, colour = "black"), legend.text = element_text(size=16, colour= "black"))+
theme_classic()
plot(p)
Annual forbs and grasses
afg <- dplyr::filter(cover_long, FXGroup == "AFG")
a <- ggplot(afg, aes(x=method, y=cover, fill = method))+
geom_boxplot()+
scale_fill_brewer(palette="Dark2")+ ylab("Cover") +
scale_x_discrete("Method", labels = c("FieldCover" = "Field","RAPCover" = "RAP"))+ theme_classic()+
theme(legend.position = "none",
axis.text = element_text(size=16, colour = "black"),
axis.title = element_text(size=16, colour = "black"))+
ggtitle("Annual forbs and grasses")+
geom_point(alpha = 0.5)+
geom_line(aes(group = PLOTID), color = 'azure4')+
stat_compare_means(paired = TRUE, method = "t.test", label.x =1.5, label.y=47 )
# geom_dotplot(binaxis='y', stackdir='center', position = position_dodge(.8),
# dotsize = .5, fill = white)
plot(a)
head(afg)## # A tibble: 6 x 4
## PLOTID FXGroup method cover
## <chr> <chr> <chr> <dbl>
## 1 101C AFG FieldCover 4.44
## 2 101C AFG RAPCover 13.4
## 3 101T AFG FieldCover 1.67
## 4 101T AFG RAPCover 8.42
## 5 102T AFG FieldCover 5.56
## 6 102T AFG RAPCover 6.63res <- t.test(cover ~ method, data = afg, paired = TRUE)
res##
## Paired t-test
##
## data: cover by method
## t = -1.1991, df = 28, p-value = 0.2405
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.362152 1.402334
## sample estimates:
## mean of the differences
## -1.979909Perennial forbs and grasses
pfg <- dplyr::filter(cover_long, FXGroup == "PFG")
p <- ggplot(pfg, aes(x=method, y=cover, fill = method))+
geom_boxplot()+
scale_fill_brewer(palette="Dark2")+ ylab("Cover") +
scale_x_discrete("Method", labels = c("FieldCover" = "Field","RAPCover" = "RAP"))+ theme_classic()+
theme(legend.position = "none",
axis.text = element_text(size=16, colour = "black"),
axis.title = element_text(size=16, colour = "black"))+
geom_point(alpha = 0.5)+
ggtitle("Perennial forbs and grasses")+
geom_line(aes(group = PLOTID), color = 'azure4')+
stat_compare_means(paired = TRUE, method = "t.test", label.x =1.5, label.y=200 )
# geom_dotplot(binaxis='y', stackdir='center', position = position_dodge(.8),
# dotsize = .5, fill = white)
plot(p)
res <- t.test(cover ~ method, data = pfg, paired = TRUE)
res##
## Paired t-test
##
## data: cover by method
## t = 27.602, df = 28, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 101.1565 117.3743
## sample estimates:
## mean of the differences
## 109.2654Shrubs
shr <- dplyr::filter(cover_long, FXGroup == "SHR")
s <- ggplot(shr, aes(x=method, y=cover, fill = method))+
geom_boxplot()+
scale_fill_brewer(palette="Dark2")+ ylab("Cover") +
scale_x_discrete("Method", labels = c("FieldCover" = "Field","RAPCover" = "RAP"))+ theme_classic()+
theme(legend.position = "none",
axis.text = element_text(size=16, colour = "black"),
axis.title = element_text(size=16, colour = "black"))+
ggtitle("Shrubs")+
geom_point(alpha = 0.5)+
geom_line(aes(group = PLOTID), color = 'azure4') +
stat_compare_means(paired = TRUE, method = "t.test", label.x =1.5, label.y=70 )
# geom_dotplot(binaxis='y', stackdir='center', position = position_dodge(.8),
# dotsize = .5, fill = white)
plot(s)
res <- t.test(cover ~ method, data = shr, paired = TRUE)
res##
## Paired t-test
##
## data: cover by method
## t = 4.0375, df = 28, p-value = 0.00038
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 6.64688 20.33703
## sample estimates:
## mean of the differences
## 13.49196Merge all the taxa groups together onto one figure
all_taxa_summary <- ggarrange(a, p, s, ncol=3)
all_taxa_summary 