多层模型及其交互效应入门教程

xxxspy 2023-09-02 15:24:07
Categories: Tags:

本教程介绍了如何使用多层模型来分析嵌套数据,以及调节效应在跨层结构中的分析方法。
我们的案例数据是日记数据(嵌套在个人中的重复事件),意思是每个人都会重复测量多次,但也适用于其他类型的嵌套数据。

教程目标

教程用到的库(请自行安装)

1
2
3
4
5
6
7
library(ggplot2)       # for 数据可视化
library(lme4) # for 混合效应模型(多层模型)的构建
library(lmerTest) # for 计算P值
library(psych) # for 描述统计
library(plyr) # for 数据整理
library(effects) # for 检验调节效应
library(interactions) # for 调节效应可视化

数据介绍

我们的数据来自于一个重复测量的调查研究,虽然我们没有找到关于数据的详细介绍,但是我们至少清楚我们用到的这几个变量是什么意义。
数据包含两个文件”daily-data.csv”和”person-data.csv”, 这两个文件的意义就是文件名所揭示的, 第一个文件是重复测量的数据,
被试的每天汇报的数据, 第二个文件是个体的特征数据,它只采集了一次。

我们用到的变量:

下面我们加载两个数据,取出用到的变量:

1
2
3
4
5
6
pdata = read.csv("person-data.csv")
ddata = read.csv("daily-data.csv")
pdata <- pdata[ ,c("id","bfi_n")]
ddata <- ddata[ ,c("id","day","negaff","pss")]
head(pdata)

输出(html):
A data.frame: 6 × 2
idbfi_n
<int><dbl>
11012.0
21022.0
31032.5
41042.5
51053.5
61061.5

我们有必要介绍一下日测数据, 变量day记录的是第几日, 每一行数据不是被试样本, 是被试每天的数据,你可以看下面的数据,
这种格式的数据叫做长格式, 用于存储重测数据的常用格式。

1
head(ddata)
输出(html):
A data.frame: 6 × 4
iddaynegaffpss
<int><int><dbl><dbl>
110103.02.50
210112.32.75
310121.03.50
410131.33.00
510141.12.75
610151.02.75

数据整理

感知压力的反向编码

pss是感知压力, 但是数据中是分数越大压力越小, 所以我们反向编码一下,
使得分数越大压力越大。

1
2
3
ddata$stress <- 4 - ddata$pss

psych::describe(ddata$stress)
输出(html):
A psych: 1 × 13
varsnmeansdmediantrimmedmadminmaxrangeskewkurtosisse
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
X1114451.3855250.68433771.251.363440.74130440.35492760.12663230.01800266

所有的变量都建议检查一下它的分布,我们可以使用直方图来看一下:

1
2
3
ggplot(data=ddata, aes(x=stress)) +
geom_histogram(fill="white", color="black",bins=19) +
labs(x = "压力")
输出(stream):
Warning message: "Removed 13 rows containing non-finite values (`stat_bin()`)."

整理变量,划分为个体间和个体内变量

现在我们分清楚个体间变量和个体内变量, 个体间变量就是不随时间变化的,同一个个体相同的值, 不同的个体不同的值,这些变量可以称为人的特质;
个体内变量是随时间变化的, 同一个被试随着时间的变化量,这种变量可以称为状态, 因为状态是随时间变化的。

我们先要要处理的是压力(stress)这个变量, 显然它是状态变量, 但是每个人的状态都是围绕的自己的均值变化的,
所以我们可以从状态中计算个体均值, 代表他的特质, 所以新生成的变量 stress_trait 就是压力特质, 不随时间变化的变量。 我们可以这样计算这个变量:

