郭迎暄 陳 達 譚志軍 米白冰 黃昌可 廉恒麗△
【提 要】 目的 基于C#語言和R語言開發(fā)一款便于臨床醫(yī)生使用的軟件,用于實現(xiàn)重復(fù)測量設(shè)計定性資料統(tǒng)計分析的自動化,確保結(jié)果準(zhǔn)確、完整和規(guī)范,提高科研效率。方法 首先用C#語言將統(tǒng)計分析數(shù)據(jù)導(dǎo)入到軟件中,然后調(diào)用R語言命令腳本,完成重復(fù)測量設(shè)計定性資料的統(tǒng)計分析與結(jié)果輸出。最后通過實例,驗證本自動化實現(xiàn)軟件運行的效果。結(jié)果 基于C#語言和R語言相結(jié)合開發(fā)出的統(tǒng)計軟件運行結(jié)果與SPSS操作結(jié)果完全一致,統(tǒng)計分析結(jié)果包括了模型選擇、模型模擬結(jié)果和模擬結(jié)果預(yù)測三部分。其操作簡便,結(jié)果自動化呈現(xiàn),減少了結(jié)果判斷與模型選擇的錯誤,確保了結(jié)果的準(zhǔn)確性和規(guī)范性。結(jié)論 該軟件能夠自動實現(xiàn)重復(fù)測量設(shè)計定性資料的統(tǒng)計分析,值得在臨床上推廣使用。
在臨床研究中,常常收集多次重復(fù)測量設(shè)計的結(jié)果變量(因變量)及其影響結(jié)果變量的數(shù)據(jù)(自變量,如性別、年齡、治療方法、疾病等)。重復(fù)測量資料是對同一受試者的同一觀察指標(biāo)在不同時間點上進行多次測量(≥3次)所得的資料,通常用于分析該項觀察指標(biāo)在不同時間點上的變化特點[1]。數(shù)據(jù)類型可分為定量數(shù)據(jù)和定性數(shù)據(jù),其中,定性數(shù)據(jù)是指測得的僅反映某一方面性質(zhì)的指標(biāo),并不能用具體的數(shù)值表示。定性重復(fù)測量資料根據(jù)因變量類型可以分為離散型變量、分類變量和等級變量,分析模型可包含固定效應(yīng)、隨機效應(yīng)或者混合效應(yīng)。對于臨床醫(yī)生而言,選擇正確的模型和統(tǒng)計方法進行數(shù)據(jù)分析,得出正確的結(jié)果和結(jié)論,具有一定難度且需要花費很大精力。因此,迫切需要開發(fā)一款能夠?qū)崿F(xiàn)定性重復(fù)測量數(shù)據(jù)統(tǒng)計分析自動化的軟件。鑒于此,本研究將結(jié)合C#語言和R語言進行軟件設(shè)計和開發(fā),并通過實例展示軟件在臨床研究中的應(yīng)用。
C#是微軟公司推出的一種面向?qū)ο蟮木幊陶Z言,具有可視化操作和高效率運行的特點,其支持快速地編寫各種基于Microsoft.NET平臺的應(yīng)用程序[2]。R是用于統(tǒng)計分析和統(tǒng)計繪圖的語言和操作環(huán)境,是一個免費軟件,擁有各種各樣的R統(tǒng)計分析包,通過這些R語言包,可以進行教育、醫(yī)療、可視化、統(tǒng)計學(xué)、人工智能等方面應(yīng)用。
本研究以C#開發(fā)平臺Microsoft Visual Studio Enter Prise 2019為基礎(chǔ),結(jié)合R軟件lme4程序包和geepack程序包實現(xiàn)重復(fù)測量設(shè)計定性資料的統(tǒng)計分析,其中,R軟件環(huán)境采用R(v.3.6.1)、R studio(1.2.5001)版本。C#對R語言的調(diào)用方法有兩種,一種是通過R語言的COM接口,直接和R語言進行交互;一種是通過RDotNet.dll與R語言進行交互。本軟件通過后者與R語言進行交互,首先,開發(fā)環(huán)境需要先安裝.NET Framework4和R.dll;然后,在C#程序中添加對RDotNet.dll項目的引用;最后,利用REngine對象的方法Evaluate、CreateNumericVector和CreateCharacterMatrix等創(chuàng)建R向量和矩陣,實現(xiàn)C#對R語言函數(shù)的調(diào)用?;趦煞N語言的特點,使得開發(fā)界面友好、操作簡便、自動化運行的統(tǒng)計軟件成為可能。
軟件根據(jù)研究目的和研究設(shè)計、因變量類型與分布、因變量與自變量關(guān)系,選擇合適的統(tǒng)計分析方法進行數(shù)據(jù)分析。不同的統(tǒng)計分析方法涉及的參數(shù)不同,對應(yīng)的界面也會有略微調(diào)整。本文用定性資料二分類GEE模型和無序多分類GLMMs模型為實例,來說明軟件如何實現(xiàn)數(shù)據(jù)導(dǎo)入和自動化輸出統(tǒng)計分析的結(jié)果。統(tǒng)計分析與模型選擇流程見圖1。
圖1 統(tǒng)計分析與模型選擇流程圖
該軟件界面(如圖2)左側(cè)紅框區(qū)為菜單欄,根據(jù)不同的資料類型選擇適用的統(tǒng)計分析方法。右側(cè)藍框區(qū)為數(shù)據(jù)導(dǎo)入格式示例區(qū),可進行數(shù)據(jù)導(dǎo)入。右側(cè)中部綠框區(qū)為參數(shù)設(shè)置區(qū),分別選擇因變量和自變量。右側(cè)下方紫框區(qū)為結(jié)果顯示區(qū),根據(jù)不同的統(tǒng)計分析方法,顯示相應(yīng)的結(jié)果。
圖2 軟件界面圖
為了測試該軟件的可靠性與有效性,本文介紹二分類重復(fù)測量資料廣義估計方程與多分類重復(fù)測量資料廣義線性混合效應(yīng)模型在該軟件的自動化實現(xiàn)。
廣義估計方程(generalized estimating equation,GEE)是Liang和Zeger在廣義線性模型和擬似然方法的基礎(chǔ)上提出的一種分析縱向數(shù)據(jù)的方法。GEE可以處理有缺失值的資料,允許每個觀察對象的觀察次數(shù)不同,觀察時間間隔亦可不同。廣義估計方程應(yīng)用條件較寬,除了正態(tài)分布,可以利用連接函數(shù)將高斯分布、二項分布、多項分布、Poisson分布、Gamma分布等多種分布的因變量擬合為相應(yīng)的統(tǒng)計模型,解決了重復(fù)測量數(shù)據(jù)非獨立性問題,可得到穩(wěn)健的參數(shù),最大程度減少測量數(shù)據(jù)的有效信息損失。
假設(shè)yij為第i個觀察對象的第j個觀察值(i=1,…,n;j=l,…,p),Xij(Xij1,Xij2,…,Xijm)為相應(yīng)的自變量向量。各觀察對象是獨立的,但同一觀察對象內(nèi)的觀察值間存在相關(guān)。模型的基本構(gòu)成如下:
(1)建立yij與各自變量Xij(Xij1,Xij2,…,Xijm)之間的函數(shù)關(guān)系
E(yij)=uijg(uij)=β0+β1Xij1+β2Xij2+…+βmXijm
(1)
其中g(shù)(uij)為聯(lián)結(jié)函數(shù),可根據(jù)數(shù)據(jù)類型選取合適的聯(lián)結(jié)函數(shù)。
(2)建立yij的方差與平均值之間的函數(shù)關(guān)系
Var(yij)=v(uij)·φ
(2)
v(uij)為已知方差函數(shù),φ為尺度參數(shù),表示y的方差不能被v(uij)解釋的部分。
(3)對yi=(yi1,…,yip)′選擇一個p×p維作業(yè)相關(guān)矩陣Ri(α),構(gòu)造廣義估計方程如下:
(3)
GEE的特點是采用實際計算得到的殘差函數(shù),作簡單回歸從而獲得作業(yè)相關(guān)矩陣。相關(guān)矩陣存在多種結(jié)構(gòu)(等相關(guān)結(jié)構(gòu)、相鄰相關(guān)結(jié)構(gòu)、自相關(guān)結(jié)構(gòu)、不確定型相關(guān)結(jié)構(gòu)、獨立相關(guān)結(jié)構(gòu)),模型擬合的好壞可以通過QIC判別準(zhǔn)則做出判斷[3]。通過QIC大小決定合適的大小相關(guān)矩陣,在同一模型中QIC值越小模型越合適[4]。對于GEE算法而言,即使對相關(guān)矩陣的結(jié)構(gòu)選擇不當(dāng),也能得到有關(guān)結(jié)果變量的回歸系數(shù)及其方差的一致性估計值[5]。當(dāng)樣本含量較大時,因?qū)ψ鳂I(yè)相關(guān)矩陣的選擇不當(dāng)而引起的效率損失可以忽略不計。
(1)背景資料
本研究為一項單中心、前瞻性干預(yù)性研究,觀察兩組不同治療方案的治療效果。研究因素為組別,即單純西醫(yī)治療組(90例)和中西醫(yī)結(jié)合治療組(90例),分別于治療后1周、1月、3月共3個時間點觀測記錄治療效果。
表1 研究變量說明
(2)R程序代碼
#原始excel數(shù)據(jù)導(dǎo)入
library(readxl)
#數(shù)據(jù)讀取操作
data<- read_excel(file.choose())
data$GROUP<- factor(data$GROUP)
data$TIME<-factor(data$TIME)
data$AGE<-factor(data$AGE)
data$ID<-factor(data$ID)
#模型適配
library(geepack)
fit1<- geeglm(EFFECT ~ GROUP + AGE + TIME,id=ID,data=data,corstr=“ar1”,family=‘binomial’)
fit2<-geeglm(EFFECT ~ GROUP + AGE + TIME,id=ID,data=data,corstr=“exchangeable”,family=‘binomial’)
fit3<- geeglm(EFFECT ~ GROUP + AGE + TIME,id=ID,data=data,corstr=“independence”,family=‘binomial’)
sapply(list(fit1,fit2,fit3),QIC)
#比較幾種模型的QIC值,選擇QIC最小值模型進行統(tǒng)計分析與結(jié)果輸出
coef(summary(fit3))
#編寫GEE95%可信區(qū)間函數(shù)
confint.geeglm<- function(object,parm,level=0.95,…){
cc<- coef(summary(object))
mult<- qnorm((1+level)/2)
citab<- with(as.data.frame(cc),
cbind(lwr=Estimate-mult*Std.err,
upr=Estimate+mult*Std.err))
rownames(citab)<- rownames(cc)
citab[parm,]
}
confint.geeglm(fit3)
#結(jié)果預(yù)測
pred=predict.glm(fit3,type=“response”,newdata=data)
predict=ifelse(pred>0.5,1,0)
data$predict=predict
library(vcd)
addmargins(table(data$PREDICTEDVALUE,data$EFFECT))
(3)結(jié)果展示與表達
本研究采用廣義估計方程研究本案例的二分類重復(fù)測量的數(shù)據(jù),運算結(jié)果如下:
a.模型選擇
根據(jù)擬似然信息準(zhǔn)則(QIC)統(tǒng)計量進行模型選擇,結(jié)果表明,independence模型QIC值最小,若遇到ra1、exchangeable指標(biāo)QIC值與independence指標(biāo)QIC值相同時,以ra1為最優(yōu)。
b.模型擬合結(jié)果
圖3 廣義估計方程參數(shù)估計結(jié)果
①從圖3可以清晰看到,組間比較結(jié)果,單純西醫(yī)的療效顯著低于中西醫(yī)結(jié)合,Wald卡方=9.701,B=-1.735<0且P<0.01,更進一步,單純西醫(yī)的有效率是中西醫(yī)結(jié)合的exp(-1.735)=17.63%;
②基線數(shù)據(jù)影響結(jié)果:年齡不能顯著影響有效率,P值均大于0.05;
③重復(fù)測量時間比較結(jié)果:治療后1月的有效率顯著高于治療后1周,Wald卡方=3.894,B=0.818>0且P=0.048<0.05,更進一步,1月的有效率是1周的2.265倍;治療后3月的有效率顯著高于治療后1周,Wald卡方=8.819,B=1.502>0且P=0.003<0.05,更進一步,3月的有效率是1周的4.490倍。
c.模型預(yù)測準(zhǔn)確率
更進一步,我們需要繼續(xù)考察以上模型的準(zhǔn)確率。模型預(yù)測的準(zhǔn)確率為:
廣義線性混合效應(yīng)模型是廣義線性模型和一般線性混合效應(yīng)模型的擴展,是在廣義線性固定效應(yīng)模型的基礎(chǔ)上引入隨機效應(yīng),在隨機效應(yīng)滿足正態(tài)分布的前提下,因變量可以是指數(shù)家族中的任一分布,指數(shù)家族可有許多基本的離散分布(包括二項分布、泊松分布和負(fù)二項式正態(tài)分布等)和連續(xù)分布(正態(tài)分布、beta分布和χ2分布等)組成,當(dāng)隨機效應(yīng)不存在時,廣義線性混合效應(yīng)模型就退化為廣義線性模型[6]。廣義線性混合效應(yīng)模型的自變量可以是分類或連續(xù)的,可以處理多個隨機效應(yīng),建模靈活,且同樣可以用于非均衡數(shù)據(jù),能較好處理含有缺失值的資料。
(1)模型框架:GLMMS利用逆連接函數(shù)來構(gòu)建線性預(yù)測值與條件均數(shù)關(guān)系的基本模型:
Y=μ+ε
μ=g-1(Xβ+Zγ)
式中,Y:n×l維觀測向量;μ:觀測的期望(均數(shù))向量;g(·):可微單調(diào)連接函數(shù),g-1(·):g(·)的轉(zhuǎn)置;X和Z分別是固定效應(yīng)和隨機效應(yīng)的設(shè)計矩陣,X:n×p維矩陣,Z:n×r維矩陣;β和γ分別是模型的固定效應(yīng)和隨機效應(yīng)的參數(shù)向量,隨機效應(yīng)γ應(yīng)滿足均數(shù)為0,方差矩陣為G的正態(tài)分布,γ~N(0,G),Var(Y)=G;殘差ε~N(0,R),var(ε)=R,R為殘差協(xié)方差矩陣[7]。
對于有序多分類結(jié)局測量,其連接函數(shù)為累積logit函數(shù)(cumulative logit function),采用多層累積logistic回歸模型來擬合數(shù)據(jù),模型可表達為:
Y=μ+ε
γ~N(0,G)var(ε)=R
其中,μ:多項式概率分布期望向量,有n個延伸的觀測。假設(shè)有4個分類,可以記作:μ=(μ11,μ12,μ13,…,μn1,μn2,μn3),μij:觀測i在分類j的概率。
(2)參數(shù)估計:GLMMS估計的最大似然目的是將如下求積似然函數(shù)(integrated likelihood function)最大化:
其中β為固定效應(yīng),θ為未知的方差/協(xié)方差參數(shù),f(Y|u)為隨機效應(yīng)u條件下的結(jié)局測量分布函數(shù),p(u)為隨機效應(yīng)的分布函數(shù)。此積分似然函數(shù)必須近似估計[8]。
(3)背景資料
本研究為一項隨訪調(diào)查研究,觀察醫(yī)學(xué)本科畢業(yè)生在剛畢業(yè)、畢業(yè)后3年和畢業(yè)后6年的去向選擇。研究因素為生源地、性別、學(xué)習(xí)成績和畢業(yè)時間,因變量為去向選擇。
表2 研究中變量說明
(4)R語言程序代碼
library(readxl)
data<- read_excel(file.choose())
data$TIME<-factor(data$TIME)
data$ID<-factor(data$ID)
data$SEX<-factor(data$SEX)
data$SCORE<-factor(data$SCORE)
data$ADDRESS<-factor(data$ADDRESS)
library(lme4)
glmms1<- glmer(DIRECTION~ SEX + SCORE + ADDRESS + TIME +(1|ID),
data=data,family=‘Gamma’)
glmms2<- glmer(DIRECTION~ SEX + SCORE + ADDRESS + TIME +(1|ID),
data=data,family=‘inverse.gaussian’)
sapply(list(glmms1,glmms2),AIC)
sapply(list(glmms1,glmms2),BIC)
coef(summary(glmms2))
confint.glmer(glmms2)
#結(jié)果預(yù)測
pred=fitted(glmms2)
pred=ifelse(pred>2.5,3,pred)
pred=ifelse(pred<2.5 & pred>1.5,2,pred)
pred=ifelse(pred<1.5,1,pred)
data$PREDICTEDVALUE=pred
library(vcd)
table(data$PREDICTEDVALUE,data$DIRECTION)
addmargins(table(data$PREDICTEDVALUE,data$DIRECTION))
(5)結(jié)果展示與表達
本研究采用廣義線性混合效應(yīng)模型研究本案例的無序多分類重復(fù)測量的數(shù)據(jù),運算結(jié)果如下。
a.模型選擇
根據(jù)赤池信息準(zhǔn)則(Akaike information criterion,AIC)和貝葉斯信息準(zhǔn)則(Bayesian information criterion,BIC)選擇最優(yōu)模型,結(jié)果表明,正態(tài)反高斯先驗?zāi)P偷腁IC、BIC值均最小,選擇正態(tài)反高斯先驗?zāi)P蜑樽顑?yōu)模型。
b.模型擬合結(jié)果
圖4 廣義線性混合效應(yīng)模型參數(shù)估計結(jié)果
從圖4可以清晰看到:
①基線數(shù)據(jù)影響結(jié)果:性別、生源地不會影響醫(yī)學(xué)生本科畢業(yè)后的選擇,P值全部大于0.05。
②考試成績的影響:學(xué)習(xí)成績會影響醫(yī)學(xué)生本科畢業(yè)后的選擇,學(xué)習(xí)成績靠后的畢業(yè)生選擇繼續(xù)深造而不選擇醫(yī)生和醫(yī)藥公司的可能性明顯低于學(xué)習(xí)成績靠前的畢業(yè)生;學(xué)習(xí)成績靠后的畢業(yè)生選擇繼續(xù)深造的可能性僅僅只有學(xué)習(xí)成績靠前的53.08%(P<0.01);學(xué)習(xí)成績中等的畢業(yè)生選擇繼續(xù)深造而不選擇醫(yī)生和醫(yī)藥公司的可能性明顯低于學(xué)習(xí)成績靠后的畢業(yè)生;學(xué)習(xí)成績中等的畢業(yè)生選擇繼續(xù)深造的可能性僅僅只有學(xué)習(xí)成績靠后的59.87%(P<0.01);基于此,學(xué)習(xí)成績好的傾向于繼續(xù)深造,學(xué)習(xí)成績中等的傾向于醫(yī)生,而成績較差的傾向于醫(yī)藥公司。
③畢業(yè)時間的影響:畢業(yè)后3年選擇繼續(xù)深造而不是醫(yī)生和醫(yī)藥公司的可能性僅僅只有剛畢業(yè)的82.53%(P<0.01),畢業(yè)后6年與剛畢業(yè)的去向比較更傾向于醫(yī)生和醫(yī)藥公司(P<0.01)。
c.模型預(yù)測準(zhǔn)確率
通過前面R語言程序與SPSS的統(tǒng)計分析結(jié)果比較,可以得出兩種統(tǒng)計方式結(jié)果一致。本軟件僅是調(diào)用R語言程序包,未做統(tǒng)計方法代碼的修改,所以本軟件結(jié)果即是R語言統(tǒng)計分析結(jié)果。因此,本軟件的運行結(jié)果準(zhǔn)確有效。
目前,國內(nèi)有公司推出了在線數(shù)據(jù)科學(xué)分析平臺[9](SPSSAU)和易侕軟件[10](EmpoertStats),能夠?qū)崿F(xiàn)自動化統(tǒng)計分析與結(jié)果輸出功能,同樣具有操作簡便、結(jié)果顯示清晰的優(yōu)點;但不足之處是費用比較高(2588元/年,200元/月等),且目前無法實現(xiàn)廣義估計方程和廣義線性混合模型的統(tǒng)計分析。
本研究基于C#語言和R語言開發(fā)了一套針對廣義估計方程和廣義線性混合模型的統(tǒng)計分析軟件,該軟件具有數(shù)據(jù)導(dǎo)入、統(tǒng)計分析、模型選擇、分析結(jié)果和結(jié)果預(yù)測等功能,實現(xiàn)對臨床重復(fù)測量定性資料的自動化統(tǒng)計分析。實際使用中,只需選擇因變量和自變量就能自動獲取統(tǒng)計分析的相關(guān)結(jié)果,且結(jié)果與SPSS軟件統(tǒng)計分析結(jié)果一致[11-12]。該軟件完全免費,安裝后醫(yī)生可以根據(jù)收集資料的性質(zhì)與分析目的,選擇適合的統(tǒng)計分析方法和統(tǒng)計圖表,只需簡單的可視化步驟,便輸出統(tǒng)計分析結(jié)果與表達,減輕了醫(yī)生的科研統(tǒng)計壓力。不足之處是,前期進行了重復(fù)測量設(shè)計定量資料的自動化實現(xiàn),尚未進行重復(fù)測量設(shè)計生存資料的統(tǒng)計分析,這部分將在以后的研究中進一步探討。