---
title: "HGEN612 - Klimek - Project 02"
output: 
  flexdashboard::flex_dashboard:
    source_code: embed
#    orientation: rows
    vertical_layout: fill 
    theme: yeti #readable - try out different themes (pleasing and aesthetic)
    logo: ../img/beetle.png
    favicon: ../img/beetle.png #Find beetle favicon
#    navbar: maybe ghost the nav bar to just blank it out
---
```{r setup, include=FALSE}
library(flexdashboard)

# Tables
library(DT) #data table

# Plotting
library(ggplot2)
library(plotly)
library(GGally)
library(ggfortify)

# Analysis
library(corrr) #correlation matrix
library(vip) 
library(performance) 

# Other
library(tidyverse)
library(readxl)
library(broom)
library(dplyr)
library(tidymodels)
library(knitr)
library(kableExtra)

```

Readme 
======

### Readme

This dashboard applies two different statistical modelling methods, Lasso and Ridge regression, to data characterizing a Jeffrey pine beetle outbreak in the Lake Tahoe Basin from 1991 to 1996. Around 10,000 trees were assessed annually over a 60 acre study area. 

The outcome of interest investigated in all models of this dashboard was "DeadDist." The predictors "TreeDiam" and "Infest_sever1" remained constant in all models while the basal area of the surrounding neighborhood was adjusted from 1/20th of an acre up to 1.5 acres. 

To perform the statistical models chosen, a feature engineering approach was performed using a tidymodels workflow while employing recipe step_x functions. All models utilized a square root of the outcome "DeadDist" and predictors were normalized to a SD of 1 and mean of 0. 

***Please see the Data Library in the "Data Overview" tab for further information about variables.***

***See recipes for further details on variable transformation.***

Egan, J.M., Sloughter, J.M., Cardoso, T. et al. Multi-temporal ecological analysis of Jeffrey pine beetle outbreak dynamics within the Lake Tahoe Basin. Popul Ecol 58, 441–462 (2016). https://doi.org/10.1007/s10144-016-0545-2











```{r loading data}

pine_tbl <- read_excel("../data/Data_1993.xlsx", sheet = 1)
#remember the "." indicates relative path
#the Rproj helps with using relative paths rather than setting hard working directories

```


Data Overview
======

Column {.tabset .tabset_fade}
----

### Spatial Plot of Healthy Pine Trees and Pine Trees Attacked By Jeffrey Pine Beetles.
```{r data overview plot}

#adding another column for ggplotly to pull names for legend
pine_tbl_explicit <- pine_tbl %>% 
  mutate(dead = ifelse(Response == 1, "JPB-attacked", "Alive"))

#mapping trees dead/alive (Smirnova Lecture)
p.all <- ggplot(data = pine_tbl_explicit, aes(x = Easting, y = Northing)) +
  geom_point(aes(color = factor(dead), 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))

#converting it to plotly
p.all_plotly <- ggplotly(p.all) 

#adjusting legend position
p.all_plotly_legend <- p.all_plotly %>% 
  layout(legend = list(x = 0, y = 1), legend.text = element_text(size = 10))

#data overview plot
p.all_plotly_legend


```

### Data Library 
```{r Variable Library}
knitr::include_graphics("../img/pine-beetle-highlighted-variables.png")

```

### Overview Table
```{r}
datatable(pine_tbl %>% 
  select(TreeNum, Response, TreeDiam, Infest_Serv1, DeadDist, SDI_20th, BA_20th, `Neigh_1/2th`, Neigh_1.5))

```

### GGPAIRS BA_20th
```{r ggpair BA_20th}
ggpairs(pine_tbl[,c("DeadDist", "TreeDiam", "BA_20th")],
        diag = list(continuous ="barDiag", discrete = "barDiag", na = "naDiag") )

```

### GGPAIRS Neigh_1/2th
```{r ggpair Neigh_1/2th}
ggpairs(pine_tbl[,c("DeadDist", "TreeDiam", "Neigh_1/2th")],
        diag = list(continuous ="barDiag", discrete = "barDiag", na = "naDiag") )
```

### GGPAIRS Neigh_1.5

