These chunks of code are used to produce all the figures in the main and supplementary manuscript.

#LOAD LIBRARIES

library(ggplot2)
library(reshape2)
library(dplyr)
library(tidyr)
library(tidyverse)
library(stringr)
library(highcharter)
library(networkD3)
library(wesanderson)
library(circlize)

#options(repos = list(CRAN="http://cran.rstudio.com/"))
#install.packages("knitr")
library(knitr)
#install.packages("tinytex")
#tinytex::install_tinytex()

#LOAD DATA

#Load Final dataset coded by team  
coded_data = read.csv("ToF-Review-Checked-Code-V3.csv")

#Load Baell's list of predatory journals downloaded from the online list
baells = read.csv("ToF-Review-Baells-List.csv")

#FIGURE 3B: YEAR WISE PUBLICATIONS IN EACH COUNTRY This chunk of code produces Figure 3B in the manuscript which shows the number of publications from each of the countries over the years 1980 to mid-2023.

#function to split the outcomes data into separate columns (of 0 and 1, where 1 is if the outcome was experienced and if not -0)
split_and_encode <- function(data, column_name) {
  data %>%
    separate_rows(column_name, sep = ";") %>%
    mutate(value = 1) %>%
    pivot_wider(names_from = column_name, values_from = value, values_fill = list(value = 0))
}

#NUMBER OF STUDIES IN EACH COUNTRY
#one study was carried out in multiple countries, so first, we split this column 
coded_data_split_countries <- split_and_encode(coded_data, "X3_1")

#the above code makes a 0-1 matrix for each country
#we select only the necessary columns for each country only
country_summary = select(coded_data_split_countries, c("X1_16_a_UID", "X1_4_Published_Year", "India", "Bangladesh", "Sri Lanka", "Pakistan", "Nepal"))

#pivot this above table from wide to long
country_summary_long = country_summary%>%tidyr::pivot_longer(!c(X1_16_a_UID, "X1_4_Published_Year"), names_to = "Country", values_to = c("Presence"))

#remove the zeros from the table 
country_summary_long_ = filter(country_summary_long, Presence==1)

#aggregate the country by year
country_aggregated = country_summary_long_%>%group_by(Country, X1_4_Published_Year)%>%count()

#aggregate the data to get the total number of publications per year (across all countries)
country_aggregated_total <- country_aggregated %>%
  group_by(X1_4_Published_Year) %>%
  summarise(total_n = sum(n))

library(wesanderson)

# Plot of publications in every year in each South Asian country
country_year_plot <- ggplot(country_aggregated, aes(x = X1_4_Published_Year, y = n, fill = Country)) + 
  geom_bar(stat = "identity", position = "stack") +  # Stacked bars using the 'n' value
  xlab("Year") + 
  ylab("Number of Publications") + 
  labs(fill = "Country") + 
  theme_classic() + 
  theme(legend.position = "bottom") + 
  scale_fill_manual(name = "Country", values = wes_palette("Darjeeling1", n = 5)) + 
  scale_x_continuous(breaks = seq(1983, 2023, by = 5)) +
  scale_y_continuous(breaks = seq(0, max(country_aggregated_total$total_n), by = 3),  
                     limits = c(0, max(country_aggregated_total$total_n) * 1.1)) +  
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

#Plot
country_year_plot

#FIGURE 4: WHICH TOF SYSTEMS EXIST in SOUTH ASIA AND WHAT IS THEIR FREQEUNCY? This chunk of code produces Figure 4 in the manuscript. The colors represent the country in which the system was studied.

#function to split the outcomes data into separate columns (of 0 and 1, where 1 is if the outcome was experienced and if not -0)
split_and_encode <- function(data, column_name) {
  data %>%
    separate_rows(column_name, sep = ",") %>%
    mutate(value = 1) %>%
    pivot_wider(names_from = column_name, values_from = value, values_fill = list(value = 0))
}

#TOF SYSTEM STUDIED IN EACH COUNTRY
#split up the different kinds of ToF systems studied 
coded_data_split_systems <- split_and_encode(coded_data, "X3_4")

#function to split the outcomes data into separate columns (of 0 and 1, where 1 is if the outcome was experienced and if not -0)
split_and_encode <- function(data, column_name) {
  data %>%
    separate_rows(column_name, sep = ";") %>%
    mutate(value = 1) %>%
    pivot_wider(names_from = column_name, values_from = value, values_fill = list(value = 0))
}
#one study was carried out in multiple countries, so first, we split this column 
coded_data_split_systems_countries <- split_and_encode(coded_data_split_systems , "X3_1")

#the above code makes a 0-1 matrix for each system and for each country
#we select only the necessary columns for each system and country
systems_country_summary = select(coded_data_split_systems_countries, c("X1_16_a_UID", "India",
                                                     "Bangladesh", "Nepal", 
                                                     "Pakistan", "Sri Lanka", 
                                                     "Monoculture plantation", "Mixed plantation", "Multipurpose trees on farms or lots", 
                                                     "shelterbelt/ hedgerow/ windbreaks", "tree/ forest (home)garden", 
                                                     "Mixed plantation crops", "Alley cropping", "Orchard", 
                                                     "Integrated production of animals (meat and dairy)/ crops and wood", 
                                                     "Aqua-silvo-fishery", "Tree or shrub pasturelands",  "Entomoforestry", 
                                                     "Afforestation", "Reforestation", "Improved or fallow rotation", 
                                                     "Meadow orchards", "Riparian buffer strip", "Wooded pasture products", "Unknown"))

#to pivot columns, first specify which are countries and which are systems 
country_columns <- c("India","Bangladesh", "Nepal", "Pakistan", "Sri Lanka")
tof_columns <- c("Monoculture plantation", "Mixed plantation", "Multipurpose trees on farms or lots", 
                                                     "shelterbelt/ hedgerow/ windbreaks", "tree/ forest (home)garden", 
                                                     "Mixed plantation crops", "Alley cropping", "Orchard", 
                                                     "Integrated production of animals (meat and dairy)/ crops and wood", 
                                                     "Aqua-silvo-fishery", "Tree or shrub pasturelands",  "Entomoforestry", 
                                                     "Afforestation", "Reforestation", "Improved or fallow rotation", 
                                                     "Meadow orchards", "Riparian buffer strip", "Wooded pasture products", "Unknown")

# Pivot longer for countries
long_country <- systems_country_summary %>%
  pivot_longer(cols = all_of(country_columns), names_to = "Country", values_to = "Presence_Country") %>%
  filter(Presence_Country == 1)  # Keep only rows where country is present

#next, pivot longer for tof systems 
long_tof_systems <- systems_country_summary %>%
  pivot_longer(cols = all_of(tof_columns), names_to = "TOF_Type", values_to = "Presence_TOF_Type") %>%
  filter(Presence_TOF_Type == 1)  # Keep only rows where forest type is present

#combine these two long tables by UID
df_combined <- full_join(long_country, long_tof_systems, by = "X1_16_a_UID")

#aggregate TOF system by country
df_aggrgated = df_combined%>%group_by(Country, TOF_Type)%>%count()

#to plot recode the ToF system column 
df_aggrgated <- df_aggrgated %>%
  mutate(TOF_Type = recode(TOF_Type, 
                             "Multipurpose trees on farms or lots" = "Multipurpose trees on farms", 
                           "shelterbelt/ hedgerow/ windbreaks" = "Shelterbelts", 
                           "tree/ forest (home)garden" = "Tree/Forest garden",
                           "Integrated production of animals (meat and dairy)/ crops and wood" = "Integrated production"))

#arrange in descending order
df_aggrgated_ordered <- df_aggrgated %>%
  group_by(TOF_Type) %>%
  summarise(total_count = sum(n)) %>%
  arrange(desc(total_count)) %>%
  mutate(TOF_Type = factor(TOF_Type, levels = TOF_Type)) %>%
  left_join(df_aggrgated, by = "TOF_Type")

#recode everything under 10 observations 
#see TOF_systems_aggregated_ordered. Recode entomoforestry, meadow orchards, unknown, wooded pasture products, Riparian buffer strip, Aqua-silvo-fishery
df_aggrgated_ordered <- df_aggrgated_ordered %>%
  mutate(TOF_Type = recode(TOF_Type, 
                             "Unknown" = "Other", "Entomoforestry" = "Other", "Wooded pasture products" = "Other", "Riparian buffer strip" = "Other", "Aqua-silvo-fishery" = "Other", "Meadow orchards" = "Other"))

#Plot bars with countries as colours
tof_aggregated_plot = ggplot(df_aggrgated_ordered, aes(x= reorder(TOF_Type, -total_count), y= n, fill = Country)) + geom_col() + theme_classic()+  labs(
  x = "Trees outside Forest System",
  y = "Frequency") +
  scale_fill_manual(name = "Country", values = wes_palette("Darjeeling1", n = 5))+
  scale_x_discrete(labels = function(x) str_wrap(x, width = 14)) + 
  scale_y_continuous(breaks = seq(0, max(df_aggrgated_ordered$total_count), by = 10)) +# Wrap text with a width of 10 characters
  theme_classic() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 12),  # Adjust legend text size
    legend.title = element_text(size = 14),  # Adjust legend title size
    axis.title.x = element_text(size = 14),   # Adjust x-axis title size
    axis.title.y = element_text(size = 14),   # Adjust y-axis title size
    axis.text.x = element_text(angle = 45, hjust = 1, size =14),
    axis.text.y = element_text(hjust = 1, size =14)
  )

tof_aggregated_plot

#FIGURE 5A: SANKEY PLOT OF OUTCOMES AND THE SYSTEMS THEY ARE STUDIED IN This chunk of code produces Figure 5B in the manuscript. The sankey plot shows the number of publications that examined the seven human wellbeing outcomes in the different trees outside forest systems.

