---
title: "FlexTrees"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    vertical_layout: fill
    theme:
      bootswatch: pulse
    source_code: embed
    favicon: tree-logo-resized.png
---

```{r setup, include=FALSE}
#Libraries
library(flexdashboard)
library(htmltools)
library(tidyverse)
library(broom)
library(car)
library(ggfortify)
library(knitr)
library(tidymodels)
library(vip)
library(performance)
library(dplyr)
library(DT)
library(ggplot2)
library(plotly)
library(corrr)
library(readxl)

# read in data
pine_tbl <- read_excel("Data_1993.xlsx", sheet = 1)

# RIDGE REGRESSION -------------------------------------------------------------

# Create training/testing data
pine_split <- initial_split(pine_tbl)
pine_train <- training(pine_split)
pine_test <- testing(pine_split)

# picking Dr. Smirnova's best lambda estimate; can estimate with tune() - see below
ridge_mod <-
  linear_reg(mixture = 0, penalty = 0.1629751) %>%  #validation sample or resampling can estimate this
  set_engine("glmnet")

# verify what we are doing
ridge_mod %>% 
  translate()

# create a new recipe; could use `add_step()` to recipe created above
pine_rec <- pine_train %>% 
  recipe(DeadDist ~ TreeDiam + Infest_Serv1 + SDI_20th + BA_20th) %>% 
  step_sqrt(all_outcomes()) %>% 
  step_corr(all_predictors()) %>% 
  step_normalize(all_numeric(), -all_outcomes()) %>% 
  step_zv(all_numeric(), -all_outcomes()) #%>% 
  #prep()

pine_ridge_wflow <- 
  workflow() %>% 
  add_model(ridge_mod) %>% 
  add_recipe(pine_rec)

pine_ridge_fit <- 
  pine_ridge_wflow %>% 
  fit(data = pine_train)

pine_ridge_fit %>% 
  extract_fit_parsnip() %>% 
  tidy()

pine_ridge_fit %>% 
  extract_preprocessor()

pine_ridge_fit %>% 
  extract_spec_parsnip()

# refit best model on training and evaluate on testing
last_fit(
  pine_ridge_wflow,
  pine_split
) %>%
  collect_metrics() #get stats


# LASSO REGRESSION
# ----------------------------------------------------

# create bootstrap samples for resampling and tuning the penalty parameter
set.seed(1234)
pine_boot <- bootstraps(pine_train) #resampling technique, 10 bootstrap samples 
# list of tibbles

# create a grid of tuning parameters
lambda_grid <- grid_regular(penalty(), levels = 50)
#created 50 lambda values, and want to see which performs the best over the different bootstraps

lasso_mod <-
  linear_reg(mixture = 1, penalty = tune()) %>% 
  set_engine("glmnet")
#do lasso -> mixture = 1, use tune in penalty

# verify what we are doing
lasso_mod %>% 
  translate()

# create workflow
pine_lasso_wflow <- 
  workflow() %>% 
  add_model(lasso_mod) %>% 
  add_recipe(pine_rec)

set.seed(2020)
lasso_grid <- tune_grid(
  pine_lasso_wflow,
  resamples = pine_boot,
  grid = lambda_grid
)

lowest_rmse <- lasso_grid %>%
  select_best(metric = "rmse")

# update our final model with lowest rmse
final_lasso <- finalize_workflow(
  pine_lasso_wflow,
  lowest_rmse
)

final_lasso %>% 
  fit(pine_train) %>%
  extract_fit_parsnip() %>% 
  tidy()
```


What's the Problem?
====================================================================
Column {data-width=350} 
-----------------------------------------------------------------------

```{r}
datatable(pine_tbl, caption = "Datatable for Pine Beetle Dataset", options = list(scrollX = FALSE, autoWidth = TRUE))
```

Column {data-width=650}
-----------------------------------------------------------------------

### The Culprit

```{r}
knitr::include_graphics("pine-beetle-trees.png")
```

### Attacked Trees by Jeffrey Pine Beetles

```{r, fig.width = 8}
p.all <- ggplot(data = pine_tbl, aes(x = Easting, y = Northing)) +
  geom_point(aes(color = factor(Response), alpha = factor(Response))) +
  scale_alpha_discrete(range = c(0.5, 1), guide = FALSE) +
  scale_color_manual("", values = c("yellowgreen", "red"),
  labels = c("Alive", "JPB-Attacked")) +
  theme_bw() +
  xlab("utm x") +
  ylab("utm y") +
  theme(legend.position = c(0, 1), legend.text = element_text(size = 15))

p.all
```