```{r ggpair Neigh1.5}
ggpairs(pine_tbl[,c("DeadDist", "TreeDiam", "Neigh_1.5")],
        diag = list(continuous ="barDiag", discrete = "barDiag", na = "naDiag") )
```

Recipes
=======

Column {.tabset .tabset_fade}
----

### Predictors and Outcomes
```{r, echo = TRUE}
# * SPLITTING DATA ----
set.seed(123)
pine_split <- initial_split(pine_tbl)
pine_train <- training(pine_split)
pine_test <- testing(pine_split)
#default split is 75/25

# * PREDICTORS AND OUTCOMES ----

#recipe with BA_20th
rec_BA20 <- recipe(DeadDist ~ TreeDiam + Infest_Serv1 + BA_20th, data = pine_train)

#recipe with Neigh 1/2th
rec_N05 <- recipe(DeadDist ~ TreeDiam + Infest_Serv1 + `Neigh_1/2th`, data = pine_train)

#recipe with Neigh 1.5
rec_N15 <- recipe(DeadDist ~ TreeDiam + Infest_Serv1 + `Neigh_1.5`, data = pine_train)

```

### PREPROCESS
```{r, echo = TRUE}

# * RECIPE PREPROCESS STEPS ----

#recipe for BA20th
pine_rec_BA20 <- rec_BA20 %>%
  step_sqrt(all_outcomes(), skip = TRUE) %>% #square root outcome
  step_corr(all_predictors()) %>% #VIF assessment and variable removal
  step_normalize(all_numeric(), -all_outcomes()) %>% #center and scale numeric data
  step_zv(all_numeric(), -all_outcomes()) #removes variables containing only a single value
#minus sign is indicating not to do it to the outcomes

#recipe with Neigh 1/2th
pine_rec_N05 <- rec_N05 %>% 
  step_sqrt(all_outcomes(), skip = TRUE) %>% 
  step_corr(all_predictors()) %>% 
  step_normalize(all_numeric(), -all_outcomes()) %>% #minus sign is indicating not to do it to the outcomes
  step_zv(all_numeric(), -all_outcomes())

#recipe with Neigh 1.5
pine_rec_N15 <- rec_N15 %>% 
  step_sqrt(all_outcomes(), skip = TRUE) %>% 
  step_corr(all_predictors()) %>% 
  step_normalize(all_numeric(), -all_outcomes()) %>% #minus sign is indicating not to do it to the outcomes
  step_zv(all_numeric(), -all_outcomes())



```

### ISSUES
```{r, echo=TRUE}
# pine_rec_BA20 <- rec_BA20 %>% 
#   step_sqrt(all_outcomes()) %>% 
#   step_corr(all_predictors()) %>% 
#   step_normalize(all_numeric(), -all_outcomes()) %>% #minus sign is indicating not to do it to the outcomes
#   step_zv(all_numeric(), -all_outcomes())

# ** ISSUE: ----  
# > prediction_test_BA20 <- predict(ridge_pine_fit_BA20, new_data = pine_test)
# Error in `step_sqrt()`:
#   ! The following required column is missing from `new_data` in step 'sqrt_AN0PC':
#   DeadDist.
# Run `rlang::last_trace()` to see where the error occurred.

# Not sure how to resolve, but including skip = TRUE allows it to work 
# however the model then compares feature engineered results to original results that are not transformed
```

Column {.tabset .tabset_fade}
----

### BA_20th Engineered Data Table
```{r}
feature_data_BA20 <- pine_rec_BA20 %>% 
  prep() %>% 
  bake(new_data = pine_train)

pine_test_processed_BA20 <- pine_rec_BA20 %>% 
  prep() %>% 
  bake(new_data = pine_test)

datatable(feature_data_BA20, fillContainer = TRUE)
```

### Neigh_1/2th Engineered Data Table
```{r}
feature_data_N05 <- pine_rec_N05 %>% 
  prep() %>% 
  bake(new_data = pine_train)

pine_test_processed_N05 <- pine_rec_N05 %>% 
  prep() %>% 
  bake(new_data = pine_test)


datatable(feature_data_N05)
```

### Neigh1.5 Engineered Data Table
```{r}
feature_data_N15 <- pine_rec_N15 %>% 
  prep() %>% 
  bake(new_data = pine_train)

pine_test_processed_N15 <- pine_rec_N15 %>% 
  prep() %>% 
  bake(new_data = pine_test)

datatable(feature_data_N15)

```