#First split the outcomes
#function to split the outcomes data into separate columns (of 0 and 1, where 1 is if the outcome was experienced and if not -0)
split_and_encode <- function(data, column_name) {
  data %>%
    separate_rows(column_name, sep = ",") %>%
    mutate(value = 1) %>%
    pivot_wider(names_from = column_name, values_from = value, values_fill = list(value = 0))
}

#WELLBEING OUTCOME STUDIED
#use the split code above to split all the buckets of outcomes we studied
coded_data_split_outcomes = split_and_encode(coded_data, "X4_1")

#the above code makes a 0-1 matrix for each human wellbeing outcome
#we select only the necessary columns for each outcome
outcomes_summary = select(coded_data_split_outcomes, c("X1_16_a_UID","X3_1", "Material and living standards related outcomes","Economic security related outcomes", "Social relationships related outcomes", "Education related outcomes", "Work and Leisure related outcomes","Agency and Political voice related outcomes", "Health related outcomes"))

#pivot this above table from wide to long
outcomes_summary_long = outcomes_summary%>%tidyr::pivot_longer(!c(X1_16_a_UID, X3_1), names_to = "Wellbeing_Outcome", values_to = c("Presence"))

#remove the zeros from the table 
outcomes_summary_long = filter(outcomes_summary_long, Presence==1)

#TOF SYSTEM STUDIED IN EACH COUNTRY
#split up the different kinds of ToF systems studied 
coded_data_split_systems <- split_and_encode(coded_data, "X3_4")

#the above code makes a 0-1 matrix for each system
#we select only the necessary columns for each system
systems_summary = select(coded_data_split_systems, c("X1_16_a_UID", "X3_1", "Monoculture plantation", "Mixed plantation", "Multipurpose trees on farms or lots", "shelterbelt/ hedgerow/ windbreaks", "tree/ forest (home)garden", "Mixed plantation crops", "Alley cropping", "Orchard", "Integrated production of animals (meat and dairy)/ crops and wood", "Aqua-silvo-fishery", "Tree or shrub pasturelands",  "Entomoforestry", "Afforestation", "Reforestation", "Improved or fallow rotation", "Meadow orchards", "Riparian buffer strip", "Wooded pasture products", "Unknown"))

#pivot this above table from wide to long
systems_summary_long = systems_summary%>%tidyr::pivot_longer(!c(X1_16_a_UID, X3_1), names_to = "ToF_System", values_to = c("Presence"))

#Combine the systems and outcomes tables made above
systems_outcomes = full_join(outcomes_summary_long, systems_summary_long)
## Joining, by = c("X1_16_a_UID", "X3_1", "Presence")
#summarise outcomes studied in each system in each country
systems_outcomes_summary = systems_outcomes%>%group_by(ToF_System, Wellbeing_Outcome, X3_1)%>%summarise(across(everything(), sum))
## `summarise()` has grouped output by 'ToF_System', 'Wellbeing_Outcome'. You can
## override using the `.groups` argument.
#remove the zeros from the table 
systems_outcomes_summary = filter(systems_outcomes_summary, Presence>=1)

#Use this above file to make a sankey diagram
#systems_outcomes_sankey = select(systems_outcomes, -c(X1_16_a_UID, #Presence))
                                 
#select only three columns
for_sankey = dplyr::select(systems_outcomes_summary, ToF_System, Wellbeing_Outcome, Presence)

#recode the systems with fewer than 10 pubs
for_sankey = for_sankey%>%mutate(ToF_System = recode(ToF_System, 
                             "Unknown" = "Other", "Entomoforestry" = "Other", "Wooded pasture products" = "Other", "Riparian buffer strip" = "Other", "Aqua-silvo-fishery" = "Other", "Meadow orchards" = "Other"))

for_sankey = for_sankey%>%group_by(ToF_System, Wellbeing_Outcome)%>%summarise(across(everything(), sum))
## `summarise()` has grouped output by 'ToF_System'. You can override using the
## `.groups` argument.
# Create a color palette of 21 color-blind friendly colors
color_blind_friendly_palette <- c(
  "Afforestation" = "#E69F00",  # Orange
  "Agency and Political voice related outcomes" = "#56B4E9",  # Sky Blue
  "Economic security related outcomes" = "#009E73",  # Teal
  "Alley cropping" = "#F0E442",  # Yellow
  "Integrated production of animals (meat and dairy)/ crops and wood" = "#0072B2",  # Blue
  "Education related outcomes" = "#D55E00",  # Vermilion
  "Improved or fallow rotation" = "#CC79A7",  # Pink
  "Mixed plantation" = "#999999",  # Grey
  "Mixed plantation crops" = "#F2A900",  # Gold
  "Health related outcomes" = "#007A33",  # Dark Green
  "Material and living standards related outcomes" = "#A30000",  # Dark Red
  "Monoculture plantation" = "#FFB6C1",  # Light Pink
  "Orchard" = "#8C8C8C",  # Medium Grey
  "shelterbelt/ hedgerow/ windbreaks" ="#C2C2C2",  # Light Grey
  "Multipurpose trees on farms or lots" = "#FF5733",  # Bright Red
  "Reforestation" ="#FFC300",  # Bright Yellow
  "Tree or shrub pasturelands" = "#DAF7A6",  # Light Green
  "Social relationships related outcomes" = "#581845",  # Dark Purple
  "tree/ forest (home)garden" ="#900C3F",  # Maroon
  "Work and Leisure related outcomes" = "#FF5733",  # Bright Orange
  #"#28B463",  # Bright Green
  "Other" = "#4B0082",   # Indigo
  "NA" = "#C70039" # Cherry Red"
)

nodes <- for_sankey %>%
  pivot_longer(cols = c(ToF_System, Wellbeing_Outcome), names_to = "outcome_type", values_to = "outcome") %>%
  select(outcome) %>%
  distinct() %>%
  mutate(id = row_number() - 1) %>%
  mutate(color = color_blind_friendly_palette[outcome])

links <- for_sankey %>%
  left_join(nodes %>% rename(outcome1_id = id, ToF_System = outcome), by = "ToF_System") %>%
  left_join(nodes %>% rename(outcome2_id = id, Wellbeing_Outcome = outcome, target_color = color), by = "Wellbeing_Outcome")%>%
  select(source = outcome1_id, target = outcome2_id, Presence, color = target_color)
## Adding missing grouping variables: `ToF_System`
#Sankey
sankey <- sankeyNetwork(
  Links = links,
  Nodes = nodes,
  Source = "source",
  Target = "target",
  Value = "Presence",
  NodeID = "outcome",
  units = "Units",
  NodeGroup = "outcome",
  fontSize = 14,
  iterations = 0,
  LinkGroup = "color") # Pass colors for links
## Links is a tbl_df. Converting to a plain data frame.
## Nodes is a tbl_df. Converting to a plain data frame.
sankey$x$nodes <-
  sankey$x$nodes %>% 
  mutate(is_source_node = name %in% links$source)

#Remove the labels, which we will manually add in PPT for better legibility
htmlwidgets::onRender(
  sankey,
  '
  function(el, x) {
    // Hide node labels by setting their opacity to 0
    d3.select(el).selectAll(".node text").style("opacity", 0);
  }
  '
)
# Plot Sankey diagram
print(sankey)
sankey 

#FIGURE 5B: CIRCOS PLOT OF THE CO-OCCURENCE OF CATEGORIES OF WELLBEING OUTCOMES This chunk of code produces Figure 5B, which is a Circos plot showing the co-occurence of the human wellbeing outcomes in the reviewed literature.

#function to split the outcomes data into separate columns (of 0 and 1, where 1 is if the outcome was experienced and if not -0)
split_and_encode <- function(data, column_name) {
  data %>%
    separate_rows(column_name, sep = ",") %>%
    mutate(value = 1) %>%
    pivot_wider(names_from = column_name, values_from = value, values_fill = list(value = 0))
}

#use the split code above to split all the buckets of outcomes we studied
coded_data_split_outcomes = split_and_encode(coded_data, "X4_1")

#the above code makes a 0-1 matrix for each human wellbeing outcome
#we select only the necessary columns for each outcome
outcomes_summary = select(coded_data_split_outcomes, c("X1_16_a_UID","X3_1", "Material and living standards related outcomes","Economic security related outcomes", "Social relationships related outcomes", "Education related outcomes", "Work and Leisure related outcomes","Agency and Political voice related outcomes", "Health related outcomes"))


# Prepare data for Circos plot
# Create a matrix of presence/absence
presence_matrix <- as.matrix(outcomes_summary[-c(1, 2)])  # Exclude ID and Country
rownames(presence_matrix) <- outcomes_summary$X3_1  # Set country names as row names

# Convert presence/absence to a more suitable format for Circos
presence_matrix[presence_matrix == 0] <- NA  # Replace 0s with NA

colors <- c(
  "Agency and Political voice related outcomes" = "#1f77b4",
  "Economic security related outcomes" = "#ff7f0e",
  "Health related outcomes" = "#2ca02c",
  "Material and living standards related outcomes" = "#d62728",
  "Social relationships related outcomes" = "#9467bd",
  "Work and Leisure related outcomes" = "#8c564b",
  "Education related outcomes" = "#FF0000")

# Create a presence-absence matrix for the outcomes
outcomes_summary <- outcomes_summary[-c(1, 2)]  # Exclude ID and Country
presence_matrix <- as.matrix(outcomes_summary)

# Calculate co-occurrences
co_occurrence <- t(presence_matrix) %*% presence_matrix
diag(co_occurrence) <- 0  # Set diagonal to 0 to ignore self-loops

# Set row and column names
rownames(co_occurrence) <- colnames(co_occurrence) <- colnames(outcomes_summary)

# Create the Circos plot
circos.clear()  # Clear previous settings if any
chordDiagram(co_occurrence, transparency = 0.5, grid.col = colors,
              annotationTrack = "grid", preAllocateTracks = 1)

