Análise Florística Chaco

Conjunto de dados

Este é um documento anexo que pertence ao artigo intitulado “Spatial structuring of woody biodiversity in the South-Central South America”, de autoria de Ângela Lúcia Bagnatori Sartori, Geraldo Alves Damasceno Junior, Emerson Pereira da Silva, Giovana Amaral Umar, Isabela Oliveira Ferreira, Mateus César Araújo Pestana, Maxwell da Rosa Oliveira, Thomaz Ricardo Favreto Sinani e Flávio Macedo Alves

Estas análises são baseadas em 112 referencias bibliográficas e base de dados dos autores realizados em remanescentes de Chaco e ecorregiões adjacentes. Nós compilamos dados de ocorrência de 1.623 táxons da flora em 307 localidades, que podem ser acessados no Material Suplementar “matrix_original.xlsx” e as análises podem ser replicadas seguindo este roteiro.

1. Carregamento de pacotes

Nesta etapa, são carregados os pacotes necessários para a manipulação de dados, análises espaciais e ordenação multivariada. A semente aleatória (set.seed) é definida para garantir a reprodutibilidade das análises que envolvem permutações.

library(readxl)
library(tidyverse)
library(vegan)
library(betapart)
library(venn)
library(ggrepel)
library(dendextend)
library(geosphere)
library(indicspecies)
library(rgbif)
library(sf)
library(terra)
library(circlize)
library(ggalluvial)
library(pheatmap)
library(writexl)
library(ape)
library(ggnewscale)
library(svglite)

set.seed(123)

Anotações: A correta execução desta etapa permite o fluxo ininterrupto das análises ecológicas, carregando todas as funções necessárias na memória do R. Se houver erro, os pacotes devem ser previamente instalados.

2. Carregamento e estruturação dos dados

2.1. Localidades alvo para caracterização

Definição das localidades alvo (áreas de interesse primário) que guiarão as análises e a identificação do domínio florístico correspondente.

  # Inserir no vetor as áreas que gostaria de obter a caracterização
  local_target <- c("CHA_PANT_01",
                  "CHA_PANT_02",
                  "CHA_PANT_03",
                  "CHA_PANT_04",
                  "CHA_PANT_07",
                  "CHA_PANT_08",
                  "CHA_PANT_09",
                  "CHA_PANT_10")

Anotações: Estas áreas servirão como marcadores nas ordenações futuras, permitindo rastrear a qual agrupamento estatístico elas pertencem e com quais regiões elas se assemelham.

2.2. Área geográfica dos domínios na delimitação dos dados

Definição das coordenadas geográficas (Bounding Box) que delimitam a área de estudo.

# Inserir coordenadas de abrangência dos dados (Bounding Box)
# c(Long_Min, Long_Max, Lat_Min, Lat_Max)
boundary_area <- c(-67.720090, -41.072987, -41.011410, -4.714250)

Anotações: Este recorte espacial atua como um filtro de segurança para garantir que a busca de dados em bases abertas (como o GBIF) retorne apenas registros georreferenciados dentro da região continental de interesse.

2.3. Base de dados

Importação e estruturação da matriz de dados de presença/ausência. Aplica-se uma filtragem iterativa matemática para remover locais com menos de 5 espécies e espécies que ocorrem em apenas um local (singletons).

# 2.1. Importar a planilha original
input_matrix <- read_excel("matrix_original.xlsx", col_types = "text")

# 2.2. Processar e pivotar os dados
pivot_matrix <- input_matrix %>%
  select(locality, taxa) %>% 
  distinct(locality, taxa) %>% 
  mutate(presence = 1) %>% 
  pivot_wider(names_from = taxa, values_from = presence, values_fill = 0)

# 2.3. Ajustar os nomes das linhas (rownames) e converter para matriz numérica
  matriz_df <- as.data.frame(pivot_matrix)
  rownames(matriz_df) <- matriz_df$locality
  matriz_df <- matriz_df %>% select(-locality)
  matrix_PA <- as.matrix(matriz_df)

# 2.4 Limpeza Recursiva da Matriz (manter localidades >= 5 e espécies > 1)
  matrix_PA_clean <- matrix_PA
  initial_row <- nrow(matrix_PA_clean)
  initial_col <- ncol(matrix_PA_clean)
  iteration <- 0

  while(any(rowSums(matrix_PA_clean) < 5) | any(colSums(matrix_PA_clean) <= 1)) {
    iteration <- iteration + 1
    
    valid_rows <- rowSums(matrix_PA_clean) >= 5
    matrix_PA_clean <- matrix_PA_clean[valid_rows, , drop = FALSE]
    
    valid_cols <- colSums(matrix_PA_clean) > 1
    matrix_PA_clean <- matrix_PA_clean[, valid_cols, drop = FALSE]
    
    if(nrow(matrix_PA_clean) == 0 | ncol(matrix_PA_clean) == 0) stop("Filtro muito rigoroso. Matriz vazia!")
  }

  cat("--- RELATÓRIO DE ESTABILIZAÇÃO DA MATRIZ ---\n")

## --- RELATÓRIO DE ESTABILIZAÇÃO DA MATRIZ ---

  cat("O R atingiu a estabilidade perfeita em", iteration, "ciclos de limpeza.\n\n")

## O R atingiu a estabilidade perfeita em 2 ciclos de limpeza.

  cat("LOCALIDADES (Linhas):\n - Originais:", initial_row, "\n - Mantidas: ", nrow(matrix_PA_clean), "\n - Removidas:", initial_row - nrow(matrix_PA_clean), "\n\n")

## LOCALIDADES (Linhas):
##  - Originais: 308 
##  - Mantidas:  305 
##  - Removidas: 3

  cat("ESPÉCIES (Colunas):\n - Originais:", initial_col, "\n - Mantidas: ", ncol(matrix_PA_clean), "\n - Removidas:", initial_col - ncol(matrix_PA_clean), "\n")

## ESPÉCIES (Colunas):
##  - Originais: 1587 
##  - Mantidas:  1061 
##  - Removidas: 526

Anotações: A filtragem iterativa estabilizou a matriz de dados, removendo vieses causados por subamostragem e espécies raras pontuais. O resultado é uma matriz ecologicamente robusta para as análises multivariadas.

2.4. Variáveis Ambientais

Importação da planilha contendo os metadados abióticos e as variáveis edafoclimáticas.

  raw_data_env <- read_excel("env_variables.xlsx")

Anotações: Estes dados ambientais ficam armazenados na memória para subsidiar etapas futuras de ordenação restrita (como RDA ou CCA).

3. Análise Exploratória de Dados

3.1. Distribuição Espacial dos Registros Originais

Visualização do esforço amostral bruto (pré-curadoria) plotado sobre as ecorregiões originais. A técnica de alpha blending (transparência cumulativa) é aplicada para evidenciar as coordenadas com alta densidade de amostragem sem a necessidade de clusterização numérica.

cat("\nPreparando dados espaciais originais e carregando ecorregiões...\n")

## 
## Preparando dados espaciais originais e carregando ecorregiões...

df_pontos_originais <- input_matrix %>%
  select(locality, longitud, latitud) %>%
  distinct() %>%
  mutate(long = as.numeric(gsub(",", ".", longitud)), lat = as.numeric(gsub(",", ".", latitud)))

ecorregioes_originais <- input_matrix %>%
  select(Ecorregiao_Olson = TNC_Ecoregion_Olson) %>% filter(!is.na(Ecorregiao_Olson)) %>% distinct() %>% pull()

url_wwf <- "[https://c402277.ssl.cf1.rackcdn.com/publications/15/files/original/official_teow.zip](https://c402277.ssl.cf1.rackcdn.com/publications/15/files/original/official_teow.zip)"
zip_dest <- "wwf_teow.zip"
dir_dest <- "wwf_teow_extracted"

if(!dir.exists(dir_dest)) {
  download.file(url_wwf, destfile = zip_dest, mode = "wb")
  unzip(zip_dest, exdir = dir_dest)
}
caminho_shp <- list.files(dir_dest, pattern = "\\.shp$", recursive = TRUE, full.names = TRUE)

sf::sf_use_s2(FALSE) 

shape_olson <- st_read(caminho_shp[1], quiet = TRUE) %>% st_transform(4326)

shape_eco_originais <- shape_olson %>%
  filter(ECO_NAME %in% ecorregioes_originais) %>%
  group_by(ECO_NAME) %>% summarise(geometry = st_union(geometry))

sf::sf_use_s2(TRUE)

# --- PADRONIZAÇÃO ESTÉTICA DO MAPA 1 (Baseado no MST) ---
cores_atuais <- c("#35978f", "#bf812d", "#8c510a", "#dfc27d", "#f6e8c3", "#f5f5f5", "#80cdc1", "#01665e", "#c7eae5")
ecos_alfabetica <- sort(unique(shape_eco_originais$ECO_NAME))
names(cores_atuais) <- ecos_alfabetica[1:length(cores_atuais)]
ordem_legenda_brbg <- names(cores_atuais)[c(3, 2, 4, 5, 6, 9, 7, 1, 8)]

mapa_amostras_iniciais <- ggplot() +
  # Ecorregiões sem borda (color = NA) e com transparência (alpha = 0.7)
  geom_sf(data = shape_eco_originais, aes(fill = ECO_NAME), alpha = 0.7, color = NA) +
  # Pontos estruturados com borda branca (shape = 21) para melhor contraste
  geom_point(data = df_pontos_originais, aes(x = long, y = lat), size = 2.5, shape = 21, fill = "black", color = "white", stroke = 0.4, alpha = 0.6) +
  scale_fill_manual(values = cores_atuais, breaks = ordem_legenda_brbg) +
  coord_sf(xlim = boundary_area[1:2], ylim = boundary_area[3:4], expand = FALSE) +
  theme_minimal() +
  labs(title = "Distribuição Espacial das Localidades Amostradas", subtitle = "Esforço amostral bruto", fill = "Ecorregião (Olson):", x = "Longitude", y = "Latitude") +
  theme(legend.position = "right", plot.title = element_text(face = "bold", size = 14), panel.grid.minor = element_blank())

print(mapa_amostras_iniciais)

3.2. Rank-Frequency

Cálculo da frequência de ocorrência das espécies para visualizar a distribuição e a dominância de táxons comuns e raros na região de estudo.

frequency_sp <- colSums(matrix_PA_clean)
frequency_order <- sort(frequency_sp, decreasing = TRUE)
percent_ocur <- (frequency_order / nrow(matrix_PA_clean)) * 100

plot(percent_ocur, type = "l", lwd = 2, col = "blue", main = "Curva de Frequencia-Ordem das Angiospermas", xlab = "Ordenacao das Especies", ylab = "Frequencia de Ocorrencia (%)")
abline(h = 5, col = "red", lty = 2)

cat("A espécie mais comum ocorre em", round(percent_ocur[1], 1), "% das áreas.\n")

## A espécie mais comum ocorre em 36.1 % das áreas.

cat("Número de espécies hiper-raras (<5% de ocorrência):", sum(percent_ocur < 5), "\n")

## Número de espécies hiper-raras (<5% de ocorrência): 792

Anotações: O gráfico em forma de “J invertido” é típico de comunidades biológicas regionais. A enorme proporção de espécies raras (abaixo da linha vermelha de 5%) sugere um alto número de especialistas locais e uma forte substituição de táxons (turnover) na América do Sul.

3.3. Heatmap

Visualização da matriz de presença/ausência da flora através de um mapa de calor (Heatmap).

matrix_image <- as.matrix(matrix_PA_clean)
matrix_image <- t(matrix_image)[, nrow(matrix_image):1]

image(x = 1:ncol(matrix_PA_clean), y = 1:nrow(matrix_PA_clean), z = matrix_image, col = c("white", "#cb181d"), main = "Mapa de Calor da Matriz", xlab = "Espécies", ylab = "Localidades")

# 1. Calcular as distâncias florísticas (Jaccard)
dist_loc <- vegdist(matrix_PA_clean, method = "jaccard") 
dist_sp <- vegdist(t(matrix_PA_clean), method = "jaccard")

tree_loc <- hclust(dist_loc, method = "ward.D2")
tree_sp <- hclust(dist_sp, method = "ward.D2")

pheatmap(as.matrix(matrix_PA_clean), color = c("#f7f7f7", "#cb181d"), cluster_rows = tree_loc, cluster_cols = tree_sp, show_rownames = FALSE, show_colnames = FALSE, border_color = NA, legend_breaks = c(0, 1), legend_labels = c("Ausência", "Presença"), main = "Heatmap Florístico Clusterizado")

Anotações: O Heatmap permite uma visualização rápida da estrutura bruta dos dados, evidenciando lacunas e agrupamentos prévios de localidades antes mesmo da aplicação de testes estatísticos avançados. Os pontos vermelhos no heatmap representam presença (1) e a parte branca representa ausência (0); é possível detectar agrupamentos prévios das espécies nas localidades.

3.4. Rarefaction

Curva de acumulação de espécies calculada por meio de permutações aleatórias das localidades amostradas baseadas em incidência.

regional_rare <- specaccum(matrix_PA_clean, method = "random", permutations = 100)
plot(regional_rare, ci.type = "poly", col = "darkgreen", lwd = 2, ci.lty = 0, ci.col = "lightgreen", main = "Curva de Acumulacao Regional", xlab = "Esforco Amostral", ylab = "Riqueza Acumulada")
grid()

Anotações: A estabilização (achatamento) da curva confirma que o esforço amostral bibliográfico compilado neste estudo foi suficiente para capturar e representar com segurança a diversidade florística regional.

3.5. Diagrama de Venn/Euler

Criação de um diagrama de conjuntos baseado nas classificações de biomas atribuídas originalmente às áreas na literatura base.

df_locality <- tibble(locality = rownames(matrix_PA_clean))
df_biome_align <- df_locality %>% left_join(input_matrix %>% select(locality, domain_biome) %>% distinct(), by = "locality")

biome_vec <- df_biome_align$domain_biome
biome_uniq <- unique(na.omit(biome_vec))
nomes_biomas_puros <- biome_uniq 

list_venn_dinamic <- lapply(biome_uniq, function(biome_actual) {
  row_biome <- matrix_PA_clean[biome_vec == biome_actual, , drop = FALSE]
  present_sp <- colnames(row_biome)[colSums(row_biome) > 0]
  return(present_sp)
})
names(list_venn_dinamic) <- biome_uniq
lista_limpa_venn <- list_venn_dinamic
biome_name <- names(list_venn_dinamic)
biome_rich <- sapply(list_venn_dinamic, length)
names(list_venn_dinamic) <- paste0(biome_name, "\n(", biome_rich, " spp.)")

# AQUI ESTÁ A CORREÇÃO: Inserção do ilabels = "counts"
venn(list_venn_dinamic, zcolor = "style", opacity = 0.4, ilabels = "counts", ilcs = 0.5, sncs = 0.8, box = FALSE, borders = FALSE)

todas_spp <- unique(unlist(lista_limpa_venn))
df_sp_venn <- data.frame(Especie = todas_spp)
for (b in nomes_biomas_puros) { df_sp_venn[[b]] <- df_sp_venn$Especie %in% lista_limpa_venn[[b]] }

df_combinacoes <- df_sp_venn %>% select(-Especie) %>% group_by(across(everything())) %>% summarise(Numero_Especies = n(), .groups = "drop")

nomear_categoria <- function(linha) {
  biomas_presentes <- nomes_biomas_puros[as.logical(linha)]
  n_pres <- length(biomas_presentes)
  n_tot <- length(nomes_biomas_puros)
  if(n_pres == 1) return(paste("Exclusiva:", biomas_presentes))
  else if (n_pres == n_tot) return("Intersecção: TODOS os Biomas")
  else return(paste("Intersecção:", paste(biomas_presentes, collapse = " + ")))
}

df_combinacoes$Categoria <- apply(df_combinacoes[, nomes_biomas_puros], 1, nomear_categoria)
df_tabela_venn <- df_combinacoes %>% arrange(desc(Numero_Especies)) %>% select(Categoria, Numero_Especies, everything())

cat("\n")

knitr::kable(df_tabela_venn, caption = "Tabela Resumo: Compartilhamento e Exclusividade.")

Tabela Resumo: Compartilhamento e Exclusividade.