Explore the Data
====================================================================
Column {data-width=650}
-----------------------------------------------------------------------

### Corr/GGPair Plot - Exploratory Analysis

```{r}
#Corr/GGPair Plot--------------------
library(tidymodels)
library(GGally)
# Extracting predictor variables from the training set
predictor_data <- select(pine_train, TreeDiam, Infest_Serv1, SDI_20th, BA_20th)

# Compute the correlation matrix
correlation_matrix <- cor(predictor_data)

# Plot correlation matrix using ggpair
ggpairs(predictor_data)
```

### LASSO Grid

```{r}
lasso_grid %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, color = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
  ),
  alpha = 0.5
  ) +
  geom_line(size = 1.5) +
  facet_wrap(~.metric, scales = "free", nrow = 2) +
  scale_x_log10() +
  theme(legend.position = "none")
```

Column {data-width=350} 
-----------------------------------------------------------------------
### Variable Importance Plot - LASSO

```{r}
# variable importance plot
final_lasso %>%
  fit(pine_train) %>%
  extract_fit_parsnip() %>%
  vi(lambda = lowest_rmse$penalty) %>%
  mutate(
    Importance = abs(Importance),
    Variable = fct_reorder(Variable, Importance)
  ) %>%
  ggplot(aes(x = Importance, y = Variable, fill = Sign)) +
  geom_col() +
  scale_x_continuous(expand = c(0, 0)) +
  labs(y = NULL)
```

### Variable Importance Plot - Ridge

```{r, fig.width = 5}
#Variable Importance Plot--------------------
# Extract fitted Ridge Regression model object
ridge_model <- extract_fit_parsnip(pine_ridge_fit)

# Get the coefficients from the Ridge Regression model
coefficients <- tidy(ridge_model)

# Select top features based on their coefficient magnitudes
top_features <- coefficients %>%
  arrange(desc(abs(estimate)))

# Plot feature importance
top_features %>%
  ggplot(aes(x = reorder(term, estimate), y = abs(estimate))) +
  geom_point(size = 4, color = "#7700CF") +
  coord_flip() +  # Horizontal bars
  theme_minimal(base_size = 18) +
  labs(title = "Feature Importance")
```

How'd We Get There?
====================================================================
Column {data-width=650}
-----------------------------------------------------------------------

### Gauge: LASSO Regression R-Squared

```{r}
output_lasso <- last_fit(
  pine_ridge_wflow,
  pine_split
) %>%
  collect_metrics()#get stats

rsq_value <- output_lasso %>%
  filter(.metric == "rsq") %>%
  pull(.estimate)

gauge(rsq_value, min = 0, max = 1, gaugeSectors(
  success = c(0.7, 1), warning = c(0.4, 0.6), danger = c(0, 0.3)
))
```

### Gauge: Ridge Regression R-Squared

```{r}
output_ridge <- last_fit(
  pine_ridge_wflow,
  pine_split
) %>%
  collect_metrics()#get stats

rsq_value <- output_ridge %>%
  filter(.metric == "rsq") %>%
  pull(.estimate)

gauge(rsq_value, min = 0, max = 1, gaugeSectors(
  success = c(0.7, 1), warning = c(0.4, 0.6), danger = c(0, 0.3)
))
```


Column {data-width=350} 
-----------------------------------------------------------------------
### Predictor Variables - VIF

```{r}
# Fit linear regression model
lm.fit <- lm(DeadDist ~ TreeDiam + Infest_Serv1 + SDI_20th + BA_20th, data = pine_tbl)

#VIF Table
vif <- vif(lm.fit)
knitr::kable(vif, digits = 3)
```

### Recipe

```{r, echo = TRUE}
pine_rec <- pine_train %>% 
  recipe(DeadDist ~ TreeDiam + Infest_Serv1
         + SDI_20th + BA_20th) %>% 
  step_sqrt(all_outcomes()) %>% 
  step_corr(all_predictors()) %>% 
  step_normalize(all_numeric(), -all_outcomes()) %>% 
  step_zv(all_numeric(), -all_outcomes()) #%>% 
  #prep()
```

CodeBook
====================================================================

```{r picture, echo = F, fig.cap = "CodeBook", out.width = '100%'}
knitr::include_graphics("pine-beetle-variables.png")
```