1 processR 简单介绍

processr是R中一个可以执行中介、调节以及有调节的中介模型的包。它和SPSS process的功能一样

1.1 支持哪些模型

看看 processR 包都可以支持哪些模型?

library(processR)
## Warning: package 'processR' was built under R version 4.0.3
## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
sort(pmacro$no)
##  [1]  0.0  1.0  2.0  3.0  4.0  4.2  4.3  5.0  6.0  6.3  6.4  7.0  8.0  9.0 10.0
## [16] 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 25.0
## [31] 26.0 27.0 28.0 29.0 30.0 31.0 35.0 36.0 40.0 41.0 45.0 49.0 50.0 58.0 58.2
## [46] 59.0 60.0 61.0 62.0 63.0 64.0 65.0 66.0 67.0 74.0 75.0 76.0

1.2 模型方程进行分析

1.2.1 查看模型的概念图与统计图

Moderated Mediation (PROCESS macro model 8)

你可以很容易地画出概念图和统计图。例如,您可以为进程宏模型8绘制概念图。

pmacroModel(8)

您可以绘制这个模型的统计图。

statisticalDiagram(8)

这个包使用lavaan和中介包进行分析。lavaan包是可定制的,如果您的测量模型需要,它还可以建模潜在变量。

用到的数据:

head(disaster)
##   id frame donate justify skeptic
## 1  1     1    5.6    2.95     1.8
## 2  2     1    4.2    2.85     5.2
## 3  3     1    4.2    3.00     3.2
## 4  4     1    4.6    3.30     1.0
## 5  5     1    3.0    5.00     7.6
## 6  6     0    5.0    3.20     4.2

1.2.2 列出分析方程

在processR中,可以看到每个模型对应的方程,此处使用processR自带的数据集disaster为例,该数据集包含5个变量:

  • id:编号

  • frame:实验条件,0为自然灾害,1为气候变化导致的灾害

  • donate:积极的捐助意愿

  • justify:消极理由negative justifications

  • skeptic:气候变化怀疑论

数据来源:

Chapman, D. A., & Little, B. (2016). Climate change and disasters: How framing affects justifications for giving or withholding aid to disaster victims. Social Psychological and Personality Science, 7, 13-20.

查看包含变量名称的模型图,这里只需开启pmacroModel()函数中的labels参数,并将之前储存变量名称的列表(list)赋值给它

labels=list(X="frame",M="justify",Y="donate",W="skeptic")
pmacroModel(8,labels=labels)

该模型包含一个调节变量skeptic,它调节的位置在第一阶段,也就是frame影响justify这个阶段,该阶段的系数称之为“a”;而且还对frame影响donate有调节,该阶段的系数称之为“c”。(无法判断系数名称时,可调出模型统计图查看)

moderator=list(name="skeptic",site=list(c("a","c")))

modmedEquation()用来生成有调节中介模型的方程

model=tripleEquation(X="frame",M="justify",Y="donate",moderator=moderator)
cat(model)
## justify~a1*frame+a2*skeptic+a3*frame:skeptic
## donate~c1*frame+c2*skeptic+c3*frame:skeptic+b*justify
## skeptic ~ skeptic.mean*1
## skeptic ~~ skeptic.var*skeptic
## CE.XonM :=a1+a3*skeptic.mean
## indirect :=(a1+a3*skeptic.mean)*(b)
## index.mod.med :=a3*b
## direct :=c1+c3*skeptic.mean
## total := direct + indirect
## prop.mediated := indirect / total
## CE.XonM.below :=a1+a3*(skeptic.mean-sqrt(skeptic.var))
## indirect.below :=(a1+a3*(skeptic.mean-sqrt(skeptic.var)))*(b)
## CE.XonM.above :=a1+a3*(skeptic.mean+sqrt(skeptic.var))
## indirect.above :=(a1+a3*(skeptic.mean+sqrt(skeptic.var)))*(b)
## direct.below:=c1+c3*(skeptic.mean-sqrt(skeptic.var))
## direct.above:=c1+c3*(skeptic.mean+sqrt(skeptic.var))
## total.below := direct.below + indirect.below
## total.above := direct.above + indirect.above
## prop.mediated.below := indirect.below / total.below
## prop.mediated.above := indirect.above / total.above

使用cat()列出方程