Categoria	Numero_Especies	Cerrado	Pantanal	Chaco	Mata Atlântica	Chiquitano
Exclusiva: Chaco	230	FALSE	FALSE	TRUE	FALSE	FALSE
Intersecção: Cerrado + Pantanal	193	TRUE	TRUE	FALSE	FALSE	FALSE
Exclusiva: Cerrado	116	TRUE	FALSE	FALSE	FALSE	FALSE
Intersecção: Cerrado + Pantanal + Mata Atlântica	82	TRUE	TRUE	FALSE	TRUE	FALSE
Exclusiva: Pantanal	79	FALSE	TRUE	FALSE	FALSE	FALSE
Intersecção: Pantanal + Chaco	79	FALSE	TRUE	TRUE	FALSE	FALSE
Intersecção: Cerrado + Pantanal + Chaco + Mata Atlântica	74	TRUE	TRUE	TRUE	TRUE	FALSE
Intersecção: Cerrado + Mata Atlântica	40	TRUE	FALSE	FALSE	TRUE	FALSE
Intersecção: Cerrado + Pantanal + Chaco	34	TRUE	TRUE	TRUE	FALSE	FALSE
Intersecção: Pantanal + Chaco + Mata Atlântica	18	FALSE	TRUE	TRUE	TRUE	FALSE
Intersecção: Cerrado + Chaco	17	TRUE	FALSE	TRUE	FALSE	FALSE
Intersecção: Cerrado + Chaco + Mata Atlântica	16	TRUE	FALSE	TRUE	TRUE	FALSE
Intersecção: TODOS os Biomas	15	TRUE	TRUE	TRUE	TRUE	TRUE
Intersecção: Chaco + Mata Atlântica	14	FALSE	FALSE	TRUE	TRUE	FALSE
Intersecção: Pantanal + Mata Atlântica	12	FALSE	TRUE	FALSE	TRUE	FALSE
Intersecção: Pantanal + Chaco + Chiquitano	8	FALSE	TRUE	TRUE	FALSE	TRUE
Intersecção: Chaco + Chiquitano	7	FALSE	FALSE	TRUE	FALSE	TRUE
Intersecção: Cerrado + Pantanal + Mata Atlântica + Chiquitano	7	TRUE	TRUE	FALSE	TRUE	TRUE
Intersecção: Cerrado + Pantanal + Chaco + Chiquitano	6	TRUE	TRUE	TRUE	FALSE	TRUE
Exclusiva: Mata Atlântica	5	FALSE	FALSE	FALSE	TRUE	FALSE
Intersecção: Cerrado + Chaco + Mata Atlântica + Chiquitano	3	TRUE	FALSE	TRUE	TRUE	TRUE
Intersecção: Cerrado + Chiquitano	2	TRUE	FALSE	FALSE	FALSE	TRUE
Intersecção: Cerrado + Pantanal + Chiquitano	2	TRUE	TRUE	FALSE	FALSE	TRUE
Intersecção: Pantanal + Chaco + Mata Atlântica + Chiquitano	1	FALSE	TRUE	TRUE	TRUE	TRUE
Intersecção: Cerrado + Mata Atlântica + Chiquitano	1	TRUE	FALSE	FALSE	TRUE	TRUE

Anotações: O diagrama ilustra o número de espécies consideradas exclusivas de cada bioma e dimensiona o elevado volume de táxons compartilhados nas zonas de transição (ecótonos).

3.6. Detecção de Outliers via Z-Score

Cálculo do Z-Score para identificar localidades com desvio multivariado extremo perante a comunidade regional, seguindo os critérios de McCune & Grace.

square_matrix_raw <- as.matrix(dist_loc) 
distm_raw <- rowMeans(square_matrix_raw)
z_scores_raw <- (distm_raw - mean(distm_raw)) / sd(distm_raw)
raw_rich <- rowSums(matrix_PA_clean)

df_diagnostic <- data.frame(Localidade = rownames(matrix_PA_clean), Riqueza = raw_rich, Z_Score = z_scores_raw)
df_diagnostic <- df_diagnostic %>% mutate(Status = ifelse(Z_Score > 2.0 | Riqueza < 5, "Excluída", "Mantida"))

ggplot(df_diagnostic, aes(x = Riqueza, y = Z_Score, color = Status)) +
  geom_point(alpha = 0.7, size = 2.5) +
  geom_hline(yintercept = 2.0, linetype = "dashed", linewidth = 1, color = "darkred") +
  annotate("text", x = max(df_diagnostic$Riqueza) * 0.5, y = 2.25, label = "Corte Multivariado (Z > 2.0)", color = "darkred", fontface = "bold") +
  geom_vline(xintercept = 4.5, linetype = "dotted", linewidth = 1, color = "blue") +
  theme_bw() + labs(title = "Diagnostico de Curadoria da Matriz", x = "Riqueza Floristica", y = "Z-Score") +
  scale_color_manual(values = c("Excluída" = "red", "Mantida" = "grey40")) + theme(legend.position = "bottom")

Anotações: O gráfico cruza a riqueza de espécies com o distanciamento florístico (Z-score). Localidades pontuadas acima da linha de corte vermelha (\(Z > 2.0\)) são anômalas à tendência central dos dados e precisam ser tratadas.

4. Tratamento dos Dados

4.1. Retirada de Outliers e Singletons

Remoção matemática das localidades identificadas como outliers estatísticos com base nas suas distâncias médias de Jaccard.

outliers_multi <- which(z_scores_raw > 2.0)
  if(length(outliers_multi) > 0) {
    info_outliers <- data.frame(Localidade = names(outliers_multi), Z_Score = round(z_scores_raw[outliers_multi], 2))
    print(knitr::kable(info_outliers[order(-info_outliers$Z_Score), ]))
  }

## 
## 
## |              |Localidade    | Z_Score|
## |:-------------|:-------------|-------:|
## |CHA_SEC_11    |CHA_SEC_11    |    2.65|
## |CERDAO_CER_33 |CERDAO_CER_33 |    2.47|
## |CHA_SEC_15    |CHA_SEC_15    |    2.18|
## |CERDAO_CER_35 |CERDAO_CER_35 |    2.14|
## |CERDAO_CER_26 |CERDAO_CER_26 |    2.12|
## |CHA_HYGRO_03  |CHA_HYGRO_03  |    2.08|
## |CILIAR_CER_17 |CILIAR_CER_17 |    2.02|

  matrix_PA_clean_inter <- matrix_PA_clean[z_scores_raw <= 2.0, ]

Anotações: A remoção destas localidades atípicas é fundamental para estabilizar os dados, garantindo que o estresse do NMDS seja aceitável e que as ordenações ecológicas não sofram distorções de escala.

4.2. Conferência final de baixa riqueza e táxons raros

Nova rodada iterativa de limpeza biológica para garantir que a remoção dos outliers não tenha gerado novos “singletons” secundários na matriz de dados.

cat("Limpeza Biológica Recursiva (Riqueza < 5 e Frequência <= 1) ---\n")

## Limpeza Biológica Recursiva (Riqueza < 5 e Frequência <= 1) ---

  iteration_bio <- 1
  continue_clean_bio <- TRUE
  matrix_PA_clean_recur <- matrix_PA_clean_inter
  
  while(continue_clean_bio) {
    cat("[Rodada", iteration_bio, "] ")
    curren_rich <- rowSums(matrix_PA_clean_recur)
    curren_freq <- colSums(matrix_PA_clean_recur)
    
    loc_remove <- which(curren_rich < 5)
    spp_remove <- which(curren_freq <= 1)
    
    if(length(loc_remove) == 0 && length(spp_remove) == 0) {
      cat("A matriz estabilizou! Nenhuma localidade < 5 spp ou espécie única detectada.\n")
      continue_clean_bio <- FALSE
    } else {
      cat("Removendo", length(loc_remove), "localidades e", length(spp_remove), "espécies...\n")
      matrix_PA_clean_recur <- matrix_PA_clean_recur[curren_rich >= 5, curren_freq > 1]
      cat("   -> Novo tamanho parcial:", nrow(matrix_PA_clean_recur), "localidades e", ncol(matrix_PA_clean_recur), "espécies.\n")
      iteration_bio <- iteration_bio + 1
      
      if(iteration_bio > 4) {
        cat("\n[!] ALERTA: Limite de segurança atingido. Parando.\n")
        continue_clean_bio <- FALSE
      }
    }
  }

## [Rodada 1 ] Removendo 0 localidades e 16 espécies...
##    -> Novo tamanho parcial: 298 localidades e 1045 espécies.
## [Rodada 2 ] A matriz estabilizou! Nenhuma localidade < 5 spp ou espécie única detectada.

  matrix_run <- matrix_PA_clean_recur
  
  cat("\n=== RESUMO DA CURADORIA DE DADOS ===\n")

## 
## === RESUMO DA CURADORIA DE DADOS ===

  cat("Matriz Definitiva: ", nrow(matrix_run), "localidades e", ncol(matrix_run), "espécies.\n")

## Matriz Definitiva:  298 localidades e 1045 espécies.

Anotações: Esta etapa de controle de qualidade garante uma matriz definitiva padronizada, totalmente livre de artefatos amostrais que poderiam interferir nos agrupamentos ecológicos das próximas seções.

5. Análises Multivariadas e Delimitação Florística

5.1. Turnover vs. Nestedness

Partição da diversidade beta utilizando o índice de Jaccard, separando a dissimilaridade florística em componentes de substituição (Turnover) e aninhamento (Nestedness).

part_beta <- beta.multi(matrix_run, index.family = "jaccard")

cat("Dissimilaridade Jaccard Total (Regiao):", round(part_beta$beta.JAC, 3), "\n")

## Dissimilaridade Jaccard Total (Regiao): 0.996

cat("Fracao por Turnover (Substituicao):", round(part_beta$beta.JTU, 3), "\n")

## Fracao por Turnover (Substituicao): 0.994

cat("Fracao por Aninhamento (Perda de especies):", round(part_beta$beta.JNE, 3), "\n")

## Fracao por Aninhamento (Perda de especies): 0.002

Anotações: Os resultados indicam que quase a totalidade da diversidade beta regional é explicada pela substituição de espécies (Turnover), evidenciando fortes descontinuidades florísticas, refutando a hipótese de que a diversidade local é apenas um subconjunto empobrecido de uma flora dominante (Nestedness).

5.2. Escalonamento Multidimensional Não-Métrico (NMDS) com Jaccard

Ordenação não supervisionada das localidades (NMDS) e agrupamento hierárquico (Ward.D2) com base exclusivamente na dissimilaridade florística de Jaccard, ignorando classificações a priori.

dissimilarity_matrix <- vegdist(matrix_run, method = "jaccard")

tree_biome <- hclust(dissimilarity_matrix, method = "ward.D2")
k_excell <- 2
my_data_group <- cutree(tree_biome, k = k_excell)
my_data_group <- factor(my_data_group, levels = c(1, 2), labels = c("Floristic Domain A", "Floristic Domain B"))
dendro_obj <- as.dendrogram(tree_biome)

grupos_na_ordem_do_plot <- my_data_group[tree_biome$order]
left_branc <- as.character(grupos_na_ordem_do_plot[1])
cores_ramos <- if (left_branc == "Floristic Domain A") c("#d95f02", "#1b9e77") else c("#1b9e77", "#d95f02")

dendro_colorido <- color_branches(dendro_obj, k = k_excell, col = cores_ramos)
dendro_colorido <- dendrapply(dendro_colorido, function(node) {
  rotulos_no <- labels(node)
  if (length(rotulos_no) > 0 && all(rotulos_no %in% local_target)) {
    ep <- attr(node, "edgePar")
    if (is.null(ep)) attr(node, "edgePar") <- list(lwd = 3) else { ep$lwd <- 3; attr(node, "edgePar") <- ep }
  }
  if (is.leaf(node)) {
    if (length(rotulos_no) > 0 && rotulos_no %in% local_target) {
      np <- attr(node, "nodePar")
      if (is.null(np)) attr(node, "nodePar") <- list(pch = "|", col = "black", cex = 1.5) else { np$pch <- "|"; np$col <- "black"; np$cex <- 1.5; attr(node, "nodePar") <- np }
    }
    attr(node, "label") <- "" 
  }
  return(node)
})

par(mar = c(2, 4, 3, 1))
plot(dendro_colorido, main = "Dendrograma Floristico: Separacao dos Dominios (k=2)", ylab = "Distancia de Ward", xlab = "Localidades")
legend("topright", legend = c("Floristic Domain A", "Floristic Domain B", "Target Areas"), fill = c("#d95f02", "#1b9e77", "white"), border = c("black", "black", "white"), pch = c(NA, NA, "|"), col = c(NA, NA, "black"), pt.cex = 1.2, bty = "n", inset = c(0.02, 0.02))

nmds_result <- metaMDS(dissimilarity_matrix, k = 2, trymax = 100, trace = FALSE)

df_nmds_blind <- data.frame(Locality = rownames(matrix_run), NMDS1 = nmds_result$points[, 1], NMDS2 = nmds_result$points[, 2], Floristic_Group = my_data_group)

val_calc <- df_nmds_blind %>% group_by(Floristic_Group) %>% summarise(Outliers_NMDS1 = list(boxplot.stats(NMDS1)$out), .groups = "drop")
nmds_outliers_id <- df_nmds_blind %>% filter(NMDS1 %in% unlist(val_calc$Outliers_NMDS1))

df_nmds_blind <- df_nmds_blind %>% mutate(Is_Outlier = ifelse(Locality %in% nmds_outliers_id$Locality, "Sim", "Não"), Point_Color = ifelse(Locality %in% local_target, "Target Areas", as.character(Floristic_Group)))

# Centroides NMDS
centroides_nmds <- df_nmds_blind %>% group_by(Floristic_Group) %>% summarise(Centro_X = mean(NMDS1), Centro_Y = mean(NMDS2))
df_spider <- df_nmds_blind %>% left_join(centroides_nmds, by = "Floristic_Group")

# Scatter NMDS
ggplot(df_spider, aes(x = NMDS1, y = NMDS2)) +
  geom_segment(aes(xend = Centro_X, yend = Centro_Y, color = Floristic_Group), alpha = 0.2, linewidth = 0.3) +
  geom_point(aes(color = Point_Color, alpha = Is_Outlier), size = 3) +
  scale_alpha_manual(values = c("Sim" = 1, "Não" = 0.4), guide = "none") +
  geom_point(data = centroides_nmds, aes(x = Centro_X, y = Centro_Y, fill = Floristic_Group), size = 4.5, shape = 23, color = "black", stroke = 0.8) +
  geom_text_repel(data = subset(df_spider, Is_Outlier == "Sim"), aes(label = Locality), size = 3.5, fontface = "bold", box.padding = 0.6) +
  theme_bw() + labs(title = "NMDS: Dispersão Florística", x = "NMDS 1", y = "NMDS 2") +
  scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) +
  scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")

# Boxplot NMDS
ggplot(df_nmds_blind, aes(x = Floristic_Group, y = NMDS1, fill = Floristic_Group)) + geom_boxplot(alpha = 0.8) + theme_bw() + scale_fill_manual(values = c("#d95f02", "#1b9e77")) + theme(legend.position = "none")

Anotações: O dendrograma sugere a quebra dos dados em dois grandes domínios florísticos centrais (\(k=2\)). O NMDS confirma visualmente essa forte separação no espaço bidimensional, atrelado a um nível de estresse estatisticamente interpretável (< 0.20).

5.3. Análise de Coordenadas Principais (PCoA)

Análise de Coordenadas Principais (PCoA) aplicada como método de ordenação métrico complementar ao NMDS, de modo a preservar as distâncias lineares originais da matriz.

pcoa_result <- cmdscale(dissimilarity_matrix, k = 2, eig = TRUE)
percent_axis1 <- round(pcoa_result$eig[1] / sum(pcoa_result$eig) * 100, 1)
percent_axis2 <- round(pcoa_result$eig[2] / sum(pcoa_result$eig) * 100, 1)

df_pcoa_blind <- data.frame(Locality = rownames(matrix_run), PCoA1 = pcoa_result$points[, 1], PCoA2 = pcoa_result$points[, 2], Floristic_Group = my_data_group)

centroides_pcoa <- df_pcoa_blind %>% group_by(Floristic_Group) %>% summarise(Centro_X = mean(PCoA1), Centro_Y = mean(PCoA2))
val_calc_pcoa <- df_pcoa_blind %>% group_by(Floristic_Group) %>% summarise(Outliers_PCoA1 = list(boxplot.stats(PCoA1)$out), .groups = "drop")
pcoa_outliers_id <- df_pcoa_blind %>% filter(PCoA1 %in% unlist(val_calc_pcoa$Outliers_PCoA1))

df_pcoa_spider <- df_pcoa_blind %>% left_join(centroides_pcoa, by = "Floristic_Group") %>% mutate(Is_Outlier = ifelse(Locality %in% pcoa_outliers_id$Locality, "Sim", "Não"), Point_Color = ifelse(Locality %in% local_target, "Target Areas", as.character(Floristic_Group)))

