05-FunctionalAnalyses.Rmd

```{r analyses2, include=FALSE, eval=T}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(tidyverse)
library(googlesheets4)
library(caret)
library(gbm)
theme_set(bayesplot::theme_default())
opts_chunk$set(
  echo = F, message = F, warning = F, fig.height = 6, fig.width = 8,
  cache = T, cache.lazy = F, eval=T)
```

# Functional analyses

In this chapter, I quickly  investigated effects of functional traits and weather on individual growth.

## Methods

I used model to explore individual growth potential relation to functional traits at the individual or species level.
I used either: (1) linear model with a step procedure, (2) boosted regression trees (BRT), and multimodel regressions.

```{r dataan2}
growth <- vroom::vroom("save/ecoevo/growthfull.tsv")
phd <- bind_rows(
  vroom::vroom("data/Measures_Symphonia - AllTraits.tsv") %>% 
    mutate(Genus = "Symphonia") %>%
    rename(Species = Morphotype) %>%
    mutate(Species = ifelse(Species == "Indet.",
                            c("globulifera", "sp.1", "sp.1")[fct_recode(Bark, "globulifera" = "G",
                                                                        "sp.1" =  "S")], Species)),
  vroom::vroom("data/Measures_Eschweilera - AllTraits.tsv") %>%
    filter(!(Plot == 14 & SubPlot == 1 & TreeFieldNum == 760))
) %>% 
  dplyr::select(idTree, SLA, LDMC, LT, LA, CC, brBT, brWD, brBD) %>% 
  left_join(dplyr::select(growth, idTree, species, gmax)) %>% 
  filter(!is.na(gmax)) %>% 
  group_by(species, idTree) %>% 
  summarise_all(mean, na.rm =T)
hydro <- vroom::vroom("data/211124_TraitDatabase_species_means_updatedNames.csv") %>% 
  mutate(species = gsub("_", " ", spName.new)) %>% 
  dplyr::rename(LSWC = LSWC_E.per_corr_mean, RWC = RWC_E.per_corr_mean, Ptlp = ptlp_mean, d13C = d13C.per_mean, CN = C.N_mean,
                stomataD = stomata_density.mm2_mean, gmin = gmin.slope_mean, LWC = LWC_mean, LDMC = LDMC.mg.g_mean,
                N = N.per_mean, C = C.per_mean, P12stem = P12_stem.Mpa_mean, P50stem = P50_stem.MPa_mean, P88stem = P88_stem.MPa_mean,
                P12leaf = P12_leaf.MPa_mean, P50leaf = P50_leaf.MPa_mean, P88leaf = P88_leaf.MPa_mean) %>% 
  dplyr::select(species, LSWC, RWC, Ptlp, d13C, CN, stomataD, gmin, LWC, LDMC, N, C, P12stem, P50stem, P88stem, P12leaf, P50leaf, P88leaf) %>% 
  left_join(group_by(growth, species) %>% summarise(gmax  = median(gmax, na.rm = T)))
fg <- readxl::read_xlsx("data/Appendix_S2-6___S8-13.xlsx", "App.S6-ok", skip = 3) %>% 
  mutate(species = gsub("_", " ", Species)) %>% 
  left_join(group_by(growth, species) %>% summarise(gmax  = median(gmax, na.rm = T))) %>% 
  filter(!is.na(gmax)) %>% 
  dplyr::select(species, gmax, Chlorophyll_content, Thickness, Toughness, Leaf_Area, SLA, C, N, `13C`, Ca, P, K, 
                Trunk_bark_thickness, Sapwood_WSG, WSG, Diameter, SRL, Tissue_density, SRTA, branchiness)
cam <- read_tsv("data/Ziegler2019 - All.tsv") %>% 
  dplyr::select(-Code) %>% 
  reshape2::melt("Species", variable.name = "Trait") %>% 
  mutate(value = gsub("− ", "-", value)) %>% 
  separate(value, c("mean", "sd"), "±", convert = T) %>% 
  mutate(mean = as.numeric(mean)) %>% 
  reshape2::dcast(Species ~ Trait, value.var = "mean") %>% 
  separate(Species, c("Genus", "Species", "Author")) %>% 
  select(-Author) %>% 
  mutate(species = paste(Genus, Species)) %>% 
  left_join(group_by(growth, species) %>% summarise(gmax  = median(gmax, na.rm = T)) %>% separate(species, c("Genus", "Species")))
santi <- read_delim("data/Santiago2018_S1.tsv", delim = " ", skip = 1,
           col_names = c("Genus", "Species", "WD", "X1", "X2",
                         "sapwood_saturated_water_content", "X3", "X4",
                         "sapwood_osmotic_potential_at_full_turgor", "X5", "X6",
                         "sapwood_water_potential_at_turgor_loss_point", "X7", "X8",
                         "total_sapwood_bulk_elastic_modulus", "X9", "X10",
                         "total_sapwood_relative_water_content_at_turgor_loss_point", "X11", "X12",
                         "sapwood_capacitance_at_full_turgor", "X13", "X14")) %>% 
  select(Genus, Species, sapwood_saturated_water_content, sapwood_osmotic_potential_at_full_turgor,
         sapwood_water_potential_at_turgor_loss_point, total_sapwood_bulk_elastic_modulus,
         total_sapwood_relative_water_content_at_turgor_loss_point, sapwood_capacitance_at_full_turgor) %>% 
  left_join(read_delim("data/Santiago2018_S2.tsv", delim = " ", skip = 1,
                       col_names = c("Genus", "Species", "WD", "X1", "X2",
                                     "P50", "X3", "X4",
                                     "slope_P50", "X5", "X6")) %>% 
              select(Genus, Species, P50, slope_P50)) %>% 
  mutate(Genus = recode(Genus, "B." = "Bocoa", "D." = "Dicorynia", "E." = "Eperua",
                        "J." = "Jacaranda", "L." = "Licania", "P." = "Pradosia", "S." = "Sextonia",
                        "T." = "Tachigali", "V." = "Vouacapoua")) %>% 
  mutate(Genus = ifelse(Species == "globulifera", "Symphonia", Genus)) %>% 
  mutate(Genus = ifelse(Species == "sagotiana", "Eschweilera", Genus)) %>% 
  mutate(Genus = ifelse(Species == "persistens", "Lecythis", Genus)) %>% 
  left_join(group_by(growth, species) %>% summarise(gmax  = median(gmax, na.rm = T)) %>% separate(species, c("Genus", "Species"))) %>% 
  filter(!is.na(gmax)) 
isa <- read_tsv("data/fec12452-sup-0002-appendixs1.csv") %>% 
  dplyr::select(Binomial, Pi_tlp) %>% 
  group_by(Binomial) %>% 
  summarise(Pi_tlp = mean(Pi_tlp)) %>% 
  separate(Binomial, c("Genus", "Species")) %>% 
  left_join(group_by(growth, species) %>% summarise(gmax  = median(gmax, na.rm = T)) %>% separate(species, c("Genus", "Species"))) %>% 
  filter(!is.na(gmax)) 
isa2 <- read_tsv("data/fec12452-sup-0002-appendixs1.csv") %>% 
  separate(Binomial, c("Genus", "Species")) %>% 
  group_by(Genus) %>% 
  summarise(Pi_tlp = mean(Pi_tlp)) %>% 
  left_join(separate(growth, species, c("Genus", "Species"))  %>% group_by(Genus) %>% summarise(gmax  = median(gmax, na.rm = T))) %>% 
  filter(!is.na(gmax)) 
guillemot <- read_tsv("data/TAB_FINAL_GUILLEMOTetal.csv") %>% 
 left_join(group_by(growth, species) %>% summarise(gmax  = median(gmax, na.rm = T))) %>% 
    filter(!is.na(gmax))  
```

