Full versus partial marginal effect using fixest package

Question

I have would like to know the full marginal effect of the continuous variable provtariff given the interaction term Female * provtariff on the outcome variable log(totalinc) as well as the coefficient of the interaction term.

Using the code:

  feols(log(totalinc) ~ i(Female, provtariff) | hhid02 + year,
        data = inc0402_p,
        weights = ~hhwt, 
        vcov = ~tinh)

I got the following results

OLS estimation, Dep. Var.: log(totalinc)
Observations: 24,966 
Weights: hhwt 
Fixed-effects: hhid02: 11,018,  year: 2
Standard-errors: Clustered (tinh) 
                     Estimate Std. Error t value  Pr(>|t|)    
Female::0:provtariff  5.79524    1.84811 3.13577 0.0026542 ** 
Female::1:provtariff  2.66994    2.09540 1.27419 0.2075088    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 7.61702     Adj. R2: 0.670289
                Within R2: 0.045238s

However, when I implement the following code

feols(log(totalinc) ~ Female*provtariff | hhid02 + year,
        data = inc0402_p,
        weights = ~hhwt, 
        vcov = ~tinh)

I get the following results

OLS estimation, Dep. Var.: log(totalinc)
Observations: 24,966 
Weights: hhwt 
Fixed-effects: hhid02: 11,018,  year: 2
Standard-errors: Clustered (tinh) 
                   Estimate Std. Error  t value   Pr(>|t|)    
Female            -0.290019   0.029894 -9.70142 6.6491e-14 ***
provtariff         4.499561   1.884625  2.38751 2.0130e-02 *  
Female:provtariff -0.433963   0.170505 -2.54516 1.3512e-02 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 7.52022     Adj. R2: 0.678592
                Within R2: 0.069349

Should the provtariff coefficient in the latter model not be the same as the coefficient for Female::0:provtariff in the first model?

Answer 1

No, clearly the two models are different because one includes two parameters and the other one includes 3. They won’t produce equivalent results. More specifically, one of your models includes only the interactions, but no “constitutive” term, whereas the other model includes both.

Here is a reproducible example with a 3rd model that reproduces your model with the * asterisk, but uses the fixest interaction syntax with i(). You’ll see that some of the coefficients and standard errors are exactly identical to those in the 2nd model, and that R2 are the same. This suggests that m2 and m3 are equivalent:

library(fixest)
library(modelsummary)
library(marginaleffects)

# Your 
m1 <- feols(mpg ~ i(am, hp) | gear, data = mtcars)
m2 <- feols(mpg ~ am * hp | gear, data = mtcars)
m3 <- feols(mpg ~ am + i(am, hp) | gear, data = mtcars)

models <- list(m1, m2, m3)
modelsummary(models)

	(1)	(2)	(3)
am = 0 × hp	-0.076		-0.056
	(0.025)		(0.006)
am = 1 × hp	-0.059		-0.071
	(0.009)		(0.021)
am		5.568	5.568
		(1.575)	(1.575)
hp		-0.056
		(0.006)
am × hp		-0.015
		(0.019)
Num.Obs.	32	32	32
R2	0.763	0.797	0.797
Std.Errors	by: gear	by: gear	by: gear
FE: gear	X	X	X

We can further check the equivalence between models 2 and 3 by computing the partial derivative of the outcome with respect to one of the predictors. In economics they call this slope a “marginal effect”, although the terminology changes across disciplines, and I am not sure if that is the quantity you are interested in when you say “marginal effects”:

marginaleffects(m2, variables = "hp") |> summary()
#>   Term    Contrast Estimate Std. Error      z   Pr(>|z|)   2.5 %  97.5 %
#> 1   hp mean(dY/dX)   -0.062    0.01087 -5.705 1.1665e-08 -0.0833 -0.0407
#> 
#> Model type:  fixest 
#> Prediction type:  response

marginaleffects(m3, variables = "hp") |> summary()
#>   Term    Contrast Estimate Std. Error      z   Pr(>|z|)   2.5 %  97.5 %
#> 1   hp mean(dY/dX)   -0.062    0.01087 -5.705 1.1665e-08 -0.0833 -0.0407
#> 
#> Model type:  fixest 
#> Prediction type:  response