library(regsim)

data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)

model <- lm(mpg ~ wt + disp + cyl, mtcars)
summary(model)

## 
## Call:
## lm(formula = mpg ~ wt + disp + cyl, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5965 -1.2361 -0.4855  1.4740  5.8043 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.041673   1.963039  17.341 3.66e-16 ***
## wt          -3.306751   1.105083  -2.992  0.00586 ** 
## disp         0.001715   0.013481   0.127  0.89972    
## cyl6        -4.305559   1.464760  -2.939  0.00666 ** 
## cyl8        -6.322786   2.598416  -2.433  0.02186 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.603 on 27 degrees of freedom
## Multiple R-squared:  0.8375, Adjusted R-squared:  0.8135 
## F-statistic:  34.8 on 4 and 27 DF,  p-value: 2.726e-10

x <- list(wt = seq(1, 5), cyl = levels(mtcars$cyl))

sim <- regsim(model, x)
summary(sim)

##      Profile 1
## wt   1        
## cyl  4        
## disp 230.7219 
## 
##        2.5%      50%    97.5%
## p1 25.56391 31.10842 36.89383
## -----------------------------
## 
##      Profile 2
## wt   2        
## cyl  4        
## disp 230.7219 
## 
##        2.5%      50%    97.5%
## p2 23.99886 27.78313 31.90591
## -----------------------------
## 
##      Profile 3
## wt   3        
## cyl  4        
## disp 230.7219 
## 
##        2.5%      50%   97.5%
## p3 21.82099 24.51221 27.3655
## ----------------------------
## 
##      Profile 4
## wt   4        
## cyl  4        
## disp 230.7219 
## 
##        2.5%      50%    97.5%
## p4 18.22641 21.28332 24.18166
## -----------------------------
## 
##      Profile 5
## wt   5        
## cyl  4        
## disp 230.7219 
## 
##       2.5%      50%    97.5%
## p5 13.5847 17.92686 22.46591
## ----------------------------
## 
##      Profile 6
## wt   1        
## cyl  6        
## disp 230.7219 
## 
##        2.5%      50%   97.5%
## p6 20.67071 26.83756 32.2676
## ----------------------------
## 
##      Profile 7
## wt   2        
## cyl  6        
## disp 230.7219 
## 
##        2.5%      50%    97.5%
## p7 19.31662 23.51637 27.17672
## -----------------------------
## 
##      Profile 8
## wt   3        
## cyl  6        
## disp 230.7219 
## 
##        2.5%      50%    97.5%
## p8 17.89272 20.16748 22.53182
## -----------------------------
## 
##      Profile 9
## wt   4        
## cyl  6        
## disp 230.7219 
## 
##        2.5%      50%    97.5%
## p9 14.54709 16.90482 19.14244
## -----------------------------
## 
##      Profile 10
## wt   5         
## cyl  6         
## disp 230.7219  
## 
##         2.5%      50%   97.5%
## p10 9.776385 13.58769 17.3142
## -----------------------------
## 
##      Profile 11
## wt   1         
## cyl  8         
## disp 230.7219  
## 
##         2.5%      50%    97.5%
## p11 19.92288 24.85967 29.66922
## ------------------------------
## 
##      Profile 12
## wt   2         
## cyl  8         
## disp 230.7219  
## 
##         2.5%     50%    97.5%
## p12 18.24777 21.4852 24.70022
## -----------------------------
## 
##      Profile 13
## wt   3         
## cyl  8         
## disp 230.7219  
## 
##         2.5%      50%    97.5%
## p13 15.58305 18.20112 20.76528
## ------------------------------
## 
##      Profile 14
## wt   4         
## cyl  8         
## disp 230.7219  
## 
##         2.5%      50%    97.5%
## p14 11.47673 14.93413 18.60458
## ------------------------------
## 
##      Profile 15
## wt   5         
## cyl  8         
## disp 230.7219  
## 
##        2.5%      50%  97.5%
## p15 6.39727 11.63853 17.192
## ---------------------------

plot(sim, ~ wt + cyl, lines.col = "Set1", lines.lwd = 1.5)