cat(model)
## justify~a1*frame+a2*skeptic+a3*frame:skeptic
## donate~c1*frame+c2*skeptic+c3*frame:skeptic+b*justify
## skeptic ~ skeptic.mean*1
## skeptic ~~ skeptic.var*skeptic
## CE.XonM :=a1+a3*skeptic.mean
## indirect :=(a1+a3*skeptic.mean)*(b)
## index.mod.med :=a3*b
## direct :=c1+c3*skeptic.mean
## total := direct + indirect
## prop.mediated := indirect / total
## CE.XonM.below :=a1+a3*(skeptic.mean-sqrt(skeptic.var))
## indirect.below :=(a1+a3*(skeptic.mean-sqrt(skeptic.var)))*(b)
## CE.XonM.above :=a1+a3*(skeptic.mean+sqrt(skeptic.var))
## indirect.above :=(a1+a3*(skeptic.mean+sqrt(skeptic.var)))*(b)
## direct.below:=c1+c3*(skeptic.mean-sqrt(skeptic.var))
## direct.above:=c1+c3*(skeptic.mean+sqrt(skeptic.var))
## total.below := direct.below + indirect.below
## total.above := direct.above + indirect.above
## prop.mediated.below := indirect.below / total.below
## prop.mediated.above := indirect.above / total.above

1.2.3 模型分析

使用这个模型语法,您可以使用lavaan包的sem()函数分析被调节的中介。

library(lavaan)
semfit=sem(model=model,data=disaster)
summary(semfit)
## lavaan 0.6-7 ended normally after 33 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         18
##                                                       
##   Number of observations                           211
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               136.428
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   justify ~                                           
##     frame     (a1)   -0.562    0.175   -3.211    0.001
##     skeptic   (a2)    0.105    0.027    3.844    0.000
##     frm:skptc (a3)    0.201    0.040    5.077    0.000
##   donate ~                                            
##     frame     (c1)    0.160    0.216    0.741    0.459
##     skeptic   (c2)   -0.043    0.034   -1.248    0.212
##     frm:skptc (c3)    0.015    0.051    0.295    0.768
##     justify    (b)   -0.923    0.083  -11.113    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   frame ~~                                            
##     frame:skeptic     0.854    0.096    8.890    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     skeptic (skp.)    3.378    0.140   24.196    0.000
##    .justify           2.452    0.120   20.416    0.000
##    .donate            7.291    0.250   29.189    0.000
##     frame             0.479    0.034   13.919    0.000
##     frm:skp           1.637    0.152   10.771    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     skeptic (skp.)    4.113    0.400   10.271    0.000
##    .justify           0.648    0.063   10.271    0.000
##    .donate            0.943    0.092   10.271    0.000
##     frame             0.250    0.024   10.271    0.000
##     frm:skp           4.877    0.475   10.271    0.000
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     CE.XonM           0.117    0.114    1.023    0.306
##     indirect         -0.108    0.106   -1.019    0.308
##     index.mod.med    -0.186    0.040   -4.618    0.000
##     direct            0.211    0.134    1.570    0.116
##     total             0.103    0.170    0.603    0.546
##     prop.mediated    -1.053    2.516   -0.418    0.676
##     CE.XonM.below    -0.291    0.142   -2.045    0.041
##     indirect.below    0.268    0.133    2.011    0.044
##     CE.XonM.above     0.525    0.140    3.743    0.000
##     indirect.above   -0.485    0.137   -3.547    0.000
##     direct.below      0.180    0.168    1.073    0.283
##     direct.above      0.241    0.170    1.420    0.156
##     total.below       0.449    0.212    2.121    0.034
##     total.above      -0.244    0.209   -1.167    0.243
##     prop.medtd.blw    0.598    0.257    2.326    0.020
##     prop.meditd.bv    1.990    1.449    1.373    0.170

1.2.4 整理结果

使用estimatesTable()将结果整理成较为简洁的格式。

library(lavaan)
estimatesTable(semfit)
##   Variables    Predictors label     B   SE      z       p     β
## 1   justify         frame    a1 -0.56 0.18  -3.21   0.001 -0.32
## 2   justify       skeptic    a2  0.11 0.03   3.84 < 0.001  0.24
## 3   justify frame:skeptic    a3  0.20 0.04   5.08 < 0.001  0.50
## 4    donate         frame    c1  0.16 0.22   0.74   0.459  0.06
## 5    donate       skeptic    c2 -0.04 0.03  -1.25   0.212 -0.07
## 6    donate frame:skeptic    c3  0.01 0.05   0.29   0.768  0.03
## 7    donate       justify     b -0.92 0.08 -11.11 < 0.001 -0.64

还可以使用estimatesTable2()将结果整理成表格。

estimatesTable2(semfit)

Variables

Predictors

label

B

SE

z

p

β

justify

frame

a1

-0.562

0.175

-3.211

0.001

-0.319

justify

skeptic

a2

0.105

0.027