## Results

### All

*Individual regressions for each trait and each dataset.*

```{r, eval=F}
phylo <- ape::read.tree("save/ecoevo/phylogeny_full.tree")
p4d <- phylo4d(phylo,
               data.frame(species = gsub("_", " ", phylo$tip.label)) %>% 
                 left_join(select(fg, species) %>% mutate(Vleminckx = 1)) %>% 
                 mutate(Vleminckx = ifelse(is.na(Vleminckx), "no", "yes")) %>% 
                 dplyr::select(-species))
data.frame(species = gsub("_", " ", phylo$tip.label)) %>% 
                 left_join(select(fg, species) %>% mutate(Vleminckx = 1)) %>% 
                 mutate(Vleminckx = ifelse(is.na(Vleminckx), "no", "yes")) %>% 
                 dplyr::select(-species) %>% 
  group_by(Vleminckx) %>% 
  summarise(N = n())
fortify(p4d) %>% 
  mutate(species = gsub("_", " ", label)) %>%
  mutate(label = species) %>% 
  ggtree(aes(color = Vleminckx), layout="circular") + 
  geom_tiplab2(size = 2) +
  theme_tree(legend.position='right', legend.text = element_text(face = "italic")) +
  scale_alpha_manual("taxon", values = c(0.2, 1)) +
  scale_size_manual("taxon", values = c(1, 2)) +
  theme(legend.position = "bottom") + xlim(NA, 200)
```

