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?
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"))