r plm marginal-effects modelsummary r-marginaleffects

Display plm Model Results: Marginal Effects, p-Values for Interaction Terms, and Plots (Fixed Effects Method)

I'm working with panel data using the plm package in R to estimate a fixed effects model that includes interaction terms. I'm currently able to get my model results using plm and display them with stargazer, but I'm looking for a more refined way to present these results, specifically focusing on the interaction terms. Additionally, I need to calculate and display the marginal effects of these interaction terms, including their p-values, in the same table. I'm also interested in plotting these marginal effects to better understand their impact.

Here's an example of my plm model setup:

fe_model_interaction <- plm(life_satisfaction ~ 
                                  employment_level_Full_Time*care_low +
                                  employment_level_Full_Time*care_high +

                                  employment_level_Part_Time*care_low +
                                  employment_level_Part_Time*care_high+
                             
                                  other_control_variables,
                                  data = data_analyse_mother, 
                                  index = c("pid", "syear"), 
                                  model = "within")

My questions are:

What is the best way to display the plm model results, especially for interaction terms, in a more detailed and customizable format than what stargazer offers? I'm looking for a method to include standard errors, coefficients, and p-values in a clean and professional-looking table.
How can I calculate and include the marginal effects of these interaction terms in the results table, along with their p-values?
What are the recommended approaches or packages in R for plotting the marginal effects of interaction terms in a fixed effects model estimated with plm?

Solution

The marginaleffects and modelsummary both support plm objects. They both have extensive customization options and detailed tutorials. Please visit the package websites to access each package’s “Get Started” page:

(Disclaimer: I am the maintainer.)

Your questions are not specific enough for me to solve your exact problem, but here is an example. First, we fit a model and display the results in a regression table using modelsummary:

library(plm)
library(modelsummary)
library(marginaleffects)

data("Produc", package = "plm")

mod <- plm(log(gsp) ~ log(pcap) * log(pc) + log(emp) + unemp,
    data = Produc, index = c("state", "year"), model = "within")

modelsummary(mod)

	(1)
log(pcap)	0.033
	(0.074)
log(pc)	0.348
	(0.069)
log(emp)	0.763
	(0.031)
unemp	-0.005
	(0.001)
log(pcap) × log(pc)	-0.005
	(0.006)
Num.Obs.	816
R2	0.941
R2 Adj.	0.937
AIC	14092.6
BIC	14120.8
RMSE	0.04

Now, we use marginaleffects to compute quantities of interest. Note that the terminology around “marginal effects” is highly inconsistent across field, so I can’t know exactly what quantity you want based only on your post. Regardless, the package is extremely flexible, and you should be able to compute the quantities of interest if your read the documentation and tutorials. Here are some examples: average slopes, average contrasts, and average predictions by region:

avg_slopes(mod)
# 
#   Term  Estimate Std. Error      z Pr(>|z|)     S     2.5 %    97.5 %
#  emp    1.16e-03   4.66e-05 24.889   <0.001 451.8  1.07e-03  1.25e-03
#  pc     1.19e-05   1.05e-06 11.319   <0.001  96.3  9.81e-06  1.39e-05
#  pcap  -1.98e-06   2.77e-06 -0.712    0.476   1.1 -7.41e-06  3.46e-06
#  unemp -5.34e-03   9.89e-04 -5.393   <0.001  23.8 -7.27e-03 -3.40e-03
# 
# Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
# Type:  response

avg_comparisons(mod)
# 
#   Term Contrast  Estimate Std. Error      z Pr(>|z|)     S     2.5 %    97.5 %
#  emp         +1  1.16e-03   4.66e-05 24.890   <0.001 451.9  1.07e-03  1.25e-03
#  pc          +1  1.19e-05   1.05e-06 11.324   <0.001  96.3  9.81e-06  1.39e-05
#  pcap        +1 -1.98e-06   2.77e-06 -0.712    0.476   1.1 -7.41e-06  3.46e-06
#  unemp       +1 -5.34e-03   9.90e-04 -5.390   <0.001  23.8 -7.28e-03 -3.40e-03
# 
# Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
# Type:  response

avg_predictions(mod, by = "region")
# 
#  region Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
#       1     9.72    0.01720  565   <0.001 Inf  9.68   9.75
#       2    11.90    0.02465  483   <0.001 Inf 11.85  11.95
#       3    11.53    0.01413  816   <0.001 Inf 11.50  11.55
#       4    10.12    0.00863 1173   <0.001 Inf 10.10  10.14
#       5    10.65    0.00480 2217   <0.001 Inf 10.64  10.66
#       6    10.58    0.00502 2109   <0.001 Inf 10.57  10.59
#       7    10.97    0.01094 1003   <0.001 Inf 10.95  10.99
#       8     9.61    0.01582  607   <0.001 Inf  9.58   9.64
#       9    11.23    0.01623  692   <0.001 Inf 11.19  11.26
# 
# Columns: region, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
# Type:  response

You can feed the resulting object back to modelsummary to display the results in a regression table. Here is an example:

modelsummary(avg_slopes(mod))

	(1)
emp	0.001
	(0.000)
pc	0.000
	(0.000)
pcap	0.000
	(0.000)
unemp	-0.005
	(0.001)
Num.Obs.	816
R2	0.941
R2 Adj.	0.937
AIC	14092.6
BIC	14120.8
RMSE	0.04

And you can use the modelsummary::modelplot() function or one of the plotting functions from marginaleffects to visualize the results easily.

modelplot(avg_slopes(mod))

plot_predictions(mod, condition = c("pcap", "pc"))