---
title: "Lake Tahoe Basin Pine Beetle Outbreak"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    vertical_layout: scroll
    source_code: embed
 
---

```{r}
library(thematic)
thematic::thematic_rmd(font = "auto")
#https://pkgs.rstudio.com/flexdashboard/articles/theme.html
```


```{r setup, include=FALSE}
library(tidymodels)
library(vip)
library(tidyverse)
library(readxl)
library(ggplot2)
library(dplyr)
library(hrbrthemes)
library(GGally)
library(tidyr)
library(generics)
library(workflows)
library(plotly)
library(dotwhisker)
library(glmnet)
library(broom)
library(kableExtra)
library(gt)


# Analysis
library(corrr)
library(flexdashboard)
library(DT)
library(rpivotTable)
library(ggcorrplot)
library(recipes)
library(parsnip)
library(performance)
library(vip)
library(MASS)
#load data

setwd("~/Downloads/data_science_2/project_1")
#overview of data
pine_beetles_dat <-readxl::read_excel("~/Downloads/data_science_2/project_1/Data_1993.xlsx")
```

Overview
===================================================================================================

Row 
----------------------------------------------------------------------------------------------------------------


```{r, include=FALSE}
#bin the severity (3 groups)

#Bins of infested trees within 1.5 acre
breaks <- c(0,20, 40, 60, 110)
#create buckets
tags <- c("[0-20)","[20-40)", "[40-60)", "[60-110)")
#complete cases only
pine_beetles_dat_complete=pine_beetles_dat[complete.cases(pine_beetles_dat), ]

# bucketing values into bins
group_tags <- cut(pine_beetles_dat$Infest_Serv1, 
                  breaks=breaks, 
                  include.lowest=TRUE, 
                  right=FALSE, 
                  labels=tags) 

summary(group_tags)

factor(group_tags, 
       levels = tags,
       ordered = TRUE)
#https://www.jdatalab.com/data_science_and_data_mining/2017/01/30/data-binning-plot.html

r <-  cor(pine_beetles_dat[, c(7:13, 17,18, 21, 22)], use="complete.obs")
round(r,2)
```

### Distribution of infestation
```{r}
#map with severity 

inf_sev<- ggplot(pine_beetles_dat, aes(Northing, Easting)) + geom_point(aes(color = Infest_Serv1)) +  scale_colour_gradientn(colours =  c("#049D73", "white", "#e59f00"))+ coord_equal() + labs(x = "Easting (m)", y = "Northing (m)", fill="Number of Trees Infested")
inf_sev+ labs(
         colour = "Infestation Severity"
        )
```

### Correlation
```{r}
ggcorrplot(r, ggtheme = ggplot2::theme_gray,
  colors = c("#049D73", "white", "#e59f00")
)
```


```{r, include=FALSE}
ggplot(data = as_tibble(group_tags), mapping = aes(x=value)) + 
  geom_bar(fill="bisque",color="white",alpha=0.7) + 
  stat_count(geom="text", aes(label=sprintf("%.4f",..count../length(group_tags))), vjust=-0.5) +
  labs(x='Infestation Severity \n to nearest response tree') +
  theme_minimal() 
```


```{r, include=FALSE}
#Ind_DeadDist=indicator if nearest brood tree is within 50m effective distance
pine_beetles_dat %>% dplyr::select('IND_BA_Infest_1/4th':IND_BA_Infest_1.5, DeadDist, Response, Ind_DeadDist) %>% #get dsm scores in long format
  tidyr::pivot_longer('IND_BA_Infest_1/4th':IND_BA_Infest_1.5,
                      names_to="neighboorhood_tree_radius",
                      values_to = "infection_response") %>%
  group_by(neighboorhood_tree_radius, infection_response, Ind_DeadDist)  %>%  
  summarise(N=n())  %>% #sum of n to find out the occurences of infection at each infected tree radius
mutate(neighboorhood_tree_radius= factor(neighboorhood_tree_radius),
       N= as.numeric(N),
       infection_response= recode_factor(infection_response, "0"= "No tree in neighboorhood infected",
                            "1"= "Tree has been infected in neighboorhood"))  %>%
  filter(infection_response != "NA") %>% ungroup()    %>%
  ggplot(aes(x=fct_reorder(neighboorhood_tree_radius, N), y=N, fill=factor(Ind_DeadDist)))+ geom_bar(position="dodge", stat="identity")+ 
  facet_wrap(infection_response~., scales="free_y", ncol=1)+ labs(x= "Neighboorhood Size (acres)", fill= "Nearest brood tree within 50m") + scale_fill_manual(values=c("#69b3a2", "#404080")) +
  scale_x_discrete(labels=c("IND_BA_Infest_20th"="All sizes", "IND_BA_Infest_1/4th"= ".25", "IND_BA_Infest_1/2th"=".5", "IND_BA_Infest_1"= "1", "IND_BA_Infest_1.5"="1.5")) + theme_ipsum()


```



