基于C#與R語言的重復(fù)測量設(shè)計(jì)資料方差分析的自動(dòng)化實(shí)現(xiàn)

郭迎暄 陳 達(dá) 張雙喆 黃昌可 廉恒麗

1 溫州醫(yī)科大學(xué)附屬眼視光醫(yī)院杭州院區(qū)，330020 浙江 杭州;2 溫州醫(yī)科大學(xué)附屬眼視光醫(yī)院，325027 浙江 溫州

醫(yī)學(xué)臨床研究中常收集多次重復(fù)測量數(shù)據(jù)，一是可以減少個(gè)體間的誤差，二是為了節(jié)約樣本量。當(dāng)重復(fù)測量次數(shù)m≥3時(shí),稱為重復(fù)測量設(shè)計(jì)(repeated measurement design)或重復(fù)測量數(shù)據(jù)(repeated measurement data)。該數(shù)據(jù)可用于分析被觀察指標(biāo)在不同時(shí)間點(diǎn)上的變化趨勢(shì)與治療效果[1]。由于多次重復(fù)測量的數(shù)據(jù)來自于同一個(gè)體，這些數(shù)據(jù)之間存在相關(guān)性，與方差分析的基本假設(shè)——數(shù)據(jù)獨(dú)立性相違背，若采用常規(guī)的方差分析，則會(huì)增加犯I類錯(cuò)誤的概率[2]，因此需要采用重復(fù)測量方差分析。重復(fù)測量方差分析可用于以下目的：比較組間有無差異；比較各時(shí)間點(diǎn)間有無差異；比較組間的時(shí)間變化趨勢(shì)有無差異(即隨時(shí)間變化的多條曲線是否平行)[3]。數(shù)據(jù)需滿足以下條件：因變量沒有明顯異常值；因變量需服從近似正態(tài)分布；因變量的方差協(xié)方差矩陣相等，也稱滿足球形假設(shè)。近年來，醫(yī)學(xué)研究越來越離不開統(tǒng)計(jì)方法的使用，而每種統(tǒng)計(jì)方法都有其適用條件,醫(yī)學(xué)雜志發(fā)表的論文存在不同程度的統(tǒng)計(jì)錯(cuò)誤,影響研究結(jié)論。為了減少這一現(xiàn)象,提高論文水平, 提高醫(yī)生科研效率，迫切需要一款界面更加友好，操作更加簡單，結(jié)果更加直觀的統(tǒng)計(jì)自動(dòng)化軟件。

R軟件是目前醫(yī)學(xué)研究領(lǐng)域廣泛使用的免費(fèi)開源軟件，具有全面完善的統(tǒng)計(jì)分析功能和數(shù)量龐大的拓展包，具備跨平臺(tái)運(yùn)行和兼容各類統(tǒng)計(jì)軟件等優(yōu)點(diǎn)[4]。近年來，在醫(yī)療衛(wèi)生各領(lǐng)域受到越來越多統(tǒng)計(jì)學(xué)者的使用和認(rèn)可[5-7]。C#是一種面向?qū)ο蟮木幊陶Z言，利用C#.NET 開發(fā)的系統(tǒng)具有界面友好、執(zhí)行速度快、易維護(hù)和升級(jí)等優(yōu)點(diǎn)[8]。通過RDotNet組件，C#可以直接訪問R語言程序包。2種語言的這些特點(diǎn)，使得C#語言和R語言聯(lián)合編程具有可行性。為了讓臨床醫(yī)生輕松又快捷地完成臨床科研中重復(fù)測量設(shè)計(jì)資料的方差分析，本研究基于2種語言開發(fā)了一款簡單易用的軟件，可自動(dòng)完成重復(fù)測量設(shè)計(jì)資料的單、雙因素兩/多水平方差分析。

1 C#語言與R語言聯(lián)合編程方法

1.1 技術(shù)路線

1)安裝3.6.3版本的R語言運(yùn)行平臺(tái)RStudio和16.7.6版本的Microsoft Visual Studio Enterprise 2019 C#開發(fā)平臺(tái)用來做調(diào)試使用。在本研究中，主要使用readxl、ez、psych、RVAideMemoire、agricolae、PairedData、ggplot 2共7個(gè)分析包編寫R語言統(tǒng)計(jì)分析腳本，在C#項(xiàng)目中引用RDotNet.dll，調(diào)用編寫好的R語言統(tǒng)計(jì)腳本，實(shí)現(xiàn)重復(fù)測量設(shè)計(jì)資料方差分析。

