read in TE std stk, carbon and summmary data. can be found on google drive [local copy used here]

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.1
stdstk <- '/Users/eyackulic/Downloads/SP-TE-StdStk.rds' |> 
  readRDS() |>
  dplyr::filter(
                substr(Species,0,2) == 'QU'#subset for oaks only
                  ) 

carbon <- '/Users/eyackulic/Downloads/SP-TE-carbon.rds' |>
  readRDS() |>
  dplyr::filter(CaseID %in% stdstk$CaseID) 

summary <- '/Users/eyackulic/Downloads/SP-TE-FVSSummary2.rds' |>
  readRDS() |>
  dplyr::filter(CaseID %in% stdstk$CaseID) |>
  dplyr::select(-c(MgmtID, RunTitle,StandID)) |>
  dplyr::right_join(carbon, by = c('CaseID','Year'))

Identify what stands are oak woodlands pre-FVS prescription

oak_woodland_stands <-
  summary |>
  dplyr::filter(Year < 2029) |>
  dplyr::filter(ForTyp %in% c(401, 403, 404, 405, 501, 502:506,
                            509, 510, 514:519, 601:605, 922:941)) |>
  dplyr::select(c('StandID','MgmtID','CaseID'))

filter the std stk table for ids that match the oak woodlands definition and experience a mortality event during 2029:2031 that is greater than any other point in the 20 year period from 2020:2040. This tells us that management occurred.

lite_treatments <- #these treatments have a muted response and need to be filtered differently
  stdstk |>
  dplyr::filter(MgmtID %in% c('MTH2', 'RMTF', 'HERB') &
                  CaseID %in% oak_woodland_stands$CaseID) |>
  dplyr::group_by(MgmtID, StandID) |>
  dplyr::select('StandID','MgmtID','CaseID') |>
  dplyr::distinct()

mort_removal <- 
  stdstk |>
  dplyr::filter(!MgmtID %in% c('MTH2', 'RMTF', 'HERB') &
                 Year %in% c(2020:2040) &
                  CaseID %in% oak_woodland_stands$CaseID) |>
  dplyr::group_by(MgmtID, StandID) |>
  dplyr::filter(MrtTPA == max(MrtTPA)) |>
  dplyr::filter(Year %in% c(2029,2030,2031)) |>
  dplyr::select('StandID','MgmtID','CaseID') |>
  dplyr::distinct() |>
  dplyr::bind_rows(lite_treatments)

calculate responses.

responses <-
summary |>
  dplyr::arrange(StandID, MgmtID,Year) |>
  dplyr::distinct() |>
  dplyr::filter(CaseID %in% mort_removal$CaseID) |>
  dplyr::mutate(
    response = (SDI - dplyr::first(SDI)) / dplyr::first(SDI),
    tpa_response = (Tpa - dplyr::first(Tpa)) / dplyr::first(Tpa)
  ) |>
  dplyr::filter(response != 0) |>
  dplyr::as_tibble() |>
  dplyr::mutate(
    response = dplyr::if_else(response >= 1, NA, response),
    tpa_response = dplyr::if_else(tpa_response >= 1, NA, tpa_response)
  )  

ggplot(data = responses, aes( x= MgmtID, y = response)) + geom_boxplot()
## Warning: Removed 1192 rows containing non-finite values (`stat_boxplot()`).

responses |> dplyr::group_by(MgmtID) |>
  dplyr::select(response)|> rstatix::get_summary_stats() |>
  dplyr::select(MgmtID, median) |>
  dplyr::bind_cols(name = c('Clearcut','Targeted Herbicide', 'Herbivory',
                     'Commercial Thin','Manual Thin' ,'Mechanical Thinning',
                     'Thinning', 'Grapple Pile', 'Mastication', 'Mowing',
                     'Aerial Rx', 'Ground Rx'))
## Adding missing grouping variables: `MgmtID`
## # A tibble: 12 × 3
##    MgmtID median name               
##    <chr>   <dbl> <chr>              
##  1 CMCC   -0.84  Clearcut           
##  2 HCTA   -0.271 Targeted Herbicide 
##  3 HERB   -0.309 Herbivory          
##  4 MRCT    0.269 Commercial Thin    
##  5 MTH2   -0.307 Manual Thin        
##  6 MTH3    0.24  Mechanical Thinning
##  7 MTTH    0.64  Thinning           
##  8 RMGP   -0.036 Grapple Pile       
##  9 RMMA   -0.529 Mastication        
## 10 RMTF   -0.196 Mowing             
## 11 RXAI   -0.458 Aerial Rx          
## 12 RXGF   -0.329 Ground Rx