Row {.tabset .tabset-fade, data-height=400}
-------------------------------------------------------------------------------------------------------
### Table

```{r pressure, echo=FALSE}
DT::datatable(pine_beetles_dat,options = list(
  pageLength=50, scrollX='400px'), filter = 'top')


```



Lasso Regression
==========================================================================================================

Column {.sidebar data-width=500}
-------------------------------------
### LASSO
#### Least Absolute Shrinkage and Selection Operation

- Lasso regression can be used when multicolinearity is present to select and shrink the parameters. It keeps all the predictor variables in the model but constrains or regularizes the regression coefficients by shrinking them toward zero.
- In the first run of lasso, I decided to identify the lambda value that produced the lowest mean squared error. This produced a lower r sqaured value compared to the pre-established lambda value.

- In the second run of lasso, I decided to pick the best model based off the lowest RMSE, however you can see from the penalty and RMSE plot that the best model was only able to explain 16% of the variance.

- In the final lasso results, we can see that as the dead distance goes down, the infestation severity goes up.


```{r}
knitr::include_graphics("lake-tahoe.png")

```



Row {data-height=350, .tabset}
-------------------------------------
### Test MSE by lambda value
#### Identify the lambda value that produces the lowest test mean squared error (MSE)
```{r, include=FALSE}
#https://tahoequarterly.com/environment/how-healthy-are-tahoes-forests
#https://bookdown.org/tpinto_home/Regularisation/lasso-regression.html
#Lasso

#we need to determine the best lambda value so we will perform k folds cross validation
#identify the lambda value that produces the lowest test mean squared error (MSE)
y <- pine_beetles_dat$DeadDist

#define matrix of predictor variables
x <- data.matrix(pine_beetles_dat[, c('TreeDiam', 'Infest_Serv1', 'SDI_20th', 'BA_20th')])

#perform k-fold cross-validation to find optimal lambda value
cv_model <- cv.glmnet(x, y, alpha = 1)
cv_model
#minimum mse is achieved when lambda 0.0191.

#find optimal lambda value that minimizes test MSE
best_lambda <- cv_model$lambda.min
best_lambda

best_model <- glmnet(x, y, alpha = 1, lambda = best_lambda)
coef(best_model)


#define new observation
new = matrix(c(24, 2.5, 3.5, 18.5), nrow=1, ncol=4) 

#use lasso regression model to predict response value
predict(best_model, s = best_lambda, newx = new)

#use fitted best model to make predictions
y_predicted <- predict(best_model, s = best_lambda, newx = x)

#find SST and SSE
sst <- sum((y - mean(y))^2)
sse <- sum((y_predicted - y)^2)

#find R-Squared
rsq <- 1 - sse/sst
rsq
#variation in the response values of the training data is 13%

```


```{r, fig.width=5,fig.height=3}

plot(cv_model) 

```



### Impact of lambdas on est. coefficients
```{r}
plot(cv_model$glmnet.fit, 
     "lambda", label=FALSE)
```


```{r, include=FALSE}
set.seed(1234)
# Create training/testing data
pine_split <- initial_split(pine_beetles_dat)
pine_train <- training(pine_split)
pine_test <- testing(pine_split)

#create bootstrap samples for resampling and tuning the penalty parameters
pine_boot <- bootstraps(pine_train)

# create a new recipe; could use `add_step()` 
pine_rec <- pine_train %>% 
  recipe(DeadDist ~ TreeDiam + Infest_Serv1 + SDI_20th + BA_20th) %>% 
  step_log(all_outcomes()) %>% 
  step_corr(all_predictors()) %>% 
  step_normalize(all_numeric(), -all_outcomes()) %>% 
  step_zv(all_numeric(), -all_outcomes()) %>% prep()
#prep function estimates the required means and standard deviations from the data in the training argument. The transformations specified by each step are also sequentially executed on the data set

#
#take old data and simulate new column and test recipe on it
pine_beetles_dat_sim<- pine_beetles_dat %>% mutate(Neigh_2=rnorm(9687))
test_ex <- bake(pine_rec, new_data = pine_beetles_dat_sim)
names(test_ex) %>% head()


# create a grid of tuning parameters
lambda_grid <- grid_regular(penalty(), levels = 50)

lasso_mod <-
  linear_reg(mixture = 1, penalty = tune()) %>% 
  set_engine("glmnet")

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

# let's look at bootstrap results
lasso_grid %>%
  collect_metrics()
```