ggplot(df_pcoa_spider, aes(x = PCoA1, y = PCoA2)) + geom_segment(aes(xend = Centro_X, yend = Centro_Y, color = Floristic_Group), alpha = 0.2, linewidth = 0.3) + geom_point(aes(color = Point_Color, alpha = Is_Outlier), size = 3) + scale_alpha_manual(values = c("Sim" = 1, "Não" = 0.4), guide = "none") + geom_point(data = centroides_pcoa, aes(x = Centro_X, y = Centro_Y, fill = Floristic_Group), size = 4, shape = 23, color = "black", stroke = 0.5) + geom_text_repel(data = subset(df_pcoa_spider, Is_Outlier == "Sim"), aes(label = Locality), size = 3.5, fontface = "bold", box.padding = 0.6, point.padding = 0.3, max.overlaps = Inf) + theme_bw() + labs(title = "PCoA: Dispersão Florística a partir do Núcleo", x = paste0("PCoA 1 (", percent_axis1, "%)"), y = paste0("PCoA 2 (", percent_axis2, "%)")) + scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")

ggplot(df_pcoa_blind, aes(x = Floristic_Group, y = PCoA1, fill = Floristic_Group)) + geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 2, outlier.color = "red") + theme_bw() + labs(title = "Separacao Floristica ao Longo do Gradiente", x = "Dominio", y = "PCoA1") + scale_fill_manual(values = c("#d95f02", "#1b9e77")) + theme(legend.position = "none")

Anotações: A PCoA corrobora perfeitamente a dicotomia observada no NMDS. A quantificação da variância retida nos Eixos 1 e 2 confirma a existência de um gradiente ambiental e florístico dominante moldando esta província intertropical.

5.4. Análise de Variância Permutacional Multivariada (PERMANOVA)

Aplicação do teste não paramétrico PERMANOVA para validar estatisticamente a diferença florística entre os agrupamentos derivados do dendrograma.

df_permanova_blind <- data.frame(Floristic_Group = my_data_group)
permanova_test <- adonis2(dissimilarity_matrix ~ Floristic_Group, data = df_permanova_blind, permutations = 999)
print(permanova_test)

## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dissimilarity_matrix ~ Floristic_Group, data = df_permanova_blind, permutations = 999)
##           Df SumOfSqs      R2      F Pr(>F)    
## Model      1   13.777 0.10437 34.492  0.001 ***
## Residual 296  118.229 0.89563                  
## Total    297  132.006 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dispers_grupos <- betadisper(dissimilarity_matrix, df_permanova_blind$Floristic_Group)
df_perm_scatter <- data.frame(Locality = rownames(matrix_run), PCoA1 = dispers_grupos$vectors[, 1], PCoA2 = dispers_grupos$vectors[, 2], Floristic_Group = df_permanova_blind$Floristic_Group, Centroid_Distance = dispers_grupos$distances)

df_centroids <- data.frame(Floristic_Group = rownames(dispers_grupos$centroids), Centroid1 = dispers_grupos$centroids[, 1], Centroid2 = dispers_grupos$centroids[, 2])
val_calc_perm <- df_perm_scatter %>% group_by(Floristic_Group) %>% summarise(Outliers_Dist = list(boxplot.stats(Centroid_Distance)$out), .groups = "drop")
perm_outliers_id <- df_perm_scatter %>% filter(Centroid_Distance %in% unlist(val_calc_perm$Outliers_Dist))

df_perm_scatter <- df_perm_scatter %>% left_join(df_centroids, by = "Floristic_Group") %>% mutate(Is_Outlier = ifelse(Locality %in% perm_outliers_id$Locality, "Sim", "Não"), Point_Color = ifelse(Locality %in% local_target, "Target Areas", as.character(Floristic_Group)))

ggplot(df_perm_scatter) + geom_segment(aes(x = Centroid1, y = Centroid2, xend = PCoA1, yend = PCoA2, color = Floristic_Group), alpha = 0.2, linewidth = 0.3) + geom_point(aes(x = PCoA1, y = PCoA2, color = Point_Color, alpha = Is_Outlier), size = 3) + scale_alpha_manual(values = c("Sim" = 1, "Não" = 0.4), guide = "none") + geom_point(aes(x = Centroid1, y = Centroid2, fill = Floristic_Group), size = 4, shape = 23, color = "black", stroke = 0.5) + geom_text_repel(data = subset(df_perm_scatter, Is_Outlier == "Sim"), aes(x = PCoA1, y = PCoA2, label = Locality), size = 3.5, fontface = "bold", box.padding = 0.6, point.padding = 0.3, max.overlaps = Inf) + theme_bw() + labs(title = "Dispersão da PERMANOVA", x = "Eixo 1", y = "Eixo 2") + scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")

df_permanova_blind$Centroid_Distance <- dispers_grupos$distances
ggplot(df_permanova_blind, aes(x = Floristic_Group, y = Centroid_Distance, fill = Floristic_Group)) + geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 2, outlier.color = "red") + theme_bw() + labs(title = "Variância Interna dos Grupos", x = "Domínio", y = "Distância Florística") + scale_fill_manual(values = c("#d95f02", "#1b9e77")) + theme(legend.position = "none")

Anotações: O teste PERMANOVA confirma que a divisão nos Domínios A e B possui significância estatística robusta (\(p \leq 0.05\)). O gráfico de Teia (Spider Plot) ilustra a distância e o confinamento de cada localidade em relação ao centro do seu próprio grupo ecológico.

5.5. Teste de Mantel

Teste de Mantel para avaliar o efeito prático do Isolamento por Distância (IBD), correlacionando a matriz de dissimilaridade florística com a matriz de distância geográfica real (em km).

df_coords_align <- df_locality %>% left_join(input_matrix %>% select(locality, longitud, latitud) %>% distinct(), by = "locality") %>% mutate(longitud = as.numeric(gsub(",", ".", longitud)), latitud = as.numeric(gsub(",", ".", latitud)))

surv_locality <- rownames(matrix_run)
df_coords_clean <- df_coords_align[match(surv_locality, df_coords_align$locality), ]

coords_numeric <- df_coords_clean %>% select(longitud, latitud)
coords_matrix <- as.matrix(coords_numeric)
geo_dist_matrix_km <- distm(coords_matrix, fun = distHaversine) / 1000
geo_distance_clean <- as.dist(geo_dist_matrix_km)

# Usa a dissimilarity_matrix original do 5.2
mantel_test <- mantel(dissimilarity_matrix, geo_distance_clean, method = "spearman", permutations = 999)
print(mantel_test)

## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dissimilarity_matrix, ydis = geo_distance_clean,      method = "spearman", permutations = 999) 
## 
## Mantel statistic r: 0.459 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0225 0.0304 0.0368 0.0425 
## Permutation: free
## Number of permutations: 999

df_mantel_visual <- data.frame(Dist_Geo = as.vector(geo_distance_clean), Dist_Flor = as.vector(dissimilarity_matrix))

ggplot(df_mantel_visual, aes(x = Dist_Geo, y = Dist_Flor)) + geom_point(alpha = 0.08, color = "#2c3e50", size = 0.4) + geom_smooth(method = "lm", formula = y ~ x, color = "red", linetype = "dashed", linewidth = 1, se = FALSE) + theme_bw() + labs(title = "Isolamento por Distancia", x = "Distancia Geografica Real (km)", y = "Dissimilaridade Floristica")

Anotações: A correlação linear positiva e significativa suporta o princípio do Isolamento por Distância: as diferenças na composição de espécies entre duas áreas estudadas se acentuam fortemente à medida que a distância geográfica geodésica entre elas aumenta.

5.6. Conectividade Florística Espacial (Minimum Spanning Tree)

Projeção de uma rede de ligação mínima (MST) baseada nas distâncias florísticas de Jaccard, colapsadas em nível de ecorregiões para revelar as principais rotas de similaridade biogeográfica.

# ==============================================================================
# 1. COLAPSAR A MATRIZ FLORÍSTICA PARA O NÍVEL DE ECORREGIÃO
# ==============================================================================
df_ecorregioes <- input_matrix %>% select(locality, Ecorregiao_Olson = TNC_Ecoregion_Olson) %>% distinct()

df_matriz <- as.data.frame(matrix_run)
df_matriz$locality <- rownames(df_matriz)

df_matriz_eco <- df_matriz %>% left_join(df_ecorregioes, by = "locality") %>% filter(!is.na(Ecorregiao_Olson))

matriz_eco_agrupada <- df_matriz_eco %>% select(-locality) %>% group_by(Ecorregiao_Olson) %>% summarise(across(everything(), sum)) %>% as.data.frame()

rownames(matriz_eco_agrupada) <- matriz_eco_agrupada$Ecorregiao_Olson
matriz_eco_agrupada <- matriz_eco_agrupada[, -1]
matriz_eco_binaria <- ifelse(matriz_eco_agrupada > 0, 1, 0)

# CÁLCULO DA RIQUEZA: Soma das ocorrências na matriz binária
riqueza_por_eco <- rowSums(matriz_eco_binaria)

dist_eco <- vegdist(matriz_eco_binaria, method = "jaccard")
arvore_espacial <- spantree(dist_eco)

# ==============================================================================
# 2. PREPARAÇÃO ESPACIAL, CENTROIDES E UNIÃO COM NMDS (k=2)
# ==============================================================================
cat("\nUnificando polígonos das ecorregiões a partir do shapefile base...\n")

## 
## Unificando polígonos das ecorregiões a partir do shapefile base...

ecorregioes_presentes <- rownames(matriz_eco_binaria)

sf::sf_use_s2(FALSE)

## Spherical geometry (s2) switched off

# Reutiliza o shape_olson carregado no chunk map_originais (3.1)
shape_eco_unificado <- shape_olson %>%
  filter(ECO_NAME %in% ecorregioes_presentes) %>%
  group_by(ECO_NAME) %>%
  summarise(geometry = st_union(geometry))

## although coordinates are longitude/latitude, st_union assumes that they are
## planar

## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar
## although coordinates are longitude/latitude, st_union assumes that they are
## planar

centroides <- st_centroid(shape_eco_unificado)
coords_centroides <- st_coordinates(centroides)

df_centroides <- data.frame(Ecorregiao = shape_eco_unificado$ECO_NAME, X = coords_centroides[, 1], Y = coords_centroides[, 2])
df_centroides <- df_centroides[match(rownames(matriz_eco_binaria), df_centroides$Ecorregiao), ]

# Adiciona a Riqueza
df_centroides$Riqueza <- riqueza_por_eco[df_centroides$Ecorregiao]

# --- A MÁGICA DOS DATAFRAMES PARA PUXAR A COR DO NMDS ---
# 1. Pega o grupo (A ou B) que foi definido no NMDS para cada localidade
df_grupos_loc <- df_nmds_blind %>% select(locality = Locality, Floristic_Group)

# 2. Cruza isso com a informação de qual ecorregião cada localidade pertence
df_eco_grupos <- df_grupos_loc %>% 
  left_join(df_ecorregioes, by = "locality") %>% 
  filter(!is.na(Ecorregiao_Olson))

# 3. Calcula qual é o domínio predominante dentro da ecorregião
df_grupos_nmds <- df_eco_grupos %>%
  group_by(Ecorregiao_Olson) %>%
  summarise(Grupo_k2 = names(which.max(table(Floristic_Group)))) %>%
  rename(Ecorregiao = Ecorregiao_Olson)

# 4. Junta essa informação aos centroides do mapa
df_centroides <- df_centroides %>% left_join(df_grupos_nmds, by = "Ecorregiao")
# --------------------------------------------------------

sf::sf_use_s2(TRUE) 

## Spherical geometry (s2) switched on

# ==============================================================================
# 3. CRIAR OS "RAMOS ESPACIAIS" DO MST E PREPARAR CORES BrBG
# ==============================================================================
df_ramos <- data.frame(
  X_Origem  = df_centroides$X[2:nrow(df_centroides)], 
  Y_Origem  = df_centroides$Y[2:nrow(df_centroides)], 
  X_Destino = df_centroides$X[arvore_espacial$kid], 
  Y_Destino = df_centroides$Y[arvore_espacial$kid],
  Similaridade_Jaccard = 1 - arvore_espacial$dist
)

df_ramos$X_Mid <- (df_ramos$X_Origem + df_ramos$X_Destino) / 2
df_ramos$Y_Mid <- (df_ramos$Y_Origem + df_ramos$Y_Destino) / 2
df_ramos$Jaccard_Rotulo <- sprintf("%.2f", df_ramos$Similaridade_Jaccard)

# Configuração da Paleta BrBG (preservando cores no mapa, alterando legenda)
cores_atuais <- c("#35978f", "#bf812d", "#8c510a", "#dfc27d", "#f6e8c3", "#f5f5f5", "#80cdc1", "#01665e", "#c7eae5")
ecos_alfabetica <- sort(unique(shape_eco_unificado$ECO_NAME))
names(cores_atuais) <- ecos_alfabetica
ordem_legenda_brbg <- names(cores_atuais)[c(3, 2, 4, 5, 6, 9, 7, 1, 8)]

# ==============================================================================
# 4. PLOTAGEM DO MAPA DE REDE (COM ggnewscale)
# ==============================================================================
mapa_ramos_colapsados <- ggplot() +
  
  # A. Polígonos das Ecorregiões (color = NA para tirar a linha branca)
  geom_sf(data = shape_eco_unificado, aes(fill = ECO_NAME), alpha = 0.7, color = NA) +
  scale_fill_manual(values = cores_atuais, breaks = ordem_legenda_brbg) +
  
  # B. Ramos da MST (Espessura dinâmica baseada na similaridade)
  geom_segment(
    data = df_ramos, 
    aes(x = X_Origem, y = Y_Origem, xend = X_Destino, yend = Y_Destino, linewidth = Similaridade_Jaccard), 
    color = "gray15", alpha = 0.85
  ) +
  
  # C. Rótulos numéricos da Similaridade nos Ramos
  geom_label(
    data = df_ramos, 
    aes(x = X_Mid, y = Y_Mid, label = Jaccard_Rotulo),
    size = 2.8, fill = "white", color = "black", fontface = "bold", label.padding = unit(0.15, "lines"), label.size = 0.2
  ) +
  
  # --- HABILITANDO A SEGUNDA ESCALA DE CORES PARA OS PONTOS ---
  new_scale_fill() +
  
  # D. Nós/Centroides (Tamanho pela Riqueza, Cor pelo NMDS k=2)
  geom_point(
    data = df_centroides, 
    aes(x = X, y = Y, size = Riqueza, fill = Grupo_k2), 
    shape = 21, 
    color = "black", 
    stroke = 1.3
  ) +
  
  # E. Escalas dos Centroides e Ramos
  scale_fill_manual(
    values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), 
    name = "Complexo Fitogeográfico (NMDS)",
    na.value = "red",
    guide = guide_legend(override.aes = list(size = 6)) 
  ) +
  scale_size_continuous(range = c(3.5, 9.5), name = "Riqueza de Espécies") +
  scale_linewidth_continuous(range = c(0.7, 2.2), guide = "none") + # Retirada da legenda do Jaccard
  
  # F. Tema e Títulos
  theme_minimal() + 
  labs(
    title = "Rede de Conexões Florísticas Inter-Ecorregiões", 
    subtitle = "Ramos: Similaridade de Jaccard | Nós: Riqueza e Domínio (NMDS)", 
    fill = "Ecorregião (Olson):", 
    x = "Longitude", 
    y = "Latitude"
  ) + 
  theme(
    legend.position = "right", 
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "#f0f0f0", linewidth = 0.3)
  )

print(mapa_ramos_colapsados)

5.7. Detrended Correspondence Analysis

Execução da DCA (Detrended Correspondence Analysis) para aferir o comprimento do gradiente florístico e determinar se o modelo de resposta das plantas é predominantemente linear ou unimodal.

dca_model <- decorana(matrix_run)
  print(dca_model)

## 
## Call:
## decorana(veg = matrix_run) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## Total inertia (scaled Chi-square): 19.1239 
## 
##                        DCA1   DCA2   DCA3   DCA4
## Eigenvalues          0.8076 0.4328 0.3257 0.2722
## Additive Eigenvalues 0.8076 0.4129 0.2933 0.2647
## Decorana values      0.8380 0.4078 0.3697 0.2478
## Axis lengths         7.5949 4.1883 4.7307 3.2118

df_dca <- as.data.frame(scores(dca_model, display = "sites"))
df_dca$Locality <- rownames(matrix_run)
df_dca$Floristic_Group <- my_data_group