```{r alldata}
all_data <- bind_rows(
  phd %>% 
    mutate(dataset = "Schmitt et al., (2020)", response = "Gmax_i", response_value = gmax) %>% 
    dplyr::select(dataset, species, response, response_value, SLA, LDMC, LT, LA, CC) %>% 
    reshape2::melt(c("dataset", "species", "response", "response_value"), variable.name = "trait", value.name = "trait_value"),
  fg %>% 
    mutate(dataset = "Vleminckx et al., (2021)", response = "Gmax_s", response_value = gmax) %>% 
    reshape2::melt(c("dataset", "species", "response", "response_value"), variable.name = "trait", value.name = "trait_value"),
  cam %>% 
    mutate(dataset = "Ziegler et al., (2019)", response = "Gmax_s", response_value = gmax) %>% 
    dplyr::select(-Species, -Genus) %>% 
    reshape2::melt(c("dataset", "species", "response", "response_value"), variable.name = "trait", value.name = "trait_value"),
  santi %>% 
    mutate(dataset = "Santiago et al., (2018)", response = "Gmax_s", response_value = gmax, species = paste(Species, Genus)) %>% 
    dplyr::select(-Species, -Genus) %>% 
    reshape2::melt(c("dataset", "species", "response", "response_value"), variable.name = "trait", value.name = "trait_value"),
  isa %>% 
    mutate(dataset = "Maréchaux et al., (2015)", response = "Gmax_s", response_value = gmax, species = paste(Species, Genus)) %>% 
    dplyr::select(-Species, -Genus) %>% 
    reshape2::melt(c("dataset", "species", "response", "response_value"), variable.name = "trait", value.name = "trait_value"),
  guillemot %>% 
    mutate(dataset = "Guillemot et al., (2022)", response = "Gmax_s", response_value = gmax) %>% 
    dplyr::select(-genus, -family, -leaf_habit) %>% 
    reshape2::melt(c("dataset", "species", "response", "response_value"), variable.name = "trait", value.name = "trait_value")
) %>% na.omit() %>% 
  filter(!(trait %in% c("slope_P50", "branchiness", "branch_Psi12", "branch_Psi88", "HSM_Psimd_Psi88",
                        "branch_vulnerability_slope", "SMleaf", "HSM_Psimd2010_Psi12", "HSM_Psimd2010_Psi50",
                        "HSM_Psimd_Psi12", "HSM_PiTLP_Psi12", "total_sapwood_bulk_elastic_modulus"))) %>% 
  mutate(trait = recode(trait,
                        "P50"	= "xylem pressure inducing 50% loss of hydraulic conductance",			
                        "TLP"	= "leaf water potential at turgor loss point",				
                        "LMA"	= "leaf mass per area",				
                        "leaf_size"	= "leaf area",				
                        "Leaf_N"	= "leaf nitrogen content",				
                        "Leaf_P"	= "leaf phosphorus content",				
                        "Wood_density"	= "trunk wood specific gravity",					
                        "Seed_mass"	= "seed mass",				
                        "max_height"	= "maximum tree height",				
                        "Pi_tlp"	= "leaf water potential at turgor loss point",
                        "sapwood_saturated_water_content"	= "sapwood saturated water content",				
                        "sapwood_osmotic_potential_at_full_turgor"	= "sapwood osmotic potential at full turgor",				
                        "sapwood_water_potential_at_turgor_loss_point"	= "sapwood water potential at turgor loss point",				
                        "total_sapwood_relative_water_content_at_turgor_loss_point"	= "total sapwood relative water content at turgor loss_point",				
                        "sapwood_capacitance_at_full_turgor"	= "sapwood capacitance at full turgor",				
                        "SLA"	= "specific leaf area",				
                        "LDMC"	= "leaf dry matter content",				
                        "LT"	= "leaf thickness",
                        "LA"	= "leaf area",				
                        "CC"	= "chlorophyll content",				
                        "Chlorophyll_content"	= "chlorophyll content",				
                        "Thickness"	= "leaf thickness",				
                        "Toughness"	= "leaf thoughness",				
                        "Leaf_Area"	= "leaf area",				
                        "C"	= "leaf carbon content",				
                        "N"	= "leaf nitrogen content",				
                        "13C"	= "leaf delta 13C",				
                        "Ca"	= "leaf calcium content",	
                        "P"	= "leaf phosphorus content",				
                        "K"	= "leaf potassium content",				
                        "Trunk_bark_thickness"	= "trunk bark thickness",				
                        "Sapwood_WSG"	= "sapwood specific gravity",				
                        "WSG"	= "roots wood specific gravity",				
                        "Diameter"	= "fine roots diameter",				
                        "SRL"	= "specific root length",				
                        "Tissue_density"	= "fine roots tissue density",				
                        "SRTA"	= "specific root tip abundance",				
                        "HSM_PiTLP_Psi88"	= "hydraulic safety margin"))
```

