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.

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

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

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