1
2
3
4
5
6
7
# 计算个体均值
# negaff_trait虽然不是模型中的变量, 但是我们后期可能会用到它,所以也在这里一并计算了
personmeans <- ddply(ddata, "id", summarize,
stress_trait = mean(stress, na.rm=TRUE),
negaff_trait = mean(negaff, na.rm=TRUE))
head(personmeans)
describe(personmeans)
输出(html):
A data.frame: 6 × 3
idstress_traitnegaff_trait
<int><dbl><dbl>
11011.062501.500000
21020.781252.218750
31031.250002.416667
41041.812501.550000
51051.750002.612500
61061.125002.071875
输出(html):
A psych: 3 × 13
varsnmeansdmediantrimmedmadminmaxrangeskewkurtosisse
<int><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
id1190318.294737130.4413245321.50000318.993421151.2252000101.0000532.0000431.000-0.04393582-1.09453989.46320829
stress_trait2190 1.395533 0.4788391 1.40625 1.394792 0.5096437 0.1875 2.5625 2.375-0.04026841-0.23375930.03473864
negaff_trait3190 2.478309 0.7335710 2.41250 2.429841 0.7227675 1.1125 5.0875 3.975 0.67593800 0.45059050.05321883
1
2
3
4
5
6
7
# 将新生成的个体特质变量合并到pdata数据框, pdata 是被试样本
pdata <- merge(pdata, personmeans, by="id")
# 后面我们会分析调节效应, 对变量中心化是必然的步骤,这一步就是变量中心化
pdata$bfi_n_c <- scale(pdata$bfi_n,center=TRUE,scale=FALSE)
pdata$stress_trait_c <- scale(pdata$stress_trait,center=TRUE,scale=FALSE)
# 个体数据的描述性统计
describe(pdata)
输出(html):
A psych: 6 × 13
varsnmeansdmediantrimmedmadminmaxrangeskewkurtosisse
<int><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
id11903.182947e+02130.4413245321.50000000 3.189934e+02151.2252000101.000000532.000000431.000-0.04393582-1.09453989.46320829
bfi_n21902.981579e+00 0.9558661 3.00000000 2.996711e+00 1.4826000 1.000000 5.000000 4.000-0.09238813-0.81730500.06934582
stress_trait31901.395533e+00 0.4788391 1.40625000 1.394792e+00 0.5096437 0.187500 2.562500 2.375-0.04026841-0.23375930.03473864
negaff_trait41902.478309e+00 0.7335710 2.41250000 2.429841e+00 0.7227675 1.112500 5.087500 3.975 0.67593800 0.45059050.05321883
bfi_n_c51901.683047e-16 0.9558661 0.01842105 1.513158e-02 1.4826000 -1.981579 2.018421 4.000-0.09238813-0.81730500.06934582
stress_trait_c61901.927149e-17 0.4788391 0.01071742-7.409148e-04 0.5096437 -1.208033 1.166967 2.375-0.04026841-0.23375930.03473864

将个体数据合并到重复测量的日测数据

1
2
3
4
5
6
7
8
9

ddata_long <- merge(ddata,pdata,by="id")

# 计算状态变量,即减去个体均值
ddata_long$stress_state <- ddata_long$stress - ddata_long$stress_trait
ddata_long$negaff_state <- ddata_long$negaff - ddata_long$negaff_trait

# 查看数据
head(ddata_long)
输出(html):
A data.frame: 6 × 12
iddaynegaffpssstressbfi_nstress_traitnegaff_traitbfi_n_cstress_trait_cstress_statenegaff_state
<int><int><dbl><dbl><dbl><dbl><dbl><dbl><dbl[,1]><dbl[,1]><dbl><dbl>
110103.02.501.5021.06251.5-0.9815789-0.3330326 0.4375 1.5
210112.32.751.2521.06251.5-0.9815789-0.3330326 0.1875 0.8
310121.03.500.5021.06251.5-0.9815789-0.3330326-0.5625-0.5
410131.33.001.0021.06251.5-0.9815789-0.3330326-0.0625-0.2
510141.12.751.2521.06251.5-0.9815789-0.3330326 0.1875-0.4
610151.02.751.2521.06251.5-0.9815789-0.3330326 0.1875-0.5
1

我们取出前25个被试的数据,通过绘制每个被试的日次数据, 我们可以大概了解到被试之间的差异有多大,
比如看下面的图, 每个被试的拟合回归线的斜率都有很大差异, 着预示着压力和消极情绪的关系是因人而异的。

1
2
3
4
5
6
7
8
9
#faceted plot
ggplot(data=ddata_long[which(ddata_long$id <= 125),], aes(x=stress_state,y=negaff)) +
geom_point() +
stat_smooth(method="lm", fullrange=TRUE) +
xlab("Stress State") + ylab("Negative Affect (Continuous)") +
facet_wrap( ~ id) +
theme(axis.title=element_text(size=16),
axis.text=element_text(size=14),
strip.text=element_text(size=14))
输出(stream):
`geom_smooth()` using formula = 'y ~ x' Warning message: "Removed 4 rows containing non-finite values (`stat_smooth()`)." Warning message: "Removed 4 rows containing missing values (`geom_point()`)."

