Comparison of EVI, NPP, ABI and ABI:NPP ratio across species and elevation

Figure

Code 
library(tidyverse)
library(emmeans)
library(multcompView)
library(kableExtra)
source("../scripts/aux.R")

abi <- read_csv("../data/abi.csv") |> 
  mutate(se = NA, sd = NA, variable = "abi") |> 
  rename(mean = abi)

npp <- read_csv("../data/npp_modis.csv") |> 
  rename(mean = npp) |> 
  mutate(se = NA, sd = NA, variable = "npp") |> 
  dplyr::select(year, sp_code, elev_code, sp_elev, mean, sd, se, variable) 

evi_landsat <- read_csv("../data/iv_landsat.csv") |> 
  filter(iv == "evi") |> 
  dplyr::select(year, sp_code, elev_code, mean, sd, se) |> 
  mutate(variable = "evi_landsat")

ratio <- read_csv("../data/ratio_abinpp.csv") |> 
  dplyr::select(year, sp_code, elev_code, sp_elev, mean = ratio) |> 
  mutate(se = NA, sd = NA, variable = "ratio") 


df <- bind_rows(
  abi, evi_landsat, npp, ratio) |> 
  mutate(elev_code = fct_relevel(elev_code, "low-Dec","low", "medium", "high")) |> 
  mutate(Specie = paste0("P. ", sp_code)) |> 
  rename(mean_y = mean, y_variable = variable)


## Summary statistics 
# Function to generate comparisons and perform letter-based grouping
generate_comparisons <- function(model, sp_code = NULL, elev_code = NULL, alpha = 0.001) {
  emm <- emmeans(model, pairwise ~ sp_code * elev_code, adjust = "bonferroni", alpha = alpha,
                 at = list(sp_code = sp_code, elev_code = elev_code))
  
  
  cld_df <- multcomp::cld(emm, Letters = letters, level = 1 - alpha)
  p_values <- summary(emm, infer = TRUE)$contrasts$'p.value'
  result_df <- cld_df |> 
    mutate(p_value = p_values[match(.group, cld_df$.group)])
  return(result_df)
}


df_avg <- df |> 
  group_by(y_variable, sp_code, elev_code, sp_elev) |> 
  summarise(mean = mean(mean_y), 
            sd = sd(mean_y), 
            se = sd/sqrt(length(mean_y))) |> ungroup()

variable <- "evi_landsat"

x <- df |> filter(y_variable == variable)  

lm_evi <- lm(mean_y ~ sp_code*elev_code, data = x) 

means_letters_lm_evi <- bind_rows(
  generate_comparisons(lm_evi, sp_code = "halepensis", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_evi, sp_code = "pinaster", elev_code = c("low", "low-Dec", "medium", "high")),
  generate_comparisons(lm_evi, sp_code = "nigra", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_evi, sp_code = "sylvestris", elev_code = c("low", "medium", "high"))
) |> 
  mutate(y_variable = variable)



variable <- "evi_landsat"

x <- df |> filter(y_variable == variable)  

lm_evi <- lm(mean_y ~ sp_code*elev_code, data = x) 

means_letters_lm_evi <- bind_rows(
  generate_comparisons(lm_evi, sp_code = "halepensis", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_evi, sp_code = "pinaster", elev_code = c("low", "low-Dec", "medium", "high")),
  generate_comparisons(lm_evi, sp_code = "nigra", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_evi, sp_code = "sylvestris", elev_code = c("low", "medium", "high"))
) |> 
  mutate(y_variable = variable)


variable <- "abi"

x <- df |> filter(y_variable == variable)  

lm_abi <- lm(mean_y ~ sp_code*elev_code, data = x) 

means_letters_lm_abi <- bind_rows(
  generate_comparisons(lm_abi, sp_code = "halepensis", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_abi, sp_code = "pinaster", elev_code = c("low", "low-Dec", "medium", "high")),
  generate_comparisons(lm_abi, sp_code = "nigra", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_abi, sp_code = "sylvestris", elev_code = c("low", "medium", "high"))
) |> 
  mutate(y_variable = variable)


variable <- "npp"

x <- df |> filter(y_variable == variable)  

lm_npp <- lm(mean_y ~ sp_code*elev_code, data = x) 

