Regression in first differences with packages `plm` and `broom`

I'm estimating a linear regression in first-differences straight out of the textbook:

• James H. Stock and Mark W. Watson, Introduction to Econometrics, Pearson, 4th Edition.
• Christoph Hanck, Martin Arnold, Alexander Gerber, and Martin Schmelzer, Introduction to Econometrics with R, https://www.econometrics-with-r.org/10-rwpd.html
• Data from package AER, https://cran.r-project.org/web/packages/AER/AER.pdf

I'm estimating the model in first-differences with the plm() function from the plm package and extracting residuals with the augment() function from the broom package. I'm getting an error message and suspect I may not be using the "fd" option correctly and/or misusing augment(). A similar attempt for model="pooling" appears to work. Help appreciated!

library(AER)
data(Fatalities)

Fatalities$fatality <- Fatalities$fatal / Fatalities$pop * 10000

library(broom)

plm1 <- plm(fatality ~ beertax, data=Fatalities, model="pooling")
tidy(plm1)  # ok
augment(plm1)  # ok


plm2 <- plm(fatality ~ beertax, data=Fatalities, index=c("state", "year"), model="fd")
tidy(plm2)  # looks ok
# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept) -0.00314    0.0119   -0.263    0.792
2 beertax      0.0137     0.285     0.0480   0.962


augment(plm2)  # not ok

Error in `$<-.data.frame`(`*tmp*`, ".resid", value = c(`2` = 0.219840293582125,  : 
  replacement has 288 rows, data has 336
In addition: Warning message:
In get(.Generic)(e1, e2) :
  longer object length is not a multiple of shorter object length

EDIT: A WORKAROUND

So I suspect the problem has something to do with the fact that the model returned by plm and the residuals do not have the same number of rows:

length(row.names(plm.diff$model)) is 96

length(names(plm.diff$residuals)) is exactly half that many at 48.

Can someone tell me if the following is the correct way to get residuals and fitted values from the first-difference estimation?

data.frame(".rownames" = names(resid(plm.diff)), plm.diff$model) %>%
  left_join(data.frame(".rownames" = names(resid(plm.diff)), 
                       ".fitted" = fitted(plm.diff),
                       ".resid" = resid(plm.diff)
                       )) -> Fatalities.augmented
head(Fatalities.augmented)
  .rownames fatality beertax    .fitted    .resid
1         2   2.1284 1.53938  0.0394902  0.398097
2         4   2.4939 1.50144  0.1309701  0.103987
3         6   2.3841 0.65036 -0.1370856  0.415633

References:

• https://cran.r-project.org/web/packages/plm/plm.pdf
• https://cran.r-project.org/web/packages/broom/broom.pdf

Reference for Edit:

• Using broom::augment Panel data models