2)自動(dòng)化輸出統(tǒng)計(jì)分析結(jié)果。首先判斷受試者內(nèi)因素的各個(gè)水平，因變量有沒有極端異常值；因變量是否服從近似正態(tài)分布；因變量的方差協(xié)方差矩陣相等(根據(jù)ezANOVA包統(tǒng)計(jì)分析結(jié)果中的Mauchly's球性檢驗(yàn)結(jié)果進(jìn)行判斷)。其中，若滿足球形假定(P>0.05)，直接輸出方差分析的結(jié)果；若不滿足球形假定(P<0.05)，需要參考GGe(Greenhouse-Geisser epsilon)或者HFe(Huynh-Feldt epsilon)對(duì)df進(jìn)行校正。當(dāng)GGe和HFe對(duì)應(yīng)的Epsilon值均>0.75,采用HFe對(duì)應(yīng)的Epsilon值對(duì)df進(jìn)行校正；如果GGe和HFe均<0.75，采用GGe對(duì)應(yīng)的Epsilon值對(duì)df進(jìn)行校正。其他情況默認(rèn)GGe的Epsilon值進(jìn)行校正。接著分析各因素間是否存在交互作用，若存在交互作用(P<0.05)，需要固定某一因素分析另一因素的單獨(dú)效應(yīng)；若不存在交互作用(P>0.05)，則直接分析各因素的主效應(yīng)。統(tǒng)計(jì)分析與結(jié)果選擇流程見圖1。

圖1 統(tǒng)計(jì)分析與結(jié)果選擇流程圖

1.2 R語言重復(fù)測量方差分析程序腳本

#1 將數(shù)據(jù)導(dǎo)入分析系統(tǒng)

library(readxl)

data <-read_excel(file.choose())

#2 正態(tài)性檢驗(yàn)(Shapiro-Wilk檢驗(yàn)，小樣本使用的正態(tài)性檢驗(yàn))、異常值檢驗(yàn)

library(RVAideMemoire)

byf.shapiro(VALUE～TIME+GROUP, data=data)

#3 基本描述,加載psych包，輸出組別和時(shí)間的均數(shù)與標(biāo)準(zhǔn)差，繪制QQ圖

library(psych)

library(car)

fit_lm <- lm(VALUE～GROUP*TIME,data= data)

qqPlot(fit_lm,simulate=TRUE,main=“Q-Q PLot”,lables=FALSE)

#4 加載R語言統(tǒng)計(jì)分析ez包，按照實(shí)驗(yàn)設(shè)計(jì)執(zhí)行統(tǒng)計(jì)模型腳本

library(ez)

model <-ezANOVA(data, dv = VALUE, wid = SUBJECT, within = .(TIME), between= GROUP,type = 3, detailed = T)

#4.1 滿足球形檢驗(yàn)直接輸出結(jié)果

model

#4.2 不滿足球形檢驗(yàn)，輸出校正后的結(jié)果

model

table(Int_DFn_TIME,Int_DFd_TIME)

table(Int_DFn_GT,Int_DFd_GT)

#5 若GROUP和TIME不存在交互作用，進(jìn)行GROUP、TIME的主效應(yīng)分析

#5.1 GROUP的主效應(yīng)分析，以及不同組別間兩兩配對(duì)比較

library (agricolae)

model <- aov (VALUE～GROUP,data=data)

Res <- LSD.test(model,"GROUP", p.adj ="bonferroni")

plot(Res)

#5.2 進(jìn)行TIME的主效應(yīng)分析，以及不同組別間兩兩配對(duì)比較

aov_T <-aov(VALUE～TIME + Error(SUBJECT/TIME), data = data)

summary(aov_T)

#6.若GROUP和TIME存在交互作用，分別進(jìn)行GROUP、TIME的單獨(dú)效應(yīng)分析

#6.1 GROUP的單獨(dú)效應(yīng)分析，各分組下不同時(shí)間點(diǎn)的兩兩配對(duì)比較

#6.2 TIME的單獨(dú)效應(yīng)分析，各分組下不同時(shí)間點(diǎn)的兩兩配對(duì)比較

aov_A <-aov(VALUE～TIME + Error(SUBJECT/TIME), data = g_d[["A"]])

summary(aov_A)

aov_B <-aov(VALUE～TIME + Error(SUBJECT/TIME), data = g_d[["B"]])

summary(aov_B)

aov_C <-aov(VALUE～TIME + Error(SUBJECT/TIME), data = g_d[["C"]])

summary(aov_C)

#6.3 不同時(shí)間點(diǎn)的兩兩配對(duì)比較

#7 繪制各因素的數(shù)據(jù)輪廓圖

fit<-aov(VALUE～GROUP*TIME + Error(SUBJECT/(TIME)),data)

summary(fit)

par(las=1)

par(mar=c(10,4,4,2))

with(data, interaction.plot( TIME,GROUP,VALUE,type="b", col=c("red","blue","green"), pch=c(16,18,20),main="Interaction Plot for GROUP and TIME"))

1.3 軟件安裝與界面

安裝時(shí)，需將包含可執(zhí)行程序的文件包拷貝至操作系統(tǒng)為Windows 7及以上版本的計(jì)算機(jī)上。雙擊文件夾中的CallRLanguage.exe文件，即可打開應(yīng)用程序。該軟件界面左側(cè)為菜單欄，右側(cè)上方為數(shù)據(jù)導(dǎo)入格式示例區(qū)，右側(cè)中部為參數(shù)設(shè)置區(qū)，右側(cè)下方為結(jié)果顯示區(qū)。見圖2。

圖2 方差分析結(jié)果顯示界面

1.4 數(shù)據(jù)格式

本軟件適用于重復(fù)測量設(shè)計(jì)資料的單、雙因素兩/多水平方差分析。對(duì)于單因素重復(fù)測量設(shè)計(jì)來說，時(shí)間(time)是唯一的一個(gè)重復(fù)測量的因素。受試對(duì)象(subject)個(gè)體是一個(gè)自變量因素。對(duì)雙因素設(shè)計(jì)來說， 時(shí)間(time)是一個(gè)重復(fù)測量的因素，分組(group)是另一個(gè)因素。單、雙因素重復(fù)測量設(shè)計(jì)下，各因素的水平可以是兩水平，也可以是多水平。在收集數(shù)據(jù)資料時(shí)，單因素?cái)?shù)據(jù)無分組(group)變量，所以該字段可以為空。2種設(shè)計(jì)類型資料的數(shù)據(jù)格式見表1。

