使用R语言实现Johnson-Neyman分析

摘  要： 为利用R语言实现Johnson-Neyman分析，本文首先利用car程序包对回归模型进行诊断，利用ggplot2程序包绘制散点图并添加拟合直线，然后利用interactions程序包进行Johnson-Neyman分析，结果显示 interactions程序包能够正确计算两组连续型因变量间差异有统计学意义时的协变量的取值范围并画出Johnson-Neyman图。R语言是实现Johnson-Neyman分析的有效工具，并且在绘图上具有优势，功能强大。

关键词： R语言;Johnson-Neyman

中图分类号： O212.4    文献标识码： A    DOI：10.3969/j.issn.1003-6970.2019.12.005

本文著录格式：高启胜. 使用R语言实现Johnson-Neyman分析[J]. 软件，2019，40（12）：2124

Johnson-Neyman Analysis in R

GAO Qi-sheng

（School of Public Health， HangZhou Medical College， Hangzhou， Zhejiang 310053 China）

【Abstract】： To implement the Johnson-Neyman analysis in R， the car package was used to diagnose the regression model firstly， the ggplot2 package was used to draw scatter diagram and add fitting straight-line， the interactions package was used to conduct Johnson-Neyman analysis. The results showed that the interactions package can correctly calculate the value range of covariate when two continuous dependent variables have a significant difference and draw Johnson-Neyman plot. R programmer is a powerful tool to implement Johnson-Neyman analysis， and has advantages in plotting.

【Key words】： R language; Johnson-Neyman

0  引言

在实验性研究中，通常对受试对象进行随机化，以使各种非处理因素如病种、病期、病型、年龄、性别、生活和心理等因素在各组间保持均衡，从而减少非处理因素的影响，使所考察的实验效应真实地显露出来。然而，对于无法对受试对象进行随机化或即使随机化仍存在组间基线不均衡的研究，通常需要应用分层分析、协方差分析、回归分析和PSM（倾向评分匹配）等方法来对重要的非处理因素予以控制。其中，协方差分析（ANCOVA）是将线性回归分析与方差分析结合起来的一种分析方法，其基本思想是将定量变量X（不可控或未控因素）对Y的影响看做协变量，建立因变量Y随X变化的线性回归关系，并利用这种回归关系把X值化为相等后再进行各组Y的修正均数的比较。其实质是从Y的总离均差平方和中扣除X对Y的回归平方和，对残差平方和作进一步分解后再进行方差分析，以更准确的评价各种处理的效应[1]。

协方差分析的应用条件包括：（1）观察对象相关独立;（2）各组因变量服从正态分布且方差相等;（3）协变量与因变量之间的关系是线性关系;   （4）在各组中协变量的回归系数（即各回归线的斜率）必须是相等的，即各组的回归线是平行[2]。正是由于假定条件的复杂性，导致协方差分析在实际中出现误用或使用不普遍等问题[3]。目前，国内对于协方差分析违背斜率平行性假定时的分析方法介绍较少。近年来，Johnson-Neyman（JN）方法作为违背平行假定时恰当的分析方法受到越来越多的关注[4-6]。JN方法在SPSS和SAS中主要通过加载PROCESS插件来实现，但是在可视化图形上存在不足[7]，本文拟通过实例展示如何应用R语言逐步实现JN法。

1  Johnson-Neyman方法基本原理

Johnson-Neyman方法在1936年由Johnson和Neyman最早提出[8]，主要用于协方差分析中不满足斜率平行假定时，计算两组间因变量差异有统计学意义时（Simultaneous regions of significance，SROS）的协变量的取值范围。2005年Bauer和Curran将JN方法拓展到回归模型中，适用于调节变量为定量变量的调节效应分析[9]。如果因变量Y与自变量X 的關系受到第三个变量M的影响，就称M为调节变量。调节变量可以是定性的，也可以是定量的。通常应用回归模型进行调节效应分析[10]，简要模型如下：

其中，Y与X的关系由回归系数b1 + b3M来刻画，（b1 + b3M）表示简单斜率（simple slope），b3衡量了调节效应（moderating effect）的大小。如果b3显著，说明M的调节效应显著。b3也代表了X与M的交互效应，调节效应与交互效应分析从统计学角度而言是一致的，然而调节效应中，调节变量是确定的，不能互换，而交互效应分析中变量的地位是等价的[11]。

JN法就是对简单斜率进行检验，以确定简单斜率显著与否的分界点[10]。

当进行两组间比较时，调节变量M的分界点R计算如下[12]：

其中，和是两组协变量平均数，和是两组协变量观察值的平方和，为整体回归分析的残差平方和，和分别为两组回归线的截距，和分别为两组回归线的回归系数，N为总样本例数，n1和n2分别是两组的样本例数，是自由度为（2，N-4）时检验水准为α（通常为0.05）时的F临界值。

2  在R语言上实现Johnson-Neyman分析

