基于C#與R語言的重復測量設計定性資料統計分析的自動化實現研究*

郭迎暄 陳 達 譚志軍 米白冰 黃昌可 廉恒麗△

【提 要】 目的 基于C#語言和R語言開發一款便于臨床醫生使用的軟件，用于實現重復測量設計定性資料統計分析的自動化，確保結果準確、完整和規范，提高科研效率。方法 首先用C#語言將統計分析數據導入到軟件中，然后調用R語言命令腳本，完成重復測量設計定性資料的統計分析與結果輸出。最后通過實例，驗證本自動化實現軟件運行的效果。結果 基于C#語言和R語言相結合開發出的統計軟件運行結果與SPSS操作結果完全一致，統計分析結果包括了模型選擇、模型模擬結果和模擬結果預測三部分。其操作簡便，結果自動化呈現，減少了結果判斷與模型選擇的錯誤，確保了結果的準確性和規范性。結論 該軟件能夠自動實現重復測量設計定性資料的統計分析，值得在臨床上推廣使用。

在臨床研究中，常常收集多次重復測量設計的結果變量(因變量)及其影響結果變量的數據(自變量，如性別、年齡、治療方法、疾病等)。重復測量資料是對同一受試者的同一觀察指標在不同時間點上進行多次測量(≥3次)所得的資料，通常用于分析該項觀察指標在不同時間點上的變化特點[1]。數據類型可分為定量數據和定性數據，其中，定性數據是指測得的僅反映某一方面性質的指標，并不能用具體的數值表示。定性重復測量資料根據因變量類型可以分為離散型變量、分類變量和等級變量，分析模型可包含固定效應、隨機效應或者混合效應。對于臨床醫生而言，選擇正確的模型和統計方法進行數據分析，得出正確的結果和結論，具有一定難度且需要花費很大精力。因此，迫切需要開發一款能夠實現定性重復測量數據統計分析自動化的軟件。鑒于此，本研究將結合C#語言和R語言進行軟件設計和開發，并通過實例展示軟件在臨床研究中的應用。

軟件設計與實現

1.軟件設計原理

C#是微軟公司推出的一種面向對象的編程語言，具有可視化操作和高效率運行的特點，其支持快速地編寫各種基于Microsoft.NET平臺的應用程序[2]。R是用于統計分析和統計繪圖的語言和操作環境，是一個免費軟件，擁有各種各樣的R統計分析包，通過這些R語言包，可以進行教育、醫療、可視化、統計學、人工智能等方面應用。

本研究以C#開發平臺Microsoft Visual Studio Enter Prise 2019為基礎，結合R軟件lme4程序包和geepack程序包實現重復測量設計定性資料的統計分析，其中，R軟件環境采用R(v.3.6.1)、R studio(1.2.5001)版本。C#對R語言的調用方法有兩種，一種是通過R語言的COM接口，直接和R語言進行交互；一種是通過RDotNet.dll與R語言進行交互。本軟件通過后者與R語言進行交互，首先，開發環境需要先安裝.NET Framework4和R.dll；然后，在C#程序中添加對RDotNet.dll項目的引用；最后，利用REngine對象的方法Evaluate、CreateNumericVector和CreateCharacterMatrix等創建R向量和矩陣，實現C#對R語言函數的調用。基于兩種語言的特點，使得開發界面友好、操作簡便、自動化運行的統計軟件成為可能。

2.軟件統計分析流程

軟件根據研究目的和研究設計、因變量類型與分布、因變量與自變量關系，選擇合適的統計分析方法進行數據分析。不同的統計分析方法涉及的參數不同，對應的界面也會有略微調整。本文用定性資料二分類GEE模型和無序多分類GLMMs模型為實例，來說明軟件如何實現數據導入和自動化輸出統計分析的結果。統計分析與模型選擇流程見圖1。

圖1 統計分析與模型選擇流程圖

3.軟件界面設計與功能介紹

該軟件界面(如圖2)左側紅框區為菜單欄，根據不同的資料類型選擇適用的統計分析方法。右側藍框區為數據導入格式示例區，可進行數據導入。右側中部綠框區為參數設置區，分別選擇因變量和自變量。右側下方紫框區為結果顯示區，根據不同的統計分析方法，顯示相應的結果。

圖2 軟件界面圖

軟件驗證

為了測試該軟件的可靠性與有效性，本文介紹二分類重復測量資料廣義估計方程與多分類重復測量資料廣義線性混合效應模型在該軟件的自動化實現。

1.廣義估計方程