表1 重復(fù)測量單、雙因素兩/多水平方差分析數(shù)據(jù)格式

2 實(shí)例分析與結(jié)果比較

2.1 數(shù)據(jù)來源

數(shù)據(jù)采用由顏虹、王彤主編的第5版《醫(yī)學(xué)統(tǒng)計(jì)學(xué)》教材[9]中的表12～13數(shù)據(jù)。某研究者欲比較治療厭食癥的3種不同成分藥物的效果,在已構(gòu)建的厭食癥模型大鼠中抽取15只成年雄性大鼠,隨機(jī)分為3組,每組飼料中添加一種藥物,藥物有效成分含量相同,連續(xù)記錄藥物治療前(T0)和治療后1 d、3 d、5 d和7 d的小鼠體質(zhì)量。數(shù)據(jù)中“觀測時(shí)間”(time)是一個(gè)重復(fù)測量因素，“實(shí)驗(yàn)分組”(group)因素為“藥物”，“大鼠編號(hào)”(subject)為個(gè)體標(biāo)記，因變量“體質(zhì)量”(value)是連續(xù)型的結(jié)果變量。將數(shù)據(jù)錄入表1對(duì)應(yīng)的格式化表格中，單擊界面中的“數(shù)據(jù)導(dǎo)入”按鈕，將整理好的數(shù)據(jù)導(dǎo)入到本軟件中。

2.2 數(shù)據(jù)分析