Models & Workflows
=======

Column {.tabset .tabset_fade}
----

### Lasso Model
```{r, echo=TRUE}
# LASSO ----
# Will get rid of variables (can be zero for slope)

# * BOOTSTRAP AND LAMBDA ----
#for re-sampling and tuning the penalty parameter 
set.seed(456)
pine_boot <- bootstraps(pine_train)

# create a grid of tuning parameters
lambda_grid <- grid_regular(penalty(), levels = 50)
#creating a grid of values from 0 to 1 for lambda values

# * CREATE MODEL ----

#model
lasso_mod <-
  linear_reg(mixture = 1, penalty = tune()) %>%  #instead of giving a penalty we are using tune()
  set_engine("glmnet")
```

### Lasso Workflow BA_20th
```{r, echo=TRUE}
# ** WORKFLOW ----
lasso_pine_wflow_BA20 <- 
  workflow() %>% 
  add_model(lasso_mod) %>% 
  add_recipe(pine_rec_BA20)

# ** BUILDING LAMBDA AND RMSE ----
set.seed(789)
lasso_grid_BA20 <- tune_grid( #building our sampling grid for what we indicated for tune() in this case penalty
  lasso_pine_wflow_BA20,
  resamples = pine_boot, 
  grid = lambda_grid # how many values to try for the tune() value we indicated (penalty)
)

# finalize our workflow, i.e. update it with this value.
se_rmse_BA20 <- lasso_grid_BA20 %>%
  select_by_one_std_err("rmse") #where from our grid do we want to select for best RMSE

# ** FINALIZE MODEL----
# update our final model with rmse
final_lasso_BA20 <- finalize_workflow(
  lasso_pine_wflow_BA20,
  se_rmse_BA20
)

final_lasso_fit_BA20 <- fit(final_lasso_BA20, data = pine_train)
```

### Lasso Workflow Neigh_1/2th
```{r, echo=TRUE}
# ** WORKFLOW ----
lasso_pine_wflow_N05 <- 
  workflow() %>% 
  add_model(lasso_mod) %>% 
  add_recipe(pine_rec_N05)

# ** BUILDING LAMBDA AND RMSE ----
set.seed(789)
lasso_grid_N05 <- tune_grid( #building our sampling grid for what we indicated for tune() in this case penalty
  lasso_pine_wflow_N05,
  resamples = pine_boot, 
  grid = lambda_grid # how many values to try for the tune() value we indicated (penalty)
)

# finalize our workflow, i.e. update it with rmse value.
se_rmse_N05 <- lasso_grid_N05 %>%
  select_by_one_std_err("rmse") #where from our grid do we want to select for best RMSE

# ** FINALIZE MODEL----
# update our final model with rmse
final_lasso_N05 <- finalize_workflow(
  lasso_pine_wflow_N05,
  se_rmse_N05
)

final_lasso_fit_N05 <- fit(final_lasso_N05, data = pine_train)
```

### Lasso Workflow Neigh_1.5
```{r, echo=TRUE}
# ** WORKFLOW ----
lasso_pine_wflow_N15 <- 
  workflow() %>% 
  add_model(lasso_mod) %>% 
  add_recipe(pine_rec_N15)

# ** BUILDING LAMBDA AND RMSE ----
set.seed(789)
lasso_grid_N15 <- tune_grid( #building our sampling grid for what we indicated for tune() in this case penalty
  lasso_pine_wflow_N15,
  resamples = pine_boot, 
  grid = lambda_grid # how many values to try for the tune() value we indicated (penalty)
)

# finalize our workflow, i.e. update it with rmse value.
se_rmse_N15 <- lasso_grid_N15 %>%
  select_by_one_std_err("rmse") #where from our grid do we want to select for best RMSE

# ** FINALIZE MODEL----
# update our final model with rmse
final_lasso_N15 <- finalize_workflow(
  lasso_pine_wflow_N15,
  se_rmse_N15
)

final_lasso_fit_N15 <- fit(final_lasso_N15, data = pine_train)

```

Column {.tabset .tabset_fade}
----