df_dca <- df_dca %>% mutate(Point_Color = ifelse(Locality %in% local_target, "Target Areas", as.character(Floristic_Group)))

ggplot(df_dca, aes(x = DCA1, y = DCA2)) + stat_ellipse(aes(color = Floristic_Group, fill = Floristic_Group), geom = "polygon", alpha = 0.15, linetype = "solid", linewidth = 0.8) + geom_point(aes(color = Point_Color), size = 3, alpha = 0.7) + theme_bw() + labs(title = "Detrended Correspondence Analysis (DCA)", x = "Eixo DCA 1", y = "Eixo DCA 2") + scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom", plot.title = element_text(face = "bold"))

Anotações: O comprimento calculado do Eixo 1 (DCA1) dita o método de ordenação ambiental subsequente (Fluxograma). Valores acentuados (ex: > 4.0 SD) justificam análises constrangidas unimodais, como a CCA, para avaliar as condicionantes edafoclimáticas. Verificar o fluxograma abaixo para a tomada de decisão em qual análise utilizar para avaliação da influência de um conjunto de variáveis independentes sobre múltiplos resultados. [ Matriz de Espécies ] | v [ Executar DCA (Função: decorana) ] | v [ Avaliar Comprimento do Eixo 1 (DCA1) ] | +——————–+——————–+ | | | v v v [ < 3.0 SD ] [ 3.0 a 4.0 SD ] [ > 4.0 SD ] | | | v v v [ Resposta Linear ] [ Gradiente Intermediário ] [ Resposta Unimodal ] (Baixo Turnover) (Zona de Transição) (Alto Turnover) | | | v v v [ tb-RDA ] [ tb-RDA ] [ CCA ]

SE (DCA1 < 3.0) ENTÃO: Hipótese: As plantas respondem de forma linear ao solo/clima. Ação prévia: Transformar matrix_run com o método de Hellinger. Análise Final: RDA (Análise de Redundância).

SE (DCA1 >= 3.0 E DCA1 <= 4.0) ENTÃO: Hipótese: Gradiente moderado. Ambas as respostas existem. Decisão: Escolher tb-RDA (preferência atual da literatura por robustez estatística).

SE (DCA1 > 4.0) ENTÃO: Hipótese: As plantas têm uma resposta em forma de sino (unimodal). Ocorrem substituições completas de espécies nos extremos do ambiente. Ação prévia: Nenhuma transformação necessária na matriz. Análise Final: CCA (Análise Canônica de Correspondências).

6. Análise dos Grupos (Domínios)

6.1. Caracterização dos Grupos

Criação de um gráfico aluvial (Sankey) e de tabelas de contingência cruzando a classificação de biomas originais da literatura com os domínios fitogeográficos estatísticos obtidos (Data-Driven).

df_lista_grupos <- data.frame(Localidade = names(my_data_group), Dominio_Floristico = as.character(my_data_group))
df_lado_a_lado <- df_lista_grupos %>% group_by(Dominio_Floristico) %>% mutate(Linha = row_number()) %>% pivot_wider(names_from = Dominio_Floristico, values_from = Localidade) %>% select(-Linha) 
df_lado_a_lado[is.na(df_lado_a_lado)] <- ""
knitr::kable(df_lado_a_lado, caption = "Divisao das localidades entre dominios.")

Divisao das localidades entre dominios.

Floristic Domain A	Floristic Domain B
CERDAO_CER_01	CHA_UM_CHA_01
CER.S.ST_CER_01	CHA_UM_CHA_02
CERDAO_CER_02	CHA_UM_CHA_03
CERDAO_CER_03	CHA_SEC_01
CERDAO_CER_04	CHA_SEC_02
CERDAO_CER_05	CHA_SEC_03
CERDAO_CER_06	CHA_SEC_04
CERDAO_CER_07	CHA_SEC_05
CERDAO_CER_08	CHA_SEC_06
CERDAO_CER_09	CHA_SEC_07
CERDAO_CER_10	CHA_ARID_01
CERDAO_CER_11	CHA_SERRA_01
CERDAO_CER_23	CHA_BOL_01
CERDAO_CER_24	CHA_PAMPA_01
CER.S.ST_PANT_01	CHA_PANT_01
CERDAO_CER_12	CHA_PANT_02
CERDAO_CER_13	CHA_PANT_03
CILIAR_CER_01	CHA_PANT_04
CILIAR_CER_02	CHA_UM_CHA_04
CILIAR_CER_03	CHA_UM_CHA_05
CILIAR_CER_04	CHA_UM_CHA_06
CILIAR_CER_05	CHA_PANT_06
CAMP.MUR_CER_01	CHA_UM_CHA_07
CAMP.MUR_CER_02	CHA_BOL_02
CERDAO_CER_14	CHA_SEM_SAL_01
CER.RUPES_CER_01	CHA_SEM_ARID_01
CILIAR_CER_06	CHA_SERRA_02
CER.S.ST_CER_02	CHA_PANT_07
CILIAR_CER_07	CHA_PANT_08
CILIAR_CER_08	CHA_PANT_09
CILIAR_CER_09	CHA_SEM_ARID_02
CERDAO_CER_15	CHA_SEM_ARID_03
CERDAO_CER_25	CHA_SEM_ARID_04
CERDAO_CER_27	CHA_SEM_ARID_05
FES_PANT_17	CHA_SEM_ARID_06
CILIAR_CER_10	CHA_SEM_ARID_07
CILIAR_CER_11	CHA_SEM_ARID_08
CERDAO_CER_28	CHA_PANT_10
CERDAO_CER_29	CHA_AREN_01
CERDAO_CER_30	CHA_AREN_02
CERDAO_CER_31	CHA_ENC_AREN_02
CILIAR_CER_12	CHA_INTERDUN_01
CERDAO_CER_16	CHA_INTERDUN_02
CERDAO_CER_17	CHA_INTERDUN_03
CERDAO_CER_18	CHA_LEQ_ALU_01
TRANS_CER_01	CHA_LEQ_ALU_02
FED_PANT_01	CHA_LEQ_ALU_03
FED_CER_01	CHA_LEQ_ALU_04
FED_PANT_02	CHA_LEQ_ALU_05
FED_PANT_03	CHA_LEQ_ALU_06
CILIAR_CER_13	CHA_LEQ_ALU_07
CILIAR_MA_01	CHA_LEQ_ALU_08
CILIAR_CER_14	CHA_LEQ_ALU_09
FES_CER_01	CHA_LEQ_ALU_10
FES_CER_02	CHA_LEQ_ALU_11
FES_CER_03	CHA_LEQ_ALU_12
FES_MA_01	CHA_LEQ_ALU_13
FES_MA_02	CHA_LEQ_ALU_14
FES_MA_03	CHA_LEQ_ALU_15
FES_MA_04	CHA_LEQ_ALU_16
FES_PANT_01	CHA_LEQ_ALU_17
FES_CHI_01	CHA_LEQ_ALU_18
CILIAR_MA_02	CHA_LEQ_ALU_19
FES_CER_04	CHA_CHIQ_01
FES_CER_05	CHA_CHIQ_02
FES_CER_06	CHA_CHIQ_03
FES_CER_07	CHA_CHIQ_04
FES_PANT_02	CHA_CHIQ_05
FES_PANT_03	CHA_CHIQ_06
FES_PANT_04	CHA_CHIQ_07
FES_PANT_05	CHA_CHIQ_08
FES_PANT_06	CHA_CHIQ_09
FES_CER_08	CHA_CHIQ_10
FES_CER_09	CHA_CHIQ_11
FES_CER_10	CHA_CHIQ_12
FES_CER_11	CHA_CHIQ_13
FES_MA_05	CHA_EASTERN_01
FES_MA_06	CHA_EASTERN_02
FES_MA_07	CHA_EASTERN_03
FES_MA_08	CHA_EASTERN_04
FED_PANT_04	CHA_LEQ_ALU_20
PARAT_PANT_03	CHA_LEQ_ALU_21
CERDAO_PANT_01	CHA_LEQ_ALU_22
CERDAO_PANT_02	CHA_LEQ_ALU_23
CERDAO_PANT_03	CHA_LEQ_ALU_24
CILIAR_PANT_01	CHA_LEQ_ALU_25
TRANS_CER_02	CHA_LEQ_ALU_26
CAPAO_PANT_13	CHA_LEQ_ALU_27
CERDAO_PANT_23	CHA_LEQ_ALU_28
CAPAO_PANT_07	CHA_SEM_AR_LEQ_ALU_01
CAPAO_PANT_06	CHA_EASTERN_05
CILIAR_PANT_02	CHA_EASTERN_06
CILIAR_PANT_03	CHA_EASTERN_07
CILIAR_PANT_04	CHA_HYGRO_01
CILIAR_PANT_05	CHA_HYGRO_02
CILIAR_PANT_06	CHA_SEC_08
CAMP.CERR_PANT_05	CHA_UM_CHA_08
FES_PANT_18	CHA_SEC_09
FED_PANT_05	CHA_SEC_10
CAPAO_PANT_01	CHA_ESPI_01
CAPAO_PANT_08	CHA_UM_CHA_09
CERDAO_PANT_04	CHA_UM_CHA_10
CERDAO_PANT_05	CHA_SEC_12
CERDAO_PANT_06	CHA_SEC_13
CERDAO_PANT_07	CHA_SEC_14
CERDAO_PANT_08	CHA_UM_CHA_11
CERDAO_PANT_09	CHA_SEC_16
CERDAO_PANT_10	CHA_BOL_05
CERDAO_PANT_11
CAPAO_PANT_02
CERDAO_PANT_20
CAPAO_PANT_03
CAMP.CERR_PANT_01
CERDAO_PANT_12
CILIAR_PANT_07
CERDAO_PANT_19
CERDAO_PANT_21
FED_PANT_07
FED_PANT_08
FED_PANT_09
FED_PANT_10
CAMP.CERR_PANT_02
CAPAO_PANT_14
PARAT_PANT_02
CAMP.MUR_PANT_01
CAPAO_PANT_12
CAPAO_PANT_09
CAMP.MUR_PANT_02
CERDAO_PANT_13
CAPAO_PANT_04
CAPAO_PANT_05
FES_PANT_07
CAMP.CERR_PANT_03
CERDAO_PANT_14
CAPAO_PANT_11
CILIAR_PANT_08
FES_PANT_08
CILIAR_PANT_09
CILIAR_PANT_10
CILIAR_PANT_11
CERDAO_PANT_22
MAT_ACURI_PANT_01
CAPAO_PANT_10
CAMB_PANT_01
CAPAO_PANT_15
CERDAO_PANT_15
CERDAO_CER_19
CILIAR_PANT_12
CERDAO_PANT_16
FES_PANT_09
FES_PANT_10
CILIAR_PANT_13
CILIAR_PANT_14
FESM_PANT_01
CILIAR_PANT_15
FES_PANT_11
CAMP.CERR_PANT_04
FES_PANT_12
CAMP.CERR_CER_01
CAMP.MUR_CER_03
CER.S.ST_CER_03
CERDAO_CER_20
BABA_CER_01
CAPAO_CER_01
CILIAR_CER_18
CILIAR_MA_03
CERDAO_CER_21
FED_PANT_06
CERDAO_CER_22
CERDAO_PANT_17
CERDAO_PANT_18
FES_PANT_13
FES_PANT_14
CER.S.ST_PANT_02
FES_PANT_15
FES_PANT_16
CAMB_PANT_02
CAMP.MUR_PANT_03
CHA_CHIQ_14
CHA_CHIQ_15
CHA_CHIQ_16
CHA_CHIQ_17
CHA_CHIQ_18
CHA_CHIQ_19
CERDAO_CER_32
CERDAO_CER_34
CERDAO_CER_36
CER.S.ST_CER_04
CHA_BOL_03
CHA_BOL_04

locs_finais <- rownames(matrix_run)
df_biome_filtrado <- df_biome_align %>% filter(locality %in% locs_finais) %>% arrange(match(locality, locs_finais))

df_flow <- data.frame(Locality = locs_finais, Dominio_Original = df_biome_filtrado$domain_biome, Dominio_Cluster = as.character(my_data_group)) %>% mutate(Flow_Color = ifelse(Locality %in% local_target, "Target Areas", as.character(Dominio_Cluster)))

ggplot(df_flow, aes(axis1 = Dominio_Original, axis2 = Dominio_Cluster)) + geom_alluvium(aes(fill = Flow_Color), width = 1/8, alpha = 0.7, knot.pos = 0.4) + geom_stratum(width = 1/4, fill = "grey80", color = "black") + geom_text_repel(stat = "stratum", aes(label = str_wrap(after_stat(stratum), width = 18)), size = 3.5, direction = "y", nudge_x = 0) + scale_x_discrete(limits = c("Bioma Original", "Grupo Floristico (k=2)"), expand = c(.05, .05)) + scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + theme_minimal() + labs(title = "Fluxo de Reclassificacao Floristica", y = "Numero de Localidades") + theme(legend.position = "bottom")

matriz_A <- matrix_run[my_data_group == "Floristic Domain A", ]
matriz_B <- matrix_run[my_data_group == "Floristic Domain B", ]

lista_venn_dominios <- list("Domain A" = colnames(matriz_A)[colSums(matriz_A) > 0], "Domain B" = colnames(matriz_B)[colSums(matriz_B) > 0])
names(lista_venn_dominios) <- paste0(names(lista_venn_dominios), "\n(", sapply(lista_venn_dominios, length), " spp.)")

# AQUI ESTÁ A CORREÇÃO: Inserção do ilabels = "counts"
venn(lista_venn_dominios, zcolor = c("#d95f02", "#1b9e77"), opacity = 0.5, ilabels = "counts", ilcs = 1.5, sncs = 1.2, box = FALSE)
title("Compartilhamento de Especies entre Dominios")

beta_A <- beta.multi(matriz_A, index.family = "jaccard")
beta_B <- beta.multi(matriz_B, index.family = "jaccard")

df_beta_longo <- data.frame(Dominio = c("Floristic Domain A", "Floristic Domain B"), Turnover = c(beta_A$beta.JTU, beta_B$beta.JTU), Nestedness = c(beta_A$beta.JNE, beta_B$beta.JNE)) %>% pivot_longer(cols = c(Turnover, Nestedness), names_to = "Componente", values_to = "Valor_Beta") %>% group_by(Dominio) %>% mutate(Total = sum(Valor_Beta), Proporcao = Valor_Beta / Total * 100, Label = paste0(round(Proporcao, 1), "%")) %>% ungroup()

ggplot(df_beta_longo, aes(x = Dominio, y = Valor_Beta, fill = Componente)) + geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7, color = "black") + geom_text(aes(label = Label), position = position_dodge(width = 0.8), vjust = -0.5, fontface = "bold", size = 4.5) + theme_bw() + labs(title = "Componentes da Diversidade Beta por Dominio", y = "Valor da Dissimilaridade") + scale_y_continuous(expand = expansion(mult = c(0, 0.15))) + scale_fill_manual(values = c("Turnover" = "#2c3e50", "Nestedness" = "#bdc3c7")) + theme(legend.position = "bottom")

Anotações: O fluxo aluvial revela a correspondência ou o espalhamento de transição dos biomas ao serem submetidos a algoritmos de clusterização imparciais. Adicionalmente, o particionamento Beta confirma a invariabilidade da predominância do Turnover dentro de cada um dos novos domínios.

6.2. Análise de espécies indicadoras

Identificação estatística das espécies indicadoras de cada Domínio através do método IndVal e subsequente validação espacial utilizando pontos de ocorrência georreferenciada resgatados do GBIF.

resultados_indval <- multipatt(matrix_run, my_data_group, func = "IndVal.g", duleg = TRUE, control = how(nperm = 999))
res_tabela <- resultados_indval$sign
res_tabela$especie <- rownames(res_tabela)
col_grupo_A <- grep("s.Floristic.Domain.A", colnames(res_tabela), value = TRUE)
col_grupo_B <- grep("s.Floristic.Domain.B", colnames(res_tabela), value = TRUE)

res_sig <- res_tabela %>% filter(p.value <= 0.05)
res_sig$dominio <- NA_character_
res_sig$dominio[res_sig[, col_grupo_A] == 1] <- "Floristic Domain A"
res_sig$dominio[res_sig[, col_grupo_B] == 1] <- "Floristic Domain B"
res_sig <- res_sig %>% filter(!is.na(dominio))

top_indicadoras <- res_sig %>% 
  group_by(dominio) %>% 
  slice_max(stat, n = 10) %>% 
  arrange(dominio, desc(stat)) %>%
  mutate(especie_clean = stringr::word(gsub("_", " ", especie), 1, 2))

