Diff in Diff Margins and Treatment Effect Heterogeneity

[Maxence Morlet]

Hi all,

I have a triple difference coefficient: the interaction between three dummy variables.These variables are post, treat_ind and treat_county:

1.post#1.treat_ind#1.treat_county

I have another variable, which is county specific. It is called "devel". It is continuous. I am certain that the treatment effect, i.e. the effect of the triple difference coefficient, varies across the levels of the continuous variable "devel". Thus I want to interact it with the triple difference coefficient in a regression. Then, I want to take the marginal effect of the triple difference coefficient at each value of "devel". I start with the following code: (reghdfe from ssc install)

levelsof devel, local(develvals)

reghdfe dep_var 1.post#1.treat_ind#1.treat_county 1.post#1.treat_ind#1.treat_county#c.devel, abs(i.countyfips#i.ind i.ind#i.time i.countyfips#i.time) cluster(countyfips)

margins, dydx(1.post#1.treat_ind#1.treat_county) at(devel=(`develvals'))

I am aware that 1.post#1.treat_ind#1.treat_county is an interaction itself.

So will the output of this margins command give me the average marginal effect of the triple difference coefficient on dep_var at each of the values of the "devel" variable?

[Maxence Morlet]

For anyone interested in the same question, I just came across the following two threads:

stata - How to run an Instrumental Variable Regression for heterogenous treatment effects - Cross Validated (stackexchange.com)

Heterogenous treatment effect with difference-in-differences - Statalist

[Dimitriy V. Masterov]

I don't think this code will run as written. Margins will not know how to calculate that derivative.

You could do something like this using nlcom:

cls
use "http://fmwww.bc.edu/repec/bocode/c/CardKrueger1994.dta", clear
drop if id == 407
xtset id t
gen chain = ""
foreach var of varlist bk kfc roys wendys {
    replace chain="`var'" if `var'==1
 }
ssc install sencode
sencode chain, replace
fvset base 0 t
fvset base 1 chain // base is Burger King
fvset base 0 treated
reghdfe fte i.treated##i.t##i.chain, abs(i.id) cluster(id)
nlcom ///
(bk:     _b[1.treated#1.t]) ///
(kfc:    _b[1.treated#1.t] + _b[1.treated#1.t#2.chain]) ///
(roys:   _b[1.treated#1.t] + _b[1.treated#1.t#3.chain]) ///
(wendys: _b[1.treated#1.t] + _b[1.treated#1.t#4.chain]), post // coefl
test _b[bk] = _b[kfc] = _b[roys] = _b[wendys]

Here you would reject the null that the effects of minimum wage are the same for all four restaurant chains.

[Dimitriy V. Masterov]

Here's a continuous DD interaction example that can be extended to DDD:

cls
use "http://fmwww.bc.edu/repec/bocode/c/CardKrueger1994.dta", clear
drop if id == 407
xtset id t
gen chain = ""
foreach var of varlist bk kfc roys wendys {
    replace chain="`var'" if `var'==1
 }
sencode chain, replace
fvset base 0 t
fvset base 1 chain // base is Burger King
fvset base 0 treated
reghdfe fte i.treated##i.t##c.chain, abs(i.id) cluster(id) // coefl
nlcom (effect:_b[1.treated#1.t] + _b[1.treated#1.t#c.chain]), post // coefl
test _b[effect] = 0

[Maxence Morlet]

Thank you!

[Maxence Morlet]

Dear Dimitriy V. Masterov,

Thank you so much for this. I have implemented your method.

Out of curiosity, is the following correct:

Suppose I were to dichotomise the "devel" variable at the median, and call the resulting binary variable "above_median_devel".

Then, forgetting the margins command, I just run:

 
 reghdfe dep_var 1.post#1.treat_ind#1.treat_county 1.post#1.treat_ind#1.treat_county#1.above_median_devel, abs(i.countyfips#i.ind i.ind#i.time i.countyfips#i.time) cluster(countyfips)

Is the following interpretation correct:

1.post#1.treat_ind#1.treat_county

is the triple difference treatment effect estimate for below median level of devel.

1.post#1.treat_ind#1.treat_county#1.above_median_devel

is the triple difference treatment effect estimate for above median level of devel.

If

1.post#1.treat_ind#1.treat_county#1.above_median_devel

is statistically significant and positive, then the triple difference treatment effect is significantly stronger for above median devel counties.

Many thanks in advance!

[Dimitriy V. Masterov]

I believe you should add
_b[1.post#1.treat_ind#1.treat_county] and
_b[1.post#1.treat_ind#1.treat_county#1.above_median_d evel] to get the effect for the above median units. Your second statement is correct.

[Maxence Morlet]

Dear Dimitriy V. Masterov,

Thank you so much for your answer!

Apologies, I have simply misunderstood something:

Is there something to add to the following equation to achieve the goal stated in #6?

 
 reghdfe dep_var 1.post#1.treat_ind#1.treat_county 1.post#1.treat_ind#1.treat_county#1.above_median_devel, abs(i.countyfips#i.ind i.ind#i.time i.countyfips#i.time) cluster(countyfips)

[Dimitriy V. Masterov]

The best way to make sure you are doing things correctly is to calculate all the means by hand, put them into the DDD formula, and make sure the regression version returns the same numbers as the formula.

[Dimitriy V. Masterov]

I suspect your code is possibly quite wrong. Let's use the C&K data to answer the question: Is the effect of a higher minimum wage greater for chicken restaurants if they had an initial FTE above the median for their chain? That can be calculated in two ways:

1. Calculate DDD for above-median chicken restaurants. Calculate the DDD for below-median chicken restaurants. Difference them to get the answer, which is a quad-diff (DDDD).
2. Estimate a DDDD using regression.

The advantage of approach 1 is that it demonstrates that a DDDD is really a difference between two DDDs. You can even go further. A DDD is a difference of two DDs, and a DD is a difference of two Ds. You could do all this by hand since it's just averages and division. Yay! The disadvantage is that it does not return any standard errors. Boo! But you should do (1) to verify (2) is doing what you expect.

Here's an example using the same data as above:

. ssc install diff
checking diff consistency and verifying not already installed...
all files already exist and are up to date.

. net from https://www.sealedenvelope.com/
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
https://www.sealedenvelope.com/
Programs from Sealed Envelope Ltd
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

These programs are provided by Sealed Envelope Ltd.
For more details on their use visit www.sealedenvelope.com

PACKAGES you could -net describe-:
    hl                Hosmer-Lemeshow goodness of fit test
    mlogplot2         Extension to mlogplot to handle two-way interaction terms in mlogit
    reformat          Reformat regression output
    rescale           Rescale categorical variables using numbers found in value labels
    time              Convert strings in HH:MM:SS format to elapsed times and back again
    xcount            Count longitudinal data
    xfill             Fill in static variables in longitudinal data
    xtab              Tabulate longitudinal data
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

. net install xfill
checking xfill consistency and verifying not already installed...
all files already exist and are up to date.

. use "http://fmwww.bc.edu/repec/bocode/c/CardKrueger1994.dta", clear
(Dataset from Card&Krueger (1994))

. drop if id == 407
(4 observations deleted)

. xtset id t

Panel variable: id (strongly balanced)
 Time variable: t, 0 to 1
         Delta: 1 unit

. gen chain = ""
(816 missing values generated)

. foreach var of varlist bk kfc roys wendys {
  2.         replace chain="`var'" if `var'==1
  3.  }
variable chain was str1 now str2
(342 real changes made)
variable chain was str2 now str3
(158 real changes made)
variable chain was str3 now str4
(198 real changes made)
variable chain was str4 now str6
(118 real changes made)

. sencode chain, replace

. gen chix = cond(kfc == 1, 1, 0)

. lab define chix 0 "Burger Joint" 1 "Fried Chicken"

. lab val chix chix

. bys chain: egen chain_median_fte = median(fte) if t == 0
(408 missing values generated)

. gen above_pre_chain_median = cond(fte > chain_median_fte, 1, 0) if !missing(chain_median_fte)
(408 missing values generated)

. xfill above_pre_chain_median, i(id)

. fvset base 0 t

. fvset base 0 chix // base is Burger Chain

. fvset base 0 treated

. fvset base 0 above_pre_chain_median

. diff fte if above_pre_chain_median == 1, t(treated) p(t) ddd(chix) cluster(chain)

TRIPLE DIFFERENCE (DDD) ESTIMATION RESULTS
Notation of DDD:
   Control (A)     treated = 0 and chix = 1
   Control (B)     treated = 0 and chix = 0
   Treated (A)     treated = 1 and chix = 1
   Treated (B)     treated = 1 and chix = 0

Number of observations in the DDD: 379
               Before      After    
   Control (A):3           3           6
   Control (B):35          35          70
   Treated (A):34          34          68
   Treated (B):118         117         235
               190         189
--------------------------------------------------------
 Outcome var.   | fte     | S. Err. |   |t|   |  P>|t|
----------------+---------+---------+---------+---------
Before          |         |         |         | 
   Control (A)  | 11.500  |         |         | 
   Control (B)  | 29.207  |         |         | 
   Treated (A)  | 12.949  |         |         | 
   Treated (B)  | 25.619  |         |         | 
   Diff (T-C)   | 5.037   | 2.100   | 2.40    | 0.096*
After           |         |         |         | 
   Control (A)  | 12.833  |         |         | 
   Control (B)  | 21.700  |         |         | 
   Treated (A)  | 12.037  |         |         | 
   Treated (B)  | 22.949  |         |         | 
   Diff (T-C)   | -2.045  | 1.617   | 1.27    | 0.295
                |         |         |         | 
DDD             | -7.082  | 1.149   | -6.16   | 0.009***
--------------------------------------------------------
R-square:    0.31
* Means and Standard Errors are estimated by linear regression
**Clustered Std. Errors
**Inference: *** p<0.01; ** p<0.05; * p<0.1

. scalar ddd1 = r(ddd)

. diff fte if above_pre_chain_median == 0, t(treated) p(t) ddd(chix) cluster(chain)

TRIPLE DIFFERENCE (DDD) ESTIMATION RESULTS
Notation of DDD:
   Control (A)     treated = 0 and chix = 1
   Control (B)     treated = 0 and chix = 0
   Treated (A)     treated = 1 and chix = 1
   Treated (B)     treated = 1 and chix = 0

Number of observations in the DDD: 418
               Before      After    
   Control (A):8           8           16
   Control (B):31          30          61
   Treated (A):34          34          68
   Treated (B):139         134         273
               212         206
--------------------------------------------------------
 Outcome var.   | fte     | S. Err. |   |t|   |  P>|t|
----------------+---------+---------+---------+---------
Before          |         |         |         | 
   Control (A)  | 6.375   |         |         | 
   Control (B)  | 14.202  |         |         | 
   Treated (A)  | 6.765   |         |         | 
   Treated (B)  | 13.381  |         |         | 
   Diff (T-C)   | 1.210   | 0.463   | 2.62    | 0.079*
After           |         |         |         | 
   Control (A)  | 8.938   |         |         | 
   Control (B)  | 15.775  |         |         | 
   Treated (A)  | 9.213   |         |         | 
   Treated (B)  | 16.450  |         |         | 
   Diff (T-C)   | -0.399  | 1.288   | 0.31    | 0.777
                |         |         |         | 
DDD             | -1.609  | 1.569   | -1.03   | 0.381
--------------------------------------------------------
R-square:    0.24
* Means and Standard Errors are estimated by linear regression
**Clustered Std. Errors
**Inference: *** p<0.01; ** p<0.05; * p<0.1

. scalar ddd0= r(ddd)

. di scalar(ddd1) - scalar(ddd0)
-5.4733994

. reghdfe fte i.treated##i.t##i.chix##i.above_pre_chain_median, cluster(chain) //coefl
(MWFE estimator converged in 1 iterations)
warning: missing F statistic; dropped variables due to collinearity or too few clusters

HDFE Linear regression                            Number of obs   =        797
Absorbing 1 HDFE group                            F(  15,      3) =          .
Statistics robust to heteroskedasticity           Prob > F        =          .
                                                  R-squared       =     0.4519
                                                  Adj R-squared   =     0.4414
                                                  Within R-sq.    =     0.4519
Number of clusters (chain)   =          4         Root MSE        =     6.7457

                                                           (Std. err. adjusted for 4 clusters in chain)
-------------------------------------------------------------------------------------------------------
                                      |               Robust
                                  fte | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
--------------------------------------+----------------------------------------------------------------
                              treated |
                                  NJ  |  -.8203179   .4630548    -1.77   0.175    -2.293965     .653329
                                  1.t |   1.573387    2.60247     0.60   0.588    -6.708833    9.855607
                                      |
                            treated#t |
                                NJ#1  |   1.494945   1.570513     0.95   0.411    -3.503129    6.493018
                                      |
                                 chix |
                       Fried Chicken  |  -7.826613   .3275663   -23.89   0.000    -8.869075   -6.784151
                                      |
                         treated#chix |
                    NJ#Fried Chicken  |   1.210024   .4630548     2.61   0.079    -.2636231    2.683671
                                      |
                               t#chix |
                     1#Fried Chicken  |   .9891129    2.60247     0.38   0.729    -7.293107    9.271333
                                      |
                       treated#t#chix |
                  NJ#1#Fried Chicken  |  -1.608915   1.570513    -1.02   0.381    -6.606989    3.389158
                                      |
             1.above_pre_chain_median |   15.00553   1.941979     7.73   0.005     8.825286    21.18577
                                      |
       treated#above_pre_chain_median |
                                NJ#1  |  -2.768181   2.250544    -1.23   0.306    -9.930417    4.394055
                                      |
             t#above_pre_chain_median |
                                 1 1  |   -9.08053    2.64601    -3.43   0.041    -17.50131   -.6597449
                                      |
     treated#t#above_pre_chain_median |
                              NJ#1#1  |   3.342272   2.083684     1.60   0.207    -3.288942    9.973485
                                      |
          chix#above_pre_chain_median |
                     Fried Chicken#1  |   -9.88053   1.941979    -5.09   0.015    -16.06077   -3.700286
                                      |
  treated#chix#above_pre_chain_median |
                  NJ#Fried Chicken#1  |   3.827004   2.250544     1.70   0.188    -3.335232    10.98924
                                      |
        t#chix#above_pre_chain_median |
                   1#Fried Chicken#1  |   7.851363    2.64601     2.97   0.059    -.5694218    16.27215
                                      |
treated#t#chix#above_pre_chain_median |
                NJ#1#Fried Chicken#1  |  -5.473399   2.083684    -2.63   0.079    -12.10461    1.157814
                                      |
                                _cons |   14.20161   .3275663    43.35   0.000     13.15915    15.24408
-------------------------------------------------------------------------------------------------------

. nlcom ///
> (DDD_chix_above:_b[1.treated#1.t#1.chix] + _b[1.treated#1.t#1.chix#1.above_pre_chain_median]) ///
> (DDD_chix_below:_b[1.treated#1.t#1.chix]) ///
> (DDDD:_b[1.treated#1.t#1.chix#1.above_pre_chain_median]), post

DDD_chix_a~e: _b[1.treated#1.t#1.chix] + _b[1.treated#1.t#1.chix#1.above_pre_chain_median]
DDD_chix_b~w: _b[1.treated#1.t#1.chix]
        DDDD: _b[1.treated#1.t#1.chix#1.above_pre_chain_median]

--------------------------------------------------------------------------------
           fte | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
DDD_chix_above |  -7.082315   1.149446    -6.16   0.000    -9.335187   -4.829443
DDD_chix_below |  -1.608915   1.570513    -1.02   0.306    -4.687064    1.469234
          DDDD |  -5.473399   2.083684    -2.63   0.009    -9.557346   -1.389453
--------------------------------------------------------------------------------

. coefplot, vertical ylab(#20) yline(0) xtitle("DDD_chix_above - DDD_chix_below = DDDD") ytitle("Effects (FTEs)") name("DDDD", replace)

. reghdfe fte 1.t#1.treated#1.chix 1.t#1.treated#1.chix#1.above_pre_chain_median, cluster(id)
(MWFE estimator converged in 1 iterations)

HDFE Linear regression                            Number of obs   =        797
Absorbing 1 HDFE group                            F(   2,    407) =      57.67
Statistics robust to heteroskedasticity           Prob > F        =     0.0000
                                                  R-squared       =     0.0584
                                                  Adj R-squared   =     0.0561
                                                  Within R-sq.    =     0.0584
Number of clusters (id)      =        408         Root MSE        =     8.7690

                                                            (Std. err. adjusted for 408 clusters in id)
-------------------------------------------------------------------------------------------------------
                                      |               Robust
                                  fte | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
--------------------------------------+----------------------------------------------------------------
                       t#treated#chix |
                  1#NJ#Fried Chicken  |  -9.076545   .9120556    -9.95   0.000    -10.86947   -7.283617
                                      |
t#treated#chix#above_pre_chain_median |
                1#NJ#Fried Chicken#1  |   2.823529   1.366058     2.07   0.039     .1381199    5.508939
                                      |
                                _cons |   18.28978   .3967238    46.10   0.000      17.5099    19.06966
-------------------------------------------------------------------------------------------------------

Note that (1) gets me -5.5 employees. Note that (2) also gets me -5.5 employees. My version of your code gets me +2.8. The issue is your code only includes the two final interactions and omits the lower-level ones and the dummies. It's possible that the variables that go into your absorb() take care of those terms, but looking at the names (since no other details were provided), suggests that may not be the case.

Now for your second question: How to estimate a DDD that is moderated by a continuous variable? Here is an example where I treat chain as if it were continuous and calculate a DDD(chain):

. reghdfe fte i.treated##i.t##i.above_pre_chain_median##c.chain, cluster(chain) //coefl
(MWFE estimator converged in 1 iterations)
warning: missing F statistic; dropped variables due to collinearity or too few clusters

HDFE Linear regression                            Number of obs   =        797
Absorbing 1 HDFE group                            F(  15,      3) =          .
Statistics robust to heteroskedasticity           Prob > F        =          .
                                                  R-squared       =     0.2817
                                                  Adj R-squared   =     0.2679
                                                  Within R-sq.    =     0.2817
Number of clusters (chain)   =          4         Root MSE        =     7.7225

                                                              (Std. err. adjusted for 4 clusters in chain)
----------------------------------------------------------------------------------------------------------
                                         |               Robust
                                     fte | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-----------------------------------------+----------------------------------------------------------------
                                 treated |
                                     NJ  |   .8645841   1.660515     0.52   0.639    -4.419916    6.149085
                                     1.t |   7.469999   1.931494     3.87   0.031     1.323122    13.61688
                                         |
                               treated#t |
                                   NJ#1  |  -3.500525   .6592664    -5.31   0.013    -5.598605   -1.402445
                                         |
                1.above_pre_chain_median |   20.48461   1.631085    12.56   0.001     15.29377    25.67545
                                         |
          treated#above_pre_chain_median |
                                   NJ#1  |  -11.24226   4.705066    -2.39   0.097    -26.21588    3.731357
                                         |
                t#above_pre_chain_median |
                                    1 1  |  -15.64305   2.294632    -6.82   0.006    -22.94559   -8.340508
                                         |
        treated#t#above_pre_chain_median |
                                 NJ#1#1  |   9.937227   .4226838    23.51   0.000     8.592059     11.2824
                                         |
                                   chain |   .4870996   1.102913     0.44   0.689    -3.022862    3.997061
                                         |
                         treated#c.chain |
                                     NJ  |  -.5937672   .6863989    -0.87   0.451    -2.778195     1.59066
                                         |
                               t#c.chain |
                                      1  |  -2.437943    1.07156    -2.28   0.107    -5.848124    .9722377
                                         |
                       treated#t#c.chain |
                                   NJ#1  |   1.917762   .2365613     8.11   0.004     1.164918    2.670606
                                         |
          above_pre_chain_median#c.chain |
                                      1  |  -2.597644   .9461896    -2.75   0.071    -5.608842    .4135539
                                         |
  treated#above_pre_chain_median#c.chain |
                                   NJ#1  |   3.278585   1.745245     1.88   0.157    -2.275563    8.832732
                                         |
        t#above_pre_chain_median#c.chain |
                                    1 1  |   3.078708   .8008859     3.84   0.031     .5299322    5.627485
                                         |
treated#t#above_pre_chain_median#c.chain |
                                 NJ#1#1  |  -2.810856   .4227301    -6.65   0.007    -4.156172    -1.46554
                                         |
                                   _cons |   11.43461    3.88128     2.95   0.060    -.9173562    23.78657
----------------------------------------------------------------------------------------------------------

. levelsof chain 
1 2 3 4

. nlcom ///
> (DDD_chain_1:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*1) ///
> (DDD_chain_2:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*2) ///
> (DDD_chain_3:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*3) ///
> (DDD_chain_4:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*4) ///
> , post // coefl

 DDD_chain_1: _b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*1
 DDD_chain_2: _b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*2
 DDD_chain_3: _b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*3
 DDD_chain_4: _b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*4

------------------------------------------------------------------------------
         fte | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
 DDD_chain_1 |   7.126371   .3909049    18.23   0.000     6.360212    7.892531
 DDD_chain_2 |   4.315516    .695955     6.20   0.000     2.951469    5.679562
 DDD_chain_3 |    1.50466    1.08319     1.39   0.165    -.6183545    3.627674
 DDD_chain_4 |  -1.306196   1.489849    -0.88   0.381    -4.226247    1.613855
------------------------------------------------------------------------------

. coefplot, vertical ylab(#20) xlab(1 2 3 4) yline(0) xtitle("Chain") ytitle("DDD Effect (FTEs)") name("DDD", replace)

. // should match _b[1.treated#1.t#1.above_pre_chain_median#c.chain] exactly
. nlcom ///
> (_b[DDD_chain_4] - _b[DDD_chain_3]) ///
> (_b[DDD_chain_3] - _b[DDD_chain_2]) ///
> (_b[DDD_chain_2] - _b[DDD_chain_1])

       _nl_1: _b[DDD_chain_4] - _b[DDD_chain_3]
       _nl_2: _b[DDD_chain_3] - _b[DDD_chain_2]
       _nl_3: _b[DDD_chain_2] - _b[DDD_chain_1]

------------------------------------------------------------------------------
         fte | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |  -2.810856   .4227301    -6.65   0.000    -3.639392    -1.98232
       _nl_2 |  -2.810856   .4227301    -6.65   0.000    -3.639392    -1.98232
       _nl_3 |  -2.810856   .4227301    -6.65   0.000    -3.639392    -1.98232
------------------------------------------------------------------------------

This can be interpreted as "the DDD effect of the minimum wage for an above-median restaurant declines by -2.8 employees each time the chain variable increments by 1." If you like calculus, this is the dDDD/dchain. There is substantial heterogeneity. The DDD effect at chain = 1 is +7.1. The effect at chain=4 is -1.3 employees and no longer significantly different from zero.

[Dimitriy V. Masterov]

Here's the code for the two examples above:

cls
ssc install diff
net from https://www.sealedenvelope.com/
net install xfill
use "http://fmwww.bc.edu/repec/bocode/c/CardKrueger1994.dta", clear
drop if id == 407
xtset id t
gen chain = ""
foreach var of varlist bk kfc roys wendys {
    replace chain="`var'" if `var'==1
 }
sencode chain, replace
gen chix = cond(kfc == 1, 1, 0)
lab define chix 0 "Burger Joint" 1 "Fried Chicken"
lab val chix chix
bys chain: egen chain_median_fte = median(fte) if t == 0
gen above_pre_chain_median = cond(fte > chain_median_fte, 1, 0) if !missing(chain_median_fte)
xfill above_pre_chain_median, i(id)
fvset base 0 t
fvset base 0 chix // base is Burger Chain
fvset base 0 treated
fvset base 0 above_pre_chain_median
diff fte if above_pre_chain_median == 1, t(treated) p(t) ddd(chix) cluster(chain)
scalar ddd1 = r(ddd)
diff fte if above_pre_chain_median == 0, t(treated) p(t) ddd(chix) cluster(chain)
scalar ddd0= r(ddd)
di scalar(ddd1) - scalar(ddd0)
reghdfe fte i.treated##i.t##i.chix##i.above_pre_chain_median, cluster(chain) //coefl
nlcom ///
(DDD_chix_above:_b[1.treated#1.t#1.chix] + _b[1.treated#1.t#1.chix#1.above_pre_chain_median]) ///
(DDD_chix_below:_b[1.treated#1.t#1.chix]) ///
(DDDD:_b[1.treated#1.t#1.chix#1.above_pre_chain_median]), post
coefplot, vertical ylab(#20) yline(0) xtitle("DDD_chix_above - DDD_chix_below = DDDD") ytitle("Effects (FTEs)") name("DDDD", replace)
reghdfe fte 1.t#1.treated#1.chix 1.t#1.treated#1.chix#1.above_pre_chain_median, cluster(id)

reghdfe fte i.treated##i.t##i.above_pre_chain_median##c.chain, cluster(chain) //coefl
levelsof chain 
nlcom ///
(DDD_chain_1:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*1) ///
(DDD_chain_2:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*2) ///
(DDD_chain_3:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*3) ///
(DDD_chain_4:_b[1.treated#1.t#1.above_pre_chain_median] + _b[1.treated#1.t#1.above_pre_chain_median#c.chain]*4) ///
, post // coefl
coefplot, vertical ylab(#20) xlab(1 2 3 4) yline(0) xtitle("Chain") ytitle("DDD Effect (FTEs)") name("DDD", replace)
// should match _b[1.treated#1.t#1.above_pre_chain_median#c.chain] exactly
nlcom ///
(_b[DDD_chain_4] - _b[DDD_chain_3]) ///
(_b[DDD_chain_3] - _b[DDD_chain_2]) ///
(_b[DDD_chain_2] - _b[DDD_chain_1])

[Maxence Morlet]

Thank you so much! I will try it out! Really, thanks for your time and effort.