means_letters_lm_npp <- bind_rows(
  generate_comparisons(lm_npp, sp_code = "halepensis", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_npp, sp_code = "pinaster", elev_code = c("low", "low-Dec", "medium", "high")),
  generate_comparisons(lm_npp, sp_code = "nigra", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_npp, sp_code = "sylvestris", elev_code = c("low", "medium", "high"))
) |> 
  mutate(y_variable = variable)


variable <- "ratio"

x <- df |> filter(y_variable == variable)  

lm_ratio <- lm(mean_y ~ sp_code*elev_code, data = x) 

means_letters_lm_ratio <- bind_rows(
  generate_comparisons(lm_ratio, sp_code = "halepensis", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_ratio, sp_code = "pinaster", elev_code = c("low", "low-Dec", "medium", "high")),
  generate_comparisons(lm_ratio, sp_code = "nigra", elev_code = c("low", "medium", "high")),
  generate_comparisons(lm_ratio, sp_code = "sylvestris", elev_code = c("low", "medium", "high"))
) |> 
  mutate(y_variable = variable)


df_letters <- bind_rows(means_letters_lm_abi, means_letters_lm_evi, means_letters_lm_npp, means_letters_lm_ratio) |> ungroup()

out <- df_avg |> inner_join(df_letters) |> 
  mutate(pvalue = ifelse(p_value < 0.0001, "<0.0001", round(p_value, 3))) 


## Figure 
df_plot <- df |> 
  mutate(sp_code = fct_relevel(sp_code, "halepensis","pinaster", "nigra", "sylvestris")) |> 
  mutate(variable2 = case_when(
    y_variable == "abi" ~ "ABI~(g~C~m^2~year^{-1})", 
    y_variable == "npp" ~ "NPP[MODIS]~(g~C~m^2~year^{-1})",
    y_variable == "evi_landsat" ~ "EVI[Landsat]",
    y_variable == "ratio" ~ "ABI:NPP~ratio"
  )) |> 
  mutate(Species = paste0("P. ", sp_code)) |> 
  ungroup() 

out_plot <- out |> 
  mutate(sp_code = fct_relevel(sp_code, "halepensis","pinaster", "nigra", "sylvestris")) |> 
  mutate(variable2 = case_when(
    y_variable == "abi" ~ "ABI~(g~C~m^2~year^{-1})", 
    y_variable == "npp" ~ "NPP[MODIS]~(g~C~m^2~year^{-1})",
    y_variable == "evi_landsat" ~ "EVI[Landsat]",
    y_variable == "ratio" ~ "ABI:NPP~ratio"
  )) |> 
  ungroup() 


y_scales <- list(
  scale_y_continuous(limits = c(0, 0.5)),
  scale_y_continuous(limits = c(0,1000)),
  scale_y_continuous(limits = c(0, 600)), 
  scale_y_continuous(limits = c(0,.7))
)

letras_indicativo_plot <- data.frame(
  variable2 = c("EVI[Landsat]", 
                "NPP[MODIS]~(g~C~m^2~year^{-1})",
                "ABI~(g~C~m^2~year^{-1})",
                "ABI:NPP~ratio"), 
  label = c("a", "b", "c", "d"),
  x = c(-Inf, -Inf, -Inf, -Inf),  
  y = c(Inf, Inf, Inf, Inf)  # Inf to put at upper 
)

custom_position = position_dodge(width = 0.3)


out_plot2 <- out_plot |> 
  inner_join(letras_indicativo_plot, by="variable2") |> 
  mutate(Species = paste0("P. ", sp_code)) |> 
  mutate(Species = fct_relevel(Species, "P. halepensis","P. pinaster", "P. nigra", "P. sylvestris")) 


colours_elev <- c("low-Dec" = "#8c510a",
                  "low" ="#fc8d59",
                  "medium" = "#2166ac",
                  "high" = "#72A55A")

