使用R語言實現(xiàn)Johnson-Neyman分析

摘 ?要： 為利用R語言實現(xiàn)Johnson-Neyman分析，本文首先利用car程序包對回歸模型進(jìn)行診斷，利用ggplot2程序包繪制散點圖并添加擬合直線，然后利用interactions程序包進(jìn)行Johnson-Neyman分析，結(jié)果顯示 interactions程序包能夠正確計算兩組連續(xù)型因變量間差異有統(tǒng)計學(xué)意義時的協(xié)變量的取值范圍并畫出Johnson-Neyman圖。R語言是實現(xiàn)Johnson-Neyman分析的有效工具，并且在繪圖上具有優(yōu)勢，功能強大。

關(guān)鍵詞： R語言;Johnson-Neyman

中圖分類號： O212.4 ? ?文獻(xiàn)標(biāo)識碼： A ? ?DOI：10.3969/j.issn.1003-6970.2019.12.005

本文著錄格式：高啟勝. 使用R語言實現(xiàn)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 ?引言

在實驗性研究中，通常對受試對象進(jìn)行隨機化，以使各種非處理因素如病種、病期、病型、年齡、性別、生活和心理等因素在各組間保持均衡，從而減少非處理因素的影響，使所考察的實驗效應(yīng)真實地顯露出來。然而，對于無法對受試對象進(jìn)行隨機化或即使隨機化仍存在組間基線不均衡的研究，通常需要應(yīng)用分層分析、協(xié)方差分析、回歸分析和PSM（傾向評分匹配）等方法來對重要的非處理因素予以控制。其中，協(xié)方差分析（ANCOVA）是將線性回歸分析與方差分析結(jié)合起來的一種分析方法，其基本思想是將定量變量X（不可控或未控因素）對Y的影響看做協(xié)變量，建立因變量Y隨X變化的線性回歸關(guān)系，并利用這種回歸關(guān)系把X值化為相等后再進(jìn)行各組Y的修正均數(shù)的比較。其實質(zhì)是從Y的總離均差平方和中扣除X對Y的回歸平方和，對殘差平方和作進(jìn)一步分解后再進(jìn)行方差分析，以更準(zhǔn)確的評價各種處理的效應(yīng)[1]。

協(xié)方差分析的應(yīng)用條件包括：（1）觀察對象相關(guān)獨立;（2）各組因變量服從正態(tài)分布且方差相等;（3）協(xié)變量與因變量之間的關(guān)系是線性關(guān)系; ? （4）在各組中協(xié)變量的回歸系數(shù)（即各回歸線的斜率）必須是相等的，即各組的回歸線是平行[2]。正是由于假定條件的復(fù)雜性，導(dǎo)致協(xié)方差分析在實際中出現(xiàn)誤用或使用不普遍等問題[3]。目前，國內(nèi)對于協(xié)方差分析違背斜率平行性假定時的分析方法介紹較少。近年來，Johnson-Neyman（JN）方法作為違背平行假定時恰當(dāng)?shù)姆治龇椒ㄊ艿皆絹碓蕉嗟年P(guān)注[4-6]。JN方法在SPSS和SAS中主要通過加載PROCESS插件來實現(xiàn)，但是在可視化圖形上存在不足[7]，本文擬通過實例展示如何應(yīng)用R語言逐步實現(xiàn)JN法。

1 ?Johnson-Neyman方法基本原理

Johnson-Neyman方法在1936年由Johnson和Neyman最早提出[8]，主要用于協(xié)方差分析中不滿足斜率平行假定時，計算兩組間因變量差異有統(tǒng)計學(xué)意義時（Simultaneous regions of significance，SROS）的協(xié)變量的取值范圍。2005年Bauer和Curran將JN方法拓展到回歸模型中，適用于調(diào)節(jié)變量為定量變量的調(diào)節(jié)效應(yīng)分析[9]。如果因變量Y與自變量X 的關(guān)系受到第三個變量M的影響，就稱M為調(diào)節(jié)變量。調(diào)節(jié)變量可以是定性的，也可以是定量的。通常應(yīng)用回歸模型進(jìn)行調(diào)節(jié)效應(yīng)分析[10]，簡要模型如下：