ggplot(top_indicadoras, aes(x = stat, y = reorder(especie_clean, stat), color = dominio)) + 
  geom_segment(aes(x = 0, xend = stat, y = reorder(especie_clean, stat), yend = reorder(especie_clean, stat)), color = "grey60", linetype = "dashed") + 
  geom_point(aes(shape = dominio), size = 4) +  # <-- Adicionado mapeamento de shape
  facet_wrap(~dominio, scales = "free_y", ncol = 1) + 
  theme_bw() + 
  labs(title = "Top 10 Especies Indicadoras", x = "Valor de Indicação (stat)", y = "") + 
  scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77")) + 
  scale_shape_manual(values = c("Floristic Domain A" = 16, "Floristic Domain B" = 17)) + # <-- 16 (Círculo), 17 (Triângulo)
  theme(legend.position = "none", axis.text.y = element_text(face = "italic"), strip.background = element_rect(fill = "#f0f0f0"), strip.text = element_text(face = "bold"))

df_indicadoras <- top_indicadoras %>% select(especie = especie_clean, dominio)

cat("\nBuscando ocorrencias no GBIF para as especies indicadoras...\n")

## 
## Buscando ocorrencias no GBIF para as especies indicadoras...

ocorrencias_list <- lapply(df_indicadoras$especie, function(sp) {
  dados_gbif <- occ_data(scientificName = sp, limit = 300, hasCoordinate = TRUE, hasGeospatialIssue = FALSE)$data 
  if(!is.null(dados_gbif) && is.data.frame(dados_gbif) && nrow(dados_gbif) > 0 && "decimalLongitude" %in% colnames(dados_gbif)) {
    dados_filtrados <- dados_gbif %>% select(decimalLongitude, decimalLatitude) %>% mutate(especie = sp) 
    return(dados_filtrados)
  } else { return(NULL) }
})

df_ocorrencias <- bind_rows(ocorrencias_list)