```{r allreg}
all_reg <- all_data %>% 
  filter(trait != "gmax") %>% 
  group_by(dataset, response, trait) %>% 
  do(lm(log(response_value) ~ scale(log(abs(trait_value))), data = .) %>% moderndive::get_regression_table()) %>% 
  filter(term != "intercept") %>% 
  dplyr::select(-term)
```

```{r allregfig, fig.height=10, fig.width=8}
all_reg %>% 
  mutate(p_value = ifelse(p_value == 0, 10^-4, p_value)) %>%
  ggplot(aes(x = trait, y = estimate, alpha = p_value < 0.01, fill  = log10(p_value), label = trait)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = .2) +
  geom_hline(yintercept = 0, col = "black") +
  coord_flip() +
  facet_grid(dataset ~ ., scale = "free_y", space = "free_y") +
  viridis::scale_fill_viridis(expression(log[10](p[value])), direction = -1) +
  ylab("Standard estimate with confidence interval") +
  theme(strip.text.y = element_text(angle = 0), axis.title.y = element_blank(), legend.position = "bottom") +
  scale_alpha_discrete(guide = "none")
```

```{r seltraits}
traits <- all_data %>% 
  filter(dataset == "Vleminckx et al., (2021)") %>% 
  group_by(species, trait) %>% 
  summarise(response_value = mean(response_value), trait_value = mean(trait_value)) %>% 
  reshape2::dcast(response_value + species ~ trait, value.var = "trait_value") %>% 
  na.omit() %>% 
  dplyr::select(-species, -`leaf delta 13C`, -response_value)
```

### Boosted regression trees

*Boosted regression trees for Vleminckx et al., (2021) data.*

```{r brt}
# brt_training <- train(
# 		gmax ~ .,
# 		data = traits,
# 		method = "gbm",
# 		verbose = FALSE,
# 		tuneGrid = expand.grid(
# 			n.trees = c(1000, 5000, 7500, 10000, 20000),
# 			interaction.depth = 1,
# 			shrinkage = c(0.0001, 0.0005, 0.001),
# 			n.minobsinnode = c(5, 10, 15, 20, 25)
# 			)
# 	)
# brt_training$bestTune
# n.trees: 10000, shrinkage: 0.0005, n.minobsinnode: 5
traits_brt <- gbm(
	gmax ~ .,
	distribution = "gaussian",
	data = traits,
	n.trees = 10000, # Change to value from best BRT tune
	interaction.depth = 1, # Change to value from best BRT tune
	n.minobsinnode = 5, # Change to value from best BRT tune
	shrinkage = 0.0005, # Change to value from best BRT tune
	bag.fraction = 0.5,
	cv.folds = 10
)
brt_tab <- tibble::as_tibble(gbm::summary.gbm(traits_brt, plotit = FALSE)) %>% 
  mutate(var = gsub("`", "", var)) %>% 
  rename(trait = var, `relative information` = "rel.inf")