```{r, include=FALSE}
# We have a couple of options for choosing our final parameter, such as 
# select_by_pct_loss() or select_by_one_std_err(), but for now let’s stick
# with just picking the lowest RMSE. After we have that parameter, we can 
# finalize our workflow, i.e. update it with this value.
# see https://juliasilge.com/blog/lasso-the-office/
#https://www.statology.org/lasso-regression-in-r/

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

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



```



Row {data-height=650, .tabset}
-------------------------------------


### Penalty and RSME
```{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")
```

### Variable Importance Plot
```{r}
final_lasso %>%
  fit(pine_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = lowest_rmse$penalty) %>% mutate(
    Importance = abs(Importance),
    Variable = fct_reorder(Variable, Importance)
  ) %>%
 dplyr::mutate(Variable = case_when(Variable =='Infest_Serv1' ~ "Infestation Severity",
                                     Variable == 'BA_20th' ~ "Basal area 1/20th acre",
                                Variable == 'TreeDiam' ~ "Tree Diameter"
                                )) %>%
  ggplot(aes(x = Importance, y = Variable, fill = Sign)) +
  geom_col() +
  scale_x_continuous(expand = c(0, 0)) +
  labs(y = NULL) 

```

### Final Fit
#### Note that penalty (lambda) is close to zero, so it is near equivalent to lm() solution
```{r, echo=TRUE, out.height="80%"}
final_lasso %>% 
  fit(pine_train) %>%
  pull_workflow_fit() %>% 
  tidy()
 
```   

Linear Regression
==========================================================================================


Column {.sidebar data-width=500}
-------------------------------------
### Overview
- Goal is to minimize the residuals or the differences between the actual and fitted response.
- So, we want to minimize the vertical distances to the data points and the slope.

- This model only explains 16% of the variance and the high negative log likelihood indicates this model is not a good fit.

$$
y_i = \beta_0 + \beta_1x_i + \epsilon_i, ~~\mbox{var}(\epsilon_i) = \sigma^2 
$$
```{r, include=FALSE}
# Tidymodels Version of lm() ----
#  Specify, fit and evaluate same models with tidymodels


set.seed(1234)
pine_split <- initial_split(pine_beetles_dat)
pine_train_linear <- training(pine_split)
pine_test_linear <- testing(pine_split)

# * Create Recipe ----
## specify variable relationships
## specify (training) data
## feature engineer
## process recipe on data
pine_rec <- pine_train_linear %>% 
  recipe(DeadDist ~ TreeDiam + Infest_Serv1 +  SDI_20th + BA_20th) %>% 
  step_sqrt(all_outcomes(), skip = TRUE)  %>% step_corr(all_outcomes(), skip = TRUE) %>%
 prep()


# View feature engineered data
pine_rec %>% 
  prep() %>% 
  juice()


# * Create Model ----
lm_mod <- 
  linear_reg() %>% 
  set_engine("lm")


# * Create Workflow ----
pine_wflow <- 
  workflow() %>% 
  add_model(lm_mod) %>% 
  add_recipe(pine_rec)

pine_wflow


pine_fit_linear <- 
  pine_wflow %>% 
  fit(data =pine_train_linear)


linear.fit<- pine_fit_linear %>% 
  extract_fit_parsnip() %>% 
  tidy()

pine_fit_linear %>% 
  extract_fit_parsnip() %>% 
  glance()

pine_fit_linear %>% 
  extract_preprocessor()


### use testing set to predict
pine_predictions <- predict(pine_fit_linear,
                            new_data = pine_test_linear)

pine_test_results <- pine_test_linear %>% 
  dplyr::select(DeadDist,TreeDiam, Infest_Serv1, SDI_20th,BA_20th) %>% 
  bind_cols(pine_predictions)

ggplot(pine_test_results, aes(x = DeadDist, y = .pred)) +
  geom_point(alpha = 0.5) + 
  geom_abline(color = 'blue', linetype = 2) +
  coord_obs_pred() +
  labs(x = 'Actual Dead Distance', y = 'Predicted Dead Distance')

#https://www.tidymodels.org/start/models/

```


Row {data-width=600, .tabset}
-------------------------------------
### Model Estimates
```{r}
tidy(pine_fit_linear) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline = geom_vline(xintercept = 0, colour = "grey50", linetype = 2)) + xlab("Estimate") 


```

### Linearity
#### Most fitted values are linear, except for larger fitted values.
```{r}
###REVIEW FIT
check_mod<- pine_fit_linear %>% 
  extract_fit_parsnip() %>% 
  check_model()

pine_fit_linear %>% 
  extract_fit_parsnip() %>% check_model(check = "linearity")

```


### QQ Plot
#### Residuals are mostly normally distributed.
```{r}
pine_fit_linear %>% 
  extract_fit_parsnip() %>% check_model(check = "qq")

```

