基于C#与R语言的重复测量设计资料方差分析的自动化实现

郭迎暄 陈 达 张双喆 黄昌可 廉恒丽

1 温州医科大学附属眼视光医院杭州院区，330020 浙江 杭州;2 温州医科大学附属眼视光医院，325027 浙江 温州

医学临床研究中常收集多次重复测量数据，一是可以减少个体间的误差，二是为了节约样本量。当重复测量次数m≥3时,称为重复测量设计(repeated measurement design)或重复测量数据(repeated measurement data)。该数据可用于分析被观察指标在不同时间点上的变化趋势与治疗效果[1]。由于多次重复测量的数据来自于同一个体，这些数据之间存在相关性，与方差分析的基本假设——数据独立性相违背，若采用常规的方差分析，则会增加犯I类错误的概率[2]，因此需要采用重复测量方差分析。重复测量方差分析可用于以下目的：比较组间有无差异；比较各时间点间有无差异；比较组间的时间变化趋势有无差异(即随时间变化的多条曲线是否平行)[3]。数据需满足以下条件：因变量没有明显异常值；因变量需服从近似正态分布；因变量的方差协方差矩阵相等，也称满足球形假设。近年来，医学研究越来越离不开统计方法的使用，而每种统计方法都有其适用条件,医学杂志发表的论文存在不同程度的统计错误,影响研究结论。为了减少这一现象,提高论文水平, 提高医生科研效率，迫切需要一款界面更加友好，操作更加简单，结果更加直观的统计自动化软件。

R软件是目前医学研究领域广泛使用的免费开源软件，具有全面完善的统计分析功能和数量庞大的拓展包，具备跨平台运行和兼容各类统计软件等优点[4]。近年来，在医疗卫生各领域受到越来越多统计学者的使用和认可[5-7]。C#是一种面向对象的编程语言，利用C#.NET 开发的系统具有界面友好、执行速度快、易维护和升级等优点[8]。通过RDotNet组件，C#可以直接访问R语言程序包。2种语言的这些特点，使得C#语言和R语言联合编程具有可行性。为了让临床医生轻松又快捷地完成临床科研中重复测量设计资料的方差分析，本研究基于2种语言开发了一款简单易用的软件，可自动完成重复测量设计资料的单、双因素两/多水平方差分析。

1 C#语言与R语言联合编程方法

1.1 技术路线

1)安装3.6.3版本的R语言运行平台RStudio和16.7.6版本的Microsoft Visual Studio Enterprise 2019 C#开发平台用来做调试使用。在本研究中，主要使用readxl、ez、psych、RVAideMemoire、agricolae、PairedData、ggplot 2共7个分析包编写R语言统计分析脚本，在C#项目中引用RDotNet.dll，调用编写好的R语言统计脚本，实现重复测量设计资料方差分析。

