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

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

# Mediterranean forest community
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
)

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")
)

## ----tree, fig.width=9, fig.height=5.5, fig.alt="Dendrogram showing species grouped by family"----
plot_taxonomic_tree(tax_tree, community = community,
                    color_by = "Family", label_size = 3.5,
                    title = "Mediterranean Forest - Taxonomic Tree")

## ----heatmap, fig.width=9, fig.height=8, fig.alt="Heatmap of pairwise taxonomic distances"----
plot_heatmap(tax_tree, label_size = 2.8,
             title = "Pairwise Taxonomic Distances")

## ----bubble, fig.width=10, fig.height=7, fig.alt="Bubble chart of species contributions"----
plot_bubble(community, tax_tree, color_by = "Family",
            title = "Species Contributions to Diversity")

## ----radar_data---------------------------------------------------------------
# Degraded community for comparison
dominant_community <- c(
  Quercus_coccifera   = 80, Quercus_infectoria  = 5,
  Pinus_brutia        = 3,  Pinus_nigra         = 2,
  Juniperus_excelsa   = 2,  Juniperus_oxycedrus = 1,
  Arbutus_andrachne   = 3,  Styrax_officinalis  = 1,
  Cercis_siliquastrum = 2,  Olea_europaea       = 1
)

communities <- list(
  Diverse  = community,
  Dominant = dominant_community
)

## ----radar, fig.width=8, fig.height=8, fig.alt="Radar chart comparing two communities"----
plot_radar(communities, tax_tree,
           title = "Diverse vs Dominant Community")

## ----iteration_data-----------------------------------------------------------
run2 <- ozkan_pto_resample(community, tax_tree, n_iter = 101, seed = 42)

## ----iteration, fig.width=9, fig.height=5, fig.alt="Scatter plot of TO values across iterations"----
plot_iteration(run2, component = "TO",
               title = "Run 2: TO Values Across Iterations")

## ----rarefaction--------------------------------------------------------------
rare <- rarefaction_taxonomic(community, tax_tree,
                               index = "shannon",
                               steps = 10, n_boot = 50, seed = 42)

## ----rarefaction_plot, fig.width=8, fig.height=5, fig.alt="Rarefaction curve with confidence interval"----
plot_rarefaction(rare)

## ----funnel_data--------------------------------------------------------------
data(anatolian_trees)

sim <- simulate_td(
  tax_tree = anatolian_trees,
  s_range = c(3, 15),
  n_sim = 99,
  index = "avtd",
  seed = 42
)

## ----funnel, fig.width=9, fig.height=6, fig.alt="Funnel plot with 95% confidence bands"----
spp <- names(community)
obs_avtd <- avtd(spp, tax_tree)

plot_funnel(sim,
            observed = data.frame(
              site = "Mediterranean",
              s = length(spp),
              value = obs_avtd
            ),
            index = "avtd",
            title = "AvTD Significance Test")

## ----custom, fig.width=8, fig.height=5, fig.alt="Customized rarefaction curve with modified theme"----
library(ggplot2)

plot_rarefaction(rare) +
  theme_minimal() +
  labs(subtitle = "Mediterranean forest community") +
  theme(plot.title = element_text(face = "bold", size = 14))

## ----custom_examples, eval=FALSE----------------------------------------------
# # Change theme
# p + theme_classic()
# 
# # Modify colors
# p + scale_color_brewer(palette = "Set2")
# 
# # Adjust text size
# p + theme(axis.text = element_text(size = 12))
# 
# # Save to file
# ggsave("my_plot.png", p, width = 10, height = 6, dpi = 300)