混合效应模型

我们使用 “lme4” 包来拟合混合效应模型,以及一些辅助的R包: “lmerTest” 提供了用于获取 参数检验的 p-vlaues 的工具;
“effects”包提供了用于计算和绘制基于模型的预测的工具;”interactions” 提供了绘制和探测交互效应的工具。

lme4 提供了函数 lmer , 它用于拟合多层数据模型,或者是混合效应模型。 它的第一个参数是data, 输入你的原始数据,
na.action 参数用于指定对缺失值的处理方法。

零模型

模型里面没有自变量, 或者没有我们关心的变量, 所以这种模型叫0模型。
零模型可以被称为 unconditional means model , 它用于考察因变量的方差有多少来自被试内, 有多少来自被试间。

1
2
3
4
5
#unconditional means model
model0_fit <- lmer(formula = negaff ~ 1 + (1|id),
data=ddata_long,
na.action=na.exclude)
summary(model0_fit)
输出(plain):
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: negaff ~ 1 + (1 | id)
Data: ddata_long

REML criterion at convergence: 3833.5

Scaled residuals:
Min 1Q Median 3Q Max
-3.8739 -0.6123 -0.1608 0.4658 3.9394

Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.4270 0.6535
Residual 0.6627 0.8141
Number of obs: 1441, groups: id, 190

Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.46368 0.05229 185.80793 47.12 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1
2
3

# 我们使用 VarCorr 函数提取方差 , id这一行对应的是被试内方差, Residual对应的是残差方差(不可以被被试差异所解释的方差
VarCorr(model0_fit)
输出(plain):
Groups Name Std.Dev.
id (Intercept) 0.65347
Residual 0.81408

我们来计算组间相关性(ICC), 它指的是被试内方差占总方差的比率, 如公式:

$ ICC_{between} = \frac{\sigma^{2}_{u0}}{\sigma^{2}_{u0} + \sigma^{2}_{e}} $

我们首先取药提取得到随机效应的方差, 我们将结果保存到数据框中, 方便提取数据, 毕竟数据框是我们在R中最熟悉的数据给格式。

1
2
randEffs <- as.data.frame(VarCorr(model0_fit))
randEffs
输出(html):
A data.frame: 2 × 5
grpvar1var2vcovsdcor
<chr><chr><chr><dbl><dbl>
id (Intercept)NA0.42702940.6534749
ResidualNA NA0.66272600.8140798

根据上面的公式, 可以计算得到ICC:

1
2
3
4
5
u0v = randEffs[1,4]
ev = randEffs[2,4]

icc <- u0v / (u0v+ev)
icc
输出(html):
0.391858025596418

根据无条件均值模型计算ICC,结果表明,在负面情绪的总方差中,39.19%归因于人与人之间的差异,60.81%归因于人内差异。
这意味着使用随时间变化的变量作为预测变量时, 存在很大一部分人内方差, 着意味着我们构建多层模型、混合效应模型是非常有必要的。
专家建议ICC超过0.05(5%)时就应该考虑使用混合效应模型。

混合效应模型的构建

我们在模型里纳入了很多自变量, 我们一一解释:

1
2
3
4
5
6
model1 <- lmer(formula = negaff ~ 1 + day + stress_trait_c + 
stress_state + stress_state:stress_trait_c +
(1 + stress_state|id),
data=ddata_long,
na.action=na.exclude)
summary(model1)
输出(plain):
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
negaff ~ 1 + day + stress_trait_c + stress_state + stress_state:stress_trait_c +
(1 + stress_state | id)
Data: ddata_long

REML criterion at convergence: 3162.4

Scaled residuals:
Min 1Q Median 3Q Max
-3.5368 -0.6127 -0.0729 0.5093 4.4164

Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.2135 0.4621
stress_state 0.1257 0.3546 0.53
Residual 0.4038 0.6355
Number of obs: 1438, groups: id, 190

Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.695e+00 4.583e-02 3.922e+02 58.806 <2e-16
day -6.580e-02 7.552e-03 1.250e+03 -8.713 <2e-16
stress_trait_c 1.038e+00 7.946e-02 1.859e+02 13.067 <2e-16
stress_state 7.647e-01 4.561e-02 1.664e+02 16.765 <2e-16
stress_trait_c:stress_state 1.550e-01 9.780e-02 1.584e+02 1.585 0.115