選擇dv(因變量)為value，wid(受試者編號(hào))為subject，within(受試內(nèi)自變量)為time,between(受試間自變量)為group，考慮自變量的主效應(yīng)和自變量間的交互效應(yīng)，type(類型)選III，detail選擇trpe。點(diǎn)擊“數(shù)據(jù)分析”按鈕，輸出統(tǒng)計(jì)分析結(jié)果。單因素重復(fù)測量方差分析時(shí)，因?yàn)闆]有分組因素和主體間交互作用，所以between字段要為空，同時(shí)type字段要選擇II，其余同雙因素方差分析。點(diǎn)擊結(jié)果輸出按鈕，輸出對(duì)應(yīng)結(jié)果。

2.3 軟件運(yùn)行結(jié)果描述

2.3.1 基本資料描述

統(tǒng)計(jì)資料的基本描述，可以點(diǎn)擊界面中的“基本描述”按鈕進(jìn)行結(jié)果和圖表的輸出，包括各因素的均數(shù)、標(biāo)準(zhǔn)差和QQ圖、箱式圖等。

2.3.2 重復(fù)測量方差分析

圖2界面顯示球形檢驗(yàn)、主體內(nèi)效應(yīng)及主體間效應(yīng)分析結(jié)果。結(jié)果顯示組別和時(shí)點(diǎn)2個(gè)因素對(duì)成年雄性大鼠體質(zhì)量的影響情況。首先，time主體內(nèi)作用和group*time的交互作用的球形校驗(yàn)結(jié)果均為P>0.05，服從球形假設(shè)。因此，劃線下面直接輸出方差分析的結(jié)果，即時(shí)點(diǎn)效應(yīng)可以影響體質(zhì)量的變化，F(xiàn)(4,48)=344.299,P<0.01，意味著大鼠體質(zhì)量會(huì)隨著時(shí)點(diǎn)出現(xiàn)變化；治療方式可以影響大鼠體質(zhì)量，F(xiàn)(2,12)=5.527,P<0.05；值得注意的是，時(shí)點(diǎn)和治療方式之間存在著交互效應(yīng)，F(xiàn)(8,48)=18.105,P<0.01，意味著大鼠體質(zhì)量隨著時(shí)點(diǎn)的變化趨勢(shì)會(huì)因?yàn)橹委煼绞降牟灰粯佣煌?。其次，基于以上的分析結(jié)果，時(shí)點(diǎn)和治療方式之間的交互效應(yīng)有統(tǒng)計(jì)學(xué)意義，則需要繼續(xù)考察時(shí)點(diǎn)和治療方式之間的單獨(dú)效應(yīng)。

2.3.3 事后多重比較

點(diǎn)擊“多重比較”按鈕，若存在交互作用，輸出時(shí)點(diǎn)和治療方式之間的單獨(dú)效應(yīng)；若不存在交互作用，輸出時(shí)點(diǎn)和治療方式之間的主效應(yīng)。均采用最小顯著差異法(least significant difference，LSD法)進(jìn)行兩兩比較。本實(shí)例研究顯示時(shí)點(diǎn)和治療方式之間存在交互作用，因此軟件會(huì)分別自動(dòng)輸出每個(gè)時(shí)間點(diǎn)下不同組別間兩兩比較結(jié)果，每個(gè)組別下不同時(shí)間點(diǎn)兩兩比較結(jié)果，并輸出相應(yīng)的統(tǒng)計(jì)圖表。部分顯示如圖3。

圖3 多重比較結(jié)果顯示界面

從多重比較結(jié)果可以得到：1)時(shí)點(diǎn)走勢(shì)的簡單效應(yīng)分析顯示，A組藥物治療后4個(gè)時(shí)間點(diǎn)與治療前比較，T2、T3、T4治療時(shí)間的體質(zhì)量有所上升(P<0.05)；B組治療后T1、T2、T3、T4時(shí)間點(diǎn)體質(zhì)量都增高(P<0.05)；C組治療后也是所有時(shí)間點(diǎn)體質(zhì)量都有增高,但T3、T4體質(zhì)量增高有統(tǒng)計(jì)學(xué)意義(P<0.05)。因此,A組藥物治療后體質(zhì)量的波動(dòng)小于B組和C組,B組藥物治療后體質(zhì)量穩(wěn)定上升。2)治療方式的簡單效應(yīng)分析顯示，治療前3組大鼠體質(zhì)量均無統(tǒng)計(jì)學(xué)差異；治療后1 d A組和C組之間的大鼠體質(zhì)量出現(xiàn)了統(tǒng)計(jì)學(xué)差異，A組體質(zhì)量低于C組；治療后3 d 3組大鼠體質(zhì)量均無統(tǒng)計(jì)學(xué)差異；治療后5 d A組、B組分別與C組之間的大鼠體質(zhì)量出現(xiàn)了統(tǒng)計(jì)學(xué)差異，A組體質(zhì)量低于B組和C組；治療后7 d A組分別與B組、C組之間的大鼠體質(zhì)量出現(xiàn)了統(tǒng)計(jì)學(xué)差異，A組體質(zhì)量低于B組和C組。