```

```{r brtfig}
ggplot(brt_tab, aes(reorder(trait, `relative information`), 
                    `relative information`, 
                    label = paste(round(`relative information`), "%"))) +
  geom_col() +
  geom_text(nudge_y = 1.5) +
  coord_flip() +
  theme(axis.title.y = element_blank())
```

```{r brtpartial}
treezy::gg_partial_plot(traits_brt, 
                        vars = as_tibble(summary.gbm(traits_brt, 
                                                     plotit = FALSE))$var[1:6]) +
  theme_bw() +
  xlab("") + ylab(expression(Gmax[s]~cm.year^-1))
```

### Multiple regressions

*Multiple regressions for Vleminckx et al., (2021) data.*

```{r multlm}
# mult_traits <- glmulti::glmulti(gmax ~ ., report = FALSE,
#                                 data = rename_all(traits, list(~ gsub(" ", "_", .))) %>% 
#                                                     mutate_all(log) %>% na.omit(), 
#                                 fitfunc = lm, level = 1, plotty = FALSE, chunk = 20)
# save(mult_traits, file = "save/mult_traits.Rdata")
load("save/mult_traits.Rdata")
lmmult_tab <- coef(mult_traits) %>% 
  as.data.frame() %>% 
  rownames_to_column("trait") %>% 
  mutate(trait = gsub("_", " ", trait)) %>% 
  arrange(desc(Importance)) %>% 
  rename(estimate = Estimate,
         variance = "Uncond. variance",
         `variables importance` = Importance,
         `confidence interval` = "+/- (alpha=0.05)") %>% 
  select(-`Nb models`)
```

### Comparisons

```{r multab}
lmmult_tab %>% 
  left_join(brt_tab) %>% 
  kable()
```

```{r multbest}
mult_traits@objects[[1]] %>% 
  sjPlot::tab_model()
```

```{r multbestrelaimpo}
relaimpo::calc.relimp(mult_traits@objects[[1]])@lmg %>% 
  data.frame(trait = names(.), lmg = .) %>% 
  ggplot(aes(reorder(trait, lmg), lmg, label = paste(round(lmg*100), "%"))) +
  geom_col() +
  geom_text(nudge_y = .01) +
  coord_flip() +
  ylab(expression(R^2)) +
  theme(axis.title.y = element_blank())
```

```{r seltraitsrels}
sel_traits <- (filter(left_join(lmmult_tab, brt_tab), 
       `variables importance` > 0.5 | `relative information` > 10) %>% 
  filter(trait != "(Intercept)"))$trait
all_data %>%  
  filter(dataset == "Vleminckx et al., (2021)") %>% 
  filter(trait %in% sel_traits) %>% 
  ggplot(aes(trait_value, response_value)) +
  geom_point() +
  facet_wrap(~ trait, scales = "free") +
  scale_x_log10() + scale_y_log10() + geom_smooth(method = "lm") + 
  ylab(expression(Gmax[s]~cm.year^-1)) +
  theme(axis.title.x = element_blank(), legend.position = "bottom") +
  scale_color_discrete(guide = "none")
