---
title: "trees"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    vertical_layout: scroll
    source_code: embed
    theme: flatly
---

```{r cars, include=FALSE}
library(tidymodels)
library(vip)
library(tidyverse)
library(readxl)
library(ggplot2)
library(dplyr)
library(hrbrthemes)
library(GGally)
library(tidyr)
# Analysis
library(corrr)
library(flexdashboard)
library(DT)
library(rpivotTable)
library(ggcorrplot)
library(recipes)
library(parsnip)
library(performance)
library(vip)
#load data


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

```

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

### Neighboorhood Infection

```{r}
pine_beetles_dat %>% dplyr::select(IND_BA_Infest_20th:IND_BA_Infest_1.5, DeadDist, Response) %>% #get dsm scores in long format
  tidyr::pivot_longer(IND_BA_Infest_20th:IND_BA_Infest_1,
                      names_to="neighboorhood_tree_radius",
                      values_to = "infection_response") %>%
  group_by(neighboorhood_tree_radius, infection_response, Response)  %>%  
  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(Response)))+ geom_bar(position="dodge", stat="identity")+ 
  facet_wrap(infection_response~., scales="free_y", ncol=1)+ labs(x= "Neighboorhood Size (acres)", fill= "Response") + 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()
  

```


```{r, include=FALSE}
################LOGISTIC REGRESSSION###########################
#Ind_DeadDist ~ BA_Infest_1/4th + Ind_DeadDist+ tree_diam
#(top 10% is infected change)

pine_beetles_dat_logistic<-pine_beetles_dat %>%mutate(deaddist_quantile= factor(quantile(DeadDist, 0.7)<DeadDist)) 
pine_beetles_dat_logistic<- pine_beetles_dat_logistic %>% mutate(deaddist_quantile = fct_recode(deaddist_quantile,"1" = "TRUE", "0"= "FALSE")) 
pine_beetles_dat_logistic$deaddist_quantile<- as.factor(pine_beetles_dat_logistic$deaddist_quantile)
pine_beetles_dat_logistic<- pine_beetles_dat_logistic %>% dplyr::select(-c(DeadDist, TreeNum))
splits <- initial_split(pine_beetles_dat_logistic)
train <- training(splits)
test <- testing(splits)
set.seed(123)


#tidy model framework

# MODELING ----
model_fit_glm <- logistic_reg() %>% #specify type of model you want
  set_engine("glm") %>% #set engine, what technique do we want to estimate the parameters
  fit(deaddist_quantile ~ Infest_Serv1 + TreeDiam  , data = training(splits)) # predict the year and use the training data (grab the training data from splits)

model_fit_glm

model_fit_glm_2<- logistic_reg() %>% 
  set_engine("glm")
# * Create Recipe ----
## specify variable relationships
## specify (training) data
## feature engineer
## process recipe on data
pine_rec <- pine_beetles_dat_logistic %>% 
  recipe(deaddist_quantile~ TreeDiam + Infest_Serv1 + BA_Infest_1.5  ) %>% 
  step_sqrt(.) %>% 
  step_corr(.) %>% 
  prep()

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

# * Create Model GLM
glm_mod <- 
  logistic_reg() %>% 
  set_engine("glm")


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

pine_wflow

#fit model
pine_fit <- 
  pine_wflow %>% 
  fit(data = train)


pine_fit %>% 
  extract_fit_parsnip() %>% 
  tidy()

pine_fit %>% 
  extract_fit_parsnip() %>% 
  glance()


pine_fit %>% 
  extract_preprocessor()

# PREDICTION ----
#doesnt work
prediction_class_test <- predict(pine_fit, new_data = testing(splits), type = "class") #tells you what the predicted class is

prediction_prob_test  <- predict(pine_fit, new_data = testing(splits), type = "prob") #tells you the probability of predicting each class

results_tbl <- bind_cols( #Combine two predictions and the testing data
  prediction_class_test,
  prediction_prob_test,
  testing(splits)
)
```

### Logistic regression Results--fix size

```{r}
 pine_fit %>% 
  extract_fit_parsnip() %>% 
check_model()
  
#p[[3]] <- p[[3]] + labs(y = "Variance Inflation\nFactor (VIF)")
#p[[6]] <- p[[6]] + xlim(c(-5, 6))
#https://github.com/easystats/performance/issues/318
```



```{r pressure, echo=FALSE, include=FALSE}
#graph of killed trees over dead distance
pine_beetles_dat$Response<- as.factor(pine_beetles_dat$Response)
pine_beetles_dat  %>% 
  ggplot( aes(x=DeadDist, fill=Response)) +
  geom_histogram( color="#e9ecef", alpha=0.6, position = 'identity') +labs(x= "Dead Distance", fill= "Response")+
  scale_fill_manual(values=c("#69b3a2", "#404080")) +
  theme_ipsum() +
  labs(fill="")
```


Row {.tabset .tabset-fade}
-------------------------------------------------------------------------------------------------------
### Table

```{r}
#pine_beetles_cor <- pine_beetles_dat %>% 
  #correlate() 

#pine_beetles_cor %>% 
  #rearrange() %>% 
  #fashion() %>% 
 # knitr::kable()

DT::datatable(pine_beetles_dat,options = list(
  pageLength=50, scrollX='400px'), filter = 'top')
```

-------------------------------------------------------------------------------------
### Descriptive Table

```{r}
rpivotTable(
             pine_beetles_dat,
             aggregatorName = "Integer Sum",
             vals="Total",
             cols = "Test",
             width="100%",
             height="400px",
             rendererName = "Table"

 )
```

### Need to fix this

```{r, include=FALSE}
#Bins of infested trees within 1.5 acre
breaks <- c(0,10,20,30, 40, 50, 60, 70, 80, 90, 110)
#create buckets
tags <- c("[0-10)","[10-20)", "[20-30)", "[30-40)", "[40-50)", "[50-60)","[60-70)", "[70-80)","[80-90)", "[90-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$BA_Infest_1.5, 
                  breaks=breaks, 
                  include.lowest=TRUE, 
                  right=FALSE, 
                  labels=tags) 

summary(group_tags)

factor(group_tags, 
       levels = tags,
       ordered = TRUE)
```


```{r}
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='Infested trees per 1.5 acre' ) +
  theme_minimal()

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

# corr plot 

r <-  cor(pine_beetles_dat[, 6:26], use="complete.obs")
round(r,2)
ggcorrplot(r)
```


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

#RIDGE REGRESSION MODEL
# 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") #mixture to 0 means using ridge regression, mixture to 1 is lasso method

#glmnet wants x and y in matrix format

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


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


# verify Ridge Regression performance with standard linear regression approach
lm(sqrt(DeadDist) ~ TreeDiam + Infest_Serv1 + BA_20th, data = pine_beetles_dat) %>% 
  glance() 


#(pine_ridge_fit) %>% pull_workflow_fit() %>%
  #tidy() %>%
 # filter(row != "(Intercept)") %>%
 # top_n(2, wt = abs(value)) %>%
#  ggplot(aes(value, reorder(row, value))) +
#  geom_point() +
 # ggtitle("Top 25 influential variables") +
#  xlab("Coefficient") +
#  ylab(NULL)

```



```{r, include=FALSE}



```


```{r}

```