plot_avg_eviabinpp <- ggplot(data = out_plot2, 
                             aes(y = mean, x = Species, group = elev_code, colour = elev_code)) + 
  geom_jitter(data = df_plot,
              aes(y = mean_y, x = Species, 
                  group = elev_code, colour = elev_code), 
              alpha = 0.3, size = .75, stroke = 0.5, 
              position = custom_position, show.legend = FALSE) + 
  geom_line(position = custom_position, size = .6) +
  geom_errorbar(aes(ymin=mean-1.96*se, ymax=mean+1.96*se), 
                colour = "black", width=.2, lwd = .7, position = custom_position) +  
  geom_point(size = 2.2, position = custom_position, 
             aes(fill = elev_code, group = elev_code),
             shape = 21, colour = "black", stroke = 0.5, show.legend = FALSE) +
  facet_wrap(~factor(variable2, c("EVI[Landsat]", 
                                  "NPP[MODIS]~(g~C~m^2~year^{-1})",
                                  "ABI~(g~C~m^2~year^{-1})",
                                  "ABI:NPP~ratio")),
             scales = "free_y", 
             labeller = label_parsed, 
             strip.position = "left",
             ncol = 1) +
  theme_bw() + 
  theme(
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.placement = "outside",
    legend.position = c(0.85, 0.84),
    legend.background = element_blank(),
    axis.text.x = element_text(face = "italic"),
    strip.text = element_text(face = "bold", size = 9)
  ) + 
  scale_colour_manual(values = colours_elev, name = "", guide = guide_legend(override.aes = list(linewidth = 2))) +
  scale_fill_manual(values = colours_elev, name = "") + 
  xlab("") + 
  ylab("") + 
  ggh4x::facetted_pos_scales(y = y_scales) +
  geom_text(aes(label = label, x = x, y = y), 
            vjust = 1.2, hjust = -.5, size = 6, 
            fontface = "bold", inherit.aes = FALSE)
# custom scales

plot_avg_eviabinpp