2)自动化输出统计分析结果。首先判断受试者内因素的各个水平，因变量有没有极端异常值；因变量是否服从近似正态分布；因变量的方差协方差矩阵相等(根据ezANOVA包统计分析结果中的Mauchly's球性检验结果进行判断)。其中，若满足球形假定(P>0.05)，直接输出方差分析的结果；若不满足球形假定(P<0.05)，需要参考GGe(Greenhouse-Geisser epsilon)或者HFe(Huynh-Feldt epsilon)对df进行校正。当GGe和HFe对应的Epsilon值均>0.75,采用HFe对应的Epsilon值对df进行校正；如果GGe和HFe均<0.75，采用GGe对应的Epsilon值对df进行校正。其他情况默认GGe的Epsilon值进行校正。接着分析各因素间是否存在交互作用，若存在交互作用(P<0.05)，需要固定某一因素分析另一因素的单独效应；若不存在交互作用(P>0.05)，则直接分析各因素的主效应。统计分析与结果选择流程见图1。

图1 统计分析与结果选择流程图

1.2 R语言重复测量方差分析程序脚本

#1 将数据导入分析系统

library(readxl)

data <-read_excel(file.choose())

#2 正态性检验(Shapiro-Wilk检验，小样本使用的正态性检验)、异常值检验

library(RVAideMemoire)

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

#3 基本描述,加载psych包，输出组别和时间的均数与标准差，绘制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语言统计分析ez包，按照实验设计执行统计模型脚本

library(ez)

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

#4.1 满足球形检验直接输出结果

model

#4.2 不满足球形检验，输出校正后的结果

model

table(Int_DFn_TIME,Int_DFd_TIME)

table(Int_DFn_GT,Int_DFd_GT)

#5 若GROUP和TIME不存在交互作用，进行GROUP、TIME的主效应分析

#5.1 GROUP的主效应分析，以及不同组别间两两配对比较

library (agricolae)

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

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

plot(Res)

#5.2 进行TIME的主效应分析，以及不同组别间两两配对比较

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

summary(aov_T)

#6.若GROUP和TIME存在交互作用，分别进行GROUP、TIME的单独效应分析

#6.1 GROUP的单独效应分析，各分组下不同时间点的两两配对比较

#6.2 TIME的单独效应分析，各分组下不同时间点的两两配对比较

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 不同时间点的两两配对比较

#7 绘制各因素的数据轮廓图

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 软件安装与界面

安装时，需将包含可执行程序的文件包拷贝至操作系统为Windows 7及以上版本的计算机上。双击文件夹中的CallRLanguage.exe文件，即可打开应用程序。该软件界面左侧为菜单栏，右侧上方为数据导入格式示例区，右侧中部为参数设置区，右侧下方为结果显示区。见图2。

图2 方差分析结果显示界面

1.4 数据格式

本软件适用于重复测量设计资料的单、双因素两/多水平方差分析。对于单因素重复测量设计来说，时间(time)是唯一的一个重复测量的因素。受试对象(subject)个体是一个自变量因素。对双因素设计来说， 时间(time)是一个重复测量的因素，分组(group)是另一个因素。单、双因素重复测量设计下，各因素的水平可以是两水平，也可以是多水平。在收集数据资料时，单因素数据无分组(group)变量，所以该字段可以为空。2种设计类型资料的数据格式见表1。

表1 重复测量单、双因素两/多水平方差分析数据格式

2 实例分析与结果比较

2.1 数据来源

数据采用由颜虹、王彤主编的第5版《医学统计学》教材[9]中的表12～13数据。某研究者欲比较治疗厌食症的3种不同成分药物的效果,在已构建的厌食症模型大鼠中抽取15只成年雄性大鼠,随机分为3组,每组饲料中添加一种药物,药物有效成分含量相同,连续记录药物治疗前(T0)和治疗后1 d、3 d、5 d和7 d的小鼠体质量。数据中“观测时间”(time)是一个重复测量因素，“实验分组”(group)因素为“药物”，“大鼠编号”(subject)为个体标记，因变量“体质量”(value)是连续型的结果变量。将数据录入表1对应的格式化表格中，单击界面中的“数据导入”按钮，将整理好的数据导入到本软件中。

2.2 数据分析

选择dv(因变量)为value，wid(受试者编号)为subject，within(受试内自变量)为time,between(受试间自变量)为group，考虑自变量的主效应和自变量间的交互效应，type(类型)选III，detail选择trpe。点击“数据分析”按钮，输出统计分析结果。单因素重复测量方差分析时，因为没有分组因素和主体间交互作用，所以between字段要为空，同时type字段要选择II，其余同双因素方差分析。点击结果输出按钮，输出对应结果。

2.3 软件运行结果描述

2.3.1 基本资料描述

统计资料的基本描述，可以点击界面中的“基本描述”按钮进行结果和图表的输出，包括各因素的均数、标准差和QQ图、箱式图等。

2.3.2 重复测量方差分析

图2界面显示球形检验、主体内效应及主体间效应分析结果。结果显示组别和时点2个因素对成年雄性大鼠体质量的影响情况。首先，time主体内作用和group*time的交互作用的球形校验结果均为P>0.05，服从球形假设。因此，划线下面直接输出方差分析的结果，即时点效应可以影响体质量的变化，F(4,48)=344.299,P<0.01，意味着大鼠体质量会随着时点出现变化；治疗方式可以影响大鼠体质量，F(2,12)=5.527,P<0.05；值得注意的是，时点和治疗方式之间存在着交互效应，F(8,48)=18.105,P<0.01，意味着大鼠体质量随着时点的变化趋势会因为治疗方式的不一样而不同。其次，基于以上的分析结果，时点和治疗方式之间的交互效应有统计学意义，则需要继续考察时点和治疗方式之间的单独效应。

2.3.3 事后多重比较

点击“多重比较”按钮，若存在交互作用，输出时点和治疗方式之间的单独效应；若不存在交互作用，输出时点和治疗方式之间的主效应。均采用最小显著差异法(least significant difference，LSD法)进行两两比较。本实例研究显示时点和治疗方式之间存在交互作用，因此软件会分别自动输出每个时间点下不同组别间两两比较结果，每个组别下不同时间点两两比较结果，并输出相应的统计图表。部分显示如图3。

图3 多重比较结果显示界面

从多重比较结果可以得到：1)时点走势的简单效应分析显示，A组药物治疗后4个时间点与治疗前比较，T2、T3、T4治疗时间的体质量有所上升(P<0.05)；B组治疗后T1、T2、T3、T4时间点体质量都增高(P<0.05)；C组治疗后也是所有时间点体质量都有增高,但T3、T4体质量增高有统计学意义(P<0.05)。因此,A组药物治疗后体质量的波动小于B组和C组,B组药物治疗后体质量稳定上升。2)治疗方式的简单效应分析显示，治疗前3组大鼠体质量均无统计学差异；治疗后1 d A组和C组之间的大鼠体质量出现了统计学差异，A组体质量低于C组；治疗后3 d 3组大鼠体质量均无统计学差异；治疗后5 d A组、B组分别与C组之间的大鼠体质量出现了统计学差异，A组体质量低于B组和C组；治疗后7 d A组分别与B组、C组之间的大鼠体质量出现了统计学差异，A组体质量低于B组和C组。