# Remove labels while keeping the bands
circos.trackPlotRegion(track.index = 1, bg.border = NA)  # Just to ensure the track is set without labels

#FIGURE 6: DIRECTION OF WELLBEING OUTCOMES FOR ENTIRE DATASET (ALL LAND TENURES) This chunk of code produces Figure 6, which is a bar plots showing the direction of change in the human wellbeing outcomes in the literature reviewed. We have combined agency and political voice outcomes and social relationship outcomes for the purpose of this figure.

#FOR ALL DATA COMBINED

#MATERIAL OUTCOMES
fuelwood_outcomes_sum = coded_data%>%group_by(X4_2_a)%>%count()
fuelwood_outcomes_sum$outcome = "Fuelwood"
fuelwood_outcomes_sum = fuelwood_outcomes_sum%>%rename(Direction = X4_2_a)

biomass_outcomes_sum = coded_data%>%group_by(X4_2_b)%>%count()
biomass_outcomes_sum$outcome = "Biomass"
biomass_outcomes_sum = biomass_outcomes_sum%>%rename(Direction = X4_2_b)

fiber_outcome_sum = coded_data%>%group_by(X4_2_c)%>%count()
fiber_outcome_sum$outcome = "Fiber"
fiber_outcome_sum = fiber_outcome_sum%>%rename(Direction = X4_2_c)

housingmat_outcome_sum = 
  coded_data%>%group_by(X4_2_d)%>%count()
housingmat_outcome_sum$outcome = "Housing materials"
housingmat_outcome_sum = housingmat_outcome_sum%>%rename(Direction = X4_2_d)

water_outcome_sum = 
  coded_data%>%group_by(X4_2_e)%>%count()
water_outcome_sum$outcome = "Water availability"
water_outcome_sum = water_outcome_sum%>%rename(Direction = X4_2_e)

assets_outcome_sum = 
  coded_data%>%group_by(X4_2_g)%>%count()
assets_outcome_sum$outcome = "Assets"
assets_outcome_sum = assets_outcome_sum%>%rename(Direction = X4_2_g)

shelter_outcome_sum = 
  coded_data%>%group_by(X4_2_h)%>%count()
shelter_outcome_sum$outcome = "Adequete shelter"
shelter_outcome_sum = shelter_outcome_sum%>%rename(Direction = X4_2_h)

material_outcomes = rbind(shelter_outcome_sum, assets_outcome_sum, water_outcome_sum, housingmat_outcome_sum, fiber_outcome_sum, biomass_outcomes_sum, fuelwood_outcomes_sum)

material_outcomes$Direction = as.factor(material_outcomes$Direction)

#recode factors with spelling mistakes
material_outcomes = material_outcomes%>%mutate(Direction = fct_recode(Direction, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))

#Remove the factor which indicates no code
to_remove = material_outcomes%>%filter(Direction == "")
to_remove_2 = material_outcomes%>%filter(Direction == "Unreported")