2.1  软件安装与程序加载

R语言是一个自由、免费、开源代码开放的软件，是用于统计计算和统计制图的优秀工具。截止至2019年3月6日，R软件最新版本为R3.5.3，已收录13775个程序包，用户可从官方网站http：//www. r-project. org下载最新的软件版本。R软件安装完毕后，需要安装和加载所需的程序包，涉及的安装包有foreign，car，ggplot2和interactions包，其中，foreign包用于导入SPSS数据文件，car包用于对回归模型进行诊断，ggplot2包用于绘制可编辑的散点图并添加拟合直线，interactions包用于进行JN分析并画出JN图。

2.2  导入数据文件

此数据为两组干预前后定量变量得分的比较，每组各20名患者，数据文件格式为sav，文件名为“jndata”，保存路径：C：＼Users＼mooshaa＼Desktop，导入数据代码为：

install.packages（"foreign"）

library（foreign）

jndata<- read.spss（"C：/Users/mooshaa/Desktop/jndata. sav"，to.data.frame=T，use.value.lables=F）

上述命令执行后，在命令栏中输入命令“jndata”，即可显示导入的数据。其中，“post”代表干预后得分，“pre”代表干预前得分，“group”代表分组变量，group1=0，group2=1。

2.3  建立线性回归模型并进行回归诊断

建立因变量为post，自变量为pre和group的多重线性回归方程，并引入pre和group的交互项。输入代码如下：

jntest<- lm（post ～ pre* group， data = jndata）

summary（jntest）

结果见表1，其中模型整体检验的F=170.60，P<0.001，调整R2为0.929，表明所建立的回归模型有统计学意义。各回归系数的假设检验结果显示，pre、group以及两者的交互项均存在统计学意义。交互项有统计学意义即表明调节效应显著，两组的回归斜率不平行。建立的回归方程为： 。对于group1而言，pre得分每增加一个单位，post得分平均增加1.04个单位，对于group2而言，pre得分每增加一个单位，post得分平均改变的单位数等于group1改變的单位数与交互作用项的偏回归系数之和，即1.04-0.95=0.09个单位。

表1  线性回归系数检验

Tab.1  Linear regression coefficient test

Coefficient Standard Error t p-value

Constant 24.26 2.57 9.42 <0.001

pre 1.04 0.05 21.67 <0.001

group 55.11 3.71 14.87 <0.001

Pre：group –0.95 0.07 –14.12 <0.001

进一步对回归模型的线性（linear）、独立（independent）、正态（normal）和等方差（equal variance）等假定条件进行验证，输入命令如下：

install.packages（"car"）

library（car）

shapiro.test（jntest$residuals）

durbinWatsonTest（jntest）

outlierTest（jntest）

ncvTest（jntest）

结果显示：正态性检验的W=1，P=0.5，残差服从正态分布;独立性检验的D-W值=2.1，P=0.84，表示各观察值相互独立;极端值检验中，没有学化生残差检验的P<0.05，不存在极端值;等方差检验中，χ2=1.6，P=0.2，表明残差的方差不随预测值的变化而改变，满足等方差条件。

2.4  交互作用可视化图形

由于pre与group存在交互作用，分别拟合两组pre和post的回归直线，输入代码如下：

install.packages（"ggplot2"）

library（ggplot2）

jndata$group <- factor（jndata$group）

regpl<- ggplot（jndata，aes（x=pre，y=post，group= group））+ geom_point（aes（shape=group））+

geom_smooth（method = "lm"，color="black"）+

theme（axis.text=element_text（face="bold"，size=10），axis.title.x=element_text（size=12，face="bold"），axis.title.y=element_text（size=12，face="bold"））+theme（panel.background =

element_rect（fill = "white"）， legend.key= element_rect（fill= "white"），axis.line.x =

element_line（colour = "black"， size = 0.5）， axis. line.y = element_line（colour = "black"， size = 0.5）） + theme（legend.position =c（0.95，0.9））

regpl

图1可见，两组间pre得分和post得分均存在线性关系，并且两条直线相交，进一步说明违背斜率平行假定，不能使用协方差分析。

2.5  JN分析

当协方差分析不满足平行假定时，通过JN方法可以进一步确定两组因变量间差异有统计学意义时的协变量的取值范围（SROS），也即简单斜率显著时调节变量的分界点。输入代码如下：

install.packages（"interactions"）

library（interactions）

johnson_neyman（jntest，pred=group，modx=pre， plot=T）

图1  分组线性回归图

Fig.1  Grouped linear regression

结果见图2，统计显示当pre得分的取值范围在[55.82，60.69]之外时，简单斜率p<0.05，即在图2中两条竖直虚线与阴影部分的重合区域范围简单斜率的置信区间包含0，此区域外，简单斜率显著不为0。其中，加粗横线的长度代表pre得分的取值范围：[20，82]。

图2  Johnson-Neyman图

