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](https://www.annualreviews.org/doi/pdf/10.1146/annurev-publhealth-040617-013507) 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} = \alpha + \gamma_{s} + \lambda_{t} + \delta \text{Prohibition}_{st} + \epsilon_{st},
$$
where you observe a continuous outcome in state $s$ in year $t$. $\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](https://stats.stackexchange.com/questions/455986/panel-data-diff-in-diff-and-the-pattern-of-the-binary-treatment-indicator) 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 this recent [working paper](http://web.mit.edu/insong/www/pdf/FEmatch-twoway.pdf) 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, `int` is interaction term (i.e., `int = treat*post`). See the toy dataset below which details the "classical" DiD 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
# int = interaction term (treat*post)
# A tibble: 18 x 5
state year treat post int
<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, `int` (i.e., `int = treat*after`) instantiates the policy variable appropriately. To be clear, `int` *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 `int`. 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:
# 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
# int = interaction term (treat*after)
# A tibble: 18 x 5
state year treat after int
<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., `int`) 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.
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., `int`). Again, `int` 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.
# 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)
# int = policy variable (prohibition)
# - treatment adoption is staggered
# - NY never adopts (always '0')
# - NJ adopts early ('1' in 2012)
# - CT adopts late ('1' in 2014)
# A tibble: 18 x 10
state year NJ CT y_2011 y_2012 y_2013 y_2014 y_2015 int
<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
Note, `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 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.
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 + int + ..., data = ...)
# preferred approach
did_2 <- lm(outcome ~ as.factor(state) + as.factor(year) + int + ..., data = ...)
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**.