processr是R中一个可以执行中介、调节以及有调节的中介模型的包。它和SPSS process的功能一样
看看 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
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
在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
使用这个模型语法,您可以使用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
使用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\)值
您可以用分析结果绘制统计图表。
statisticalDiagram(8,labels=labels,fit=semfit,whatLabel="est")
可以使用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
您可以用表格总结两种模型。
#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)
你可以用表格总结条件直接和间接影响。默认情况下,方程使用调节因子的sd均值。下表总结了调节因子为mean、mean + sd和mean - sd时的直接效应和间接效应。。
x=modmedSummary(semfit,mod="skeptic")
modmedSummaryTable(x)
Indirect Effect | Direct Effect | ||||
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 | |||||
。
你可以总结一下条件直接效应和间接效应。
conditionalEffectPlot(semfit,data=disaster,mod="skeptic")
library(processR)
library(lavaan)
fit=lm(justify ~ frame*skeptic,data=disaster)
labels=list(X="frame",W="skeptic",Y="justify")
pmacroModel(1,labels=labels)
statisticalDiagram(1,labels=labels)
#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 | ||||
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))
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]
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
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"))
pmacroModel(1,labels=labels,covar=covar)
statisticalDiagram(1,labels=labels,covar=covar)
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 | ||||
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))
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).
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'