其中，Y與X的關(guān)系由回歸系數(shù)b1 + b3M來刻畫，（b1 + b3M）表示簡單斜率（simple slope），b3衡量了調(diào)節(jié)效應(yīng)（moderating effect）的大小。如果b3顯著，說明M的調(diào)節(jié)效應(yīng)顯著。b3也代表了X與M的交互效應(yīng)，調(diào)節(jié)效應(yīng)與交互效應(yīng)分析從統(tǒng)計學(xué)角度而言是一致的，然而調(diào)節(jié)效應(yīng)中，調(diào)節(jié)變量是確定的，不能互換，而交互效應(yīng)分析中變量的地位是等價的[11]。

JN法就是對簡單斜率進(jìn)行檢驗，以確定簡單斜率顯著與否的分界點[10]。

當(dāng)進(jìn)行兩組間比較時，調(diào)節(jié)變量M的分界點R計算如下[12]：

其中，和是兩組協(xié)變量平均數(shù)，和是兩組協(xié)變量觀察值的平方和，為整體回歸分析的殘差平方和，和分別為兩組回歸線的截距，和分別為兩組回歸線的回歸系數(shù)，N為總樣本例數(shù)，n1和n2分別是兩組的樣本例數(shù)，是自由度為（2，N-4）時檢驗水準(zhǔn)為α（通常為0.05）時的F臨界值。

2 ?在R語言上實現(xiàn)Johnson-Neyman分析

2.1 ?軟件安裝與程序加載

R語言是一個自由、免費、開源代碼開放的軟件，是用于統(tǒng)計計算和統(tǒng)計制圖的優(yōu)秀工具。截止至2019年3月6日，R軟件最新版本為R3.5.3，已收錄13775個程序包，用戶可從官方網(wǎng)站http：//www. r-project. org下載最新的軟件版本。R軟件安裝完畢后，需要安裝和加載所需的程序包，涉及的安裝包有foreign，car，ggplot2和interactions包，其中，foreign包用于導(dǎo)入SPSS數(shù)據(jù)文件，car包用于對回歸模型進(jìn)行診斷，ggplot2包用于繪制可編輯的散點圖并添加擬合直線，interactions包用于進(jìn)行JN分析并畫出JN圖。

2.2 ?導(dǎo)入數(shù)據(jù)文件

此數(shù)據(jù)為兩組干預(yù)前后定量變量得分的比較，每組各20名患者，數(shù)據(jù)文件格式為sav，文件名為“jndata”，保存路徑：C：＼Users＼mooshaa＼Desktop，導(dǎo)入數(shù)據(jù)代碼為：

install.packages（"foreign"）

library（foreign）

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

上述命令執(zhí)行后，在命令欄中輸入命令“jndata”，即可顯示導(dǎo)入的數(shù)據(jù)。其中，“post”代表干預(yù)后得分，“pre”代表干預(yù)前得分，“group”代表分組變量，group1=0，group2=1。

2.3 ?建立線性回歸模型并進(jìn)行回歸診斷

建立因變量為post，自變量為pre和group的多重線性回歸方程，并引入pre和group的交互項。輸入代碼如下：

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

summary（jntest）

結(jié)果見表1，其中模型整體檢驗的F=170.60，P<0.001，調(diào)整R2為0.929，表明所建立的回歸模型有統(tǒng)計學(xué)意義。各回歸系數(shù)的假設(shè)檢驗結(jié)果顯示，pre、group以及兩者的交互項均存在統(tǒng)計學(xué)意義。交互項有統(tǒng)計學(xué)意義即表明調(diào)節(jié)效應(yīng)顯著，兩組的回歸斜率不平行。建立的回歸方程為： 。對于group1而言，pre得分每增加一個單位，post得分平均增加1.04個單位，對于group2而言，pre得分每增加一個單位，post得分平均改變的單位數(shù)等于group1改變的單位數(shù)與交互作用項的偏回歸系數(shù)之和，即1.04-0.95=0.09個單位。

表1 ?線性回歸系數(shù)檢驗

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

進(jìn)一步對回歸模型的線性（linear）、獨立（independent）、正態(tài)（normal）和等方差（equal variance）等假定條件進(jìn)行驗證，輸入命令如下：

install.packages（"car"）

library（car）

shapiro.test（jntest$residuals）

durbinWatsonTest（jntest）

outlierTest（jntest）

ncvTest（jntest）