### Ridge Model
```{r,echo=TRUE}
# RIDGE REGRESSION ---- 
# keeps all but makes them very close to zero

# * CREATE MODEL ----

ridge_mod <- linear_reg(mode = "regression",
                        penalty = 0.1629751,
                        mixture = 0) %>% 
  set_engine("glmnet")
```

### Ridge Workflow BA_20th
```{r, echo=TRUE}

# ** WORKFLOW ----

pine_wflow_glm_BA20 <- workflow() %>% 
  add_recipe(pine_rec_BA20) %>% 
  add_model(ridge_mod) 
  
ridge_pine_fit_BA20 <- fit(pine_wflow_glm_BA20, data = pine_train)

```

### Ridge Workflow Neigh_1/2th
```{r, echo=TRUE}

# ** WORKFLOW ----

pine_wflow_glm_N05 <- workflow() %>% 
  add_recipe(pine_rec_N05) %>% 
  add_model(ridge_mod) 

ridge_pine_fit_N05 <- fit(pine_wflow_glm_N05, data = pine_train)

```

### Ridge Workflow Neigh_1.5
```{r, echo=TRUE}

# ** WORKFLOW ----

pine_wflow_glm_N15 <- workflow() %>% 
  add_recipe(pine_rec_N15) %>% 
  add_model(ridge_mod) 

ridge_pine_fit_N15 <- fit(pine_wflow_glm_N15, data = pine_train)

```

Lasso BA_20th
=======
```{r, echo=FALSE}

# ** PREDICT ----

lasso_prediction_test_BA20 <- predict(final_lasso_fit_BA20, new_data = pine_test)

lasso_results_tbl_BA20 <- bind_cols(lasso_prediction_test_BA20, pine_test_processed_BA20 %>% select(DeadDist) )

```


Column 
-----

### VIF
```{r}
# variable importance plot
final_lasso_BA20 %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = se_rmse_BA20$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)

```

Column 
-------

### Metrics {data-height=240}
```{r}

# Metrics

lasso_pine_metric_BA20 <- last_fit(
  final_lasso_BA20,
  pine_split
) %>%
  collect_metrics() 

table_BA20<-lasso_pine_metric_BA20 %>% kbl() %>% kable_styling()

table_BA20

```

### Adjusted Actual vs Adjusted Predicted DeadDist {data-height=780}
```{r}

ggplot(lasso_results_tbl_BA20, aes(x = DeadDist, y = .pred)) + 
  geom_point(alpha = 0.5) + 
  stat_smooth(method="lm",se=F) +
  labs(y = "Adjusted Predicted Dead Distance", x = "Adjusted Actual Dead Distance")

```

Lasso Neigh_1/2th
=======
```{r, echo=FALSE}

# ** PREDICT ----

lasso_prediction_test_N05 <- predict(final_lasso_fit_N05, new_data = pine_test)

lasso_results_tbl_N05 <- bind_cols(lasso_prediction_test_N05, pine_test_processed_N05 %>% select(DeadDist) )

```


Column 
-----

### VIF 
```{r}

final_lasso_N05 %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = se_rmse_N05$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)

```

Column 
-------

### Metrics {data-height=240}
```{r}

# Metrics

lasso_pine_metric_N05 <- last_fit(
  final_lasso_N05,
  pine_split
) %>%
  collect_metrics() 

table_N05 <- lasso_pine_metric_N05 %>% kbl() %>% kable_styling()

table_N05

```

### Adjusted Actual vs Adjusted Predicted DeadDist {data-height=780}
```{r}

ggplot(lasso_results_tbl_N05, aes(x = DeadDist, y = .pred)) + 
  geom_point(alpha = 0.5) + 
  stat_smooth(method="lm",se=F) +
  labs(y = "Adjusted Predicted Dead Distance", x = "Adjusted Actual Dead Distance")


```


Lasso Neigh_1.5
=======
```{r, echo=FALSE}

# ** PREDICT ----

lasso_prediction_test_N15 <- predict(final_lasso_fit_N15, new_data = pine_test)

lasso_results_tbl_N15 <- bind_cols(lasso_prediction_test_N15, pine_test_processed_N15 %>% select(DeadDist) )

```


Column 
-----

### VIF 
```{r}

final_lasso_N15 %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = se_rmse_N15$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)

```