廣義估計方程(generalized estimating equation，GEE)是Liang和Zeger在廣義線性模型和擬似然方法的基礎上提出的一種分析縱向數據的方法。GEE可以處理有缺失值的資料，允許每個觀察對象的觀察次數不同，觀察時間間隔亦可不同。廣義估計方程應用條件較寬，除了正態分布，可以利用連接函數將高斯分布、二項分布、多項分布、Poisson分布、Gamma分布等多種分布的因變量擬合為相應的統計模型，解決了重復測量數據非獨立性問題，可得到穩健的參數，最大程度減少測量數據的有效信息損失。

假設yij為第i個觀察對象的第j個觀察值(i=1，…，n；j=l，…，p)，Xij(Xij1，Xij2，…，Xijm)為相應的自變量向量。各觀察對象是獨立的，但同一觀察對象內的觀察值間存在相關。模型的基本構成如下：

(1)建立yij與各自變量Xij(Xij1，Xij2，…，Xijm)之間的函數關系

E(yij)=uijg(uij)=β0+β1Xij1+β2Xij2+…+βmXijm

(1)

其中g(uij)為聯結函數，可根據數據類型選取合適的聯結函數。

(2)建立yij的方差與平均值之間的函數關系

Var(yij)=v(uij)·φ

(2)

v(uij)為已知方差函數，φ為尺度參數，表示y的方差不能被v(uij)解釋的部分。

(3)對yi=(yi1，…，yip)′選擇一個p×p維作業相關矩陣Ri(α)，構造廣義估計方程如下：

(3)

GEE的特點是采用實際計算得到的殘差函數，作簡單回歸從而獲得作業相關矩陣。相關矩陣存在多種結構(等相關結構、相鄰相關結構、自相關結構、不確定型相關結構、獨立相關結構)，模型擬合的好壞可以通過QIC判別準則做出判斷[3]。通過QIC大小決定合適的大小相關矩陣，在同一模型中QIC值越小模型越合適[4]。對于GEE算法而言，即使對相關矩陣的結構選擇不當，也能得到有關結果變量的回歸系數及其方差的一致性估計值[5]。當樣本含量較大時，因對作業相關矩陣的選擇不當而引起的效率損失可以忽略不計。

(1)背景資料

本研究為一項單中心、前瞻性干預性研究，觀察兩組不同治療方案的治療效果。研究因素為組別，即單純西醫治療組(90例)和中西醫結合治療組(90例)，分別于治療后1周、1月、3月共3個時間點觀測記錄治療效果。

表1 研究變量說明

(2)R程序代碼

#原始excel數據導入

library(readxl)

#數據讀取操作

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最小值模型進行統計分析與結果輸出

coef(summary(fit3))

#編寫GEE95%可信區間函數

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)

#結果預測

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)結果展示與表達

本研究采用廣義估計方程研究本案例的二分類重復測量的數據，運算結果如下：

a.模型選擇

根據擬似然信息準則(QIC)統計量進行模型選擇，結果表明，independence模型QIC值最小，若遇到ra1、exchangeable指標QIC值與independence指標QIC值相同時，以ra1為最優。

b.模型擬合結果

圖3 廣義估計方程參數估計結果

①從圖3可以清晰看到，組間比較結果，單純西醫的療效顯著低于中西醫結合，Wald卡方=9.701，B=-1.735<0且P<0.01，更進一步，單純西醫的有效率是中西醫結合的exp(-1.735)=17.63%；

②基線數據影響結果：年齡不能顯著影響有效率，P值均大于0.05；

③重復測量時間比較結果：治療后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.模型預測準確率

更進一步，我們需要繼續考察以上模型的準確率。模型預測的準確率為：

2.廣義線性混合效應模型

廣義線性混合效應模型是廣義線性模型和一般線性混合效應模型的擴展，是在廣義線性固定效應模型的基礎上引入隨機效應，在隨機效應滿足正態分布的前提下，因變量可以是指數家族中的任一分布，指數家族可有許多基本的離散分布(包括二項分布、泊松分布和負二項式正態分布等)和連續分布(正態分布、beta分布和χ2分布等)組成，當隨機效應不存在時，廣義線性混合效應模型就退化為廣義線性模型[6]。廣義線性混合效應模型的自變量可以是分類或連續的，可以處理多個隨機效應，建模靈活，且同樣可以用于非均衡數據，能較好處理含有缺失值的資料。

(1)模型框架：GLMMS利用逆連接函數來構建線性預測值與條件均數關系的基本模型：

Y=μ+ε

μ=g-1(Xβ+Zγ)