2.3.4 繪制數(shù)據(jù)輪廓圖

點(diǎn)擊界面中“繪圖輸出”按鈕，結(jié)果顯示區(qū)呈現(xiàn)數(shù)據(jù)輪廓圖。見圖4。

圖4 數(shù)據(jù)輪廓圖

2.4 軟件運(yùn)行結(jié)果驗(yàn)證

采用SPSS軟件進(jìn)行相應(yīng)的統(tǒng)計(jì)分析，結(jié)果表明， R語言程序與SPSS的統(tǒng)計(jì)分析結(jié)果一致，并與第5版《醫(yī)學(xué)統(tǒng)計(jì)學(xué)》教材201～207頁結(jié)果一致,本軟件的運(yùn)行結(jié)果準(zhǔn)確有效。。本軟件僅是調(diào)用R語言程序包，未做統(tǒng)計(jì)方法代碼的修改，所以本軟件結(jié)果即是R語言統(tǒng)計(jì)分析結(jié)果。

3 討論

本研究將C#語言和R語言相結(jié)合，開發(fā)了一款重復(fù)測量資料的單、雙因素兩/多水平方差分析軟件，該軟件操作簡單，能為臨床研究人員提供準(zhǔn)確的統(tǒng)計(jì)分析結(jié)果，提高科研統(tǒng)計(jì)效率。

3.1 軟件準(zhǔn)確性

通過結(jié)果比較，可以得知本軟件運(yùn)行結(jié)果與SPSS軟件以及手工計(jì)算結(jié)果均一致。說明統(tǒng)計(jì)分析結(jié)果準(zhǔn)確可靠，臨床醫(yī)生可以放心使用。

3.2 軟件便捷性

本軟件無需安裝，完全免費(fèi)，可視化界面簡潔，數(shù)據(jù)處理過程簡單，操作簡便，分析結(jié)果一鍵生成，直觀清晰，統(tǒng)計(jì)結(jié)果無需二次加工整理，論文可直接引用。SPSS和R語言軟件對(duì)于無操作經(jīng)驗(yàn)和編程經(jīng)驗(yàn)的初學(xué)者來說，在短時(shí)間內(nèi)熟練使用這些方法十分困難。本軟件的一鍵生成和友好界面能夠滿足臨床研究人員

的統(tǒng)計(jì)分析需求，具有良好的便捷性。

3.3 軟件智能性

使用SPSS和R語言完成重復(fù)測量方差分析，操作復(fù)雜，步驟繁多，結(jié)果需人為解讀及判斷。而本軟件具有良好的智能性，統(tǒng)計(jì)分析無需人為干預(yù)，自動(dòng)進(jìn)行重復(fù)測量數(shù)據(jù)的正態(tài)性、異常值的檢驗(yàn)，并根據(jù)檢驗(yàn)結(jié)果，智能選擇統(tǒng)計(jì)分析步驟并輸出結(jié)果。

綜上所述，基于C#與R語言的重復(fù)測量設(shè)計(jì)資料方差分析的自動(dòng)化實(shí)現(xiàn)軟件，分析結(jié)果準(zhǔn)確，使用便捷，具有良好的智能性，值得在臨床科研上推廣使用。目前，該軟件還存在一些不足，僅能實(shí)現(xiàn)單變量重復(fù)測量定量資料的方差分析。后續(xù)計(jì)劃陸續(xù)完善重復(fù)測量定量及定性資料其他統(tǒng)計(jì)方法的自動(dòng)化統(tǒng)計(jì)分析，更好地滿足臨床研究需求。