```

```{r seltraitsvars}
all_data %>% 
  filter(dataset == "Vleminckx et al., (2021)") %>% 
  filter(trait %in% c("leaf thoughness", "leaf nitrogen content", "leaf potassium content",
                      "sapwood specific gravity", "fine roots diameter", "specific root length")) %>% 
  left_join(dplyr::select(growth, Family, Genus, species) %>% unique()) %>% 
  dplyr::rename(gmax = response_value) %>% 
  reshape2::dcast(Family + Genus + species + gmax ~ trait, value.var = "trait_value") %>% 
  reshape2::melt(c("Family", "Genus", "species")) %>%
  # split(.$variable) %>% lapply(function(x) nlme::lme(log(abs(value)) ~  1, random=~1|Family/Genus, data = x)) %>%  sjPlot::tab_model(show.icc = F)
  group_by(variable) %>% 
  do(nlme::lme(log(abs(value)) ~  1, random=~1|Family/Genus, data = .) %>% 
       ape::varcomp(scale = F, cum = F) %>% 
       as.vector() %>% 
       data.frame(level = c("Family", "Genus", "Species"), variance = as.vector(.)) %>% 
       select(-`.`)) %>% 
  mutate(pct = variance / sum(variance)*100) %>% 
  mutate(level = factor(level, levels = c("Family", "Genus", "Species"))) %>% 
  mutate(pct_text = paste0(round(pct), "%")) %>% 
  mutate(pct_text = gsub("^0%", "", pct_text)) %>% 
  ggplot(aes(x = variable, y = pct, fill = level)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(y = pct, label = pct_text), col = "white", position = position_stack(vjust = .5)) +
  scale_fill_discrete(expression(sigma^2)) +
  coord_flip() +
  theme(axis.title = element_blank(), axis.line = element_blank(), axis.text.x = element_blank(), axis.ticks = element_blank())
```


```{r, eval=FALSE}
lmmult_tab %>% 
  left_join(brt_tab) %>% 
  select(estimate, `confidence interval`, `variables importance`, `relative information`) %>% 
  mutate_all(round, 2) %>% 
  write_tsv("~/Téléchargements/tab.tsv")
```


## SI: With Phylogeny

```{r phylov, eval=F}
library(V.PhyloMaker)
t <- all_data %>% 
  filter(dataset == "Vleminckx et al., (2021)") %>% 
  filter(trait %in% c("leaf thoughness", "leaf nitrogen content", "leaf potassium content",
                      "sapwood specific gravity", "fine roots diameter", "specific root length")) %>% 
  left_join(dplyr::select(growth, Family, Genus, species) %>% unique()) %>% 
  dplyr::rename(gmax = response_value) %>% 
  reshape2::dcast(Family + Genus + species + gmax ~ trait, value.var = "trait_value")
splist <- t %>% 
  dplyr::select(Family, Genus, species) %>% 
  unique() %>% 
  mutate(genus = Genus, family = Family) %>% 
  dplyr::select(species, genus, family)
tree <- phylo.maker(sp.list = splist, tree = GBOTB.extended, nodes = nodes.info.1, scenarios = "S3")
ape::write.tree(tree$scenario.3, "save/functional/phylogeny_Vleminckx.tree")
```

```{r, eval=F}
library(adephylo)
library(ape)
library(phylobase)
library(phylosignal)
t <- all_data %>% 
  filter(dataset == "Vleminckx et al., (2021)") %>% 
  filter(trait %in% c("leaf thoughness", "leaf nitrogen content", "leaf potassium content",
                      "sapwood specific gravity", "fine roots diameter", "specific root length")) %>% 
  left_join(dplyr::select(growth, Family, Genus, species) %>% unique()) %>% 
  dplyr::rename(gmax = response_value) %>% 
  reshape2::dcast(Family + Genus + species + gmax ~ trait, value.var = "trait_value")
phylo <- ape::read.tree("save/phylogeny_Vleminckx.tree")
p4d <- phylo4d(phylo,
               data.frame(species = gsub("_", " ", phylo$tip.label)) %>% 
                 left_join(t %>% group_by(Family, Genus, species) %>% mutate_all(abs) %>% mutate_all(log)) %>% 
                 dplyr::select(-species))   
t <- lapply(names(p4d@data)[-c(1,2)], function(t) {
  c <- as.data.frame(phyloCorrelogram(p4d, trait = t, ci.bs	= 10)$res)
  names(c) <- c("phylo_dist", "ll", "hh", "m")
  return(c)
})
names(t) <- names(p4d@data)[-c(1,2)]
bind_rows(t, .id = "trait") %>% 
  mutate(trait = gsub(".", " ", trait, fixed = T)) %>% 
  ggplot(aes(x = phylo_dist, col = trait, fill = trait)) +
  geom_ribbon(aes(ymin = ll, ymax = hh), alpha = 0.2, col = NA) +
  geom_line(aes(y = m), size = 1.5) +
  geom_hline(yintercept = 0, col = "black") +
  theme(legend.position = c(0.8, 0.7)) + 
  xlab("phylogenetic distance") + ylab("Correlation")
```