3.844

< 0.001

0.242

justify

frame:skeptic

a3

0.201

0.040

5.077

< 0.001

0.504

donate

frame

c1

0.160

0.216

0.741

0.459

0.063

donate

skeptic

c2

-0.043

0.034

-1.248

0.212

-0.068

donate

frame:skeptic

c3

0.015

0.051

0.295

0.768

0.026

donate

justify

b

-0.923

0.083

-11.113

< 0.001

-0.636

开启vanilla参数,便可生成三线表。

estimatesTable2(semfit,vanilla = FALSE)

Variables

Predictors

label

B

SE

z

p

β

justify

frame

a1

-0.562

0.175

-3.211

0.001

-0.319

justify

skeptic

a2

0.105

0.027

3.844

< 0.001

0.242

justify

frame:skeptic

a3

0.201

0.040

5.077

< 0.001

0.504

donate

frame

c1

0.160

0.216

0.741

0.459

0.063

donate

skeptic

c2

-0.043

0.034

-1.248

0.212

-0.068

donate

frame:skeptic

c3

0.015

0.051

0.295

0.768

0.026

donate

justify

b

-0.923

0.083

-11.113

< 0.001

-0.636

\(B\)代表回归系数、\(SEM\)代表标准误、\(\gamma\)代表卡方值、\(P\)代表\(p\)

1.2.5 绘制包含分析结果的统计图

您可以用分析结果绘制统计图表。

statisticalDiagram(8,labels=labels,fit=semfit,whatLabel="est")

1.3 用简单回归模型进行分析

可以使用lm函数分析这个模型。你可以为调节变量和因变量建立回归方程。

equations=regEquation(X="frame",M="justify",Y="donate",moderator=moderator)
cat(equations)
## justify ~ frame
## donate ~  frame+justify+frame*skeptic

有了这个方程,您可以执行线性回归。首先,您可以使用调节因子作为因变量。

equations
## [1] "justify ~ frame\ndonate ~  frame+justify+frame*skeptic"
eq=unlist(strsplit(equations,"\n"))

eq
## [1] "justify ~ frame"                      
## [2] "donate ~  frame+justify+frame*skeptic"
fit=lapply(1:2,function(i) {
    lm(as.formula(eq[i]),data=disaster)
})



summary(fit[[1]])
## 
## Call:
## lm(formula = as.formula(eq[i]), data = disaster)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8024 -0.6867 -0.0024  0.5976  3.8633 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.80236    0.08862  31.623   <2e-16 ***
## frame        0.13437    0.12809   1.049    0.295    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9294 on 209 degrees of freedom
## Multiple R-squared:  0.005238,   Adjusted R-squared:  0.0004783 
## F-statistic:   1.1 on 1 and 209 DF,  p-value: 0.2954

第二个回归模型使用因变量。

summary(fit[[2]])
## 
## Call:
## lm(formula = as.formula(eq[i]), data = disaster)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4326 -0.5205  0.1216  0.6898  3.1508 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7.29147    0.27368  26.642   <2e-16 ***
## frame          0.16032    0.26766   0.599    0.550    
## justify       -0.92269    0.08403 -10.981   <2e-16 ***
## skeptic       -0.04257    0.04693  -0.907    0.365    
## frame:skeptic  0.01492    0.06892   0.217    0.829    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9828 on 206 degrees of freedom
## Multiple R-squared:  0.454,  Adjusted R-squared:  0.4434 
## F-statistic: 42.82 on 4 and 206 DF,  p-value: < 2.2e-16

1.4 表格总结两种模型

您可以用表格总结两种模型。

#x=modelsSummary(fit,labels=labels)
#modelsSummaryTable(x)
#x
x=modelsSummary(fit,labels=labels)
modelsSummaryTable(x)

Consequent

justify(M)

donate(Y)

Antecedent

Coef

SE

t

p

Coef

SE

t

p

frame(X)

a

0.134

0.128

1.049

.295

c'1

0.160

0.268

0.599

.550

justify(M)

b

-0.923

0.084

-10.981

<.001

skeptic(W)

c'2

-0.043

0.047

-0.907

.365

frame:skeptic(X:W)

c'3

0.015

0.069

0.217

.829

Constant

iM

2.802

0.089

31.623

<.001

iY

7.291

0.274

26.642

<.001

Observations

211

211

R2

0.005

0.454

Adjusted R2

0.000

0.443

Residual SE

0.929 ( df = 209)

0.983 ( df = 206)

F statistic

F(1,209) = 1.100, p = .295

F(4,206) = 42.816, p < .001

是否可以输出 LaTex 格式的文件?