(Intercept) ***
day ***
stress_trait_c ***
stress_state ***
stress_trait_c:stress_state
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
(Intr) day strs__ strss_
day -0.569
strss_trt_c 0.004 0.007
stress_stat 0.216 0.012 0.002
strss_tr_:_ 0.034 -0.057 0.268 -0.118

结果解读

固定效应 Fixed Effects:

随机效应 Random Effects:

预测值

基于样本数据和已有模型, 可以估计因变量的值, 也叫因变量的预测值, 如何获取因变量预测值:

1
2
3
# 保存模型的预测结果
ddata_long$pred_m1 <- predict(model1)
head(ddata_long)
输出(html):
A data.frame: 6 × 13
iddaynegaffpssstressbfi_nstress_traitnegaff_traitbfi_n_cstress_trait_cstress_statenegaff_statepred_m1
<int><int><dbl><dbl><dbl><dbl><dbl><dbl><dbl[,1]><dbl[,1]><dbl><dbl><dbl>
110103.02.501.5021.06251.5-0.9815789-0.3330326 0.4375 1.52.122456
210112.32.751.2521.06251.5-0.9815789-0.3330326 0.1875 0.81.908520
310121.03.500.5021.06251.5-0.9815789-0.3330326-0.5625-0.51.398305
410131.33.001.0021.06251.5-0.9815789-0.3330326-0.0625-0.21.628788
510141.12.751.2521.06251.5-0.9815789-0.3330326 0.1875-0.41.711131
610151.02.751.2521.06251.5-0.9815789-0.3330326 0.1875-0.51.645335

获取参数的置信区间

可以使用confint函数获取参数的置信区间, 但是结果比较难看懂, 因为里面很多参数是我们不熟悉的(sig01, sig02 等),
因为这些参数的名字其实是来自于模型的推导中用到的, 你没有深入了解过混合效应模型的参数估计方法, 你可能看不懂,
不过, 我们只能试着让你有一个感性认识:

1
confint(model1)
输出(stream):
Computing profile confidence intervals ...
输出(html):
A matrix: 9 × 2 of type dbl
2.5 %97.5 %
.sig01 0.40390297 0.52245693
.sig02 0.27635505 0.77129752
.sig03 0.25230527 0.45069358
.sigma 0.60962930 0.66244352
(Intercept) 2.60530435 2.78468853
day-0.08058845-0.05099124
stress_trait_c 0.88249219 1.19390741
stress_state 0.67339998 0.85449747
stress_trait_c:stress_state-0.03788643 0.34703329

可视化

先看被试层面的变量, negaff_trait 是消极情绪特质, 其实就是每天的消极情绪取被试内的均值, 代表了被试每天的平均消极情绪,
这个变量很可能受到压力特质(被试每天压力状态的一个被试内均值)的影响, 因此可以绘制这两个变量的关系:

1
2
3
4
5
6
7
8
9
ggplot(data=personmeans, aes(x=stress_trait, y=negaff_trait, group=factor(id)), legend=FALSE) +
geom_point(colour="gray40") +
geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
xlab("Trait Stress") + ylab("Trait Negative Affect") +
theme_classic() +
theme(axis.title=element_text(size=16),
axis.text=element_text(size=12),
plot.title=element_text(size=16, hjust=.5)) +
ggtitle("Between-Person Association Plot\nTrait Stress & Negative Affect")
输出(stream):
Warning message: "Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0. ℹ Please use `linewidth` instead." `geom_smooth()` using formula = 'y ~ x'

下面绘制negaff_state和stress_state的关系 , 因为这两个变量都是日测数据, 因此每个被试都应当是不同的。

1
2
3
4
5
6
7
8
9
ggplot(data=ddata_long, aes(x=stress_state, y=negaff_state, group=factor(id), colour="gray"), legend=FALSE) +
geom_smooth(method=lm, se=FALSE, fullrange=FALSE, lty=1, size=.5, color="gray40") +
geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
xlab("Stress State") + ylab("Predicted State Negative Affect") +
theme_classic() +
theme(axis.title=element_text(size=18),
axis.text=element_text(size=14),
plot.title=element_text(size=18, hjust=.5)) +
ggtitle("Within-Person Association Plot\nPerceived Stress & Negative Affect")
输出(stream):
`geom_smooth()` using formula = 'y ~ x' Warning message: "Removed 20 rows containing non-finite values (`stat_smooth()`)." `geom_smooth()` using formula = 'y ~ x' Warning message: "Removed 20 rows containing non-finite values (`stat_smooth()`)."

