categorical predictors

Author: Carly Lovas

R/maine/lobster_did.R

@@ -1,7 +1,6 @@
 ## Lobster DiD ####
 ###################
 
-
 # Install & load the libraries ----
 library(tidyverse) 
 library(tidygraph)
@@ -23,24 +22,26 @@ portfolio |>
 shrimpers |>
   filter(lc1 == 1 | lc2 == 1 | lc3 == 1 | lco == 1 | lc2o == 1 | lc3o == 1) -> shrimp_w_lobster
 
-portfolio |>
-  filter(landings_number %in% shrimp_w_lobster$landings_number) |> 
-  select(license_year, landings_number, csc, css, lc1, lco, lc2, lc2o, lc3, lc3o) |>
-  arrange(landings_number, license_year) |>
-  filter(license_year >= 2009 & !license_year == 2015) -> cs_portfolio
-
+# All commercial lobster harvesters ----
 portfolio |>
   filter(lc1 == 1 | lc2 == 1 | lc3 == 1 | lco == 1 | lc2o == 1 | lc3o == 1) |>
-  select(license_year, landings_number, lc1, lco, lc2, lc2o, lc3, lc3o) |> ## we need to compare lobster class changes across all harvesters t
+  select(license_year, landings_number, lc1, lco, lc2, lc2o, lc3, lc3o) |>
   pivot_longer(lc1:lc3o, names_to = "lobster_class", values_to = "count") |>
   filter(!count == 0) |>
   mutate(class = case_when(
     lobster_class %in% c("lco", "lc1")  ~ 1,
-    lobster_class %in% c("lc2", "lc2o") ~ 2, 
+    lobster_class %in% c("lc2", "lc2o") ~ 2,
     lobster_class %in% c("lc3", "lc3o") ~ 3
-  )) |> 
+  )) |>
+  
+  ## Ensure population of lobsterers had license in the year before and after the closure
   filter(license_year %in% c(2014,2016)) |>
-  group_by(landings_number, lobster_class, class) |>
+  
+  ## Remove instances of mulitple license types issued in a year...
+  mutate(n = n(), .by = c(landings_number, license_year)) |>
+  filter(!n > 1L) |>
+  select(!n) |>
+  group_by(landings_number) |>
   summarise(years = n_distinct(license_year)) |>
   filter(years == 2) -> lobsterers
 
@@ -54,24 +55,10 @@ portfolio |>
     lobster_class %in% c("lc2", "lc2o") ~ 2, 
     lobster_class %in% c("lc3", "lc3o") ~ 3
   )) |>
-  filter(license_year %in% c(2014,2016)) |> 
+  select(landings_number, license_year, class) |>
+  filter(license_year %in% c(2014,2016)) |>
   pivot_wider(names_from = license_year, values_from = class) |>
-  mutate(switch = ifelse(`2014` < `2016`, 1, 0)) -> switch # come back to the NAs later, may represent exit
-
-# portfolio |>
-#   select(license_year, landings_number, csc, css, lc1, lco, lc2, lc2o, lc3, lc3o) |>
-#   filter(license_year >= 2009 & !license_year == 2015) |>
-#   left_join(switch |> select(landings_number, switch)) |>
-#   drop_na() |>
-#   mutate(shrimp_w_lobster = ifelse(landings_number %in% shrimp_w_lobster$landings_number, 1, 0),
-#          closure = ifelse(license_year >= 2016, 1, 0),
-#          did = shrimp_w_lobster*closure) |>
-#   select(license_year, landings_number, switch, shrimp_w_lobster, closure, did) -> lob_did
-
-# lm(shrimp_w_lobster+closure+did, data = lob_did) 
-
-### Yeaaaaaah I don't think DiD is the right test here...
-
+  mutate(switch = ifelse(`2014` < `2016`, 1, 0)) -> switch 
 
 portfolio |>
   select(license_year, landings_number, lc1, lco, lc2, lc2o, lc3, lc3o) |> 
@@ -88,10 +75,14 @@ portfolio |>
          closure = ifelse(license_year >= 2016, 1, 0),
          did = shrimp_w_lobster*closure) -> lob_did
 
-lm(class ~ shrimp_w_lobster+closure+did, data = lob_did)
 
-summary(lm(license_year ~ class+shrimp_w_lobster+closure+did, data = lob_did))
+# Categorical dependent variable ----
+## Did the lobster class type switch after the closure of the shrimp fishery? 
+summary(MASS::polr(factor(class) ~ shrimp_w_lobster+closure+did, data = lob_did))
+
+
+# Category 1 --> 2 ----
+
+
 
-## Better but still not right...
-## Maybe the regression should be proportion of license classes ~ shrimp-lobster coholders and license switch...?