#library(stargazer)
#stargazer(x, align = TRUE)

1.5 条件直接效应和间接效应

你可以用表格总结条件直接和间接影响。默认情况下,方程使用调节因子的sd均值。下表总结了调节因子为mean、mean + sd和mean - sd时的直接效应和间接效应。。

x=modmedSummary(semfit,mod="skeptic")
modmedSummaryTable(x)

Indirect Effect
(a1+a3*W)*(b)

Direct Effect
c1+c3*W

skeptic(W)

estimate

95% Bootstrap CI

estimate

95% Bootstrap CI

1.350

0.268

(0.007 to 0.530)

0.180

(-0.149 to 0.510)

3.378

-0.108

(-0.316 to 0.100)

0.211

(-0.052 to 0.474)

5.406

-0.485

(-0.752 to -0.217)

0.241

(-0.092 to 0.574)

boot.ci.type = perc

1.6 有条件的直接和间接影响的图

你可以总结一下条件直接效应和间接效应。

conditionalEffectPlot(semfit,data=disaster,mod="skeptic")

2 调解效应的可视化

2.1 引入相关的包

library(processR)
library(lavaan)
fit=lm(justify ~ frame*skeptic,data=disaster)

2.2 概念模型

labels=list(X="frame",W="skeptic",Y="justify")
pmacroModel(1,labels=labels)

2.3 统计模型

statisticalDiagram(1,labels=labels)

2.4 模型统计信息

#modelsSummaryTable(list(fit),labels=labels)
x=modelsSummary(fit,labels=labels)
x
## =============================================================== 
##                                        Consequent               
##                          -------------------------------------- 
##                                       justify(Y)              
##                          -------------------------------------  
##        Antecedent               Coef     SE      t       p     
## --------------------------------------------------------------- 
##         frame(X)          c1  -0.562   0.218  -2.581    .011  
##        skeptic(W)         c2   0.105   0.038   2.756    .006  
##    frame:skeptic(X:W)     c3   0.201   0.055   3.640   <.001  
##         Constant          iY   2.452   0.149  16.449   <.001  
## --------------------------------------------------------------- 
##       Observations                        211                 
##            R2                            0.246                
##        Adjusted R2                       0.235                
##        Residual SE                 0.813 ( df = 207)          
##        F statistic            F(3,207) = 22.543, p < .001     
## ===============================================================
modelsSummaryTable(list(fit),labels=labels)

Consequent

justify(Y)

Antecedent

Coef

SE

t

p

frame(X)

c1

-0.562

0.218

-2.581

.011

skeptic(W)

c2

0.105

0.038

2.756

.006

frame:skeptic(X:W)

c3

0.201

0.055

3.640

<.001

Constant

iY

2.452

0.149

16.449

<.001

Observations

211

R2

0.246

Adjusted R2

0.235

Residual SE

0.813 ( df = 207)

F statistic

F(3,207) = 22.543, p < .001

2.5 调解效果可视化

condPlot(fit,rangemode=2,xpos=0.7,labels=c("Climate change(X=1)","Natural causes(X=0)"))
## Joining, by = "frame"
## Joining, by = "x"

condPlot(fit,mode=2,xpos=0.6)

condPlot(fit,mode=3,rangemode=2,xpos=0.5,ypos=c(0,2))

2.6 Johnson-Neyman Plot

jnPlot(fit,plot=FALSE)
## JOHNSON-NEYMAN INTERVAL 
## 
## When skeptic is OUTSIDE the interval [1.171, 3.934], the slope of frame is
## p < .05.
## 
## Note: The range of observed values of skeptic is [1.000, 9.000]

2.7 统计数据

moderator=list(name="skeptic",site=list("c"))
model=tripleEquation(labels=labels,moderator=moderator)
cat(model)
## justify~c1*frame+c2*skeptic+c3*frame:skeptic
## skeptic ~ skeptic.mean*1
## skeptic ~~ skeptic.var*skeptic
## direct :=c1+c3*skeptic.mean
## direct.below:=c1+c3*(skeptic.mean-sqrt(skeptic.var))
## direct.above:=c1+c3*(skeptic.mean+sqrt(skeptic.var))
semfit=sem(model=model,data=disaster)
## Warning in lavaan::lavaan(model = model, data = disaster, model.type = "sem", :
## lavaan WARNING: syntax contains parameters involving exogenous covariates;
## switching to fixed.x = FALSE
summary(semfit)
## lavaan 0.6-7 ended normally after 24 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         12
##                                                       
##   Number of observations                           211
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               136.428
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   justify ~                                           
##     frame     (c1)   -0.562    0.175   -3.211    0.001
##     skeptic   (c2)    0.105    0.027    3.844    0.000
##     frm:skptc (c3)    0.201    0.040    5.077    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   frame ~~                                            
##     frame:skeptic     0.854    0.096    8.890    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     skeptic (skp.)    3.378    0.140   24.196    0.000
##    .justify           2.452    0.120   20.416    0.000
##     frame             0.479    0.034   13.919    0.000
##     frm:skp           1.637    0.152   10.771    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     skeptic (skp.)    4.113    0.400   10.271    0.000
##    .justify           0.648    0.063   10.271    0.000
##     frame             0.250    0.024   10.271    0.000
##     frm:skp           4.877    0.475   10.271    0.000
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     direct            0.117    0.114    1.023    0.306
##     direct.below     -0.291    0.142   -2.045    0.041
##     direct.above      0.525    0.140    3.743    0.000