### Homogenity of Variance
#### There is no homogenity in variance.
```{r}
pine_fit_linear %>% 
  extract_fit_parsnip() %>% check_model(check = "homogeneity")

```


### Model Results
#### Strong negative log likelihood shows that this model is not a good fit.

```{r}

linear.fit %>% gt::gt()


pine_fit_linear %>% 
  extract_fit_parsnip() %>% 
  glance() %>% gt::gt() %>% tab_options(container.width=700)
```

### Prediction Accuracy
```{r, include=FALSE}
linear_mod <-   linear_reg() %>% 
  set_engine("lm") %>% set_mode("regression")

# Train linear_model with last_fit()
linear_fit <- linear_mod %>% 
  last_fit(DeadDist ~ ., split = pine_split)

# Collect predictions and view results
predictions_df <- linear_fit %>% collect_predictions()
predictions_df
```


```{r}
ggplot(predictions_df, aes(x = DeadDist, y = .pred)) +
  geom_point(alpha = 0.5) + 
  geom_abline(color = 'blue', linetype = 2) +
  coord_obs_pred() +
  labs(x = 'Actual Dead Distance', y = 'Predicted Dead Distance')
```

```

```{r, include=FALSE}
#regline <-predict(lm_fit, data.frame(SDI_20th=(c(1.8, 3.5))),interval="confidence", data = pine_beetles_dat)

#newpred <-predict(lm_fit,data.frame(SDI_20th=(c(1.8, 3.5))),interval="prediction", data =pine_beetles_dat)


pred.int <-predict(pine_fit_linear, pine_train_linear, interval = "prediction")
#mydata <-cbind(pine_beetles_dat, pred.int)

#p <-ggplot(mydata,aes(BA_20th, SDI_20th))+geom_poin()+stat_smooth(method =lm,linetype =
 #"dashed",
 #size = 0.5)
# 3. Add prediction interval
#p <-p +geom_line(aes(y = lwr),color ="red",linetype = "dashed")+geom_line(aes(y =upr),color ="red",linetype="dashed") +
#theme_bw()
```



```{r, include=FALSE }

#Ind_DeadDist=indicator if nearest brood tree is within 50m effective distance
pine_beetles_dat %>% dplyr::select('IND_BA_Infest_1/4th':IND_BA_Infest_1.5, DeadDist, Response, Ind_DeadDist) %>% #get dsm scores in long format
  tidyr::pivot_longer('IND_BA_Infest_1/4th':IND_BA_Infest_1.5,
                      names_to="neighboorhood_tree_radius",
                      values_to = "infection_response") %>%
  group_by(neighboorhood_tree_radius, infection_response, Ind_DeadDist)  %>%  
  summarise(N=n())  %>% #sum of n to find out the occurences of infection at each infected tree radius
mutate(neighboorhood_tree_radius= factor(neighboorhood_tree_radius),
       N= as.numeric(N),
       infection_response= recode_factor(infection_response, "0"= "No tree in neighboorhood infected",
                            "1"= "Tree has been infected in neighboorhood"))  %>%
  filter(infection_response != "NA") %>% ungroup()    %>%
  ggplot(aes(x=fct_reorder(neighboorhood_tree_radius, N), y=N, fill=factor(Ind_DeadDist)))+ geom_bar(position="dodge", stat="identity")+ 
  facet_wrap(infection_response~., scales="free_y", ncol=1)+ labs(x= "Neighboorhood Size (acres)", fill= "Nearest brood tree within 50m") + scale_fill_manual(values=c("#69b3a2", "#404080")) +
  scale_x_discrete(labels=c("IND_BA_Infest_20th"="All sizes", "IND_BA_Infest_1/4th"= ".25", "IND_BA_Infest_1/2th"=".5", "IND_BA_Infest_1"= "1", "IND_BA_Infest_1.5"="1.5")) + theme_ipsum()



pine_beetles_dat %>% dplyr::select('IND_BA_Infest_1/4th':IND_BA_Infest_1.5, DeadDist, Response, Ind_DeadDist, Infest_Serv1) %>% #get dsm scores in long format
  tidyr::pivot_longer('IND_BA_Infest_1/4th':IND_BA_Infest_1.5,
                      names_to="neighboorhood_tree_radius",
                      values_to = "infection_response")  %>%
ggplot(aes(x=Infest_Serv1, y=DeadDist)) + 
    geom_point(size=6, color="#69b3a2") +
    theme_ipsum()
```


```{r, include=FALSE}

hexagon<- ggplot(pine_beetles_dat, aes(TreeDiam, DeadDist)) +
  geom_hex(bins = 20, color = "white")+
  scale_fill_gradient(low =  "#00AFBB", high = "#FC4E07")+
  theme_minimal() + xlab('Tree Diameter') + ylab("Dead Distance")
```