增加神经质变量作为预测变量

bfi_n_c 是大五人格量表中的神经质变量, 它经过处理得到了中心化后的变量,
这个变量中心化的目的也是因为这个变量会参与调节效应, 构建模型如下:

1
2
3
4
5
6
7
8
9
10
# fit model
model2 <- lmer(formula = negaff ~ 1 + day + stress_trait_c +
bfi_n_c + stress_trait_c:bfi_n_c +
stress_state + stress_state:stress_trait_c +
stress_state:bfi_n_c + stress_state:stress_trait_c:bfi_n_c +
(1 + stress_state|id),
data=ddata_long,
na.action=na.exclude)
#Look at results
summary(model2)
输出(plain):
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
negaff ~ 1 + day + stress_trait_c + bfi_n_c + stress_trait_c:bfi_n_c +
stress_state + stress_state:stress_trait_c + stress_state:bfi_n_c +
stress_state:stress_trait_c:bfi_n_c + (1 + stress_state | id)
Data: ddata_long

REML criterion at convergence: 3161.8

Scaled residuals:
Min 1Q Median 3Q Max
-3.4271 -0.6011 -0.0749 0.5045 4.4732

Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.1955 0.4422
stress_state 0.1238 0.3518 0.51
Residual 0.4040 0.6356
Number of obs: 1438, groups: id, 190

Fixed effects:
Estimate Std. Error df t value
(Intercept) 2.690e+00 4.556e-02 3.944e+02 59.055
day -6.545e-02 7.572e-03 1.247e+03 -8.644
stress_trait_c 9.695e-01 7.878e-02 1.835e+02 12.307
bfi_n_c 1.543e-01 3.917e-02 1.809e+02 3.939
stress_state 7.687e-01 4.682e-02 1.673e+02 16.418
stress_trait_c:bfi_n_c 3.715e-02 7.832e-02 1.824e+02 0.474
stress_trait_c:stress_state 1.254e-01 1.015e-01 1.692e+02 1.235
bfi_n_c:stress_state 7.595e-02 4.845e-02 1.549e+02 1.568
stress_trait_c:bfi_n_c:stress_state -3.167e-02 1.024e-01 1.722e+02 -0.309
Pr(>|t|)
(Intercept) < 2e-16 ***
day < 2e-16 ***
stress_trait_c < 2e-16 ***
bfi_n_c 0.000117 ***
stress_state < 2e-16 ***
stress_trait_c:bfi_n_c 0.635819
stress_trait_c:stress_state 0.218473
bfi_n_c:stress_state 0.119006
stress_trait_c:bfi_n_c:stress_state 0.757565
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
(Intr) day strs__ bf_n_c strss_ st__:__ st__:_ bf__:_
day -0.576
strss_trt_c 0.003 0.005
bfi_n_c -0.008 0.007 -0.224
stress_stat 0.204 0.004 0.000 0.001
strss_t_:__ -0.179 0.008 0.008 0.040 -0.054
strss_tr_:_ 0.042 -0.073 0.248 -0.056 -0.087 0.003
bf_n_c:str_ -0.033 0.058 -0.058 0.258 0.016 0.009 -0.244
strs__:__:_ -0.061 0.032 0.004 0.008 -0.233 0.243 -0.103 -0.059

结果解读

固定效应

随机效应

与之前的结论一致, 在此略去

获取模型预测值

1
2
ddata_long$pred_m2 <- predict(model2)
head(ddata_long)
输出(html):
A data.frame: 6 × 14
iddaynegaffpssstressbfi_nstress_traitnegaff_traitbfi_n_cstress_trait_cstress_statenegaff_statepred_m1pred_m2
<int><int><dbl><dbl><dbl><dbl><dbl><dbl><dbl[,1]><dbl[,1]><dbl><dbl><dbl><dbl>
110103.02.501.5021.06251.5-0.9815789-0.3330326 0.4375 1.52.1224562.097006
210112.32.751.2521.06251.5-0.9815789-0.3330326 0.1875 0.81.9085201.888382
310121.03.500.5021.06251.5-0.9815789-0.3330326-0.5625-0.51.3983051.393416
410131.33.001.0021.06251.5-0.9815789-0.3330326-0.0625-0.21.6287881.614307
510141.12.751.2521.06251.5-0.9815789-0.3330326 0.1875-0.41.7111311.692027
610151.02.751.2521.06251.5-0.9815789-0.3330326 0.1875-0.51.6453351.626575