Fig.2  Johnson-Neyman plot

为进一步展现pre得分在临界点时，两组间post得分的差异情况，输入代码如下：

regpl+geom_vline（xintercept =c（55.82，60.69），line type="dotted"， ， size=1 ）+

annotate（"text"， x = 37， y=42， label = "No significant differences in post ＼n  between groups when pre are 55.82 to 60.69 "）

结果见图3，可见，当pre≤55.82时，group2中的post得分显著高于group1;当pre≥60.69时group1中的post得分显著高于group2;pre的取值范围在（55.82，60.69）时，post得分在两组的差异无统计学意义。有统计学意义并不意味着有实际意义，为进一步了解两组间在有统计学意义区域的实际效应值大小，可以分别以調节变量pre的临界点筛选数据，然后分别对post得分进行独立样本t检验，输入代码如下：

subjndata <- subset（jndata，pre<=55.82）

t.test（post～group，data=subjndata）

subjndata1 <- subset（jndata，pre>=60.69）;subjndata1

t.test（post～group，data=subjndata1）

结果显示，当pre≤55.82时，group2（=83.73，S=4.00）比group1（=67.23，S=11.14）平均高16.50（95%CI：9.44～23.54）;当pre≥60.69时，group1（=96.86，S=3.72）比group2（=84.63，S=1.77）平均高12.23（95%CI：8.71～15.75）。

图3  两组间post得分显著区域

Fig.3  Significant area of post score between the two groups

3  结论

研究表明，无论协变量是否有显著差异，在满足假定条件下，均可以采用协方差分析，以消除基线协变量的影响，提高统计检验效能[2]。但当违背应用条件时，协方差分析就不再适用。其中，当不满足回归直线平行假定时，Johnson-Neyman方法可作为协方差分析恰当的替代方法，它从调节效应分析角度，可以计算出简单斜率显著与否时连续性定量调节变量的临界值。值得注意的是，根据临界值落在协变量取值范围的不同情况，简单斜率显著或者SROS的区域范围也有多种情况[7，10]，若R1或者R2有一个落在[Mmin，Mmax]，则M≤R1/R2或者M≥R1/R2时简单斜率显著;若Mmin≤R1

参考文献

[1]金丕焕， 陈峰. 医学统计方法（第三版）. 上海： 复旦大学出版社， 2009， 131-132.

[2]张天嵩. 调整基线差异：协方差分析[J]. 临床与病理杂志， 2015， 35（12）： 2043-2048.

[3]鲍贵. 外语教学研究中的协方差分析误区[J]. 外语测试与教学， 2016（3）： 51-59.

[4]Miller J W， Stromeyer W R， Schwieterman M A. Extensions of the Johnson-Neyman Technique to Linear Models With Curvilinear Effects： Derivations and Analytical Tools. Multivariate Behavioral Research， 2013， 48（2）： 267-300.

[5]Rast P， Rush J， Piccinin A， et al. The Identification of Regions of Significance in the Effect of Multimorbidity on Depressive Symptoms Using Longitudinal Data： An Application of the Johnson-Neyman Technique. Gerontology， 2014， 60（3）： 274-281.

[6]Lazar AA， Gansky SA， Halstead DD， et al. Improving Patient Care Using the Johnson-Neyman Analysis of Heterogeneity of Treatment Effects According to Individuals Baseline Characteristics.J Dent Oral Craniofac Epidemiol. 2013，1（3）： 19-33.

[7]Carden SW， Holtzman N S， Strube MJ. CAHOST： An Excel Workbook for Facilitating the Johnson-Neyman Technique for Two-Way Interactions in Multiple Regression. Frontiers in Psychology， 2017， 8： 1293.

[8]Johnson PO， Neyman J. Tests of Certain Linear Hypotheses and their Application to Some Educational Problems. Statistical Research Memoirs， 1936， 1： 57-93.

[9]Bauer DJ， Curran PJ. Probing Interactions in Fixed and Multilevel Regression： Inferential and Graphical Techniques. Multivariate Behavioral Research， 2005， 40（3）： 373-400.

[10]方杰， 溫忠麟， 梁东梅， 等. 基于多元回归的调节效应分析[J]. 心理科学， 2015（3）： 715-720.

[11]温忠麟， 侯杰泰， 张雷. 调节效应与中介效应的比较和应用[J]. 心理学报， 2005， 37（2）： 268-274.

[12]Johnson TR. Violation of the Homogeneity of Regression Slopes Assumption in ANCOVA for Two-Group Pre-Post Designs： Tutorial on A Modified Johnson-Neyman Procedure. The Quantitative Methods for Psychology， 2016， 12（3）： 253-263.

[13]Lazar AA， Zerbe GO. Solutions for Determining the Significance Region Using the Johnson-Neyman Type Procedure in Generalized Linear （Mixed） Models. Journal of Educational and Behavioral Statistics， 2011， 36（6）： 699-719.