2.3.4 绘制数据轮廓图

点击界面中“绘图输出”按钮，结果显示区呈现数据轮廓图。见图4。

图4 数据轮廓图

2.4 软件运行结果验证

采用SPSS软件进行相应的统计分析，结果表明， R语言程序与SPSS的统计分析结果一致，并与第5版《医学统计学》教材201～207页结果一致,本软件的运行结果准确有效。。本软件仅是调用R语言程序包，未做统计方法代码的修改，所以本软件结果即是R语言统计分析结果。

3 讨论

本研究将C#语言和R语言相结合，开发了一款重复测量资料的单、双因素两/多水平方差分析软件，该软件操作简单，能为临床研究人员提供准确的统计分析结果，提高科研统计效率。

3.1 软件准确性

通过结果比较，可以得知本软件运行结果与SPSS软件以及手工计算结果均一致。说明统计分析结果准确可靠，临床医生可以放心使用。

3.2 软件便捷性

本软件无需安装，完全免费，可视化界面简洁，数据处理过程简单，操作简便，分析结果一键生成，直观清晰，统计结果无需二次加工整理，论文可直接引用。SPSS和R语言软件对于无操作经验和编程经验的初学者来说，在短时间内熟练使用这些方法十分困难。本软件的一键生成和友好界面能够满足临床研究人员

的统计分析需求，具有良好的便捷性。

3.3 软件智能性

使用SPSS和R语言完成重复测量方差分析，操作复杂，步骤繁多，结果需人为解读及判断。而本软件具有良好的智能性，统计分析无需人为干预，自动进行重复测量数据的正态性、异常值的检验，并根据检验结果，智能选择统计分析步骤并输出结果。

综上所述，基于C#与R语言的重复测量设计资料方差分析的自动化实现软件，分析结果准确，使用便捷，具有良好的智能性，值得在临床科研上推广使用。目前，该软件还存在一些不足，仅能实现单变量重复测量定量资料的方差分析。后续计划陆续完善重复测量定量及定性资料其他统计方法的自动化统计分析，更好地满足临床研究需求。