rpanel-dataplmfixest

How to Plot the "time" effects results from PLM/ Fixest?


I am running a fixed effects model with a continuous variable (say parental wealth) on another continuous variable (children's wealth). I have a few control variables as well.

I want to plot the relationship between X and Y, over "Time". Is it possible to create a year-wise plot from the coefficients of the below equation? One way to do it would be manually collecting Time-related effects, and plotting them.

summary(plm(Y~X+a+b+as.factor(Time), data=df, index=c("ID","Time"),model = "within",effect = "individual"))

I want to examine how the relationship changes over time.


Solution

  • Not sure if you meant something like this where a separate OLS model is estimated per time period.

    library(plm)
    data("Grunfeld", package = "plm")
    
    form <- inv ~ value + capital
    
    pvcm.wi.mod <- pvcm(form, data = Grunfeld, model = "within", effect = "time")
    df <- pvcm.wi.mod$coefficients
    
    library(ggplot2)
    ggplot2::ggplot() +
    geom_point(aes(x = rownames(df), y = df$value),   color="red") + 
    geom_point(aes(x = rownames(df), y = df$capital), color="blue") + 
    ylab("Values") + xlab("Year")
    

    enter image description here

    Or maybe to plot the time effects of a two-way FE model (note that in the plot the years are numbered):

    fe.mod <- plm(form, data = Grunfeld, model = "within", effect = "twoways")
    fe.mod.fix <- fixef(fe.mod)
    print(fe.mod.fix)
    #>    1935    1936    1937    1938    1939    1940    1941    1942    1943    1944 
    #>  -86.90 -106.10 -127.59 -126.13 -156.37 -131.14 -105.70 -108.04 -129.88 -130.00 
    #>    1945    1946    1947    1948    1949    1950    1951    1952    1953    1954 
    #> -142.58 -118.07 -126.29 -130.62 -160.40 -162.80 -149.38 -151.53 -154.62 -180.43
    plot(fe.mod.fix, ylab = "time FE", xlab = "year")
    

    enter image description here

    Here fixef() gives the fixed effects, you could use fixef(., type = "dfirst", effect = "time") to get what the estimation via LSDV gives you.

    fe.mod.fix.dfirst <- fixef(fe.mod, effect = "time", type = "dfirst")
    print(fe.mod.fix.dfirst)
    #>    1936    1937    1938    1939    1940    1941    1942    1943    1944    1945 
    #> -19.197 -40.690 -39.226 -69.470 -44.235 -18.804 -21.140 -42.978 -43.099 -55.683 
    #>    1946    1947    1948    1949    1950    1951    1952    1953    1954 
    #> -31.169 -39.392 -43.717 -73.495 -75.896 -62.481 -64.632 -67.718 -93.526
    lsdv <- lm(inv ~ value + capital + factor(year) + factor(firm), data = Grunfeld)
    print(lsdv)
    #> 
    #> Call:
    #> lm(formula = inv ~ value + capital + factor(year) + factor(firm), 
    #>     data = Grunfeld)
    #> 
    #> Coefficients:
    #>      (Intercept)             value           capital  factor(year)1936  
    #>         -86.9002            0.1177            0.3579          -19.1974  
    #> factor(year)1937  factor(year)1938  factor(year)1939  factor(year)1940  
    #>         -40.6900          -39.2264          -69.4703          -44.2351  
    #> factor(year)1941  factor(year)1942  factor(year)1943  factor(year)1944  
    #>         -18.8045          -21.1398          -42.9776          -43.0988  
    #> factor(year)1945  factor(year)1946  factor(year)1947  factor(year)1948  
    #>         -55.6830          -31.1693          -39.3922          -43.7165  
    #> factor(year)1949  factor(year)1950  factor(year)1951  factor(year)1952  
    #>         -73.4951          -75.8961          -62.4809          -64.6323  
    #> factor(year)1953  factor(year)1954     factor(firm)2     factor(firm)3  
    #>         -67.7180          -93.5262          207.0542         -135.2308  
    #>    factor(firm)4     factor(firm)5     factor(firm)6     factor(firm)7  
    #>          95.3538           -5.4386          102.8886           51.4666  
    #>    factor(firm)8     factor(firm)9    factor(firm)10  
    #>          67.4905           30.2176          126.8371