How do I plot marginal effects of two-way fixed effects regression with interaction term

Question

I have the following regression:

tal_lpm4 <- feols(tal ~ provtariff + provtariff*i(sex) + i(sex) + as.numeric(educ) + age + age^2| year + tinh,
                         data = employment0204,
                         vcov_cluster(~year + tinh),
                    weights = ~hhwt)

OLS estimation, Dep. Var.: tal
Observations: 332,147 
Fixed-effects: year: 2,  tinh: 61
Standard-errors: Clustered (year & tinh) 
                        Estimate Std. Error   t value Pr(>|t|)    
provtariff              0.200287   0.101019   1.98266 0.297390    
sex::Female             0.038688   0.003849  10.05185 0.063126 .  
as.numeric(educ)        0.002093   0.000290   7.22092 0.087606 .  
age                     0.001084   0.000682   1.58972 0.357463    
I(age^2)               -0.000027   0.000014  -1.90659 0.307521    
provtariff:sex::Male   -0.139618   0.017293  -8.07366 0.078452 .  
provtariff:sex::Female -0.216309   0.014045 -15.40086 0.041279 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 6.23738     Adj. R2: 0.036039
                Within R2: 0.017567

I want to plot the marginal effect of provtariff of each sex. When I used the code

iplot(tal_lpm4)

the following plot is shown:

The Male variable's coefficient is 0, even though it should be the provtariff coefficient; whereas coefficient for Female should be provtariff * sex plus provtariff.

Answer 1

Turns out placement of the interaction term matters:

tal_lpm4 <- feols(tal ~ provtariff + provtariff + age + age^2 + i(sex, provtariff) + i(urban) +i(sex)| year + tinh,
                         data = employment0204,
                         vcov_cluster(~year + tinh),
                    weights = ~hhwt)

Gave me the desired outcome