if(nrow(df_ocorrencias) > 0 && "decimalLongitude" %in% colnames(df_ocorrencias)) {
  df_ocorrencias <- df_ocorrencias %>% filter(decimalLongitude >= boundary_area[1] & decimalLongitude <= boundary_area[2], decimalLatitude >= boundary_area[3] & decimalLatitude <= boundary_area[4])
  
  if(nrow(df_ocorrencias) == 0) stop("GBIF encontrou dados, mas NENHUM PONTO caiu dentro do 'boundary_area'.")
  
  df_ocorrencias <- df_ocorrencias %>% left_join(df_indicadoras, by = "especie")
  points_sf <- st_as_sf(df_ocorrencias, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
  shape_biomas <- st_read("mapa_pontos.gpkg", quiet = TRUE) %>% st_transform(4326) 
  
  mapa_validacao <- ggplot() + 
    geom_sf(data = shape_biomas, fill = "grey90", color = "grey50", linewidth = 0.5) + 
    # 1. Adicionado shape = dominio no aes()
    geom_sf(data = points_sf, aes(color = dominio, shape = dominio), size = 1.5, alpha = 0.8) + 
    coord_sf(xlim = boundary_area[1:2], ylim = boundary_area[3:4], expand = FALSE) + 
    scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77")) + 
    # 2. Adicionado scale_shape_manual para garantir os círculos e triângulos
    scale_shape_manual(values = c("Floristic Domain A" = 16, "Floristic Domain B" = 17)) +
    theme_minimal() + 
    # 3. Adicionado shape = "Dominio:" para mesclar as legendas de cor e formato
    labs(title = "Validacao Espacial das Indicadoras", color = "Dominio:", shape = "Dominio:") + 
    theme(legend.position = "right")
  print(mapa_validacao)
} else { stop("Nenhuma especie indicadora retornou ocorrencias do GBIF.") }

Anotações: O Dot Plot destaca os táxons mais fiéis a cada agrupamento. O mapa final cruza esses resultados teóricos com a realidade, demonstrando inequivocamente que a distribuição natural das espécies indicadoras acompanha as projeções geográficas dos domínios estatísticos.

6.3. Separação e Eligibilidades do Grupo (Zoom Target)

Isolamento computacional do Domínio Fitogeográfico específico que englobou as áreas-alvo brasileiras de estudo, desvinculando-o dos atritos florísticos externos à sua ecorregião.

df_conferencia <- data.frame(Localidade = names(my_data_group), Dominio_Estatistico = as.character(my_data_group))
resultado_alvos <- df_conferencia %>% filter(Localidade %in% local_target) 

dominio_alvo <- names(which.max(table(resultado_alvos$Dominio_Estatistico)))
cat("\n[!] O domínio que engloba as suas áreas de interesse é o:", dominio_alvo, "\n")

## 
## [!] O domínio que engloba as suas áreas de interesse é o: Floristic Domain B

matrix_dominio_foco <- matrix_run[my_data_group == dominio_alvo, ]
spp_presentes <- colSums(matrix_dominio_foco) > 0
matrix_alvo_limpa <- matrix_dominio_foco[, spp_presentes]

dissimilarity_alvo <- vegdist(matrix_alvo_limpa, method = "jaccard")
tree_alvo <- hclust(dissimilarity_alvo, method = "ward.D2")
dend_alvo <- as.dendrogram(tree_alvo)

cor_dominio <- ifelse(dominio_alvo == "Floristic Domain A", "#d95f02", "#1b9e77")
dend_alvo <- color_branches(dend_alvo, k = 1, col = cor_dominio)

dend_alvo <- dendrapply(dend_alvo, function(node) {
  rotulos_no <- labels(node)
  if (length(rotulos_no) > 0 && all(rotulos_no %in% local_target)) {
    ep <- attr(node, "edgePar")
    if (is.null(ep)) attr(node, "edgePar") <- list(lwd = 3) else { ep$lwd <- 3; attr(node, "edgePar") <- ep }
  }
  if (is.leaf(node)) {
    if (length(rotulos_no) > 0 && rotulos_no %in% local_target) {
      np <- attr(node, "nodePar")
      if (is.null(np)) attr(node, "nodePar") <- list(pch = "|", col = "black", cex = 1.5) else { np$pch <- "|"; np$col <- "black"; np$cex <- 1.5; attr(node, "nodePar") <- np }
    }
    attr(node, "label") <- "" 
  }
  return(node)
})

par(mar = c(2, 4, 3, 1)) 
plot(dend_alvo, main = paste("Estrutura Interna -", dominio_alvo), ylab = "Distância de Ward", xlab = "Localidades")
title(xlab = "Localidades", line = 0.5, font.lab = 2)
legend("topright", legend = c(paste("Domínio:", dominio_alvo), "Target Areas"), fill = c(cor_dominio, "white"), border = c("black", "white"), pch = c(NA, "|"), col = c(NA, "black"), pt.cex = 1.2, bty = "n", cex = 0.9, inset = c(0.02, 0.02))

Anotações: Ao reconstruir o dendrograma utilizando apenas a matriz do Domínio Alvo, a topologia se refina e ganha precisão analítica. Isso foca o estudo estritamente na homogeneidade de escala restrita e no ecótono adjacente ao Chaco Brasileiro.

6.4. Análise do dominio B

6.4.1. Separação de Subgrupos do Domínio B

Varredura automatizada em busca do número ideal de subgrupos biológicos (\(k\)) dentro do domínio alvo, baseando-se no maior acúmulo de espécies indicadoras significativas (IndVal).

cat("Busca automatica de subgrupos para o", dominio_alvo, "...\n")

## Busca automatica de subgrupos para o Floristic Domain B ...

k_atual <- 2
quedas_consecutivas <- 0
limite_quedas <- 3
max_spp_registrado <- 0
limite_seguranca <- 20

df_resultados_sub <- data.frame(K = integer(), Spp_Significativas = integer(), Soma_IndVal = numeric())
set.seed(123)

while(quedas_consecutivas < limite_quedas && k_atual <= limite_seguranca) {
  grupos_temp <- cutree(tree_alvo, k = k_atual)
  indval_temp <- multipatt(matrix_alvo_limpa, grupos_temp, func = "IndVal.g", duleg = TRUE, control = how(nperm = 999))
  
  res <- indval_temp$sign
  spp_sig <- res[!is.na(res$p.value) & res$p.value <= 0.05, ]
  
  qtd_spp <- nrow(spp_sig)
  df_resultados_sub <- rbind(df_resultados_sub, data.frame(K = k_atual, Spp_Significativas = qtd_spp, Soma_IndVal = sum(spp_sig$stat, na.rm = TRUE)))
  
  if(qtd_spp >= max_spp_registrado) { max_spp_registrado <- qtd_spp; quedas_consecutivas <- 0 } else { quedas_consecutivas <- quedas_consecutivas + 1 }
  k_atual <- k_atual + 1
}

pico_sub_k <- df_resultados_sub$K[which.max(df_resultados_sub$Spp_Significativas)]

df_plot_sub <- df_resultados_sub %>% pivot_longer(cols = c(Spp_Significativas, Soma_IndVal), names_to = "Metrica", values_to = "Valor")
df_plot_sub$Metrica <- factor(df_plot_sub$Metrica, levels = c("Spp_Significativas", "Soma_IndVal"), labels = c("Total Espécies", "Soma IndVal"))

grafico_otimizacao_sub <- ggplot(df_plot_sub, aes(x = K, y = Valor)) + geom_line(color = "#2980b9", linewidth = 1.2) + geom_point(size = 3, color = "#2c3e50") + geom_vline(xintercept = pico_sub_k, linetype = "dashed", color = "#c0392b", linewidth = 1) + facet_wrap(~Metrica, scales = "free_y", ncol = 1) + scale_x_continuous(breaks = min(df_resultados_sub$K):max(df_resultados_sub$K)) + theme_bw() + labs(title = paste("Otimizacao de K em", dominio_alvo), subtitle = paste("Pico biológico em k =", pico_sub_k))
print(grafico_otimizacao_sub)

k_sub <- pico_sub_k 
paleta_sub <- c("#d1ece1", "#a3d8c6", "#75c5ac", "#48b191", "#147659", "#0d4f3b")[1:k_sub]

dend_base <- as.dendrogram(tree_alvo)
grupos_sub <- cutree(dend_base, k = k_sub, order_clusters_as_data = FALSE)
dend_base_colorido <- color_branches(dend_base, k = k_sub, col = paleta_sub)

dend_alvo_circ <- dendrapply(dend_base_colorido, function(node) {
  rotulos_no <- labels(node)
  if (length(rotulos_no) > 0 && all(rotulos_no %in% local_target)) {
    ep <- attr(node, "edgePar")
    if (is.null(ep)) attr(node, "edgePar") <- list(lwd = 3) else { ep$lwd <- 3; attr(node, "edgePar") <- ep }
  }
  if (is.leaf(node)) attr(node, "label") <- "" 
  return(node)
})

par(mar = c(2, 2, 2, 8))
suppressWarnings(circlize_dendrogram(dend_alvo_circ, labels = FALSE, main = paste("Subdivisões", dominio_alvo)))
legend("center", legend = paste("Agrupamento", 1:k_sub), fill = paleta_sub, border = rep("black", k_sub), bty = "n", cex = 0.75, title = "Subgrupos")

for (i in 1:k_sub) {
  locs_grupo <- names(grupos_sub[grupos_sub == i])
  locs_remover <- labels(dend_base_colorido)[!labels(dend_base_colorido) %in% locs_grupo]
  dend_individual <- prune(dend_base_colorido, locs_remover)
  dend_individual <- set(dend_individual, "labels_cex", 0.5)
  dend_individual <- set(dend_individual, "branches_lwd", 2)
  
  dend_individual <- dendrapply(dend_individual, function(node) {
    rotulos_no <- labels(node)
    if (length(rotulos_no) > 0 && all(rotulos_no %in% local_target)) {
      ep <- attr(node, "edgePar")
      if (is.null(ep)) attr(node, "edgePar") <- list(lwd = 3) else { ep$lwd <- 3; attr(node, "edgePar") <- ep }
    }
    if (is.leaf(node)) {
      rotulo <- attr(node, "label")
      if (length(rotulo) > 0 && rotulo %in% local_target) {
        np <- attr(node, "nodePar")
        if (is.null(np)) attr(node, "nodePar") <- list(lab.col = "black", lab.font = 2) else { np$lab.col <- "black"; np$lab.font <- 2; attr(node, "nodePar") <- np }
      }
    }
    return(node)
  })
  
  par(mar = c(5, 2, 4, 15))
  plot(dend_individual, horiz = TRUE, main = paste("Agrupamento", i, ":", dominio_alvo), xlab = "Distância de Ward")
}

Anotações: O modelo de parada estocástica (Early Stopping) identificou o pico biológico paramétrico da meta-comunidade antes de fragmentá-la irrealisticamente. O dendrograma circular expõe com elegância estas ramificações florísticas menores em formato de árvore filogenética.

6.4.2. Espécies Indicadoras dos Subgrupos (k ótimo)

Identificação e ranqueamento das espécies indicadoras (IndVal) que sustentam biologicamente a divisão dos subgrupos encontrados na meta-comunidade do domínio alvo.

indval_sub <- multipatt(matrix_alvo_limpa, grupos_sub, func = "IndVal.g", duleg = TRUE, control = how(nperm = 999))
res_sub <- indval_sub$sign
res_sub$especie <- rownames(res_sub)

res_sub_sig <- res_sub %>% filter(!is.na(p.value) & p.value <= 0.05)
res_sub_sig$Subgrupo <- NA_character_
for (i in 1:k_sub) {
  coluna_grupo <- paste0("s.", i)
  if (coluna_grupo %in% colnames(res_sub_sig)) res_sub_sig$Subgrupo[res_sub_sig[[coluna_grupo]] == 1] <- paste("Agrupamento", i)
}
res_sub_sig <- res_sub_sig %>% filter(!is.na(Subgrupo))

top_indicadoras_sub <- res_sub_sig %>% group_by(Subgrupo) %>% slice_max(stat, n = 10) %>% arrange(Subgrupo, desc(stat))
print(knitr::kable(top_indicadoras_sub %>% select(Subgrupo, especie, stat, p.value), digits = 3, caption = paste("Top 10 Indicadoras", dominio_alvo)))

## 
## 
## Table: Top 10 Indicadoras Floristic Domain B
## 
## |Subgrupo      |especie                                                      |  stat| p.value|
## |:-------------|:------------------------------------------------------------|-----:|-------:|
## |Agrupamento 1 |Neltuma_x_vinalillo_(Stuck.)_C.E.Hughes_&_G.P.Lewis          | 0.606|   0.001|
## |Agrupamento 1 |Neltuma_nigra_(Griseb.)_C.E.Hughes_&_G.P.Lewis               | 0.540|   0.007|
## |Agrupamento 1 |Harrisia_pomanensis_(F.A.C.Weber_ex_K.Schum.)_Britton_&_Rose | 0.535|   0.002|
## |Agrupamento 1 |Parkinsonia_praecox_(Ruiz_&_Pav.)_Hawkins                    | 0.522|   0.014|
## |Agrupamento 1 |Neltuma_ruscifolia_(Griseb.)_C.E.Hughes_&_G.P.Lewis          | 0.482|   0.024|
## |Agrupamento 1 |Annona_nutans_(R.E.Fr.)_R.E.Fr.                              | 0.474|   0.002|
## |Agrupamento 1 |Mimosa_hexandra_Micheli                                      | 0.452|   0.007|
## |Agrupamento 1 |Jodina_rhombifolia_(Hook._&_Arn.)_Hook._&_Arn._ex_Reissek    | 0.434|   0.030|
## |Agrupamento 1 |Chloroleucon_tenuiflorum_(Benth.)_Barneby_&_J.W.Grimes       | 0.429|   0.015|
## |Agrupamento 1 |Ruprechtia_exploratricis_Sandwith                            | 0.404|   0.011|
## |Agrupamento 1 |Neltuma_rubriflora_(Hassl.)_C.E.Hughes_&_G.P.Lewis           | 0.404|   0.006|
## |Agrupamento 1 |Mimosa_glutinosa_Malme                                       | 0.404|   0.014|
## |Agrupamento 2 |Ptilochaeta_nudipes_Griseb.                                  | 0.757|   0.001|
## |Agrupamento 2 |Cereus_hildmannianus_K.Schum.                                | 0.708|   0.001|
## |Agrupamento 2 |Quiabentia_verticillata_(Vaupel)_Vaupel_ex_A.Berger          | 0.664|   0.001|
## |Agrupamento 2 |Senegalia_emilioana_(Fortunato_&_Ciald.)_Seigler_&_Ebinger   | 0.645|   0.001|
## |Agrupamento 2 |Cereus_stenogonus_K.Schum.                                   | 0.643|   0.002|
## |Agrupamento 2 |Coccoloba_argentinensis_Speg.                                | 0.639|   0.001|
## |Agrupamento 2 |Ceiba_insignis_(Kunth)_P.E.Gibbs_&_Semir                     | 0.631|   0.001|
## |Agrupamento 2 |Argythamnia_breviramea_Müll.Arg.                             | 0.621|   0.001|
## |Agrupamento 2 |Chloroleucon_chacoense_(Burkart)_Barneby_&_J.W.Grimes        | 0.612|   0.002|
## |Agrupamento 2 |Mimosa_castanoclada_Barneby_&_Fortunato                      | 0.597|   0.001|
## |Agrupamento 3 |Phyllostylon_rhamnoides_(J.Poiss.)_Taub.                     | 0.628|   0.001|
## |Agrupamento 3 |Harrisia_bonplandii_(Parm._ex_Pfeiff.)_Britton_&_Rose        | 0.598|   0.001|
## |Agrupamento 3 |Bulnesia_sarmientoi_Lorentz_ex_Griseb.                       | 0.589|   0.003|
## |Agrupamento 3 |Salta_triflora_(Griseb.)_Adr.Sanchez                         | 0.581|   0.006|
## |Agrupamento 3 |Morisonia_tweedieana_(Eichler)_Christenh._&_Byng             | 0.577|   0.004|
## |Agrupamento 3 |Trithrinax_schizophylla_Drude                                | 0.517|   0.002|
## |Agrupamento 3 |Opuntia_discolor_Britton_&_Rose                              | 0.487|   0.011|
## |Agrupamento 3 |Adelia_membranifolia_(Müll.Arg.)_Chodat_&_Hassl.             | 0.427|   0.017|
## |Agrupamento 3 |Acanthocalycium_rhodotrichum_(K.Schum.)_Schlumpb.            | 0.404|   0.017|
## |Agrupamento 3 |Allophylus_pauciflorus_Radlk.                                | 0.394|   0.046|
## |Agrupamento 4 |Neltuma_alba_(Griseb.)_C.E.Hughes_&_G.P.Lewis                | 0.641|   0.001|
## |Agrupamento 4 |Celtis_iguanaea_(Jacq.)_Sarg.                                | 0.579|   0.002|
## |Agrupamento 4 |Salix_humboldtiana_Willd.                                    | 0.573|   0.001|
## |Agrupamento 4 |Aloysia_gratissima_(Gillies_&_Hook.)_Tronc.                  | 0.548|   0.006|
## |Agrupamento 4 |Ceiba_chodatii_(Hassl.)_Ravenna                              | 0.531|   0.003|
## |Agrupamento 4 |Cestrum_parqui_(Lam.)_L'Hér.                                 | 0.531|   0.002|
## |Agrupamento 4 |Sapium_haematospermum_Müll.Arg.                              | 0.514|   0.010|
## |Agrupamento 4 |Geoffroea_decorticans_(Gillies_ex_Hook._&_Arn.)_Burkart      | 0.495|   0.017|
## |Agrupamento 4 |Syagrus_romanzoffiana_(Cham.)_Glassman                       | 0.488|   0.003|
## |Agrupamento 4 |Muehlenbeckia_sagittifolia_(Ortega)_Meisn.                   | 0.488|   0.002|

if(nrow(top_indicadoras_sub) > 0) {
  # Vetor de formas geométricas para manter o padrão visual do mapa
  formas_dinamicas <- c(21, 24, 22, 23, 25, 21)[1:k_sub]

  grafico_indval_sub <- ggplot(top_indicadoras_sub, aes(x = stat, y = reorder(especie, stat))) + 
    geom_segment(aes(x = 0, xend = stat, y = reorder(especie, stat), yend = reorder(especie, stat)), color = "grey60", linetype = "dashed") + 
    # Mapeando o preenchimento (fill) e formato (shape) por subgrupo, com borda preta
    geom_point(aes(fill = Subgrupo, shape = Subgrupo), size = 4, color = "black") + 
    facet_wrap(~Subgrupo, scales = "free_y", ncol = 2) + 
    theme_bw() + 
    labs(title = "Espécies Indicadoras", x = "IndVal", y = "") + 
    # Aplicando as escalas idênticas às dos mapas
    scale_fill_manual(values = paleta_sub) + 
    scale_shape_manual(values = formas_dinamicas) + 
    theme(legend.position = "none", axis.text.y = element_text(face = "italic"), strip.background = element_rect(fill = "#f0f0f0"), strip.text = element_text(face = "bold"))
  
  print(grafico_indval_sub)
}

6.4.3. Distribuição Espacial dos Subgrupos (Domínio Alvo)

Projeção geográfica dos subgrupos delimitados no dendrograma do Domínio Alvo. O mapa utiliza a caixa delimitadora (Bounding Box) exclusiva das coordenadas amostrais, acrescida de 10% de margem, sobrepondo os subgrupos à paleta de cores original das ecorregiões de fundo.

cat("\nPreparando dados espaciais dos subgrupos e sobrepondo às ecorregiões originais...\n")

## 
## Preparando dados espaciais dos subgrupos e sobrepondo às ecorregiões originais...

# 1. Extrair os subgrupos do dendrograma e juntar com as coordenadas
df_subgrupos <- data.frame(locality = names(grupos_sub), Subgrupo = paste("Agrupamento", grupos_sub))

df_sub_coords <- df_subgrupos %>%
  left_join(input_matrix %>% select(locality, longitud, latitud) %>% distinct(), by = "locality") %>%
  mutate(long = as.numeric(gsub(",", ".", longitud)), lat = as.numeric(gsub(",", ".", latitud)))
  
# Transformar em objeto espacial
points_sub_sf <- st_as_sf(df_sub_coords, coords = c("long", "lat"), crs = 4326)

# Capturar os limites geográficos EXCLUSIVOS dos pontos dos subgrupos
bbox_pontos <- st_bbox(points_sub_sf)

# Criar uma margem de segurança (10%) para os pontos não tocarem na borda da figura
margem_x <- (bbox_pontos["xmax"] - bbox_pontos["xmin"]) * 0.1
margem_y <- (bbox_pontos["ymax"] - bbox_pontos["ymin"]) * 0.1

limites_x <- c(bbox_pontos["xmin"] - margem_x, bbox_pontos["xmax"] + margem_x)
limites_y <- c(bbox_pontos["ymin"] - margem_y, bbox_pontos["ymax"] + margem_y)

# Vetor de formas geométricas dinâmicas
formas_dinamicas <- c(21, 24, 22, 23, 25, 21)[1:k_sub]

# Resgata a configuração original de cores das Ecorregiões (Paleta BrBG)
cores_atuais <- c("#35978f", "#bf812d", "#8c510a", "#dfc27d", "#f6e8c3", "#f5f5f5", "#80cdc1", "#01665e", "#c7eae5")
ecos_alfabetica <- sort(unique(shape_eco_originais$ECO_NAME))
names(cores_atuais) <- ecos_alfabetica[1:length(cores_atuais)]
ordem_legenda_brbg <- names(cores_atuais)[c(3, 2, 4, 5, 6, 9, 7, 1, 8)]

# 3. Plotagem do Mapa com enquadramento cirúrgico e múltiplas escalas
mapa_subgrupos <- ggplot() + 
  
  # --- CAMADA 1: FUNDO DAS ECORREGIÕES ---
  geom_sf(data = shape_eco_originais, aes(fill = ECO_NAME), alpha = 0.7, color = NA) + 
  scale_fill_manual(values = cores_atuais, breaks = ordem_legenda_brbg, name = "Ecorregião (Olson):") +
  
  # Habilita a segunda escala de cores para os pontos não "brigarem" com o fundo
  new_scale_fill() +
  
  # --- CAMADA 2: PONTOS DOS SUBGRUPOS ---
  geom_sf(data = points_sub_sf, aes(fill = Subgrupo, shape = Subgrupo), color = "black", stroke = 0.4, size = 2.5, alpha = 0.9) + 
  scale_fill_manual(values = paleta_sub) + 
  scale_shape_manual(values = formas_dinamicas) + 
  
  # Aplica o zoom cravado nas coordenadas dos pontos com a margem calculada
  coord_sf(xlim = limites_x, ylim = limites_y, expand = FALSE) + 
  
  theme_minimal() + 
  labs(
    title = paste("Estruturação Espacial -", dominio_alvo), 
    subtitle = paste("Projeção focada na nuvem de dispersão dos", k_sub, "subgrupos florísticos"),
    fill = "Subgrupo:",
    shape = "Subgrupo:"
  ) + 
  theme(
    legend.position = "right", 
    plot.title = element_text(face = "bold", size = 14), 
    panel.grid.minor = element_blank()
  )
  
print(mapa_subgrupos)

6.4.4. Conectividade Florística Espacial dos Subgrupos (MST)

Projeção de uma rede de ligação mínima (Minimum Spanning Tree) baseada nas distâncias florísticas de Jaccard entre os subgrupos. Os centróides geográficos foram calculados a partir da média espacial das coordenadas amostrais de cada unidade florística.

cat("\nColapsando matriz florística e calculando centróides geográficos dos subgrupos...\n")

## 
## Colapsando matriz florística e calculando centróides geográficos dos subgrupos...

# ==============================================================================
# 1. COLAPSAR A MATRIZ PARA O NÍVEL DE SUBGRUPO
# ==============================================================================
df_matriz_alvo <- as.data.frame(matrix_alvo_limpa)
df_matriz_alvo$locality <- rownames(df_matriz_alvo)

df_grupos_match <- data.frame(locality = names(grupos_sub), Subgrupo = paste("Agrupamento", grupos_sub))

matriz_sub_agrupada <- df_matriz_alvo %>% 
  left_join(df_grupos_match, by = "locality") %>% 
  select(-locality) %>% 
  group_by(Subgrupo) %>% 
  summarise(across(everything(), sum)) %>% 
  as.data.frame()

rownames(matriz_sub_agrupada) <- matriz_sub_agrupada$Subgrupo
matriz_sub_agrupada <- matriz_sub_agrupada[, -1]

# Matriz binária de presença/ausência (1 ou 0) para o subgrupo
matriz_sub_binaria <- ifelse(matriz_sub_agrupada > 0, 1, 0)
riqueza_sub <- rowSums(matriz_sub_binaria)

# Calcular similaridade de Jaccard e a Rede (MST)
dist_sub <- vegdist(matriz_sub_binaria, method = "jaccard")
arvore_sub <- spantree(dist_sub)

# ==============================================================================
# 2. CALCULAR CENTRÓIDES GEOGRÁFICOS
# ==============================================================================
# O centróide é a média das latitudes e longitudes de todos os pontos do agrupamento
df_centroides_sub <- df_sub_coords %>%
  group_by(Subgrupo) %>%
  summarise(X = mean(long), Y = mean(lat), .groups = 'drop') %>%
  mutate(Riqueza = riqueza_sub[Subgrupo])

# Garantir a mesma ordem da árvore MST
df_centroides_sub <- df_centroides_sub[match(rownames(matriz_sub_binaria), df_centroides_sub$Subgrupo), ]

# ==============================================================================
# 3. CRIAR OS RAMOS (LINHAS) DA MST
# ==============================================================================
df_ramos_sub <- data.frame(
  X_Origem  = df_centroides_sub$X[2:nrow(df_centroides_sub)], 
  Y_Origem  = df_centroides_sub$Y[2:nrow(df_centroides_sub)], 
  X_Destino = df_centroides_sub$X[arvore_sub$kid], 
  Y_Destino = df_centroides_sub$Y[arvore_sub$kid],
  Similaridade_Jaccard = 1 - arvore_sub$dist
)

# Ponto central da linha para colocar o rótulo de texto
df_ramos_sub$X_Mid <- (df_ramos_sub$X_Origem + df_ramos_sub$X_Destino) / 2
df_ramos_sub$Y_Mid <- (df_ramos_sub$Y_Origem + df_ramos_sub$Y_Destino) / 2
df_ramos_sub$Jaccard_Rotulo <- sprintf("%.2f", df_ramos_sub$Similaridade_Jaccard)

# ==============================================================================
# 4. PLOTAGEM DO MAPA (COM ZOOM E MÚLTIPLAS ESCALAS)
# ==============================================================================
mapa_mst_sub <- ggplot() + 
  
  # A. Camada de Fundo: Ecorregiões originais (mais opacas para não poluir)
  geom_sf(data = shape_eco_originais, aes(fill = ECO_NAME), alpha = 0.5, color = NA) + 
  scale_fill_manual(values = cores_atuais, breaks = ordem_legenda_brbg, name = "Ecorregião (Olson):") +
  new_scale_fill() +
  
  # B. Pontos reais das localidades (pequenos e no fundo)
  geom_sf(data = points_sub_sf, aes(color = Subgrupo), size = 1.2, alpha = 0.4, show.legend = FALSE) +
  scale_color_manual(values = paleta_sub) +
  
  # C. Ramos (Linhas do MST) conectando os centróides
  geom_segment(data = df_ramos_sub, aes(x = X_Origem, y = Y_Origem, xend = X_Destino, yend = Y_Destino, linewidth = Similaridade_Jaccard), color = "gray15", alpha = 0.85) +
  geom_label(data = df_ramos_sub, aes(x = X_Mid, y = Y_Mid, label = Jaccard_Rotulo), size = 3, fill = "white", color = "black", fontface = "bold", label.padding = unit(0.15, "lines")) +
  
  # D. Centróides (Nós da rede)
  geom_point(data = df_centroides_sub, aes(x = X, y = Y, size = Riqueza, fill = Subgrupo), shape = 21, color = "black", stroke = 1.2) +
  scale_fill_manual(values = paleta_sub, name = "Centróide Florístico:") + 
  scale_size_continuous(range = c(5, 12), name = "Riqueza de Espécies") +
  scale_linewidth_continuous(range = c(0.8, 2.5), guide = "none") +
  
  # E. Zoom Geográfico (Mantém a Bounding Box do Chunk 6.4.3)
  coord_sf(xlim = limites_x, ylim = limites_y, expand = FALSE) + 
  
  theme_minimal() + 
  labs(
    title = paste("Conectividade Espacial (MST) -", dominio_alvo), 
    subtitle = "Rede de similaridade florística conectando os centróides geográficos dos subgrupos",
    x = "Longitude", y = "Latitude"
  ) + 
  theme(
    legend.position = "right", 
    plot.title = element_text(face = "bold", size = 14), 
    panel.grid.minor = element_blank()
  )

print(mapa_mst_sub)

6.5 Análise do domínio B em k=3

Reavaliação das subdivisões florísticas internas do Domínio B impondo manualmente a quebra em 3 subgrupos (\(k=3\)), espelhando a estruturação esperada ou discutida em literaturas correlatas.

dominio_B_name <- "Floristic Domain B"
matrix_dominio_B <- matrix_run[my_data_group == dominio_B_name, ]
matrix_B_limpa <- matrix_dominio_B[, colSums(matrix_dominio_B) > 0]

dissimilarity_B <- vegdist(matrix_B_limpa, method = "jaccard")
tree_B <- hclust(dissimilarity_B, method = "ward.D2")

k_sub_B <- 3 
paleta_sub_B <- c("#00468B", "#ED0000", "#0099B4", "#925E9F", "#FDAF91", "#AD002A")[1:k_sub_B]

dend_base_B <- as.dendrogram(tree_B)
grupos_sub_B <- cutree(dend_base_B, k = k_sub_B, order_clusters_as_data = FALSE)
dend_base_colorido_B <- color_branches(dend_base_B, k = k_sub_B, col = paleta_sub_B)

dend_B_circ <- dendrapply(dend_base_colorido_B, function(node) {
  rotulos_no <- labels(node)
  if (length(rotulos_no) > 0 && all(rotulos_no %in% local_target)) {
    ep <- attr(node, "edgePar")
    if (is.null(ep)) attr(node, "edgePar") <- list(lwd = 3) else { ep$lwd <- 3; attr(node, "edgePar") <- ep }
  }
  if (is.leaf(node)) attr(node, "label") <- "" 
  return(node)
})

par(mar = c(2, 2, 2, 8))
suppressWarnings(circlize_dendrogram(dend_B_circ, labels = FALSE, main = paste("Subdivisões", dominio_B_name)))
legend("center", legend = paste("Agrupamento", 1:k_sub_B), fill = paleta_sub_B, border = rep("black", k_sub_B), bty = "n", cex = 0.75, title = "Subgrupos")

for (i in 1:k_sub_B) {
  locs_grupo <- names(grupos_sub_B[grupos_sub_B == i])
  locs_remover <- labels(dend_base_colorido_B)[!labels(dend_base_colorido_B) %in% locs_grupo]
  dend_individual <- prune(dend_base_colorido_B, locs_remover)
  dend_individual <- set(dend_individual, "labels_cex", 0.5)
  dend_individual <- set(dend_individual, "branches_lwd", 2)
  
  dend_individual <- dendrapply(dend_individual, function(node) {
    rotulos_no <- labels(node)
    if (length(rotulos_no) > 0 && all(rotulos_no %in% local_target)) {
      ep <- attr(node, "edgePar")
      if (is.null(ep)) attr(node, "edgePar") <- list(lwd = 3) else { ep$lwd <- 3; attr(node, "edgePar") <- ep }
    }
    if (is.leaf(node)) {
      rotulo <- attr(node, "label")
      if (length(rotulo) > 0 && rotulo %in% local_target) {
        np <- attr(node, "nodePar")
        if (is.null(np)) attr(node, "nodePar") <- list(lab.col = "black", lab.font = 2) else { np$lab.col <- "black"; np$lab.font <- 2; attr(node, "nodePar") <- np }
      }
    }
    return(node)
  })
  
  par(mar = c(5, 2, 4, 15))
  plot(dend_individual, horiz = TRUE, main = paste("Agrupamento", i, ":", dominio_B_name), xlab = "Distância de Ward")
}

Anotações: O detalhamento em “lupas horizontais” isola os ramos do dendrograma em destaque, permitindo que se avalie visualmente e com extrema precisão como as unidades amostrais da região do Chaco se distribuem nestes três estratos florísticos convergentes.

7. Exportação de Resultados e Produtos Analíticos

Exportação consolidada do Material Suplementar. Reúne tabelas de interesse e salva painéis gráficos gerados em dataframes do Excel e imagens em formato PNG com alta resolução de publicação (300 DPI).

# ==============================================================================
# 7.1. PREPARAÇÃO DO DIRETÓRIO DE EXPORTAÇÃO E PACOTES
# ==============================================================================
if(!dir.exists("Resultados_Exportados")) {
  dir.create("Resultados_Exportados")
}
if(!require(writexl)) install.packages("writexl")
if(!require(svglite)) install.packages("svglite") # O melhor renderizador de vetores editáveis do R
library(writexl)
library(svglite)

# ==============================================================================
# 7.2. CATEGORIA I: TABELAS E DADOS SUPLEMENTARES
# ==============================================================================
cat("Consolidando dataframes em planilhas multi-abas...\n")

## Consolidando dataframes em planilhas multi-abas...

df_gbif_texto <- if(exists("df_ocorrencias")) df_ocorrencias else data.frame(Aviso="Rodar item 6.2")
if(inherits(df_gbif_texto, "sf")) df_gbif_texto <- st_drop_geometry(df_gbif_texto)

df_indicadoras_completas_k2 <- if(exists("res_sig")) { res_sig %>% select(dominio, especie, stat, p.value) %>% arrange(dominio, desc(stat)) } else data.frame()
df_indicadoras_completas_sub <- if(exists("res_sub_sig")) { res_sub_sig %>% select(Subgrupo, especie, stat, p.value) %>% arrange(Subgrupo, desc(stat)) } else data.frame()

lista_tabelas <- list(
  "1_Matriz_Curada" = if(exists("matrix_run")) data.frame(Localidade = rownames(matrix_run), matrix_run, row.names = NULL) else data.frame(),
  "2_Outliers" = if(exists("info_outliers")) info_outliers else data.frame(Aviso="Sem outliers"),
  "3_Dados_Venn" = if(exists("df_tabela_venn")) df_tabela_venn else data.frame(Aviso="Sem dados"),
  "4_Dominios" = if(exists("df_lista_grupos")) df_lista_grupos else data.frame(),
  "5_Indicadoras_Gerais_k2" = df_indicadoras_completas_k2,
  "6_Otimizacao_K" = if(exists("df_resultados_sub")) df_resultados_sub else data.frame(),
  "7_Indicadoras_Subgrupos" = df_indicadoras_completas_sub,
  "8_Validacao_GBIF" = df_gbif_texto
)

write_xlsx(lista_tabelas, path = "Resultados_Exportados/Material_Suplementar_Tabelas.xlsx")
cat("✓ Tabelas consolidadas com sucesso.\n")

## ✓ Tabelas consolidadas com sucesso.

# ==============================================================================
# 7.3. CATEGORIA II: MAPAS ESPACIAIS E BIOGEOGRÁFICOS
# ==============================================================================
cat("\nRenderizando gráficos e mapas em PNG (Alta Res) e SVG (Editável)...\n")

## 
## Renderizando gráficos e mapas em PNG (Alta Res) e SVG (Editável)...

if(exists("mapa_amostras_iniciais")) {
  ggsave("Resultados_Exportados/00_Mapa_Esforco_Amostral.png", plot = mapa_amostras_iniciais, width = 12, height = 8, dpi = 300, bg = "white")
  ggsave("Resultados_Exportados/00_Mapa_Esforco_Amostral.svg", plot = mapa_amostras_iniciais, width = 12, height = 8, bg = "white")
}
if(exists("mapa_ramos_colapsados")) {
  ggsave("Resultados_Exportados/15_MST_Ecorregioes.png", plot = mapa_ramos_colapsados, width = 12, height = 8, dpi = 300, bg="white")
  ggsave("Resultados_Exportados/15_MST_Ecorregioes.svg", plot = mapa_ramos_colapsados, width = 12, height = 8, bg="white")
}
if(exists("mapa_validacao")) { 
  ggsave("Resultados_Exportados/20_Mapa_GBIF.png", plot = mapa_validacao, width = 12, height = 8, dpi = 300, bg="white") 
  
  # Adicionado o svglite para o mapa também
  ggsave("Resultados_Exportados/20_Mapa_GBIF.svg", plot = mapa_validacao, width = 12, height = 8, bg="white", device = svglite::svglite) 
}

# Novos Mapas Espaciais dos Subgrupos (Domínio Alvo)
if(exists("mapa_subgrupos")) { 
  ggsave("Resultados_Exportados/22b_Mapa_Subgrupos.png", plot = mapa_subgrupos, width = 12, height = 8, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/22b_Mapa_Subgrupos.svg", plot = mapa_subgrupos, width = 12, height = 8, bg="white") 
}
if(exists("mapa_mst_sub")) { 
  ggsave("Resultados_Exportados/22c_MST_Subgrupos.png", plot = mapa_mst_sub, width = 12, height = 8, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/22c_MST_Subgrupos.svg", plot = mapa_mst_sub, width = 12, height = 8, bg="white") 
}

# ==============================================================================
# 7.4. CATEGORIA III: ORDENAÇÕES MULTIVARIADAS
# ==============================================================================
if(exists("df_spider") && exists("centroides_nmds") && exists("df_nmds_blind")) { 
  p_nmds_scat <- ggplot(df_spider, aes(x = NMDS1, y = NMDS2)) + 
    geom_segment(aes(xend = Centro_X, yend = Centro_Y, color = Floristic_Group), alpha = 0.2, linewidth = 0.3) + 
    geom_point(aes(color = Point_Color, alpha = Is_Outlier), size = 3) + 
    scale_alpha_manual(values = c("Sim" = 1, "Não" = 0.4), guide = "none") + 
    geom_point(data = centroides_nmds, aes(x = Centro_X, y = Centro_Y, fill = Floristic_Group), size = 4.5, shape = 23, color = "black", stroke = 0.8) + 
    geom_text_repel(data = subset(df_spider, Is_Outlier == "Sim"), aes(label = Locality), size = 3.5, fontface = "bold", box.padding = 0.6, point.padding = 0.3, max.overlaps = Inf) + 
    theme_bw() + labs(title = "NMDS: Dispersão Florística a partir do Núcleo", subtitle = "Associação ao complexo fitogeográfico correspondente", x = "Eixo NMDS 1", y = "Eixo NMDS 2", color = "Legenda") + 
    scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")
  ggsave("Resultados_Exportados/07_NMDS_Scatter.png", plot = p_nmds_scat, width = 10, height = 7, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/07_NMDS_Scatter.svg", plot = p_nmds_scat, width = 10, height = 7, bg="white") 
  
  p_nmds_box <- ggplot(df_nmds_blind, aes(x = Floristic_Group, y = NMDS1, fill = Floristic_Group)) + 
    geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 2, outlier.color = "red") + 
    theme_bw() + labs(title = "Separacao Floristica ao Longo do Gradiente", x = "Dominio Fitogeografico (k=2)", y = "Eixo 1 (NMDS1)") + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77")) + theme(legend.position = "none", plot.title = element_text(face = "bold"))
  ggsave("Resultados_Exportados/08_NMDS_Boxplot.png", plot = p_nmds_box, width = 8, height = 6, dpi = 300, bg="white")
  ggsave("Resultados_Exportados/08_NMDS_Boxplot.svg", plot = p_nmds_box, width = 8, height = 6, bg="white")
}

if(exists("df_pcoa_spider") && exists("df_pcoa_blind") && exists("centroides_pcoa")) { 
  p_pcoa_scat <- ggplot(df_pcoa_spider, aes(x = PCoA1, y = PCoA2)) + 
    geom_segment(aes(xend = Centro_X, yend = Centro_Y, color = Floristic_Group), alpha = 0.2, linewidth = 0.3) + 
    geom_point(aes(color = Point_Color, alpha = Is_Outlier), size = 3) + scale_alpha_manual(values = c("Sim" = 1, "Não" = 0.4), guide = "none") + 
    geom_point(data = centroides_pcoa, aes(x = Centro_X, y = Centro_Y, fill = Floristic_Group), size = 4, shape = 23, color = "black", stroke = 0.5) + 
    geom_text_repel(data = subset(df_pcoa_spider, Is_Outlier == "Sim"), aes(label = Locality), size = 3.5, fontface = "bold", box.padding = 0.6, point.padding = 0.3, max.overlaps = Inf) + 
    theme_bw() + labs(title = "PCoA: Dispersão Florística a partir do Núcleo", subtitle = "Ordenação métrica preservando distâncias originais", x = "PCoA 1", y = "PCoA 2", color = "Legenda") + 
    scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")
  ggsave("Resultados_Exportados/09_PCoA_Scatter.png", plot = p_pcoa_scat, width = 10, height = 7, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/09_PCoA_Scatter.svg", plot = p_pcoa_scat, width = 10, height = 7, bg="white") 
  
  p_pcoa_box <- ggplot(df_pcoa_blind, aes(x = Floristic_Group, y = PCoA1, fill = Floristic_Group)) + 
    geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 2, outlier.color = "red") + 
    theme_bw() + labs(title = "Separacao Floristica ao Longo do Gradiente", x = "Dominio Fitogeografico (k=2)", y = "Eixo 1 (PCoA1)") + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77")) + theme(legend.position = "none", plot.title = element_text(face = "bold"))
  ggsave("Resultados_Exportados/10_PCoA_Boxplot.png", plot = p_pcoa_box, width = 8, height = 6, dpi = 300, bg="white")
  ggsave("Resultados_Exportados/10_PCoA_Boxplot.svg", plot = p_pcoa_box, width = 8, height = 6, bg="white")
}

if(exists("df_perm_scatter") && exists("df_permanova_blind")) { 
  p_perm_spider <- ggplot(df_perm_scatter) + 
    geom_segment(aes(x = Centroid1, y = Centroid2, xend = PCoA1, yend = PCoA2, color = Floristic_Group), alpha = 0.2, linewidth = 0.3) + 
    geom_point(aes(x = PCoA1, y = PCoA2, color = Point_Color, alpha = Is_Outlier), size = 3) + 
    scale_alpha_manual(values = c("Sim" = 1, "Não" = 0.4), guide = "none") + 
    geom_point(aes(x = Centroid1, y = Centroid2, fill = Floristic_Group), size = 4, shape = 23, color = "black", stroke = 0.5) + 
    geom_text_repel(data = subset(df_perm_scatter, Is_Outlier == "Sim"), aes(x = PCoA1, y = PCoA2, label = Locality), size = 3.5, fontface = "bold", box.padding = 0.6, point.padding = 0.3, max.overlaps = Inf) + 
    theme_bw() + labs(title = "Dispersão da PERMANOVA", subtitle = "Associação baseada na variância multivariada", x = "Eixo 1", y = "Eixo 2", color = "Legenda") + 
    scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")
  ggsave("Resultados_Exportados/11_PERMANOVA_Spider.png", plot = p_perm_spider, width = 10, height = 7, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/11_PERMANOVA_Spider.svg", plot = p_perm_spider, width = 10, height = 7, bg="white") 
  
  p_perm_box <- ggplot(df_permanova_blind, aes(x = Floristic_Group, y = Centroid_Distance, fill = Floristic_Group)) + 
    geom_boxplot(alpha = 0.8, outlier.shape = 16, outlier.size = 2, outlier.color = "red") + 
    theme_bw() + labs(title = "Variância Interna dos Grupos (PERMANOVA)", x = "Domínio Fitogeográfico (k=2)", y = "Distância Florística") + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77")) + theme(legend.position = "none", plot.title = element_text(face = "bold"))
  ggsave("Resultados_Exportados/12_PERMANOVA_Boxplot.png", plot = p_perm_box, width = 8, height = 6, dpi = 300, bg="white")
  ggsave("Resultados_Exportados/12_PERMANOVA_Boxplot.svg", plot = p_perm_box, width = 8, height = 6, bg="white")
}

if(exists("df_mantel_visual")) { 
  p_mantel <- ggplot(df_mantel_visual, aes(x = Dist_Geo, y = Dist_Flor)) + 
    geom_point(alpha = 0.08, color = "#2c3e50", size = 0.4) + 
    geom_smooth(method = "lm", formula = y ~ x, color = "red", linetype = "dashed", linewidth = 1, se = FALSE) + 
    theme_bw() + labs(title = "Isolamento por Distancia", x = "Distancia Geografica Real (km)", y = "Dissimilaridade Floristica (Jaccard)")
  ggsave("Resultados_Exportados/13_Mantel_Scatter.png", plot = p_mantel, width = 8, height = 6, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/13_Mantel_Scatter.svg", plot = p_mantel, width = 8, height = 6, bg="white") 
}

if(exists("df_dca")) { 
  p_dca <- ggplot(df_dca, aes(x = DCA1, y = DCA2)) + 
    stat_ellipse(aes(color = Floristic_Group, fill = Floristic_Group), geom = "polygon", alpha = 0.15, linetype = "solid", linewidth = 0.8) + 
    geom_point(aes(color = Point_Color), size = 3, alpha = 0.7) + 
    theme_bw() + labs(title = "Detrended Correspondence Analysis (DCA)", subtitle = "Ordenacao das localidades para avaliacao do gradiente ambiental", x = "Eixo DCA 1", y = "Eixo DCA 2") + 
    scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77"), guide = "none") + theme(legend.position = "bottom")
  ggsave("Resultados_Exportados/14_DCA_Scatter.png", plot = p_dca, width = 10, height = 7, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/14_DCA_Scatter.svg", plot = p_dca, width = 10, height = 7, bg="white") 
}

# ==============================================================================
# 7.5. CATEGORIA IV: GRÁFICOS EXPLORATÓRIOS E DE COMPOSIÇÃO
# ==============================================================================
if(exists("percent_ocur")) { 
  png("Resultados_Exportados/01_Rank_Frequency.png", width=2400, height=1800, res=300)
  plot(percent_ocur, type = "l", lwd = 2, col = "blue", main = "Curva de Frequencia-Ordem das Angiospermas", xlab = "Ordenacao das Especies (da mais comum a mais rara)", ylab = "Frequencia de Ocorrencia (%)")
  abline(h = 5, col = "red", lty = 2)
  dev.off() 
  
  svglite("Resultados_Exportados/01_Rank_Frequency.svg", width=8, height=6)
  plot(percent_ocur, type = "l", lwd = 2, col = "blue", main = "Curva de Frequencia-Ordem das Angiospermas", xlab = "Ordenacao das Especies (da mais comum a mais rara)", ylab = "Frequencia de Ocorrencia (%)")
  abline(h = 5, col = "red", lty = 2)
  dev.off() 
}

## png 
##   2

if(exists("matrix_PA_clean")) { 
  matrix_image <- t(as.matrix(matrix_PA_clean))[, nrow(matrix_PA_clean):1]
  png("Resultados_Exportados/02_Heatmap.png", width=3000, height=2400, res=300)
  image(x = 1:ncol(matrix_PA_clean), y = 1:nrow(matrix_PA_clean), z = matrix_image, col = c("white", "#cb181d"), main = "Mapa de Calor da Matriz de Presença/Ausência", xlab = "Espécies", ylab = "Localidades")
  dev.off()
  
  svglite("Resultados_Exportados/02_Heatmap.svg", width=10, height=8)
  image(x = 1:ncol(matrix_PA_clean), y = 1:nrow(matrix_PA_clean), z = matrix_image, col = c("white", "#cb181d"), main = "Mapa de Calor da Matriz de Presença/Ausência", xlab = "Espécies", ylab = "Localidades")
  dev.off()
}

## png 
##   2

if(exists("matrix_PA_clean") && exists("tree_loc")) {
  pheatmap(as.matrix(matrix_PA_clean), filename = "Resultados_Exportados/02b_Heatmap_Clusterizado.png", width = 12, height = 10, color = c("#f7f7f7", "#cb181d"), cluster_rows = tree_loc, cluster_cols = tree_sp, show_rownames = FALSE, show_colnames = FALSE, border_color = NA, legend_breaks = c(0, 1), legend_labels = c("Ausência", "Presença"), main = "Heatmap Florístico Clusterizado (Localidades vs Espécies)")
  
  svglite("Resultados_Exportados/02b_Heatmap_Clusterizado.svg", width = 12, height = 10)
  pheatmap(as.matrix(matrix_PA_clean), color = c("#f7f7f7", "#cb181d"), cluster_rows = tree_loc, cluster_cols = tree_sp, show_rownames = FALSE, show_colnames = FALSE, border_color = NA, legend_breaks = c(0, 1), legend_labels = c("Ausência", "Presença"), main = "Heatmap Florístico Clusterizado (Localidades vs Espécies)")
  dev.off()
}

## pdf 
##   3

if(exists("regional_rare")) { 
  png("Resultados_Exportados/03_Rarefaction.png", width=2400, height=1800, res=300)
  plot(regional_rare, ci.type = "poly", col = "darkgreen", lwd = 2, ci.lty = 0, ci.col = "lightgreen", main = "Curva de Acumulacao Regional de Especies", xlab = "Esforco Amostral", ylab = "Riqueza Acumulada de Especies")
  grid()
  dev.off() 
  
  svglite("Resultados_Exportados/03_Rarefaction.svg", width=8, height=6)
  plot(regional_rare, ci.type = "poly", col = "darkgreen", lwd = 2, ci.lty = 0, ci.col = "lightgreen", main = "Curva de Acumulacao Regional de Especies", xlab = "Esforco Amostral", ylab = "Riqueza Acumulada de Especies")
  grid()
  dev.off() 
}

## pdf 
##   3

if(exists("list_venn_dinamic")) { 
  png("Resultados_Exportados/04_Venn_Biomas.png", width=2400, height=2400, res=300)
  venn(list_venn_dinamic, zcolor = "style", opacity = 0.4, ilabels = "counts", ilcs = 0.5, sncs = 0.8, box = FALSE, borders = FALSE)
  dev.off() 
  
  svglite("Resultados_Exportados/04_Venn_Biomas.svg", width=8, height=8)
  venn(list_venn_dinamic, zcolor = "style", opacity = 0.4, ilabels = "counts", ilcs = 0.5, sncs = 0.8, box = FALSE, borders = FALSE)
  dev.off() 
}

## pdf 
##   3

if(exists("df_diagnostic")) { 
  p_zscore <- ggplot(df_diagnostic, aes(x = Riqueza, y = Z_Score, color = Status)) + 
    geom_point(alpha = 0.7, size = 2.5) + 
    geom_hline(yintercept = 2.0, linetype = "dashed", linewidth = 1, color = "darkred") + 
    annotate("text", x = max(df_diagnostic$Riqueza) * 0.5, y = 2.25, label = "Corte Multivariado (Z > 2.0)", color = "darkred", fontface = "bold") + 
    geom_vline(xintercept = 4.5, linetype = "dotted", linewidth = 1, color = "blue") + 
    annotate("text", x = 1.5, y = max(df_diagnostic$Z_Score) * 0.85, label = "Corte Biologico (< 5 spp)", color = "blue", fontface = "bold", angle = 90, hjust = 1) + 
    theme_bw() + labs(title = "Diagnostico de Curadoria da Matriz Floristica", subtitle = "Identificacao de localidades com desvio multivariado extremo ou subamostragem", x = "Riqueza Floristica", y = "Distancia Media Multivariada (Z-Score)") + 
    scale_color_manual(values = c("Excluída (Outlier/Subamostrada)" = "red", "Mantida (Matriz Final)" = "grey40")) + 
    theme(legend.position = "bottom", plot.title = element_text(face = "bold"))
  ggsave("Resultados_Exportados/05_Diagnostico_ZScore.png", plot = p_zscore, width = 10, height = 7, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/05_Diagnostico_ZScore.svg", plot = p_zscore, width = 10, height = 7, bg="white") 
}

if(exists("df_flow")) { 
  p_alluvial <- ggplot(df_flow, aes(axis1 = Dominio_Original, axis2 = Dominio_Cluster)) + 
    geom_alluvium(aes(fill = Flow_Color), width = 1/8, alpha = 0.7, knot.pos = 0.4) + 
    geom_stratum(width = 1/4, fill = "grey80", color = "black") + 
    geom_text_repel(stat = "stratum", aes(label = str_wrap(after_stat(stratum), width = 18)), size = 3.5, direction = "y", nudge_x = 0) + 
    scale_x_discrete(limits = c("Bioma Original", "Grupo Floristico (k=2)"), expand = c(.05, .05)) + 
    scale_fill_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77", "Target Areas" = "black")) + 
    theme_minimal() + labs(title = "Fluxo de Reclassificacao Floristica", y = "Numero de Localidades") + theme(legend.position = "bottom")
  ggsave("Resultados_Exportados/16_Alluvial_Flux.png", plot = p_alluvial, width = 10, height = 8, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/16_Alluvial_Flux.svg", plot = p_alluvial, width = 10, height = 8, bg="white") 
}

# (ATUALIZADO: Segundo Diagrama de Venn com os valores de intersecção)
if(exists("lista_venn_dominios")) { 
  png("Resultados_Exportados/17_Venn_Dominios_k2.png", width=2400, height=2400, res=300)
  venn(lista_venn_dominios, zcolor = c("#d95f02", "#1b9e77"), opacity = 0.5, ilabels = "counts", ilcs = 1.5, sncs = 1.2, box = FALSE)
  title("Compartilhamento de Especies entre Dominios")
  dev.off() 
  
  svglite("Resultados_Exportados/17_Venn_Dominios_k2.svg", width=8, height=8)
  venn(lista_venn_dominios, zcolor = c("#d95f02", "#1b9e77"), opacity = 0.5, ilabels = "counts", ilcs = 1.5, sncs = 1.2, box = FALSE)
  title("Compartilhamento de Especies entre Dominios")
  dev.off()
}

## pdf 
##   3

if(exists("df_beta_longo")) { 
  p_beta <- ggplot(df_beta_longo, aes(x = Dominio, y = Valor_Beta, fill = Componente)) + 
    geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7, color = "black") + 
    geom_text(aes(label = Label), position = position_dodge(width = 0.8), vjust = -0.5, fontface = "bold", size = 4.5) + 
    theme_bw() + labs(title = "Componentes da Diversidade Beta por Dominio", y = "Valor da Dissimilaridade") + 
    scale_y_continuous(expand = expansion(mult = c(0, 0.15))) + 
    scale_fill_manual(values = c("Turnover" = "#2c3e50", "Nestedness" = "#bdc3c7")) + theme(legend.position = "bottom")
  ggsave("Resultados_Exportados/18_Beta_Diversity.png", plot = p_beta, width = 8, height = 6, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/18_Beta_Diversity.svg", plot = p_beta, width = 8, height = 6, bg="white") 
}

if(exists("top_indicadoras")) { 
  # Atualizado: usa 'especie_clean', shape = dominio e scale_shape_manual
  p_indval <- ggplot(top_indicadoras, aes(x = stat, y = reorder(especie_clean, stat), color = dominio)) + 
    geom_segment(aes(x = 0, xend = stat, y = reorder(especie_clean, stat), yend = reorder(especie_clean, stat)), color = "grey60", linetype = "dashed") + 
    geom_point(aes(shape = dominio), size = 4) + 
    facet_wrap(~dominio, scales = "free_y", ncol = 1) + 
    theme_bw() + 
    labs(title = "Top 10 Especies Indicadoras", x = "Valor de Indicação (stat)", y = "") + 
    scale_color_manual(values = c("Floristic Domain A" = "#d95f02", "Floristic Domain B" = "#1b9e77")) + 
    scale_shape_manual(values = c("Floristic Domain A" = 16, "Floristic Domain B" = 17)) +
    theme(legend.position = "none", axis.text.y = element_text(face = "italic"), strip.background = element_rect(fill = "#f0f0f0"), strip.text = element_text(face = "bold"))
  
  ggsave("Resultados_Exportados/19_IndVal_Dotplot.png", plot = p_indval, width = 8, height = 10, dpi = 300, bg="white") 
  
  # Adicionado o svglite para preservar o texto, caso o sistema permita
  ggsave("Resultados_Exportados/19_IndVal_Dotplot.svg", plot = p_indval, width = 8, height = 10, bg="white", device = svglite::svglite) 
}

if(exists("grafico_otimizacao_sub")) { 
  ggsave("Resultados_Exportados/22_Otimizacao_K.png", plot = grafico_otimizacao_sub, width = 8, height = 10, dpi = 300, bg="white") 
  ggsave("Resultados_Exportados/22_Otimizacao_K.svg", plot = grafico_otimizacao_sub, width = 8, height = 10, bg="white") 
}

# (ATUALIZADO: Exporta o novo gráfico de indicadoras com formas e cores dos subgrupos)
if(exists("grafico_indval_sub")) {
  ggsave("Resultados_Exportados/25_IndVal_Subgrupos.png", plot = grafico_indval_sub, width = 14, height = 8, dpi = 300, bg="white")
  ggsave("Resultados_Exportados/25_IndVal_Subgrupos.svg", plot = grafico_indval_sub, width = 14, height = 8, bg="white")
}

# ==============================================================================
# 7.6. CATEGORIA VI: DENDROGRAMAS E ESTRUTURAS HIERÁRQUICAS
# ==============================================================================
if(exists("dendro_colorido")) { 
  png("Resultados_Exportados/06_Dendrograma_k2.png", width=3600, height=1800, res=300)
  par(mar = c(2, 4, 3, 1))
  plot(dendro_colorido, main = "Dendrograma Floristico: Separacao dos Dominios (k=2)", ylab = "Distancia de Ward", xlab = "Localidades")
  title(xlab = "Localidades", line = 0.5, font.lab = 2)
  legend("topright", legend = c("Floristic Domain A", "Floristic Domain B", "Target Areas"), fill = c("#d95f02", "#1b9e77", "white"), border = c("black", "black", "white"), pch = c(NA, NA, "|"), col = c(NA, NA, "black"), pt.cex = 1.2, bty = "n", cex = 0.9, inset = c(0.02, 0.02))
  dev.off() 
  
  svglite("Resultados_Exportados/06_Dendrograma_k2.svg", width=12, height=6)
  par(mar = c(2, 4, 3, 1))
  plot(dendro_colorido, main = "Dendrograma Floristico: Separacao dos Dominios (k=2)", ylab = "Distancia de Ward", xlab = "Localidades")
  title(xlab = "Localidades", line = 0.5, font.lab = 2)
  legend("topright", legend = c("Floristic Domain A", "Floristic Domain B", "Target Areas"), fill = c("#d95f02", "#1b9e77", "white"), border = c("black", "black", "white"), pch = c(NA, NA, "|"), col = c(NA, NA, "black"), pt.cex = 1.2, bty = "n", cex = 0.9, inset = c(0.02, 0.02))
  dev.off() 
}

## pdf 
##   3

if(exists("dend_alvo") && exists("dominio_alvo") && exists("cor_dominio")) { 
  png("Resultados_Exportados/21_Dendrograma_Alvo.png", width=3600, height=1800, res=300)
  par(mar = c(2, 4, 3, 1))
  plot(dend_alvo, main = paste("Estrutura Interna -", dominio_alvo), ylab = "Distância de Ward", xlab = "Localidades")
  title(xlab = "Localidades", line = 0.5, font.lab = 2)
  legend("topright", legend = c(paste("Domínio:", dominio_alvo), "Target Areas"), fill = c(cor_dominio, "white"), border = c("black", "white"), pch = c(NA, "|"), col = c(NA, "black"), pt.cex = 1.2, bty = "n", cex = 0.9, inset = c(0.02, 0.02))
  dev.off() 
  
  svglite("Resultados_Exportados/21_Dendrograma_Alvo.svg", width=12, height=6)
  par(mar = c(2, 4, 3, 1))
  plot(dend_alvo, main = paste("Estrutura Interna -", dominio_alvo), ylab = "Distância de Ward", xlab = "Localidades")
  title(xlab = "Localidades", line = 0.5, font.lab = 2)
  legend("topright", legend = c(paste("Domínio:", dominio_alvo), "Target Areas"), fill = c(cor_dominio, "white"), border = c("black", "white"), pch = c(NA, "|"), col = c(NA, "black"), pt.cex = 1.2, bty = "n", cex = 0.9, inset = c(0.02, 0.02))
  dev.off() 
}

## pdf 
##   3

if(exists("dend_alvo_circ") && exists("paleta_sub")) { 
  png("Resultados_Exportados/23_Dendro_Circular_k_Optimo.png", width=3000, height=3000, res=300)
  par(mar = c(2, 2, 2, 8))
  suppressWarnings(circlize_dendrogram(dend_alvo_circ, labels = FALSE, main = paste("Subdivisões", dominio_alvo)))
  legend("center", legend = paste("Agrupamento", 1:k_sub), fill = paleta_sub, border = rep("black", k_sub), bty = "n", cex = 0.75, title = "Subgrupos")
  dev.off() 
  
  svglite("Resultados_Exportados/23_Dendro_Circular_k_Optimo.svg", width=10, height=10)
  par(mar = c(2, 2, 2, 8))
  suppressWarnings(circlize_dendrogram(dend_alvo_circ, labels = FALSE, main = paste("Subdivisões", dominio_alvo)))
  legend("center", legend = paste("Agrupamento", 1:k_sub), fill = paleta_sub, border = rep("black", k_sub), bty = "n", cex = 0.75, title = "Subgrupos")
  dev.off() 
}

## pdf 
##   3

if(exists("dend_B_circ") && exists("paleta_sub_B")) { 
  png("Resultados_Exportados/24_Dendro_Circular_k3.png", width=3000, height=3000, res=300)
  par(mar = c(2, 2, 2, 8))
  suppressWarnings(circlize_dendrogram(dend_B_circ, labels = FALSE, main = paste("Subdivisões Florísticas do", dominio_B_name)))
  legend("center", legend = paste("Agrupamento", 1:k_sub_B), fill = paleta_sub_B, border = rep("black", k_sub_B), bty = "n", cex = 0.75, title = "Subgrupos")
  dev.off() 
  
  svglite("Resultados_Exportados/24_Dendro_Circular_k3.svg", width=10, height=10)
  par(mar = c(2, 2, 2, 8))
  suppressWarnings(circlize_dendrogram(dend_B_circ, labels = FALSE, main = paste("Subdivisões Florísticas do", dominio_B_name)))
  legend("center", legend = paste("Agrupamento", 1:k_sub_B), fill = paleta_sub_B, border = rep("black", k_sub_B), bty = "n", cex = 0.75, title = "Subgrupos")
  dev.off() 
}

## pdf 
##   3

if(exists("tree_biome")) { 
  write.tree(as.phylo(tree_biome), file = "Resultados_Exportados/26_dendrograma.nwk") 
}

cat("\n✓ Exportação finalizada com sucesso! Todas as imagens foram exportadas em duplas (PNG e SVG editável).\n")

## 
## ✓ Exportação finalizada com sucesso! Todas as imagens foram exportadas em duplas (PNG e SVG editável).

Anotações: A automatização desta última etapa consolida as diretrizes de dados abertos (Open Science). A padronização da saída assegura que as matrizes de confusão, índices biológicos e scatterplots espaciais estejam diretamente aptos e adequados para anexação no periódico acadêmico.