式中，Y：n×l維觀測向量；μ：觀測的期望(均數)向量；g(·)：可微單調連接函數，g-1(·)：g(·)的轉置；X和Z分別是固定效應和隨機效應的設計矩陣，X：n×p維矩陣，Z：n×r維矩陣；β和γ分別是模型的固定效應和隨機效應的參數向量，隨機效應γ應滿足均數為0，方差矩陣為G的正態分布，γ～N(0，G)，Var(Y)=G；殘差ε～N(0，R)，var(ε)=R，R為殘差協方差矩陣[7]。

對于有序多分類結局測量，其連接函數為累積logit函數(cumulative logit function)，采用多層累積logistic回歸模型來擬合數據，模型可表達為：

Y=μ+ε

γ～N(0，G)var(ε)=R

其中，μ：多項式概率分布期望向量，有n個延伸的觀測。假設有4個分類，可以記作：μ=(μ11，μ12，μ13，…，μn1，μn2，μn3)，μij：觀測i在分類j的概率。

(2)參數估計：GLMMS估計的最大似然目的是將如下求積似然函數(integrated likelihood function)最大化：

其中β為固定效應，θ為未知的方差/協方差參數，f(Y|u)為隨機效應u條件下的結局測量分布函數，p(u)為隨機效應的分布函數。此積分似然函數必須近似估計[8]。

(3)背景資料

本研究為一項隨訪調查研究，觀察醫學本科畢業生在剛畢業、畢業后3年和畢業后6年的去向選擇。研究因素為生源地、性別、學習成績和畢業時間，因變量為去向選擇。

表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)

#結果預測

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)結果展示與表達

本研究采用廣義線性混合效應模型研究本案例的無序多分類重復測量的數據，運算結果如下。

a.模型選擇

根據赤池信息準則(Akaike information criterion，AIC)和貝葉斯信息準則(Bayesian information criterion，BIC)選擇最優模型，結果表明，正態反高斯先驗模型的AIC、BIC值均最小，選擇正態反高斯先驗模型為最優模型。

b.模型擬合結果

圖4 廣義線性混合效應模型參數估計結果

從圖4可以清晰看到：

①基線數據影響結果：性別、生源地不會影響醫學生本科畢業后的選擇，P值全部大于0.05。

②考試成績的影響：學習成績會影響醫學生本科畢業后的選擇，學習成績靠后的畢業生選擇繼續深造而不選擇醫生和醫藥公司的可能性明顯低于學習成績靠前的畢業生；學習成績靠后的畢業生選擇繼續深造的可能性僅僅只有學習成績靠前的53.08%(P<0.01)；學習成績中等的畢業生選擇繼續深造而不選擇醫生和醫藥公司的可能性明顯低于學習成績靠后的畢業生；學習成績中等的畢業生選擇繼續深造的可能性僅僅只有學習成績靠后的59.87%(P<0.01)；基于此，學習成績好的傾向于繼續深造，學習成績中等的傾向于醫生，而成績較差的傾向于醫藥公司。

③畢業時間的影響：畢業后3年選擇繼續深造而不是醫生和醫藥公司的可能性僅僅只有剛畢業的82.53%(P<0.01)，畢業后6年與剛畢業的去向比較更傾向于醫生和醫藥公司(P<0.01)。

c.模型預測準確率

3.結果比較

通過前面R語言程序與SPSS的統計分析結果比較，可以得出兩種統計方式結果一致。本軟件僅是調用R語言程序包，未做統計方法代碼的修改，所以本軟件結果即是R語言統計分析結果。因此，本軟件的運行結果準確有效。

討 論

目前，國內有公司推出了在線數據科學分析平臺[9](SPSSAU)和易侕軟件[10](EmpoertStats)，能夠實現自動化統計分析與結果輸出功能，同樣具有操作簡便、結果顯示清晰的優點；但不足之處是費用比較高(2588元/年，200元/月等)，且目前無法實現廣義估計方程和廣義線性混合模型的統計分析。

本研究基于C#語言和R語言開發了一套針對廣義估計方程和廣義線性混合模型的統計分析軟件，該軟件具有數據導入、統計分析、模型選擇、分析結果和結果預測等功能，實現對臨床重復測量定性資料的自動化統計分析。實際使用中，只需選擇因變量和自變量就能自動獲取統計分析的相關結果，且結果與SPSS軟件統計分析結果一致[11-12]。該軟件完全免費，安裝后醫生可以根據收集資料的性質與分析目的，選擇適合的統計分析方法和統計圖表，只需簡單的可視化步驟，便輸出統計分析結果與表達，減輕了醫生的科研統計壓力。不足之處是，前期進行了重復測量設計定量資料的自動化實現，尚未進行重復測量設計生存資料的統計分析，這部分將在以后的研究中進一步探討。