結(jié)果顯示：正態(tài)性檢驗的W=1，P=0.5，殘差服從正態(tài)分布;獨立性檢驗的D-W值=2.1，P=0.84，表示各觀察值相互獨立;極端值檢驗中，沒有學(xué)化生殘差檢驗的P<0.05，不存在極端值;等方差檢驗中，χ2=1.6，P=0.2，表明殘差的方差不隨預(yù)測值的變化而改變，滿足等方差條件。

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得分均存在線性關(guān)系，并且兩條直線相交，進(jìn)一步說明違背斜率平行假定，不能使用協(xié)方差分析。

2.5 ?JN分析

當(dāng)協(xié)方差分析不滿足平行假定時，通過JN方法可以進(jìn)一步確定兩組因變量間差異有統(tǒng)計學(xué)意義時的協(xié)變量的取值范圍（SROS），也即簡單斜率顯著時調(diào)節(jié)變量的分界點。輸入代碼如下：

install.packages（"interactions"）

library（interactions）

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

圖1 ?分組線性回歸圖

Fig.1 ?Grouped linear regression

結(jié)果見圖2，統(tǒng)計顯示當(dāng)pre得分的取值范圍在[55.82，60.69]之外時，簡單斜率p<0.05，即在圖2中兩條豎直虛線與陰影部分的重合區(qū)域范圍簡單斜率的置信區(qū)間包含0，此區(qū)域外，簡單斜率顯著不為0。其中，加粗橫線的長度代表pre得分的取值范圍：[20，82]。

圖2 ?Johnson-Neyman圖

Fig.2 ?Johnson-Neyman plot

為進(jìn)一步展現(xiàn)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 "）

結(jié)果見圖3，可見，當(dāng)pre≤55.82時，group2中的post得分顯著高于group1;當(dāng)pre≥60.69時group1中的post得分顯著高于group2;pre的取值范圍在（55.82，60.69）時，post得分在兩組的差異無統(tǒng)計學(xué)意義。有統(tǒng)計學(xué)意義并不意味著有實際意義，為進(jìn)一步了解兩組間在有統(tǒng)計學(xué)意義區(qū)域的實際效應(yīng)值大小，可以分別以調(diào)節(jié)變量pre的臨界點篩選數(shù)據(jù)，然后分別對post得分進(jìn)行獨立樣本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）

結(jié)果顯示，當(dāng)pre≤55.82時，group2（=83.73，S=4.00）比group1（=67.23，S=11.14）平均高16.50（95%CI：9.44～23.54）;當(dāng)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得分顯著區(qū)域

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

3 ?結(jié)論

研究表明，無論協(xié)變量是否有顯著差異，在滿足假定條件下，均可以采用協(xié)方差分析，以消除基線協(xié)變量的影響，提高統(tǒng)計檢驗效能[2]。但當(dāng)違背應(yīng)用條件時，協(xié)方差分析就不再適用。其中，當(dāng)不滿足回歸直線平行假定時，Johnson-Neyman方法可作為協(xié)方差分析恰當(dāng)?shù)奶娲椒?，它從調(diào)節(jié)效應(yīng)分析角度，可以計算出簡單斜率顯著與否時連續(xù)性定量調(diào)節(jié)變量的臨界值。值得注意的是，根據(jù)臨界值落在協(xié)變量取值范圍的不同情況，簡單斜率顯著或者SROS的區(qū)域范圍也有多種情況[7，10]，若R1或者R2有一個落在[Mmin，Mmax]，則M≤R1/R2或者M(jìn)≥R1/R2時簡單斜率顯著;若Mmin≤R1

參考文獻(xiàn)

[1]金丕煥， 陳峰. 醫(yī)學(xué)統(tǒng)計方法（第三版）. 上海： 復(fù)旦大學(xué)出版社， 2009， 131-132.

[2]張?zhí)灬? 調(diào)整基線差異：協(xié)方差分析[J]. 臨床與病理雜志， 2015， 35（12）： 2043-2048.

[3]鮑貴. 外語教學(xué)研究中的協(xié)方差分析誤區(qū)[J]. 外語測試與教學(xué)， 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]方杰， 溫忠麟， 梁東梅， 等. 基于多元回歸的調(diào)節(jié)效應(yīng)分析[J]. 心理科學(xué)， 2015（3）： 715-720.

[11]溫忠麟， 侯杰泰， 張雷. 調(diào)節(jié)效應(yīng)與中介效應(yīng)的比較和應(yīng)用[J]. 心理學(xué)報， 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.