Load Libraries and Data
## Load libraries
library(tidyverse)
library(ggthemes)
## Load functions
source("scripts//figureFormat.r")
## Load data on study information
studyData <- read.csv("data//StudyData.csv", stringsAsFactors = F)
## load all data on grazing and merge with simplified datasets
data <- read.csv("data//grazingData.csv", stringsAsFactors = F)
est <- read.csv("data//simpleCategories//estimateCategories.csv", stringsAsFactors = F)
graze <- read.csv("data//simpleCategories//grazingCategories.csv", stringsAsFactors = F)
data <- merge(data, est, by="Estimate", stringsAsFactors = F)
data <- merge(data, graze, by="grazing.estimate", stringsAsFactors = F)
A distribution of where the studies took place
## Remove sites without lat lon
gps <- studyData %>% filter(Latitude != "")
### Split multiple GPS coordinates for the same site
gpsSplit <- gps %>% transform(Longitude = strsplit(Longitude, ";")) %>% transform(Latitude = strsplit(Latitude, ";")) %>% unnest(Longitude, Latitude) %>% ## unnest lat/lon
mutate(Latitude = as.numeric(Latitude), Longitude=as.numeric(Longitude))
require(ggmap)
require(maps)
### Start with base map of world
mp <- NULL
mapWorld <- borders("world", colour="white", fill="gray75") # create a layer of borders
mp <- ggplot() + theme_classic()+ mapWorld + xlim(-180,180)
mp <- mp+ geom_point(data=gpsSplit %>% filter(Grazer.status != "") , aes(x=Longitude, y=Latitude, fill=Grazer.status), size=3, pch=21) + scale_fill_manual(values=c("#D55E00", "#009E73", "#56B4E9","grey50")) + ylab("Latitude") + xlab("Longitude")
mp
Visualize patterns among studies
## Replicates
names(studyData)
## [1] "UniqueID" "Title"
## [3] "Authors" "Source.title"
## [5] "Publication.year" "Grazer.species"
## [7] "Grazer.status" "Last.grazing.event"
## [9] "Study.duration" "Latitude"
## [11] "Longitude" "Elevation"
## [13] "Country" "N.sites"
## [15] "Survey.technique" "Wild.grazer.species"
## [17] "Domestic.grazer.species" "Plant.community"
## [19] "Ecosystem.class" "Fenced"
## [21] "Tilled" "Herbivores"
## [23] "Fertilization" "Fire"
## [25] "Ownership" "Year.initiated"
## [27] "Year.finished"
plot1 <- ggplot(studyData, aes(x = as.numeric(N.sites))) + geom_histogram(fill="Grey50", bins=10) + theme_Publication() + xlab("Number of sites") + ylab("Frequency of studies") + scale_x_continuous(trans = "log2")
plot2 <- ggplot(studyData, aes(x = as.numeric(Study.duration))) + geom_histogram(fill="Grey50",bins=10) + theme_Publication() + xlab("Duration of study in years") + ylab("Frequency of studies")+ scale_x_continuous(trans = "log2")
gridExtra::grid.arrange(plot1, plot2, ncol=2)
# studyData$n.sites <- as.numeric(studyData$n.sites)
#
# nstudies <- function(x){
# nfilt <- studyData %>% filter(study.duration > x) %>% nrow()
# return(as.vector(nfilt))
# }
#
# studiesSites <- unlist(sapply(1:30, nstudies))
# studiesYears <- unlist(sapply(1:30, nstudies))
# studiesList <- data.frame(rep=c(1:30,1:30), type= c(rep("years",30),rep("sites",30)), nstudies = c(studiesYears, studiesSites))
#
# ggplot(studiesList, aes(x=rep, y=nstudies, color=type)) +geom_point()
Pattern among Vegetation Classes
## Split by semi-colon
unnestVeg <- studyData %>% unnest(Ecosystem.class = strsplit(Ecosystem.class, ";"))
## Summarize by each class
VegClass <- unnestVeg %>% group_by(Ecosystem.class) %>% summarize(n=length(Ecosystem.class)) %>% filter(!is.na(Ecosystem.class))
## Capitalize the names of each
VegClass[,"capitalVegClass"] <- str_to_title(VegClass$Ecosystem.class)
ggplot(VegClass, aes(x=capitalVegClass, y=n)) + geom_bar(stat="identity") + theme_Publication() + xlab("") + ylab("Number of studies") + coord_flip()
Most common grazers
## Split by semi-colon
unnestGraze <- studyData %>% unnest(Herbivores = strsplit(Herbivores, ";"))
## Summarize by each class
herbsum <- unnestGraze %>% group_by(Herbivores) %>% summarize(n=length(Herbivores)) %>% filter(!is.na(Herbivores))
## Simplified name
herbSimple <- read.csv("data//simpleCategories//herbsum.csv") %>% select(-n)
herbsum <- merge(herbsum, herbSimple) %>% group_by(simpleHerbivoreName) %>% summarize(n=sum(n))
topherb <- herbsum %>% top_n(20) %>% arrange(-n) ## select 10 most commonly surveyed
topherb
## # A tibble: 21 x 2
## simpleHerbivoreName n
## <fct> <int>
## 1 cattle 165
## 2 sheep 84
## 3 rodents 43
## 4 insects 42
## 5 deer 35
## 6 lagomorphs 32
## 7 horses 28
## 8 antelopes 26
## 9 voles 22
## 10 goats 21
## # ... with 11 more rows
## Capitalize the names of each
topherb[,"simpleHerbivoreName"] <- str_to_title(topherb$simpleHerbivoreName)
ggplot(topherb, aes(x=reorder(simpleHerbivoreName, -n), y=n)) + geom_bar(stat="identity") + theme_Publication() + xlab("") + ylab("Number of studies") + coord_flip()
Summary statistics of other categories
studyData %>% group_by(Fenced) %>% summarize(n=length(UniqueID))
## # A tibble: 4 x 2
## Fenced n
## <chr> <int>
## 1 "" 214
## 2 "fenced" 12
## 3 "open" 18
## 4 "open;fenced" 1
studyData %>% group_by(Tilled) %>% summarize(n=length(UniqueID))
## # A tibble: 3 x 2
## Tilled n
## <chr> <int>
## 1 "" 170
## 2 "LA" 71
## 3 "T" 4
studyData %>% group_by(Fire) %>% summarize(n=length(UniqueID))
## # A tibble: 3 x 2
## Fire n
## <chr> <int>
## 1 "" 217
## 2 "no" 2
## 3 "yes" 26
studyData %>% group_by(Fertilization) %>% summarize(n=length(UniqueID))
## # A tibble: 4 x 2
## Fertilization n
## <chr> <int>
## 1 "" 228
## 2 "no" 8
## 3 "yes" 6
## 4 "yes;no" 3
Comparisons of binary variables
## Load packages and functions
library(reshape2)
library(metafor)
source("scripts//meta.evalSE.R")
## Reference data
meta <- data
## Create Unique identifier column
meta[,"UniqueSite"] <- paste(meta$UniqueID, meta$Higher.taxa, meta$Taxa, meta$Estimate, meta$SiteID, sep="-")
## convert se to sd
meta[meta$Stat=="StDev","Value"] <- meta[meta$Stat=="StDev","Value"] / sqrt(as.numeric(meta[meta$Stat=="StDev","replicate"] ))
meta[meta$Stat=="StDev","Stat"] <- "SE"
## drop comparisons that are not pairwise
meta <- meta %>% filter(grazing.compare == "ordinal/binary")
## Use function to extract summary statistics for comparisons
## meta.eval arguments are (meta.data, compare, ids , stats)
grazed.compare <- meta.eval(meta, grazing.level, UniqueSite, Stat)
## Combine the lists into same dataframe
## Rename Columns in second dataframe
grazed.stat <- grazed.compare [[2]] ## extracted statistics
names(grazed.stat) <- c("UniqueSite","grazed_mean","grazed_se","ungrazed_mean","ungrazed_se","grazed_n","ungrazed_n") ## rename columns to match
grazed.raw <- grazed.compare[[1]] ## calculated statistics from raw values
grazed.raw
## UniqueSite extensive_mean
## 1 --Annelida-percent occurrence- 0
## 2 --Arthropoda-biomass (mg)- 0
## 3 --Arthropoda-consumption/kg ha- 0
## 4 --Arthropoda-Density of total Collembla (#/m2)- 0
## 5 --Arthropoda-density/m2- 0
## 6 --Arthropoda-grasshopper days- 0
## 7 --Arthropoda-grasshopper days for all species- 0
## 8 --Arthropoda-percent occurrence- 0
## 9 --Aves-number of territory- 463
## 10 --Mammalia-juvenile survival- 0
## 11 --Mammalia-mean pregnancies- 0
## 12 --Mammalia-n caught/400 trap nights- 0
## 13 --Mammalia-Shannon Diversity Index- 0
## 14 --Mammalia-summer survival- 0
## 15 --Mammalia-winter survival- 0
## 16 --Microbes-Number of species soil - 0
## 17 --Microbes-Number of species vegetation and litter- 0
## 18 --Mollusca-abundance- 0
## 19 --Mollusca-percent occurrence- 0
## 20 --Nematoda-percent occurrence- 0
## 21 --Plantae-Visual obstruction (score)- 0
## grazed_mean intensive_mean non-grazed_mean rotational_mean season-long_mean
## 1 50.0000000 0 0.000000 0.0000 0.0000
## 2 6.9973796 0 0.000000 0.0000 0.0000
## 3 0.0000000 0 0.000000 46.1800 208.6200
## 4 89.5855009 0 0.000000 0.0000 0.0000
## 5 0.0000000 0 0.000000 3.0000 12.8600
## 6 0.0000000 0 0.000000 22.7258 100.6642
## 7 0.0000000 0 0.000000 1044.5834 1334.5832
## 8 17.0976331 0 0.000000 0.0000 0.0000
## 9 0.0000000 285 0.000000 0.0000 0.0000
## 10 1.6618497 0 0.000000 0.0000 0.0000
## 11 8.7829739 0 0.000000 0.0000 0.0000
## 12 4.0714286 0 0.000000 0.0000 0.0000
## 13 0.9125000 0 0.000000 0.0000 0.0000
## 14 0.7011952 0 0.000000 0.0000 0.0000
## 15 0.8120000 0 0.000000 0.0000 0.0000
## 16 5.8243801 0 0.000000 0.0000 0.0000
## 17 7.4950690 0 0.000000 0.0000 0.0000
## 18 83.0000000 0 0.000000 0.0000 0.0000
## 19 39.2631579 0 0.000000 0.0000 0.0000
## 20 65.7500000 0 0.000000 0.0000 0.0000
## 21 0.0000000 0 1.354839 0.0000 0.0000
## ungrazed_mean extensive_se grazed_se intensive_se non-grazed_se
## 1 50.0000000 0 18.08807816 0 0
## 2 13.6700301 0 1.74490614 0 0
## 3 0.0000000 0 0.00000000 0 0
## 4 35.2418464 0 30.38197552 0 0
## 5 0.0000000 0 0.00000000 0 0
## 6 0.0000000 0 0.00000000 0 0
## 7 0.0000000 0 0.00000000 0 0
## 8 22.3452381 0 1.34418703 0 0
## 9 0.0000000 1 0.00000000 1 0
## 10 0.9682081 0 0.67919075 0 0
## 11 4.3769803 0 5.68384365 0 0
## 12 4.8571429 0 1.36851266 0 0
## 13 0.8200000 0 0.09629598 0 0
## 14 0.7270916 0 0.13147410 0 0
## 15 0.8460000 0 0.04800000 0 0
## 16 6.2045349 0 0.56896033 0 0
## 17 9.1814596 0 0.71302071 0 0
## 18 171.0000000 0 0.00000000 0 0
## 19 33.4000000 0 9.28806129 0 0
## 20 63.0000000 0 6.38194067 0 0
## 21 0.0000000 0 0.00000000 0 1
## rotational_se season-long_se ungrazed_se extensive_n grazed_n intensive_n
## 1 0.0000000 0.000000 50.00000000 0 8 0
## 2 0.0000000 0.000000 1.85140894 0 191 0
## 3 10.8420662 77.029017 0.00000000 0 0 0
## 4 0.0000000 0.000000 12.90353196 0 65 0
## 5 0.7204165 4.233155 0.00000000 0 0 0
## 6 6.1950138 25.183664 0.00000000 0 0 0
## 7 334.0559528 258.889798 0.00000000 0 0 0
## 8 0.0000000 0.000000 3.02312250 0 338 0
## 9 0.0000000 0.000000 0.00000000 1 0 1
## 10 0.0000000 0.000000 0.70809249 0 2 0
## 11 0.0000000 0.000000 2.17177548 0 4 0
## 12 0.0000000 0.000000 2.46927747 0 28 0
## 13 0.0000000 0.000000 0.16000000 0 4 0
## 14 0.0000000 0.000000 0.02988048 0 2 0
## 15 0.0000000 0.000000 0.01800000 0 2 0
## 16 0.0000000 0.000000 0.52200932 0 13 0
## 17 0.0000000 0.000000 0.88888792 0 13 0
## 18 0.0000000 0.000000 0.00000000 0 3 0
## 19 0.0000000 0.000000 20.50268275 0 19 0
## 20 0.0000000 0.000000 1.00000000 0 4 0
## 21 0.0000000 0.000000 0.00000000 0 0 0
## non-grazed_n rotational_n season-long_n ungrazed_n
## 1 0 0 0 2
## 2 0 0 0 191
## 3 0 5 5 0
## 4 0 0 0 65
## 5 0 5 5 0
## 6 0 30 30 0
## 7 0 10 10 0
## 8 0 0 0 84
## 9 0 0 0 0
## 10 0 0 0 2
## 11 0 0 0 4
## 12 0 0 0 14
## 13 0 0 0 2
## 14 0 0 0 2
## 15 0 0 0 2
## 16 0 0 0 13
## 17 0 0 0 13
## 18 0 0 0 3
## 19 0 0 0 5
## 20 0 0 0 1
## 21 1 0 0 0
Which countries were most examined
library(maps)
library(raster)
gpsSplit[,"countries"] <- map.where(x = gpsSplit$Longitude, y = gpsSplit$Latitude)
## Number of unique countries
gpsSplit %>% filter(!is.na(countries)) %>% summarize(n=length(unique(countries)))
## # A tibble: 1 x 1
## n
## <int>
## 1 41