Comparisons of the annual EVI (Enhanced Vegetation Index) (1986-2022), annual NPP (Net Primary Productivity (2001-2022), ABI (Aboveground Biomass Increment) (1952-2022) and ABI:NPP ratio (2001-2022) among species and elevation. Each small dot represents one mean annual value. Big dots represent average values. Error bars indicate 95 % confidence interval.
Code 
out |> 
  ungroup() |> 
  dplyr::select(-sp_elev) |> 
  dplyr::select(-p_value) |> 
  arrange(y_variable, sp_code, elev_code) |> 
  kbl(digits = c(0,0,0,3,3,3, 3,3,0,3,3,0,5), 
      caption = "Posthoc differences within species for average ABI, EVI, NPP") |> 
  kable_paper("hover", full_width = F)
Posthoc differences within species for average ABI, EVI, NPP
y_variable sp_code elev_code mean sd se emmean SE df lower.CL upper.CL .group pvalue
abi halepensis low 37.530 36.400 5.551 37.530 12.459 694 -3.641 78.700 a 0.088
abi halepensis medium 76.225 64.156 9.900 76.225 12.606 694 34.568 117.883 a 0.088
abi halepensis high 55.754 33.091 4.136 55.754 10.212 694 22.008 89.501 a 0.088
abi nigra low 78.402 40.264 5.287 78.402 10.727 694 42.953 113.852 a 1
abi nigra medium 68.667 44.275 5.764 68.667 10.636 694 33.519 103.814 a 1
abi nigra high 142.916 96.256 13.348 142.916 11.329 694 105.478 180.355 b <0.0001
abi pinaster low-Dec 48.876 27.511 3.494 48.876 10.375 694 14.590 83.163 a 0.281
abi pinaster low 77.250 45.168 5.399 77.250 9.765 694 44.982 109.518 a 0.281
abi pinaster medium 167.705 115.563 15.045 167.705 10.636 694 132.557 202.852 b <0.0001
abi pinaster high 257.671 158.973 25.789 257.671 13.253 694 213.876 301.467 c <0.0001
abi sylvestris low 87.516 51.346 6.861 87.516 10.917 694 51.440 123.593 a <0.0001
abi sylvestris medium 179.616 102.904 14.003 179.616 11.117 694 142.877 216.355 b <0.0001
abi sylvestris high 226.002 146.879 20.772 226.002 11.554 694 187.822 264.182 c 0.012
evi_landsat halepensis low 0.186 0.029 0.005 0.186 0.005 467 0.168 0.203 a 0.04
evi_landsat halepensis medium 0.204 0.028 0.005 0.204 0.005 467 0.187 0.222 b 0.318
evi_landsat halepensis high 0.192 0.025 0.004 0.192 0.005 467 0.175 0.210 ab 1
evi_landsat nigra low 0.231 0.024 0.004 0.231 0.005 467 0.214 0.249 ab 0.074
evi_landsat nigra medium 0.241 0.022 0.004 0.241 0.005 467 0.223 0.258 b 0.001
evi_landsat nigra high 0.214 0.033 0.005 0.214 0.005 467 0.197 0.232 a 0.618
evi_landsat pinaster low-Dec 0.215 0.024 0.004 0.215 0.005 467 0.198 0.233 a 0.018
evi_landsat pinaster low 0.238 0.033 0.005 0.238 0.005 467 0.220 0.255 b <0.0001
evi_landsat pinaster medium 0.353 0.034 0.006 0.353 0.005 467 0.335 0.370 c <0.0001
evi_landsat pinaster high 0.344 0.063 0.010 0.344 0.005 467 0.327 0.362 c <0.0001
evi_landsat sylvestris low 0.358 0.026 0.004 0.358 0.005 467 0.341 0.376 a <0.0001
evi_landsat sylvestris medium 0.416 0.030 0.005 0.416 0.005 467 0.399 0.434 b <0.0001
evi_landsat sylvestris high 0.424 0.027 0.005 0.424 0.005 467 0.406 0.441 b <0.0001
npp halepensis low 291.700 56.247 11.992 291.700 13.978 273 245.203 338.197 a <0.0001
npp halepensis medium 637.305 56.993 12.151 637.305 13.978 273 590.808 683.801 c <0.0001
npp halepensis high 509.632 59.907 12.772 509.632 13.978 273 463.135 556.128 b <0.0001
npp nigra low 690.250 68.047 14.508 690.250 13.978 273 643.753 736.747 a 0.452
npp nigra medium 718.736 70.079 14.941 718.736 13.978 273 672.240 765.233 a 0.452
npp nigra high 720.136 71.862 15.321 720.136 13.978 273 673.640 766.633 a 0.452
npp pinaster low-Dec 580.368 60.299 12.856 580.368 13.978 273 533.872 626.865 a 0.003
npp pinaster low 650.586 68.523 14.609 650.586 13.978 273 604.090 697.083 b <0.0001
npp pinaster medium 759.386 77.945 16.618 759.386 13.978 273 712.890 805.883 c <0.0001
npp pinaster high 741.209 68.329 14.568 741.209 13.978 273 694.713 787.706 c <0.0001
npp sylvestris low 913.482 68.459 14.596 913.482 13.978 273 866.985 959.978 c <0.0001
npp sylvestris medium 845.050 51.023 10.878 845.050 13.978 273 798.553 891.547 b <0.0001
npp sylvestris high 636.918 69.313 14.778 636.918 13.978 273 590.422 683.415 a 0.002
ratio halepensis low 0.136 0.094 0.020 0.136 0.017 240 0.080 0.192 a 1
ratio halepensis medium 0.112 0.067 0.016 0.112 0.019 240 0.050 0.174 a 1
ratio halepensis high 0.089 0.046 0.010 0.089 0.017 240 0.033 0.145 a 1
ratio nigra low 0.113 0.057 0.012 0.113 0.017 240 0.056 0.170 a 0.12
ratio nigra medium 0.060 0.038 0.009 0.060 0.019 240 -0.003 0.124 a 0.12
ratio nigra high 0.207 0.085 0.019 0.207 0.018 240 0.148 0.265 b <0.0001
ratio pinaster low-Dec 0.086 0.038 0.008 0.086 0.017 240 0.029 0.144 a 1
ratio pinaster low 0.100 0.049 0.011 0.100 0.017 240 0.043 0.157 ab 0.079
ratio pinaster medium 0.164 0.090 0.022 0.164 0.019 240 0.100 0.227 b <0.0001
ratio pinaster high 0.296 0.155 0.035 0.296 0.018 240 0.238 0.355 c 0.017
ratio sylvestris low 0.113 0.029 0.007 0.113 0.019 240 0.051 0.175 a <0.0001
ratio sylvestris medium 0.242 0.101 0.024 0.242 0.019 240 0.180 0.304 b <0.0001
ratio sylvestris high 0.493 0.076 0.018 0.493 0.019 240 0.431 0.554 c <0.0001

Comparison table