Is this the kind of data manipulation you are looking for?
DF <- data.frame(X = 1:25, Y = 1:25 * 1.2 + 0.3 + rnorm(25))
FIT <- lm(Y ~ X, data = DF)
DF$Resid <- residuals(FIT)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
DF <- mutate(DF, Resid_1 = lag(Resid, 1), Resid_2 = lag(Resid, 2))
head(DF)
#> X Y Resid Resid_1 Resid_2
#> 1 1 0.9432692 0.01602679 NA NA
#> 2 2 1.8582220 -0.30047839 0.01602679 NA
#> 3 3 4.0887476 0.69858930 -0.30047839 0.01602679
#> 4 4 4.9185789 0.29696255 0.69858930 -0.30047839
#> 5 5 5.4193480 -0.43372624 0.29696255 0.69858930
#> 6 6 5.0194083 -2.06512400 -0.43372624 0.29696255
Created on 2020-05-26 by the reprex package (v0.3.0)