To BruceET's point, you want to do everything you can to show that the intervention is plausibly un-confounded. In my opinion, difference-in-differences (DiD) can work in this setting. See here for a gentle introduction to the DiD methodology. I will address some caveats with your study later.
One group has always had a "treatment" throughout the time period (always had a specific legal prohibition). One group never had the "treatment" within the time period (never had the prohibition). One group had the "treatment", but each member of the group had the "treatment" at a different time within the period (individual states enacted prohibition in different years).
First, we should get your empirical model correct. I should note, however, that you cannot proceed with the "classical" DiD approach. You do not have well-defined pre-/post-treatment periods. Treatment only affects a subset of U.S. states, which gives you a group of "non-adopter" states to serve as a control group. Treated state jurisdictions enter into treatment at different times, a subset of which have only one treatment history (i.e., the "always" treated). In staggered adoption settings, you are outside the realm of the "classical" approach. Instead, you must use with the "generalized" approach, which is also a two-way fixed effects estimator. If I understand your question, you wish to estimate the following:
$$
\text{Outcome}_{st} = \gamma_{s} + \lambda_{t} + \delta \text{Prohibition}_{st} + \epsilon_{st},
$$
where you observe a continuous outcome in state $s$ in year $t$. The parameters $\gamma_{s}$ and $\lambda_{t}$ denote fixed effects for states and years, respectively. The variable $\text{Prohibition}_{st}$ is your policy (i.e., treatment) dummy.
Quoting from the comments:
So I test a model like "outcome ~ year + policy + year x policy" where the "year x policy" term is the DiD term and "policy" is a binomial variable indicating each year whether or not the specific state has the specific policy?
Not quite right. You need to include state and year effects. In addition, your interaction should instantiate the policy variable appropriately. I recommend creating the policy variable manually. It is still your interaction term, you just need to define it in a different way as your "post-treatment" periods vary across states. Again, your policy variable, $\text{Prohibition}_{st}$, is a dummy equal to 1 in precisely the state-adoption years, 0 otherwise. The "never treated" states will be "always 0" for the entire observation period. Similarly, the "always treated" states will be "always 1" in all state-years. In R, your model should look something like this:
model <- lm(outcome ~ as.factor(state) + as.factor(year) + prohibition + ..., data = ...)
In general, the policy dummy is allowed to take any pattern, which is something noted in Jeffrey Wooldridge's text Introductory Econometrics: A Modern Approach (7th Edition). Your estimate of $\delta$ is your treatment effect.
It should be noted, though, that the pattern of the policy dummy matters in DiD settings. See this post where the pattern of the treatment dummy is addressed in detail. The main issue is the "always treated" jurisdictions. It will be difficult—actually impossible—to demonstrate (visually) common trends without sufficient "pre-event" data. You could circumvent this issue by excluding the "always treated" from your analysis. To Dimitriy's point (see comments), this severely limits your sample, leaving you with approximately 30 states, 8 of which will be treated as some point between 1984–2018. You will also have to look into some finite-sample adjustments to deal with your standard errors.
I don't know of any valid methods for dealing with the "always treated" (i.e., always equal to 1) states in these settings. It will be difficult to make claims regarding the inter-temporal evolution of the "always treated" state trends relative to the "never treated" state trends (or soon to be treated state trends) before the policy goes into effect—because you never observe them (i.e., the "always treated") pre-policy! A surrogate measure might help if it allows you to obtain more pre-event data and you believe treatment affects the surrogate outcome in the same way.
With respect to the staggered adoption states, you could center your states on the official adoption year and assess parallelism of the group trends. If possible, try and obtain more pre-policy data on the "always treated" states. Or, maybe you could subject the "always treated" to a separate analysis. Do these states repeal any of the prohibition policies over time? This might be something to consider.
You might also find the research by Imai and Kim (2020) quite interesting.
Update (See Comments)
So, no interactions in the model? Since you say 'dummies', I should use dummy coding and not effect coding.
This is still an interaction model, we just define it in a different way. I am going to simulate some datasets in R to help with your intuition. Let's consider a panel dataset with 3 states (i.e., New York, New Jersey, and Connecticut) observed across 6 years (i.e., 2010–2015). Treatment is at the state level. I live in the Tri-State area so the decision regarding which states to use is completely arbitrary.
Now suppose treatment goes into effect in two states (e.g., New Jersey and Connecticut) in 2013 and stays in effect. The roll-out phase is not staggered. In fact, treatment begins at precisely the same time for all treated states. Note, proh is your interaction term (i.e., proh = treat*post) denoting the prohibition policy variable. See the toy dataset below which details the "classical" DiD setup:
# --- The "classical" difference-in-differences setup
# state = state id variable (3 states)
# year = year id variable (6 years)
# treat = dummy for treated states (NY is the control state)
# post = dummy for all years after treatment in both treatment and control groups
# proh = prohibition policy variable (proh = treat*post)
# A tibble: 18 x 5
state year treat post proh
<fct> <int> <dbl> <dbl> <dbl>
1 NY 2010 0 0 0
2 NY 2011 0 0 0
3 NY 2012 0 0 0
4 NY 2013 0 1 0
5 NY 2014 0 1 0
6 NY 2015 0 1 0
7 NJ 2010 1 0 0
8 NJ 2011 1 0 0
9 NJ 2012 1 0 0
10 NJ 2013 1 1 1
11 NJ 2014 1 1 1
12 NJ 2015 1 1 1
13 CT 2010 1 0 0
14 CT 2011 1 0 0
15 CT 2012 1 0 0
16 CT 2013 1 1 1
17 CT 2014 1 1 1
18 CT 2015 1 1 1
Now suppose the new law/policy is enacted in New Jersey and Connecticut in different years. How do we define post in this setting? Clearly, we cannot use the post variable in the "classical" case, as it will vary across states. You can proceed by defining a new variable called after to index the post-treatment years specific to treated jurisdictions. In my toy example (see below), New Jersey adopts the policy early (i.e., 'turns on' in 2012), while Connecticut is a late adopter (i.e., 'turns on' in 2014). Note, proh (i.e., proh = treat*after) instantiates the policy variable appropriately. To be clear, proh is the policy variable in the empirical specification above (i.e., $\text{Prohibition}_{st}$). You can think of this as the interaction of a treatment indicator with a time dummy indexing the "after" periods. But I recommend you create this variable manually!
In the example below, imagine we first created a column of 0's and we named it proh. Starting from the top row and working our way down the column, we replace a state-year cell with a value of 1 if it meets two conditions: it is a treated state and it is in the post-treatment period. I only simulated the second dataset to help with your intuition. See below:
# --- Coding of the treatment dummy in a "generalized" difference-in-differences setting
# state = state id variable (3 states)
# year = year id variable (6 years)
# treat = dummy for treated states (NY is the control state)
# after = dummy for time periods when treatment 'switches on'
# proh = prohibition policy variable (proh = treat*after)
# A tibble: 18 x 5
state year treat after proh
<fct> <int> <dbl> <dbl> <dbl>
1 NY 2010 0 0 0
2 NY 2011 0 0 0
3 NY 2012 0 0 0
4 NY 2013 0 0 0
5 NY 2014 0 0 0
6 NY 2015 0 0 0
7 NJ 2010 1 0 0
8 NJ 2011 1 0 0
9 NJ 2012 1 1 1
10 NJ 2013 1 1 1
11 NJ 2014 1 1 1
12 NJ 2015 1 1 1
13 CT 2010 1 0 0
14 CT 2011 1 0 0
15 CT 2012 1 0 0
16 CT 2013 1 0 0
17 CT 2014 1 1 1
18 CT 2015 1 1 1
The policy variable (i.e., proh) is the variable of interest. The terms treat and after act as placeholders in this example. I used them to demonstrate how we code the interaction term in a "generalized" DiD setting. The treatment dummy (i.e., proh) can accommodate multiple treatment regimes. Suppose New Jersey repeals their policy early. In that case, the treatment dummy must reflect this new "off" period of treatment. The same formula applies: the dummy equals 1 if you are a treated state and in the post-treatment period, 0 in all other conditions. It is typical in policy analysis to observe units (i.e., states, counties, firms, etc.) move into and out of treatment status multiple times. Again, the same formula applies.
To put this all together, a "generalized" DiD equation regresses your outcome on a full set of state dummies, a full set of year dummies, and the policy variable (i.e., proh). Again, proh is equal to 1 in all state-year combinations when the policy is in effect, 0 otherwise! Here is how the final dataset should look. I excluded one state dummy (i.e., NY) and one year dummy (i.e., 2010) to avoid collinearity.
# --- The "generalized" difference-in-differences setup
# state = state id variable (3 states)
# year = year id variable (6 years)
# NJ = dummy for New Jersey (NY is dropped)
# CT = dummy for Connecticut (NY is dropped)
# y_ = a set of year indicators (2010 is dropped)
# proh = prohibition policy variable
# --- Staggered adoption summary:
# NY never adopts (always '0')
# NJ adopts early ('turns on' in 2012)
# CT adopts late ('turns on' in 2014)
# A tibble: 18 x 10
state year NJ CT y_2011 y_2012 y_2013 y_2014 y_2015 proh
<fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 NY 2010 0 0 0 0 0 0 0 0
2 NY 2011 0 0 1 0 0 0 0 0
3 NY 2012 0 0 0 1 0 0 0 0
4 NY 2013 0 0 0 0 1 0 0 0
5 NY 2014 0 0 0 0 0 1 0 0
6 NY 2015 0 0 0 0 0 0 1 0
7 NJ 2010 1 0 0 0 0 0 0 0
8 NJ 2011 1 0 1 0 0 0 0 0
9 NJ 2012 1 0 0 1 0 0 0 1
10 NJ 2013 1 0 0 0 1 0 0 1
11 NJ 2014 1 0 0 0 0 1 0 1
12 NJ 2015 1 0 0 0 0 0 1 1
13 CT 2010 0 1 0 0 0 0 0 0
14 CT 2011 0 1 1 0 0 0 0 0
15 CT 2012 0 1 0 1 0 0 0 0
16 CT 2013 0 1 0 0 1 0 0 0
17 CT 2014 0 1 0 0 0 1 0 1
18 CT 2015 0 1 0 0 0 0 1 1
To go back to my earlier point, all treated states have pre- and post-policy periods. It will be easier to defend the validity of the methodology in this setting. In the real world, however, if too many states are "always treated" (i.e., always equal to 1), then it will be difficult to disentangle a treatment effect from some underlying trend that may have already been in place in treated states.
As for your code, there are multiple ways for you to proceed. Two specifications are presented below. In my opinion, the last equation is preferred.
# explicit but messy
did_1 <- lm(outcome ~ NJ + CT + y_2011 + y_2012 + y_2013 + y_2014 + y_2015 + proh + ..., data = ...)
# preferred approach
did_2 <- lm(outcome ~ as.factor(state) + as.factor(year) + proh + ..., data = ...)
In the latter equation, as.factor(state) and as.factor(year) will 'dummy out' all of your states and years for free and with no additional work on your part. You can think of as.factor(state) (i.e., state fixed effects) and as.factor(year) (i.e., year fixed effects) as replacing treat and post in the "classical" DiD equation, respectively. You could certainly create these variables manually and include them inside of the lm() function, but this method is error-prone. Imagine if you acquired a panel with 1,000 individuals; it would be unwieldy to code these dummies by hand.
Note, each formulaic expression only includes the policy variable! Any use of a post variable is not appropriate as the post-treatment epoch varies across states.