3 第二个例子

labels=list(X="negemot",Y="govact",W="age",C1="posemot",C2="ideology",C3="sex")
moderator=list(name="age",site=list("c"))
covar=list(name=c("posemot","ideology","sex"),site=list("Y","Y","Y"))

3.1 概念模型

pmacroModel(1,labels=labels,covar=covar)

3.2 统计模型

statisticalDiagram(1,labels=labels,covar=covar)

3.3 模型统计信息

fit=lm(govact~negemot*age+ideology+sex,data=glbwarm)
modelsSummaryTable(list(fit),labels=labels)

Consequent

govact(Y)

Antecedent

Coef

SE

t

p

negemot(X)

c1

0.114

0.082

1.388

.165

age(W)

c2

-0.024

0.006

-4.044

<.001

ideology(C2)

g2

-0.211

0.027

-7.883

<.001

sex(C3)

g3

-0.016

0.076

-0.212

.832

negemot:age(X:W)

c3

0.006

0.002

4.145

<.001

Constant

iY

5.132

0.334

15.371

<.001

Observations

815

R2

0.400

Adjusted R2

0.397

Residual SE

1.057 ( df = 809)

F statistic

F(5,809) = 108.033, p < .001

3.4 调解效果可视化

condPlot(fit,xmode=2,pred.values=c(30,70),hjust=c(-0.1,-0.1,1.1))
## Joining, by = "age"
## Joining, by = "x"

condPlot(fit,xmode=2,mode=2,pred.values=c(30,50,70),xpos=0.2)

condPlot(fit,xmode=2,modx.values=2:5,mode=3,xpos=0.6,hjust=c(-0.1,-0.1,-0.1,1.1))

3.5 Johnson-Neyman Plot

jnPlot(fit,plot=FALSE)
## JOHNSON-NEYMAN INTERVAL 
## 
## When age is OUTSIDE the interval [-80.462, 5.112], the slope of negemot is
## p < .05.
## 
## Note: The range of observed values of age is [17.000, 87.000]
## Warning: Removed 1 rows containing missing values (geom_text).

3.6 Make lavaan syntax

model=tripleEquation(labels=labels,moderator=moderator,covar=covar)
cat(model)
## govact ~ c1*negemot+c2*age+c3*negemot:age + g1*posemot + g2*ideology + g3*sex
## age ~ age.mean*1
## age ~~ age.var*age
## direct :=c1+c3*age.mean
## direct.below:=c1+c3*(age.mean-sqrt(age.var))
## direct.above:=c1+c3*(age.mean+sqrt(age.var))
#用`broom::augment()`和ggplot2做出类似的残差图
#假定数据是
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.0.3
## -- Attaching packages ------------------------------------------------------------------------- tidyverse 1.3.0 --
## √ ggplot2 3.3.2     √ purrr   0.3.4
## √ tibble  3.0.3     √ dplyr   1.0.2
## √ tidyr   1.1.0     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 4.0.3
## Warning: package 'dplyr' was built under R version 4.0.3
## Warning: package 'stringr' was built under R version 4.0.3
## -- Conflicts ---------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
df <- tibble(
  x = runif(30, 2, 10),
  y = -2*x + rnorm(30, 0, 5)
  )


fitted_lm <- lm(y ~ x, data = df)



fitted_lm %>% 
  broom::augment() %>% 
  select(x, y, predicted = .fitted, residuals = .resid) %>% 

  ggplot(aes(x = x, y = y)) +
  geom_smooth(method = "lm", se = FALSE, color = "gray50") +
  geom_segment(aes(xend= x, yend = predicted), alpha = 0.2) +
  geom_point(aes(size = abs(residuals), color = abs(residuals)))
## `geom_smooth()` using formula 'y ~ x'