rplmmarginal-effectsmodelsummaryr-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:


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