## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 5
)

## ----setup--------------------------------------------------------------------
library(taxdiv)

## ----single_site--------------------------------------------------------------
community <- c(
  Quercus_coccifera    = 25,
  Quercus_infectoria   = 18,
  Pinus_brutia         = 30,
  Pinus_nigra          = 12,
  Juniperus_excelsa    = 8,
  Juniperus_oxycedrus  = 6,
  Arbutus_andrachne    = 15,
  Styrax_officinalis   = 4,
  Cercis_siliquastrum  = 3,
  Olea_europaea        = 10
)

## ----tree---------------------------------------------------------------------
tax_tree <- build_tax_tree(
  species = names(community),
  Genus   = c("Quercus", "Quercus", "Pinus", "Pinus",
              "Juniperus", "Juniperus", "Arbutus", "Styrax",
              "Cercis", "Olea"),
  Family  = c("Fagaceae", "Fagaceae", "Pinaceae", "Pinaceae",
              "Cupressaceae", "Cupressaceae", "Ericaceae", "Styracaceae",
              "Fabaceae", "Oleaceae"),
  Order   = c("Fagales", "Fagales", "Pinales", "Pinales",
              "Pinales", "Pinales", "Ericales", "Ericales",
              "Fabales", "Lamiales")
)

tax_tree

## ----multi_site---------------------------------------------------------------
multi_data <- data.frame(
  Site = rep(c("Forest_A", "Forest_B"), each = 5),
  Species = c("Quercus_coccifera", "Pinus_brutia", "Juniperus_excelsa",
              "Arbutus_andrachne", "Olea_europaea",
              "Quercus_coccifera", "Quercus_infectoria", "Pinus_brutia",
              "Pinus_nigra", "Cercis_siliquastrum"),
  Abundance = c(25, 30, 8, 15, 10,
                40, 20, 15, 10, 5),
  Genus = c("Quercus", "Pinus", "Juniperus", "Arbutus", "Olea",
            "Quercus", "Quercus", "Pinus", "Pinus", "Cercis"),
  Family = c("Fagaceae", "Pinaceae", "Cupressaceae", "Ericaceae", "Oleaceae",
             "Fagaceae", "Fagaceae", "Pinaceae", "Pinaceae", "Fabaceae"),
  Order = c("Fagales", "Pinales", "Pinales", "Ericales", "Lamiales",
            "Fagales", "Fagales", "Pinales", "Pinales", "Fabales")
)

head(multi_data)

## ----single_compare-----------------------------------------------------------
results <- compare_indices(community, tax_tree)
results

## ----batch--------------------------------------------------------------------
batch_results <- batch_analysis(multi_data)
batch_results

## ----batch_summary------------------------------------------------------------
summary(batch_results)

## ----full_pipeline------------------------------------------------------------
full <- ozkan_pto_full(community, tax_tree, n_iter = 101, seed = 42)

cat("Run 1 (deterministic):\n")
cat("  uTO+ =", round(full$run1$uTO_plus, 4), "\n")
cat("  TO+  =", round(full$run1$TO_plus, 4), "\n\n")

cat("Run 2 (stochastic maximum):\n")
cat("  uTO+ =", round(full$run2$uTO_plus_max, 4), "\n")
cat("  TO+  =", round(full$run2$TO_plus_max, 4), "\n\n")

cat("Run 3 (sensitivity maximum):\n")
cat("  uTO+ =", round(full$run3$uTO_plus_max, 4), "\n")
cat("  TO+  =", round(full$run3$TO_plus_max, 4), "\n")

## ----compare_plot, fig.width=10, fig.height=6, fig.alt="Bar chart comparing diversity indices between intact and degraded forest communities"----
degraded <- c(
  Quercus_coccifera   = 40,
  Pinus_brutia        = 35,
  Juniperus_oxycedrus = 10
)
tax_degraded <- tax_tree[tax_tree$Species %in% names(degraded), ]

communities <- list(
  "Intact Forest" = community,
  "Degraded"      = degraded
)

comp <- compare_indices(communities, tax_tree, plot = TRUE)
comp$plot

## ----tree_plot, fig.width=9, fig.height=5.5, fig.alt="Dendrogram showing community taxonomic structure colored by family"----
plot_taxonomic_tree(tax_tree, community = community,
                    color_by = "Family",
                    title = "Community Taxonomic Structure")

## ----rarefaction, fig.width=8, fig.height=5, fig.alt="Shannon rarefaction curve with bootstrap confidence interval"----
rare <- rarefaction_taxonomic(community, tax_tree,
                               index = "shannon",
                               steps = 10, n_boot = 50, seed = 42)
plot_rarefaction(rare)

## ----export_csv, eval=FALSE---------------------------------------------------
# # Single community
# df <- as.data.frame(t(compare_indices(community, tax_tree)$table))
# write.csv(df, "results.csv", row.names = FALSE)
# 
# # Multi-site batch
# write.csv(as.data.frame(batch_results$results),
#           "batch_results.csv", row.names = FALSE)

## ----export_xlsx, eval=FALSE--------------------------------------------------
# # Requires writexl package
# # install.packages("writexl")
# writexl::write_xlsx(as.data.frame(batch_results$results),
#                     "batch_results.xlsx")