调节效应的绘制

根据上面的结果,stress_state:bfi_n_c的调节效应是显著的, 意味着我们有必要进一步将调节效应可视化,
通常我们有两种可视化调节效应的方法, 1是选点法, 就是自变量和调节变量选择 M±SD 的值作为点, 绘制在不同调节变量取值下, 自变量与因变量的关系;
另一种方法是绘制JN图,它可以看到调节变量取值什么范围内, 自变量对因变量的效应是显著的。

选点法

我们先看下自变量和调节变量的描述性统计, 因为我们需要用到均值和标准差。

1
2
describe(ddata_long$bfi_n_c)
describe(ddata_long$stress_state)
输出(html):
A psych: 1 × 13
varsnmeansdmediantrimmedmadminmaxrangeskewkurtosisse
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
X111458-0.012100210.95646190.018421050.0025820121.4826-1.9815792.0184214-0.07586913-0.79088490.02504891
输出(html):
A psych: 1 × 13
varsnmeansdmediantrimmedmadminmaxrangeskewkurtosisse
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
X1114452.378342e-180.4943735-0.03125-0.015474130.4633125-1.752.1253.8750.35801080.79241720.01300533

我们使用 effect 函数来计算不同自变量和调节变量取值下因变量的值, effect 函数的term参数就是你关注的交互项,即自变量和调节变量的乘积。
mod 参数就是我们之前拟合得到的模型; xlevels 用于设置选点值, 比如自变量stress_state的正负一个标准差的值就是c(-0.49,+0.49) 。

1
2
3
#calculate effect
effects_model2 <- effect(term="bfi_n_c*stress_state", mod=model2,xlevels=list(bfi_n_c=c(-0.96, +0.96), stress_state=c(-0.49,+0.49)))
summary(effects_model2)
输出(stream):
NOTE: bfi_n_c:stress_state is not a high-order term in the model Warning message in Analyze.model(focal.predictors, mod, xlevels, default.levels, : "the predictors stress_trait_c, bfi_n_c are one-column matrices that were converted to vectors"
输出(plain):

bfi_n_c*stress_state effect
stress_state
bfi_n_c -0.49 0.49
-0.96 1.964450 2.644774
0.96 2.188157 3.011998

Lower 95 Percent Confidence Limits
stress_state
bfi_n_c -0.49 0.49
-0.96 1.858012 2.510688
0.96 2.080619 2.876440

Upper 95 Percent Confidence Limits
stress_state
bfi_n_c -0.49 0.49
-0.96 2.070889 2.778860
0.96 2.295694 3.147556

这个结果输出的是自变量和调节变量不同取值下, 因变量的值; 同时输出了因变量值的95%置信区间。有了这些数据, 我们就可以绘制简单效应图。

1
2
3
4
5
6
7
8
9
10
11
#convert to dataframe
effectsdata <- as.data.frame(effects_model2)
#plotting the effect evaluation (with standard error ribbon)
ggplot(data=effectsdata, aes(x=stress_state, y=fit, group=bfi_n_c), legend=FALSE) +
geom_point() +
geom_line() +
#geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.3) +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.15) +
xlab("Stress State") + xlim(-2,2) +
ylab("Predicted Negative Affect") + ylim(1,7) +
ggtitle("Differences in Stress Reactivity across Neuroticism")

johnson_neyman 这个函数可以用于绘制 JN图 , 关于这个图的原理可以看这篇文章《Johnson-Neyman图原理和制作Excel工具分享》, 并且这篇文章介绍了如何使用excel绘制JN图

1
johnson_neyman(model=model2, pred=stress_state, modx=bfi_n_c)
输出(plain):
JOHNSON-NEYMAN INTERVAL

When bfi_n_c is INSIDE the interval [-4.45, 40.14], the slope of
stress_state is p < .05.

Note: The range of observed values of bfi_n_c is [-1.98, 2.02]

一些有用函数

1
BIC(logLik(model2))
输出(html):
3256.29897530418
1
logLik(logLik(model2))
输出(plain):
'log Lik.' -1580.888 (df=13)
1
BIC(logLik(model2))
输出(html):
3256.29897530418

数据下载

本教程所有用到的代码和数据都可以在这里下载。

注意
统计咨询请加QQ 2726725926, 微信 shujufenxidaizuo, SPSS统计咨询是收费的, 不论什么模型都可以, 只限制于1个研究内.