Column 
-------

### Metrics {data-height = 240}
```{r}

lasso_pine_metric_N15 <- last_fit(
  final_lasso_N15,
  pine_split
) %>%
  collect_metrics() 

table_N15 <- lasso_pine_metric_N15 %>% kbl() %>% kable_styling()

table_N15

```

### Adjusted Actual vs Adjusted Predicted DeadDist {data-height=780}
```{r}

ggplot(lasso_results_tbl_N15, aes(x = DeadDist, y = .pred)) + 
  geom_point(alpha = 0.5) + 
  stat_smooth(method="lm",se=F) +
  labs(y = "Adjusted Predicted Dead Distance", x = "Adjusted Actual Dead Distance")

```


Ridge BA_20th
=======
```{r}

# ** PREDICT ----

prediction_test_BA20 <- predict(ridge_pine_fit_BA20, new_data = pine_test)

ridge_results_tbl_BA20 <- bind_cols(prediction_test_BA20, pine_test_processed_BA20 %>% select(DeadDist) )

```

Column
----

### Metrics {data-height = 240}
```{r}

ridge_pine_metric_BA20 <- last_fit(
  pine_wflow_glm_BA20,
  pine_split
) %>%
  collect_metrics()

ridge_table_BA20 <- ridge_pine_metric_BA20 %>% kbl() %>% kable_styling()

ridge_table_BA20

```

### Adjusted Actual vs Adjusted Predicted DeadDist {data-height=780}
```{r}

ggplot(ridge_results_tbl_BA20, aes(x = DeadDist, y = .pred)) + 
  geom_point(alpha = 0.5) + 
  stat_smooth(method="lm",se=F) +
  labs(y = "Adjusted Predicted Dead Distance", x = "Adjusted Actual Dead Distance")

```

Column 
----

### VIP
```{r}

ridge_pine_fit_BA20 %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = 0.1629751) %>%
  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)

```


Ridge Neigh_1/2th
=======
```{r}

# ** PREDICT ----

prediction_test_N05 <- predict(ridge_pine_fit_N05, new_data = pine_test)

ridge_results_tbl_N05 <- bind_cols(prediction_test_N05, pine_test_processed_N05 %>% select(DeadDist) )


```

Column
----

### Metrics {data-height = 240}
```{r}

ridge_pine_metric_N05 <- last_fit(
  pine_wflow_glm_N05,
  pine_split
) %>%
  collect_metrics()

ridge_table_N05 <- ridge_pine_metric_N05 %>% kbl() %>% kable_styling()

ridge_table_N05

```

### Adjusted Actual vs Adjusted Predicted DeadDist {data-height=780}
```{r}

ggplot(ridge_results_tbl_N05, aes(x = DeadDist, y = .pred)) + 
  geom_point(alpha = 0.5) + 
  stat_smooth(method="lm",se=F) +
  labs(y = "Adjusted Predicted Dead Distance", x = "Adjusted Actual Dead Distance")


```

Column 
----

### VIP
```{r}

ridge_pine_fit_BA20 %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = 0.1629751) %>%
  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)

```

Ridge Neigh_1.5
=======
```{r}

# ** PREDICT ----

prediction_test_N15 <- predict(ridge_pine_fit_N15, new_data = pine_test)

ridge_results_tbl_N15 <- bind_cols(prediction_test_N15, pine_test_processed_N15 %>% select(DeadDist) )

```

Column
----

### Metrics {data-height=240}
```{r}

ridge_pine_metric_N15 <- last_fit(
  pine_wflow_glm_N15,
  pine_split
) %>%
  collect_metrics()

ridge_table_N15 <- ridge_pine_metric_N15 %>% kbl() %>% kable_styling()

ridge_table_N15

```

### Adjusted Actual vs Adjusted Predicted DeadDist {data-height=780}
```{r}

ggplot(ridge_results_tbl_N15, aes(x = DeadDist, y = .pred)) + 
  geom_point(alpha = 0.5) + 
  stat_smooth(method="lm",se=F) +
  labs(y = "Adjusted Predicted Dead Distance", x = "Adjusted Actual Dead Distance")


```

Column 
----

### VIP
```{r}

ridge_pine_fit_N15 %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = 0.1629751) %>%
  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)

```