material_outcomes = anti_join(material_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
material_outcomes = anti_join(material_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
material_outcomes$Direction <- droplevels(material_outcomes$Direction)
material_outcomes <- material_outcomes%>%drop_na(Direction)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

material_plot = ggplot(material_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Material outcomes")+ ylab("Number of publications")+ theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#HEALTH OUTCOMES
disease_outcomes_sum = coded_data%>%group_by(X4_2_i)%>%count()
disease_outcomes_sum$outcome = "Disease"
disease_outcomes_sum = disease_outcomes_sum%>%rename(Direction = X4_2_i)

medicine_outcomes_sum = coded_data%>%group_by(X4_2_j)%>%count()
medicine_outcomes_sum$outcome = "Medical care"
medicine_outcomes_sum = medicine_outcomes_sum%>%rename(Direction = X4_2_j)

mental_outcomes_sum = coded_data%>%group_by(X4_2_k)%>%count()
mental_outcomes_sum$outcome = "Mental health"
mental_outcomes_sum = mental_outcomes_sum%>%rename(Direction = X4_2_k)

long_outcomes_sum = coded_data%>%group_by(X4_2_l)%>%count()
long_outcomes_sum$outcome = "Longevity"
long_outcomes_sum = long_outcomes_sum%>%rename(Direction = X4_2_l)

insurance_outcomes_sum = coded_data%>%group_by(X4_2_m)%>%count()
insurance_outcomes_sum$outcome = "Insurance"
insurance_outcomes_sum = insurance_outcomes_sum%>%rename(Direction = X4_2_m)

diet_outcomes_sum = coded_data%>%group_by(X4_2_n)%>%count()
diet_outcomes_sum$outcome = "Diet"
diet_outcomes_sum = diet_outcomes_sum%>%rename(Direction = X4_2_n)

health_outcomes = rbind(diet_outcomes_sum, mental_outcomes_sum, insurance_outcomes_sum, long_outcomes_sum, medicine_outcomes_sum, disease_outcomes_sum)

health_outcomes$Direction = as.factor(health_outcomes$Direction)
levels(health_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "neutral"                                       
## [6] "Neutral"                                       
## [7] "Unreported"
#recode factors with spelling mistakes
health_outcomes = health_outcomes%>%mutate(Direction = fct_recode(Direction, "Neutral" = "neutral"))

#Remove the factor which indicates no code
to_remove = health_outcomes%>%filter(Direction == "")
to_remove_2 = health_outcomes%>%filter(Direction == "Unreported")

health_outcomes = anti_join(health_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
health_outcomes = anti_join(health_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
health_outcomes$Direction <- droplevels(health_outcomes$Direction)
health_outcomes <- health_outcomes%>%drop_na(Direction)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

health_plot = ggplot(health_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Health outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#EDUCATIONAL OUTCOMES
education_outcomes_sum = coded_data%>%group_by(X4_2_o)%>%count()
education_outcomes_sum$outcome = "Education"
education_outcomes_sum = education_outcomes_sum%>%rename(Direction = X4_2_o)

vocation_outcomes_sum = coded_data%>%group_by(X4_2_p)%>%count()
vocation_outcomes_sum$outcome = "Vocational training"
vocation_outcomes_sum = vocation_outcomes_sum%>%rename(Direction = X4_2_p)

education_outcomes = rbind(education_outcomes_sum, vocation_outcomes_sum)

education_outcomes$Direction = as.factor(education_outcomes$Direction)
levels(education_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Neutral"                                       
## [6] "Unreported"
#Remove the factor which indicates no code
to_remove = education_outcomes%>%filter(Direction == "")
to_remove_2 = education_outcomes%>%filter(Direction == "Unreported")

education_outcomes = anti_join(education_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
education_outcomes = anti_join(education_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
education_outcomes$Direction <- droplevels(education_outcomes$Direction)
education_outcomes = education_outcomes%>%drop_na(Direction)

education_plot = ggplot(education_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Educational outcomes")+  ylab("Number of publications")+ theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#LEISURE OUTCOMES 
nature_outcomes_sum = coded_data%>%group_by(X4_2_q)%>%count()
nature_outcomes_sum$outcome = "Access to nature"
nature_outcomes_sum = nature_outcomes_sum%>%rename(Direction = X4_2_q)

qualtiy_outcomes_sum = coded_data%>%group_by(X4_2_r)%>%count()
qualtiy_outcomes_sum$outcome = "Quality of life"
qualtiy_outcomes_sum = qualtiy_outcomes_sum%>%rename(Direction = X4_2_r)

recreation_outcomes_sum = coded_data%>%group_by(X4_2_s)%>%count()
recreation_outcomes_sum$outcome = "Recreation"
recreation_outcomes_sum = recreation_outcomes_sum%>%rename(Direction = X4_2_s)
  
belief_outcomes_sum = coded_data%>%group_by(X4_2_t)%>%count()
belief_outcomes_sum$outcome = "Belief system"
belief_outcomes_sum = belief_outcomes_sum%>%rename(Direction = X4_2_t)

move_outcomes_sum = coded_data%>%group_by(X4_2_u)%>%count()
move_outcomes_sum$outcome = "Freedom to move"
move_outcomes_sum = move_outcomes_sum%>%rename(Direction = X4_2_u)

leisure_outcomes = rbind(move_outcomes_sum, nature_outcomes_sum, belief_outcomes_sum, recreation_outcomes_sum, qualtiy_outcomes_sum)

leisure_outcomes$Direction = as.factor(leisure_outcomes$Direction)
levels(leisure_outcomes$Direction)
## [1] ""                                              
## [2] "Increase"                                      
## [3] "Mixed - increased for some, decreased for some"
## [4] "Neutral"                                       
## [5] "Unreported"
#Remove the factor which indicates no code
to_remove = leisure_outcomes%>%filter(Direction == "")
to_remove_2 = leisure_outcomes%>%filter(Direction == "Unreported")

leisure_outcomes = anti_join(leisure_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
leisure_outcomes = anti_join(leisure_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
leisure_outcomes$Direction <- droplevels(leisure_outcomes$Direction)
leisure_outcomes = leisure_outcomes%>%drop_na(Direction)

leisure_plot = ggplot(leisure_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Leisure/Work outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  scale_y_continuous(breaks = seq(0, 10, by= 2))

#AGENCY OUTCOMES
political_outcomes_sum = coded_data%>%group_by(X4_2_v)%>%count()
political_outcomes_sum$outcome = "Political voice"
political_outcomes_sum = political_outcomes_sum%>%rename(Direction = X4_2_v)

recognition_outcomes_sum = coded_data%>%group_by(X4_2_w)%>%count()
recognition_outcomes_sum$outcome = "Local recognition"
recognition_outcomes_sum = recognition_outcomes_sum%>%rename(Direction = X4_2_w)

equity_outcomes_sum = coded_data%>%group_by(X4_2_x)%>%count()
equity_outcomes_sum$outcome = "Social equity"
equity_outcomes_sum = equity_outcomes_sum%>%rename(Direction = X4_2_x)

interactions_outcomes_sum = coded_data%>%group_by(X4_2_y)%>%count()
interactions_outcomes_sum$outcome = "Community interactions"
interactions_outcomes_sum = interactions_outcomes_sum%>%rename(Direction = X4_2_y)

trust_outcomes_sum = coded_data%>%group_by(X4_2_z)%>%count()
trust_outcomes_sum$outcome = "Trust"
trust_outcomes_sum = trust_outcomes_sum%>%rename(Direction = X4_2_z)

agency_outcomes = rbind(political_outcomes_sum, recognition_outcomes_sum, trust_outcomes_sum, interactions_outcomes_sum, equity_outcomes_sum)

agency_outcomes$Direction = as.factor(agency_outcomes$Direction)
levels(agency_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Neutral"                                       
## [6] "Unreported"                                    
## [7] "UNreported"
agency_outcomes = agency_outcomes%>%mutate(Direction = fct_recode(Direction, "Unreported" = "UNreported"))

#Remove the factor which indicates no code
to_remove = agency_outcomes%>%filter(Direction == "")
to_remove_2 = agency_outcomes%>%filter(Direction == "Unreported")

agency_outcomes = anti_join(agency_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
agency_outcomes = anti_join(agency_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
agency_outcomes$Direction <- droplevels(agency_outcomes$Direction)
agency_outcomes = agency_outcomes%>%drop_na(Direction)

agency_plot = ggplot(agency_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+ geom_bar(position="stack", stat="identity") +xlab("Agency/ Social outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#ECONOMIC OUTCOMES
income_outcomes_sum = coded_data%>%group_by(X4_2_aa)%>%count()
income_outcomes_sum$outcome = "Employment/ income"
income_outcomes_sum = income_outcomes_sum%>%rename(Direction = X4_2_aa)

loan_outcomes_sum = coded_data%>%group_by(X4_2_ab)%>%count()
loan_outcomes_sum$outcome = "Access to loans"
loan_outcomes_sum = loan_outcomes_sum%>%rename(Direction = X4_2_ab)

debt_outcomes_sum = coded_data%>%group_by(X4_2_ac)%>%count()
debt_outcomes_sum$outcome = "Lack of debt"
debt_outcomes_sum = debt_outcomes_sum%>%rename(Direction = X4_2_ac)

savings_outcomes_sum = coded_data%>%group_by(X4_2_ad)%>%count()
savings_outcomes_sum$outcome = "Savings/ provision for dependents"
savings_outcomes_sum = savings_outcomes_sum%>%rename(Direction = X4_2_ad)

economic_outcomes = rbind(savings_outcomes_sum, debt_outcomes_sum, loan_outcomes_sum, income_outcomes_sum)

economic_outcomes$Direction = as.factor(economic_outcomes$Direction)
levels(economic_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Neutral"                                       
## [6] "Unreported"
#Remove the factor which indicates no code
to_remove = economic_outcomes%>%filter(Direction == "")
to_remove_2 = economic_outcomes%>%filter(Direction == "Unreported")

economic_outcomes = anti_join(economic_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
economic_outcomes = anti_join(economic_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
economic_outcomes$Direction <- droplevels(economic_outcomes$Direction)
economic_outcomes = economic_outcomes%>%drop_na(Direction)
                                          
economic_plot = ggplot(economic_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Economic outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#PLot all together

library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.1.3
benefits_plot = ggarrange(material_plot, leisure_plot, health_plot, agency_plot, economic_plot, education_plot, ncol = 2, nrow = 3, common.legend = TRUE)

benefits_plot

#FIGURE S1: DIRECTION OF WELLBEING OUTCOMES FOR COMMUNITY-LAND TENURE ONLY This chunk of code produces Figure S1, which is a bar plots showing the direction of change in the human wellbeing outcomes for publications focused on communally held land tenure alone. We have combined agency and political voice outcomes and social relationship outcomes in this figure.

#FOR COMMUNALLY HELD LAND

community_data = coded_data%>%dplyr::filter(X3_8 == "Community-held")

#MATERIAL OUTCOMES
fuelwood_outcomes_sum = community_data%>%group_by(X4_2_a)%>%count()
fuelwood_outcomes_sum$outcome = "Fuelwood"
fuelwood_outcomes_sum = fuelwood_outcomes_sum%>%rename(Direction = X4_2_a)

biomass_outcomes_sum = community_data%>%group_by(X4_2_b)%>%count()
biomass_outcomes_sum$outcome = "Biomass"
biomass_outcomes_sum = biomass_outcomes_sum%>%rename(Direction = X4_2_b)

fiber_outcome_sum = community_data%>%group_by(X4_2_c)%>%count()
fiber_outcome_sum$outcome = "Fiber"
fiber_outcome_sum = fiber_outcome_sum%>%rename(Direction = X4_2_c)

housingmat_outcome_sum = 
  community_data%>%group_by(X4_2_d)%>%count()
housingmat_outcome_sum$outcome = "Housing materials"
housingmat_outcome_sum = housingmat_outcome_sum%>%rename(Direction = X4_2_d)

water_outcome_sum = 
  community_data%>%group_by(X4_2_e)%>%count()
water_outcome_sum$outcome = "Water availability"
water_outcome_sum = water_outcome_sum%>%rename(Direction = X4_2_e)

assets_outcome_sum = 
  community_data%>%group_by(X4_2_g)%>%count()
assets_outcome_sum$outcome = "Assets"
assets_outcome_sum = assets_outcome_sum%>%rename(Direction = X4_2_g)

shelter_outcome_sum = 
  community_data%>%group_by(X4_2_h)%>%count()
shelter_outcome_sum$outcome = "Adequete shelter"
shelter_outcome_sum = shelter_outcome_sum%>%rename(Direction = X4_2_h)

material_outcomes = rbind(shelter_outcome_sum, assets_outcome_sum, water_outcome_sum, housingmat_outcome_sum, fiber_outcome_sum, biomass_outcomes_sum, fuelwood_outcomes_sum)

material_outcomes$Direction = as.factor(material_outcomes$Direction)

#recode factors with spelling mistakes
material_outcomes = material_outcomes%>%mutate(Direction = fct_recode(Direction, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Nuetral, decrease, unreported
#Remove the factor which indicates no code
to_remove = material_outcomes%>%filter(Direction == "")
to_remove_2 = material_outcomes%>%filter(Direction == "Unreported")

material_outcomes = anti_join(material_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
material_outcomes = anti_join(material_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
material_outcomes$Direction <- droplevels(material_outcomes$Direction)
material_outcomes <- material_outcomes%>%drop_na(Direction)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

material_plot = ggplot(material_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Material outcomes")+ ylab("Number of publications")+ theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#HEALTH OUTCOMES
disease_outcomes_sum = community_data%>%group_by(X4_2_i)%>%count()
disease_outcomes_sum$outcome = "Disease"
disease_outcomes_sum = disease_outcomes_sum%>%rename(Direction = X4_2_i)

medicine_outcomes_sum = community_data%>%group_by(X4_2_j)%>%count()
medicine_outcomes_sum$outcome = "Medical care"
medicine_outcomes_sum = medicine_outcomes_sum%>%rename(Direction = X4_2_j)

mental_outcomes_sum = community_data%>%group_by(X4_2_k)%>%count()
mental_outcomes_sum$outcome = "Mental health"
mental_outcomes_sum = mental_outcomes_sum%>%rename(Direction = X4_2_k)

long_outcomes_sum = community_data%>%group_by(X4_2_l)%>%count()
long_outcomes_sum$outcome = "Longevity"
long_outcomes_sum = long_outcomes_sum%>%rename(Direction = X4_2_l)

insurance_outcomes_sum = community_data%>%group_by(X4_2_m)%>%count()
insurance_outcomes_sum$outcome = "Insurance"
insurance_outcomes_sum = insurance_outcomes_sum%>%rename(Direction = X4_2_m)

diet_outcomes_sum = community_data%>%group_by(X4_2_n)%>%count()
diet_outcomes_sum$outcome = "Diet"
diet_outcomes_sum = diet_outcomes_sum%>%rename(Direction = X4_2_n)

health_outcomes = rbind(diet_outcomes_sum, mental_outcomes_sum, insurance_outcomes_sum, long_outcomes_sum, medicine_outcomes_sum, disease_outcomes_sum)

health_outcomes$Direction = as.factor(health_outcomes$Direction)
levels(health_outcomes$Direction)
## [1] ""           "Increase"   "neutral"    "Unreported"
#recode factors with spelling mistakes
health_outcomes = health_outcomes%>%mutate(Direction = fct_recode(Direction, "Neutral" = "neutral"))

#Remove the factor which indicates no code
to_remove = health_outcomes%>%filter(Direction == "")
to_remove_2 = health_outcomes%>%filter(Direction == "Unreported")

health_outcomes = anti_join(health_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
health_outcomes = anti_join(health_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
health_outcomes$Direction <- droplevels(health_outcomes$Direction)
health_outcomes <- health_outcomes%>%drop_na(Direction)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

health_plot = ggplot(health_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Health outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#EDUCATIONAL OUTCOMES
education_outcomes_sum = community_data%>%group_by(X4_2_o)%>%count()
education_outcomes_sum$outcome = "Education"
education_outcomes_sum = education_outcomes_sum%>%rename(Direction = X4_2_o)

vocation_outcomes_sum = community_data%>%group_by(X4_2_p)%>%count()
vocation_outcomes_sum$outcome = "Vocational training"
vocation_outcomes_sum = vocation_outcomes_sum%>%rename(Direction = X4_2_p)

education_outcomes = rbind(education_outcomes_sum, vocation_outcomes_sum)

education_outcomes$Direction = as.factor(education_outcomes$Direction)
levels(education_outcomes$Direction)
## [1] ""           "Increase"   "Unreported"
#Remove the factor which indicates no code
to_remove = education_outcomes%>%filter(Direction == "")
to_remove_2 = education_outcomes%>%filter(Direction == "Unreported")

education_outcomes = anti_join(education_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
education_outcomes = anti_join(education_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
education_outcomes$Direction <- droplevels(education_outcomes$Direction)
education_outcomes = education_outcomes%>%drop_na(Direction)

education_plot = ggplot(education_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Educational outcomes")+  ylab("Number of publications")+ theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#LEISURE OUTCOMES 
nature_outcomes_sum = community_data%>%group_by(X4_2_q)%>%count()
nature_outcomes_sum$outcome = "Access to nature"
nature_outcomes_sum = nature_outcomes_sum%>%rename(Direction = X4_2_q)

qualtiy_outcomes_sum = community_data%>%group_by(X4_2_r)%>%count()
qualtiy_outcomes_sum$outcome = "Quality of life"
qualtiy_outcomes_sum = qualtiy_outcomes_sum%>%rename(Direction = X4_2_r)

recreation_outcomes_sum = community_data%>%group_by(X4_2_s)%>%count()
recreation_outcomes_sum$outcome = "Recreation"
recreation_outcomes_sum = recreation_outcomes_sum%>%rename(Direction = X4_2_s)
  
belief_outcomes_sum = community_data%>%group_by(X4_2_t)%>%count()
belief_outcomes_sum$outcome = "Belief system"
belief_outcomes_sum = belief_outcomes_sum%>%rename(Direction = X4_2_t)

move_outcomes_sum = community_data%>%group_by(X4_2_u)%>%count()
move_outcomes_sum$outcome = "Freedom to move"
move_outcomes_sum = move_outcomes_sum%>%rename(Direction = X4_2_u)

leisure_outcomes = rbind(move_outcomes_sum, nature_outcomes_sum, belief_outcomes_sum, recreation_outcomes_sum, qualtiy_outcomes_sum)

leisure_outcomes$Direction = as.factor(leisure_outcomes$Direction)
levels(leisure_outcomes$Direction)
## [1] ""                                              
## [2] "Increase"                                      
## [3] "Mixed - increased for some, decreased for some"
## [4] "Unreported"
#Remove the factor which indicates no code
to_remove = leisure_outcomes%>%filter(Direction == "")
to_remove_2 = leisure_outcomes%>%filter(Direction == "Unreported")

leisure_outcomes = anti_join(leisure_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
leisure_outcomes = anti_join(leisure_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
leisure_outcomes$Direction <- droplevels(leisure_outcomes$Direction)
leisure_outcomes = leisure_outcomes%>%drop_na(Direction)

leisure_plot = ggplot(leisure_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Leisure/Work outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  scale_y_continuous(breaks = seq(0, 10, by= 2))

#AGENCY OUTCOMES
political_outcomes_sum = community_data%>%group_by(X4_2_v)%>%count()
political_outcomes_sum$outcome = "Political voice"
political_outcomes_sum = political_outcomes_sum%>%rename(Direction = X4_2_v)

recognition_outcomes_sum = community_data%>%group_by(X4_2_w)%>%count()
recognition_outcomes_sum$outcome = "Local recognition"
recognition_outcomes_sum = recognition_outcomes_sum%>%rename(Direction = X4_2_w)

equity_outcomes_sum = community_data%>%group_by(X4_2_x)%>%count()
equity_outcomes_sum$outcome = "Social equity"
equity_outcomes_sum = equity_outcomes_sum%>%rename(Direction = X4_2_x)

interactions_outcomes_sum = community_data%>%group_by(X4_2_y)%>%count()
interactions_outcomes_sum$outcome = "Community interactions"
interactions_outcomes_sum = interactions_outcomes_sum%>%rename(Direction = X4_2_y)

trust_outcomes_sum = community_data%>%group_by(X4_2_z)%>%count()
trust_outcomes_sum$outcome = "Trust"
trust_outcomes_sum = trust_outcomes_sum%>%rename(Direction = X4_2_z)

agency_outcomes = rbind(political_outcomes_sum, recognition_outcomes_sum, trust_outcomes_sum, interactions_outcomes_sum, equity_outcomes_sum)

agency_outcomes$Direction = as.factor(agency_outcomes$Direction)
levels(agency_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Unreported"
agency_outcomes = agency_outcomes%>%mutate(Direction = fct_recode(Direction, "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: UNreported
## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported
#Remove the factor which indicates no code
to_remove = agency_outcomes%>%filter(Direction == "")
to_remove_2 = agency_outcomes%>%filter(Direction == "Unreported")

agency_outcomes = anti_join(agency_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
agency_outcomes = anti_join(agency_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
agency_outcomes$Direction <- droplevels(agency_outcomes$Direction)
agency_outcomes = agency_outcomes%>%drop_na(Direction)

agency_plot = ggplot(agency_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+ geom_bar(position="stack", stat="identity") +xlab("Agency/ Social outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#ECONOMIC OUTCOMES
income_outcomes_sum = community_data%>%group_by(X4_2_aa)%>%count()
income_outcomes_sum$outcome = "Employment/ income"
income_outcomes_sum = income_outcomes_sum%>%rename(Direction = X4_2_aa)

loan_outcomes_sum = community_data%>%group_by(X4_2_ab)%>%count()
loan_outcomes_sum$outcome = "Access to loans"
loan_outcomes_sum = loan_outcomes_sum%>%rename(Direction = X4_2_ab)

debt_outcomes_sum = community_data%>%group_by(X4_2_ac)%>%count()
debt_outcomes_sum$outcome = "Lack of debt"
debt_outcomes_sum = debt_outcomes_sum%>%rename(Direction = X4_2_ac)

savings_outcomes_sum = community_data%>%group_by(X4_2_ad)%>%count()
savings_outcomes_sum$outcome = "Savings/ provision for dependents"
savings_outcomes_sum = savings_outcomes_sum%>%rename(Direction = X4_2_ad)

economic_outcomes = rbind(savings_outcomes_sum, debt_outcomes_sum, loan_outcomes_sum, income_outcomes_sum)

economic_outcomes$Direction = as.factor(economic_outcomes$Direction)
levels(economic_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Unreported"
#Remove the factor which indicates no code
to_remove = economic_outcomes%>%filter(Direction == "")
to_remove_2 = economic_outcomes%>%filter(Direction == "Unreported")

economic_outcomes = anti_join(economic_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
economic_outcomes = anti_join(economic_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
economic_outcomes$Direction <- droplevels(economic_outcomes$Direction)
economic_outcomes = economic_outcomes%>%drop_na(Direction)
                                          
economic_plot = ggplot(economic_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Economic outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#PLot all together

library(ggpubr)

community_benefits_plot = ggarrange(material_plot, leisure_plot, health_plot, agency_plot, economic_plot, education_plot, ncol = 2, nrow = 3, common.legend = TRUE)

community_benefits_plot

#FIGURE S2: DIRECTION OF WELLBEING OUTCOMES FOR PRIVATE LAND TENURE ONLY This chunk of code produces Figure S2, which is a bar plots showing the direction of change in the human wellbeing outcomes for publications focused on private land tenure alone. We have combined agency and political voice outcomes and social relationship outcomes in this figure.

#FOR PRIVATE LAND TENURE

private_data = coded_data%>%dplyr::filter(X3_8 == "Private")

#MATERIAL OUTCOMES
fuelwood_outcomes_sum = private_data%>%group_by(X4_2_a)%>%count()
fuelwood_outcomes_sum$outcome = "Fuelwood"
fuelwood_outcomes_sum = fuelwood_outcomes_sum%>%rename(Direction = X4_2_a)

biomass_outcomes_sum = private_data%>%group_by(X4_2_b)%>%count()
biomass_outcomes_sum$outcome = "Biomass"
biomass_outcomes_sum = biomass_outcomes_sum%>%rename(Direction = X4_2_b)

fiber_outcome_sum = private_data%>%group_by(X4_2_c)%>%count()
fiber_outcome_sum$outcome = "Fiber"
fiber_outcome_sum = fiber_outcome_sum%>%rename(Direction = X4_2_c)

housingmat_outcome_sum = 
  private_data%>%group_by(X4_2_d)%>%count()
housingmat_outcome_sum$outcome = "Housing materials"
housingmat_outcome_sum = housingmat_outcome_sum%>%rename(Direction = X4_2_d)

water_outcome_sum = 
  private_data%>%group_by(X4_2_e)%>%count()
water_outcome_sum$outcome = "Water availability"
water_outcome_sum = water_outcome_sum%>%rename(Direction = X4_2_e)

assets_outcome_sum = 
  private_data%>%group_by(X4_2_g)%>%count()
assets_outcome_sum$outcome = "Assets"
assets_outcome_sum = assets_outcome_sum%>%rename(Direction = X4_2_g)

shelter_outcome_sum = 
  private_data%>%group_by(X4_2_h)%>%count()
shelter_outcome_sum$outcome = "Adequete shelter"
shelter_outcome_sum = shelter_outcome_sum%>%rename(Direction = X4_2_h)

material_outcomes = rbind(shelter_outcome_sum, assets_outcome_sum, water_outcome_sum, housingmat_outcome_sum, fiber_outcome_sum, biomass_outcomes_sum, fuelwood_outcomes_sum)

material_outcomes$Direction = as.factor(material_outcomes$Direction)

#recode factors with spelling mistakes
material_outcomes = material_outcomes%>%mutate(Direction = fct_recode(Direction, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease

## Warning: Unknown levels in `f`: decrease
#Remove the factor which indicates no code
to_remove = material_outcomes%>%filter(Direction == "")
to_remove_2 = material_outcomes%>%filter(Direction == "Unreported")

material_outcomes = anti_join(material_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
material_outcomes = anti_join(material_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
material_outcomes$Direction <- droplevels(material_outcomes$Direction)
material_outcomes <- material_outcomes%>%drop_na(Direction)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

material_plot = ggplot(material_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Material outcomes")+ ylab("Number of publications")+ theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#HEALTH OUTCOMES
disease_outcomes_sum = private_data%>%group_by(X4_2_i)%>%count()
disease_outcomes_sum$outcome = "Disease"
disease_outcomes_sum = disease_outcomes_sum%>%rename(Direction = X4_2_i)

medicine_outcomes_sum = private_data%>%group_by(X4_2_j)%>%count()
medicine_outcomes_sum$outcome = "Medical care"
medicine_outcomes_sum = medicine_outcomes_sum%>%rename(Direction = X4_2_j)

mental_outcomes_sum = private_data%>%group_by(X4_2_k)%>%count()
mental_outcomes_sum$outcome = "Mental health"
mental_outcomes_sum = mental_outcomes_sum%>%rename(Direction = X4_2_k)

long_outcomes_sum = private_data%>%group_by(X4_2_l)%>%count()
long_outcomes_sum$outcome = "Longevity"
long_outcomes_sum = long_outcomes_sum%>%rename(Direction = X4_2_l)

insurance_outcomes_sum = private_data%>%group_by(X4_2_m)%>%count()
insurance_outcomes_sum$outcome = "Insurance"
insurance_outcomes_sum = insurance_outcomes_sum%>%rename(Direction = X4_2_m)

diet_outcomes_sum = private_data%>%group_by(X4_2_n)%>%count()
diet_outcomes_sum$outcome = "Diet"
diet_outcomes_sum = diet_outcomes_sum%>%rename(Direction = X4_2_n)

health_outcomes = rbind(diet_outcomes_sum, mental_outcomes_sum, insurance_outcomes_sum, long_outcomes_sum, medicine_outcomes_sum, disease_outcomes_sum)

health_outcomes$Direction = as.factor(health_outcomes$Direction)
levels(health_outcomes$Direction)
## [1] ""                                              
## [2] "Increase"                                      
## [3] "Mixed - increased for some, decreased for some"
## [4] "neutral"                                       
## [5] "Neutral"                                       
## [6] "Unreported"
#recode factors with spelling mistakes
health_outcomes = health_outcomes%>%mutate(Direction = fct_recode(Direction, "Neutral" = "neutral"))

#Remove the factor which indicates no code
to_remove = health_outcomes%>%filter(Direction == "")
to_remove_2 = health_outcomes%>%filter(Direction == "Unreported")

health_outcomes = anti_join(health_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
health_outcomes = anti_join(health_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
health_outcomes$Direction <- droplevels(health_outcomes$Direction)
health_outcomes <- health_outcomes%>%drop_na(Direction)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

health_plot = ggplot(health_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Health outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#EDUCATIONAL OUTCOMES
education_outcomes_sum = private_data%>%group_by(X4_2_o)%>%count()
education_outcomes_sum$outcome = "Education"
education_outcomes_sum = education_outcomes_sum%>%rename(Direction = X4_2_o)

vocation_outcomes_sum = private_data%>%group_by(X4_2_p)%>%count()
vocation_outcomes_sum$outcome = "Vocational training"
vocation_outcomes_sum = vocation_outcomes_sum%>%rename(Direction = X4_2_p)

education_outcomes = rbind(education_outcomes_sum, vocation_outcomes_sum)

education_outcomes$Direction = as.factor(education_outcomes$Direction)
levels(education_outcomes$Direction)
## [1] ""                                              
## [2] "Increase"                                      
## [3] "Mixed - increased for some, decreased for some"
## [4] "Unreported"
#Remove the factor which indicates no code
to_remove = education_outcomes%>%filter(Direction == "")
to_remove_2 = education_outcomes%>%filter(Direction == "Unreported")

education_outcomes = anti_join(education_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
education_outcomes = anti_join(education_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
education_outcomes$Direction <- droplevels(education_outcomes$Direction)
education_outcomes = education_outcomes%>%drop_na(Direction)

education_plot = ggplot(education_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Educational outcomes")+  ylab("Number of publications")+ theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#LEISURE OUTCOMES 
nature_outcomes_sum = private_data%>%group_by(X4_2_q)%>%count()
nature_outcomes_sum$outcome = "Access to nature"
nature_outcomes_sum = nature_outcomes_sum%>%rename(Direction = X4_2_q)

qualtiy_outcomes_sum = private_data%>%group_by(X4_2_r)%>%count()
qualtiy_outcomes_sum$outcome = "Quality of life"
qualtiy_outcomes_sum = qualtiy_outcomes_sum%>%rename(Direction = X4_2_r)

recreation_outcomes_sum = private_data%>%group_by(X4_2_s)%>%count()
recreation_outcomes_sum$outcome = "Recreation"
recreation_outcomes_sum = recreation_outcomes_sum%>%rename(Direction = X4_2_s)
  
belief_outcomes_sum = private_data%>%group_by(X4_2_t)%>%count()
belief_outcomes_sum$outcome = "Belief system"
belief_outcomes_sum = belief_outcomes_sum%>%rename(Direction = X4_2_t)

move_outcomes_sum = private_data%>%group_by(X4_2_u)%>%count()
move_outcomes_sum$outcome = "Freedom to move"
move_outcomes_sum = move_outcomes_sum%>%rename(Direction = X4_2_u)

leisure_outcomes = rbind(move_outcomes_sum, nature_outcomes_sum, belief_outcomes_sum, recreation_outcomes_sum, qualtiy_outcomes_sum)

leisure_outcomes$Direction = as.factor(leisure_outcomes$Direction)
levels(leisure_outcomes$Direction)
## [1] ""           "Increase"   "Neutral"    "Unreported"
#Remove the factor which indicates no code
to_remove = leisure_outcomes%>%filter(Direction == "")
to_remove_2 = leisure_outcomes%>%filter(Direction == "Unreported")

leisure_outcomes = anti_join(leisure_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
leisure_outcomes = anti_join(leisure_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
leisure_outcomes$Direction <- droplevels(leisure_outcomes$Direction)
leisure_outcomes = leisure_outcomes%>%drop_na(Direction)

leisure_plot = ggplot(leisure_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Leisure/Work outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  scale_y_continuous(breaks = seq(0, 10, by= 2))

#AGENCY OUTCOMES
political_outcomes_sum = private_data%>%group_by(X4_2_v)%>%count()
political_outcomes_sum$outcome = "Political voice"
political_outcomes_sum = political_outcomes_sum%>%rename(Direction = X4_2_v)

recognition_outcomes_sum = private_data%>%group_by(X4_2_w)%>%count()
recognition_outcomes_sum$outcome = "Local recognition"
recognition_outcomes_sum = recognition_outcomes_sum%>%rename(Direction = X4_2_w)

equity_outcomes_sum = private_data%>%group_by(X4_2_x)%>%count()
equity_outcomes_sum$outcome = "Social equity"
equity_outcomes_sum = equity_outcomes_sum%>%rename(Direction = X4_2_x)

interactions_outcomes_sum = private_data%>%group_by(X4_2_y)%>%count()
interactions_outcomes_sum$outcome = "Community interactions"
interactions_outcomes_sum = interactions_outcomes_sum%>%rename(Direction = X4_2_y)

trust_outcomes_sum = private_data%>%group_by(X4_2_z)%>%count()
trust_outcomes_sum$outcome = "Trust"
trust_outcomes_sum = trust_outcomes_sum%>%rename(Direction = X4_2_z)

agency_outcomes = rbind(political_outcomes_sum, recognition_outcomes_sum, trust_outcomes_sum, interactions_outcomes_sum, equity_outcomes_sum)

agency_outcomes$Direction = as.factor(agency_outcomes$Direction)
levels(agency_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Neutral"                                       
## [6] "Unreported"
agency_outcomes = agency_outcomes%>%mutate(Direction = fct_recode(Direction, "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: UNreported
## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported

## Warning: Unknown levels in `f`: UNreported
#Remove the factor which indicates no code
to_remove = agency_outcomes%>%filter(Direction == "")
to_remove_2 = agency_outcomes%>%filter(Direction == "Unreported")

agency_outcomes = anti_join(agency_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
agency_outcomes = anti_join(agency_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
agency_outcomes$Direction <- droplevels(agency_outcomes$Direction)
agency_outcomes = agency_outcomes%>%drop_na(Direction)

agency_plot = ggplot(agency_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+ geom_bar(position="stack", stat="identity") +xlab("Agency/ Social outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#ECONOMIC OUTCOMES
income_outcomes_sum = private_data%>%group_by(X4_2_aa)%>%count()
income_outcomes_sum$outcome = "Employment/ income"
income_outcomes_sum = income_outcomes_sum%>%rename(Direction = X4_2_aa)

loan_outcomes_sum = private_data%>%group_by(X4_2_ab)%>%count()
loan_outcomes_sum$outcome = "Access to loans"
loan_outcomes_sum = loan_outcomes_sum%>%rename(Direction = X4_2_ab)

debt_outcomes_sum = private_data%>%group_by(X4_2_ac)%>%count()
debt_outcomes_sum$outcome = "Lack of debt"
debt_outcomes_sum = debt_outcomes_sum%>%rename(Direction = X4_2_ac)

savings_outcomes_sum = private_data%>%group_by(X4_2_ad)%>%count()
savings_outcomes_sum$outcome = "Savings/ provision for dependents"
savings_outcomes_sum = savings_outcomes_sum%>%rename(Direction = X4_2_ad)

economic_outcomes = rbind(savings_outcomes_sum, debt_outcomes_sum, loan_outcomes_sum, income_outcomes_sum)

economic_outcomes$Direction = as.factor(economic_outcomes$Direction)
levels(economic_outcomes$Direction)
## [1] ""                                              
## [2] "Decrease"                                      
## [3] "Increase"                                      
## [4] "Mixed - increased for some, decreased for some"
## [5] "Neutral"                                       
## [6] "Unreported"
#Remove the factor which indicates no code
to_remove = economic_outcomes%>%filter(Direction == "")
to_remove_2 = economic_outcomes%>%filter(Direction == "Unreported")

economic_outcomes = anti_join(economic_outcomes, to_remove)
## Joining, by = c("Direction", "n", "outcome")
economic_outcomes = anti_join(economic_outcomes, to_remove_2)
## Joining, by = c("Direction", "n", "outcome")
economic_outcomes$Direction <- droplevels(economic_outcomes$Direction)
economic_outcomes = economic_outcomes%>%drop_na(Direction)
                                          
economic_plot = ggplot(economic_outcomes, aes(x= factor(outcome), y= n, fill = factor(Direction)))+  geom_bar(position="stack", stat="identity") +xlab("Economic outcomes")+  ylab("Number of publications")+theme_classic()+ scale_fill_manual(name = "Direction", values = colour_palette) + coord_flip()+   scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

#PLot all together

library(ggpubr)

private_benefits_plot = ggarrange(material_plot, leisure_plot, health_plot, agency_plot, economic_plot, education_plot, ncol = 2, nrow = 3, common.legend = TRUE)

private_benefits_plot

#FIGURE S3(A-F): CLOSER EXAMINATION OF THE DIRECTION OF WELLBEING OUTCOMES AND TOF SYSTEMS THEY ARE EXPERIENCED IN This chunk of code produces figures S3(a-f), which are bar plots of the change in direction of all the categories of wellbeing outcomes of trees outside forest systems linked to the ToF systems studied

#create a function to split the outcomes data into separate columns (of 0 and 1, where 1 is if the outcome was experienced and if not -0)
split_and_encode <- function(data, column_name) {
  data %>%
    separate_rows(column_name, sep = ",") %>%
    mutate(value = 1) %>%
    pivot_wider(names_from = column_name, values_from = value, values_fill = list(value = 0))
}

#MATERIAL OUTCOMES

#Splitting the type of ToF system
types_tof <- split_and_encode(coded_data, "X3_4")

types_tof_sum = types_tof %>%
  pivot_longer(c(`Monoculture plantation`, `Mixed plantation`, `Mixed plantation crops`, Afforestation, Reforestation, `Wooded pasture products`, `Multipurpose trees on farms or lots`, `tree/ forest (home)garden`, `Tree or shrub pasturelands`, `Alley cropping`, `shelterbelt/ hedgerow/ windbreaks`, `Integrated production of animals (meat and dairy)/ crops and wood`, Orchard, `Aqua-silvo-fishery`, Entomoforestry, `Improved or fallow rotation`, `Meadow orchards`, `Riparian buffer strip`), names_to = 'ToF_type')

#to plot recode the ToF system column 
types_tof_sum <- types_tof_sum %>%
  mutate(ToF_type = recode(ToF_type, 
                             "Multipurpose trees on farms or lots" = "Multipurpose trees on farms", 
                           "shelterbelt/ hedgerow/ windbreaks" = "Shelterbelts", 
                           "tree/ forest (home)garden" = "Tree/Forest garden",
                           "Integrated production of animals (meat and dairy)/ crops and wood" = "Integrated production"))


#filter for systems with value = 1 the table 
types_tof_sum = filter(types_tof_sum, value>=1)

#Count the number of papers per human wellbeing outcome for each TOF system

#firewood
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_a, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_a'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_a)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_a == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_a", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_a == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_a", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_a = fct_recode(X4_2_a, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
## Warning: Unknown levels in `f`: Nuetral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_a <- factor(types_tof_sum_count$X4_2_a, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

firewood_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_a))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ggtitle("Access to firewood")

firewood_responses

#biomass
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_b, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_b'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_b)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_b == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_b", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_b == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_b", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_b = fct_recode(X4_2_b, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_b <- factor(types_tof_sum_count$X4_2_b, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

biomass_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_b))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ggtitle("Access to biomass")

biomass_responses

#fiber
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_c, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_c'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_c)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_c == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_c", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_c == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_c", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_c = fct_recode(X4_2_c, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_c <- factor(types_tof_sum_count$X4_2_c, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

fiber_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_c))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette) +ggtitle("Access to fiber")

fiber_responses

#housingmat
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_d, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_d'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_d)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_d == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_d", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_d == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_d", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_d = fct_recode(X4_2_d, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_d <- factor(types_tof_sum_count$X4_2_d, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

housingmat_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_d))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to housing materials")

housingmat_responses

#water
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_e, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_e'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_e)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_e == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_e", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_e == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_e", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_e = fct_recode(X4_2_e, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_e <- factor(types_tof_sum_count$X4_2_e, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

water_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_e))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette) +ggtitle("Water availability")

water_responses

#assets
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_g, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_g'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_g)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_g == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_g", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_g == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_g", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_g = fct_recode(X4_2_g, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_g <- factor(types_tof_sum_count$X4_2_g, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

assets_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_g))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to assets")

assets_responses

#shelter
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_h, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_h'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_h)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_h == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_h", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_h == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_h", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_h == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_h", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_h = fct_recode(X4_2_h, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_h <- factor(types_tof_sum_count$X4_2_h, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

shelter_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_h))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to shelter")

shelter_responses

#Plot all together
library(ggpubr)

material_plot = ggarrange(firewood_responses, biomass_responses, assets_responses, water_responses, housingmat_responses, fiber_responses, shelter_responses, ncol = 2, nrow = 4, common.legend = TRUE)

material_plot

#HEALTH OUTCOMES
#Disease
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_i, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_i'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_i)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_i == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_i", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_i == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_i", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_i = fct_recode(X4_2_i, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "neutral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_i <- factor(types_tof_sum_count$X4_2_i, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

disease_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_i))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Protection from disease")

disease_responses

#Medicine
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_j, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_j'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_j)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_j == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_j", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_j == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_j", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_j == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_j", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_j = fct_recode(X4_2_j, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "neutral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_j <- factor(types_tof_sum_count$X4_2_j, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

medicine_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_j))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to medicine")

medicine_responses

#Mental health
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_k, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_k'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_k)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_k == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_k", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_k == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_k", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_k == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_k", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_k = fct_recode(X4_2_k, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "neutral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, neutral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_k <- factor(types_tof_sum_count$X4_2_k, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

mental_health_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_k))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Mental health")

mental_health_responses

#longevity
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_l, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_l'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_l)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_l == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_l", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_l == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_l", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_l == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_l", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_l = fct_recode(X4_2_l, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_l <- factor(types_tof_sum_count$X4_2_l, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

longevity_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_l))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Longevity")

longevity_responses

#health insurance 
#there are no papers in health insurance, so no plot for it
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_m, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_m'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_m)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#diet
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_n, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_n'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_n)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_n == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_n", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_n == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_n", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_n == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_n", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_n = fct_recode(X4_2_n, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_n <- factor(types_tof_sum_count$X4_2_n, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

diet_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_n))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to diversified diet")

diet_responses

#Plot all together
library(ggpubr)

health_plot = ggarrange(diet_responses, mental_health_responses, medicine_responses, disease_responses, longevity_responses, ncol = 2, nrow = 3, common.legend = TRUE)

health_plot 

#EDUCATION OUTCOMES
#education
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_o, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_o'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_o)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_o == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_o", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_o == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_o", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_o == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_o", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_o = fct_recode(X4_2_o, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_o <- factor(types_tof_sum_count$X4_2_o, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

education_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_o))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to educational opportunities")

education_responses

#vocational training
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_p, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_p'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_p)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_p == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_p", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_p == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_p", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_p == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_p", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_p = fct_recode(X4_2_p, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_p <- factor(types_tof_sum_count$X4_2_p, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

vocational_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_p))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to vocational training")

vocational_responses

#Plot all together
library(ggpubr)

education_plot = ggarrange(education_responses, vocational_responses, ncol = 2, nrow = 1, common.legend = TRUE)

education_plot 

#LEISURE OUTCOMES 
#access to nature 

types_tof_sum_count = types_tof_sum%>%group_by(X4_2_q, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_q'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_q)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_q == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_q", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_q == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_q", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_q == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_q", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_q = fct_recode(X4_2_q, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_q <- factor(types_tof_sum_count$X4_2_q, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

nature_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_q))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to nature")

nature_responses

#quality of life 
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_r, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_r'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_r)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_r == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_r", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_r == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_r", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_r == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_r", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_r = fct_recode(X4_2_r, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_r <- factor(types_tof_sum_count$X4_2_r, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

quality_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_r))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Quality of life")

quality_responses

#recreation
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_s, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_s'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_s)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_s == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_s", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_s == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_s", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_s == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_s", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_s = fct_recode(X4_2_s, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_s <- factor(types_tof_sum_count$X4_2_s, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

recreation_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_s))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to recreation")

recreation_responses

#belief system
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_t, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_t'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_t)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_t == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_t", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_t == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_t", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_t == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_t", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_t = fct_recode(X4_2_t, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_t <- factor(types_tof_sum_count$X4_2_t, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

belief_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_t))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Ability to have belief system")

belief_responses

#freedom to move

types_tof_sum_count = types_tof_sum%>%group_by(X4_2_u, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_u'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_u)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_u == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_u", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_u == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_u", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_u == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_u", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_u = fct_recode(X4_2_u, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_u <- factor(types_tof_sum_count$X4_2_u, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

freedom_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_u))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Freedom to move")

freedom_responses

#Plot all together
library(ggpubr)

leisure_plot = ggarrange(freedom_responses, nature_responses, quality_responses, belief_responses, recreation_responses, ncol = 2, nrow = 3, common.legend = TRUE)

leisure_plot 

#AGENCY OUTCOMES 
#political voice 
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_v, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_v'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_v)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_v == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_v", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_v == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_v", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_v == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_v", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_v = fct_recode(X4_2_v, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_v <- factor(types_tof_sum_count$X4_2_v, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

political_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_v))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Ability to have political voice")

political_responses

#local recognition  
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_w, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_w'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_w)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_w == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_w", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_w == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_w", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_w == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_w", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_w = fct_recode(X4_2_w, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_w <- factor(types_tof_sum_count$X4_2_w, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

recognition_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_w))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Local recognition")

recognition_responses

#Social equity
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_x, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_x'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_x)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_x == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_x", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_x == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_x", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_x == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_x", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_x = fct_recode(X4_2_x, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_x <- factor(types_tof_sum_count$X4_2_x, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

social_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_x))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Social equity")

social_responses

#interactions 
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_y, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_y'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_y)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_y == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_y", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_y == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_y", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_y = fct_recode(X4_2_y, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "unreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, unreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_y <- factor(types_tof_sum_count$X4_2_y, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

interactions_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_y))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Community interactions")

interactions_responses

#trust
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_z, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_z'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_z)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_z == "UNreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_z", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_z == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_z", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_z == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_z", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_z = fct_recode(X4_2_z, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_z <- factor(types_tof_sum_count$X4_2_z, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

trust_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_z))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Trust within community")

trust_responses

#Plot all together
library(ggpubr)

social_agency_plot = ggarrange(trust_responses, interactions_responses, political_responses, social_responses, recognition_responses, ncol = 2, nrow = 3, common.legend = TRUE)

social_agency_plot 

#ECONOMIC OUTCOMES
#income
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_aa, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_aa'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_aa)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_aa == "unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_aa", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_aa == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_aa", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_aa == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_aa", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_aa = fct_recode(X4_2_aa, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_aa <- factor(types_tof_sum_count$X4_2_aa, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

income_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_aa))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Income")

income_responses

#loans
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_ab, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_ab'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_ab)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_ab == "UNreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ab", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_ab == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ab", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_ab == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ab", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_ab = fct_recode(X4_2_ab, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_ab <- factor(types_tof_sum_count$X4_2_ab, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

loans_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_ab))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Access to loans")

loans_responses

#debt
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_ac, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_ac'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_ac)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_ac == "UNreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ac", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_ac == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ac", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_ac == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ac", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_ac = fct_recode(X4_2_ac, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_ac <- factor(types_tof_sum_count$X4_2_ac, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

debt_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_ac))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Lack of debt")

debt_responses

#savings 
types_tof_sum_count = types_tof_sum%>%group_by(X4_2_ad, ToF_type)%>%summarise(total_papers = sum(value))
## `summarise()` has grouped output by 'X4_2_ad'. You can override using the
## `.groups` argument.
#drop nas
types_tof_sum_count = types_tof_sum_count%>%drop_na(X4_2_ad)
types_tof_sum_count = types_tof_sum_count%>%drop_na(total_papers)

#remove unreported and no responses (This is becuase of the way the survey was set up and not because there is missing data)
to_remove = types_tof_sum_count%>%filter(X4_2_ad == "UNreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ad", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_ad == "Unreported")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ad", "ToF_type", "total_papers")
to_remove = types_tof_sum_count%>%filter(X4_2_ad == "")
types_tof_sum_count = anti_join(types_tof_sum_count, to_remove)
## Joining, by = c("X4_2_ad", "ToF_type", "total_papers")
#recode factors with spelling mistakes
types_tof_sum_count = types_tof_sum_count%>%mutate(X4_2_ad = fct_recode(X4_2_ad, "Mixed - increased for some, decreased for some" = "Mixed - increased for some, decreased for others", "Neutral" = "Nuetral", "Decrease" = "decrease", "Unreported" = "UNreported"))
## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported

## Warning: Unknown levels in `f`: Mixed - increased for some, decreased for
## others, Nuetral, decrease, UNreported
# Define the desired order of levels for species richness
desired_order <- c("Increase", "Decrease", "Mixed - increased for some, decreased for some", "Neutral")

# Convert the character column to a factor with specified levels
types_tof_sum_count$X4_2_ad <- factor(types_tof_sum_count$X4_2_ad, levels = desired_order)

#convert to percentage
#Compute total counts for each agricultural system
types_tof_total <- types_tof_sum_count %>%
  group_by(ToF_type) %>%
  summarise(total_count = sum(total_papers), .groups = 'drop')

# Join total counts back to the original dataframe
types_tof_sum_count <- types_tof_sum_count %>%
  left_join(types_tof_total, by = "ToF_type") %>%
  mutate(percentage = (total_papers / total_count) * 100)

colour_palette =  c('Increase'='#009E73','Decrease'='#D55E00','Mixed - increased for some, decreased for some'='#999999','Neutral'='#56B4E9')

savings_responses = ggplot(types_tof_sum_count, aes(x= factor(ToF_type), y = percentage)) + geom_col(aes(fill = factor(X4_2_ad))) +xlab("ToF system")+ ylab("Percentage of Responses")+ theme_classic()+coord_flip()+ scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) + scale_fill_manual(name = "Direction", values = colour_palette)+ ggtitle("Savings/ provisions for self/ dependents")

savings_responses

#Plot all together
library(ggpubr)

economic_plot = ggarrange(income_responses, loans_responses, debt_responses, savings_responses, ncol = 2, nrow = 3, common.legend = TRUE)

economic_plot 

#LIMITATIONS OF THE DATASET:CONFLICT OF INTEREST, QUALITY OF STUDY, OTHER DETAILS This chunk of code produces statistics that were used in the Discussion section ‘Limitations of this dataset’ to describe the quality of the studies reviewed.

# What % of the data has an implied conflict? 
implied_conflict = coded_data%>%group_by(X1_19)%>%summarise(sum = n())
#approx 15% of the data has implied conflict of interest

# What % of the papers have replicable methos? 
replicable_methods = coded_data%>%group_by(X6_4)%>%summarise(sum = n())
#approx 23%% of the data has non replicable methods 

# what % of the papers have justified sampling strategy? 
sampling_strategy = coded_data%>%group_by(X6_5_a)%>%summarise(sum = n())

# How many of these journals are on Baell's list of predatory journals? 
baells_yes = baells%>%filter(Beall_list == "Yes")
baells_yes = baells_yes %>% distinct(X1_6_Journal)
baells_yes$On_List = "On_Bealls_List"

coded_baells = left_join(coded_data, baells_yes, by ="X1_6_Journal")

baells_list_total = coded_baells%>%group_by(On_List)%>%summarise(sum = n())

Session information

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_India.1252  LC_CTYPE=English_India.1252   
## [3] LC_MONETARY=English_India.1252 LC_NUMERIC=C                  
## [5] LC_TIME=English_India.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggpubr_0.4.0      knitr_1.40        circlize_0.4.16   wesanderson_0.3.7
##  [5] networkD3_0.4     highcharter_0.9.4 lubridate_1.9.3   forcats_1.0.0    
##  [9] stringr_1.5.1     purrr_1.0.1       readr_2.1.4       tibble_3.2.1     
## [13] tidyverse_2.0.0   tidyr_1.3.0       dplyr_1.0.10      reshape2_1.4.4   
## [17] ggplot2_3.4.4    
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.2          jsonlite_1.8.8      carData_3.0-4      
##  [4] bslib_0.4.0         assertthat_0.2.1    TTR_0.24.3         
##  [7] highr_0.11          cellranger_1.1.0    yaml_2.2.1         
## [10] pillar_1.9.0        backports_1.2.1     lattice_0.20-45    
## [13] glue_1.6.2          rlist_0.4.6.2       digest_0.6.27      
## [16] ggsignif_0.6.3      colorspace_2.0-2    cowplot_1.1.1      
## [19] htmltools_0.5.5     plyr_1.8.7          pkgconfig_2.0.3    
## [22] broom_1.0.1         haven_2.4.1         scales_1.2.1       
## [25] openxlsx_4.2.4      rio_0.5.29          tzdb_0.4.0         
## [28] timechange_0.2.0    generics_0.1.3      farver_2.1.1       
## [31] car_3.0-11          ellipsis_0.3.2      cachem_1.0.6       
## [34] withr_2.5.0         cli_3.6.1           quantmod_0.4.26    
## [37] readxl_1.4.3        magrittr_2.0.3      evaluate_0.23      
## [40] fansi_0.5.0         rstatix_0.7.2       xts_0.12.1         
## [43] foreign_0.8-81      tools_4.1.2         data.table_1.14.0  
## [46] hms_1.1.2           GlobalOptions_0.1.2 lifecycle_1.0.3    
## [49] munsell_0.5.0       zip_2.2.0           compiler_4.1.2     
## [52] jquerylib_0.1.4     rlang_1.1.0         grid_4.1.2         
## [55] rstudioapi_0.14     htmlwidgets_1.6.1   igraph_1.3.5       
## [58] labeling_0.4.2      rmarkdown_2.16      gtable_0.3.1       
## [61] abind_1.4-5         DBI_1.1.3           curl_4.3.2         
## [64] R6_2.5.1            gridExtra_2.3       zoo_1.8-11         
## [67] fastmap_1.1.0       utf8_1.2.1          shape_1.4.6        
## [70] stringi_1.7.6       Rcpp_1.0.7          vctrs_0.6.1        
## [73] tidyselect_1.2.0    xfun_0.33