两因素重复测量方差分析(Two-Way Repeated-Measures ANOVA)一(无交互作用)——Python软件实现

关键词：Python; 重复测量; 重复测量资料; 重复测量方差分析; 两因素重复测量方差分析; 球形检验; 交互作用; 主效应; 单独效应

一、案例介绍

研究A、B两种饲料对家兔的增重效果，选择20只家兔，随机分成两组，第一组用饲料A饲养，第二组用饲料B饲养，并于试验开始后第一个月(time1)、第二个月(time2)、第三个月(time3)分别测量2组家兔体重(kg)，试比较A、B两种饲料对家兔的增重效果有无差别？部分数据见图1。本案例数据可从“附件下载”处下载。

二、问题分析

本案例的分析目的是比较A、B两种饲料对家兔的增重效果有无差别。由于3个时间点的数据属于重复测量数据，且有两个组别，可以使用两因素(时间因素time和分组因素group)重复测量方差分析。但需要满足以下6个条件：

条件1：观察变量唯一，且为连续变量。本研究中观察变量只有体重，且为连续变量，该条件满足。

条件2：有两个分析因素。本研究有时间因素time和分组因素group两个因素，该条件满足。

条件3：观察变量为重复测量数据，即不满足独立性。本研究中两个组别在三个时间点时测量的体重均是针对同一批样本，因此不满足独立性，该条件满足。

条件4：观察变量不存在显著的异常值，该条件需要通过软件分析后判断。

条件5：各组、各水平(时间点)观察变量为正态(或近似正态)分布，该条件需要通过软件分析后判断。

条件6：相互比较的各处理水平(组别)的总体方差齐，即方差齐同，该条件需要通过软件分析后判断。

三、软件操作及结果解读

(一) 导入数据

import pandas as pd #导入pandas包
df = pd.read_csv('两因素重复测量方差分析一.csv')
df

在数据栏目(图1)中可以查看全部数据情况，数据集中共有4个变量和20条观察数据，4个变量分别代表饲料种类(group)、试验开始后第一个月(time1)、第二个月(time2)、第三个月(time3) 的家兔体重(kg)。

(二) 适用条件判断

1. 条件4判断(异常值判断)

(1) 软件操作

## 分组描述数据 ##

pd.set_option('display.max_columns',None)#设置展示所有列数
df.groupby("group").describe()

## 绘制箱线图 ##

#将group重新赋值
df.loc[df['group']==1,'group']='A'
df.loc[df['group']==2,'group']='B'
df

## 分组绘制箱线图 ##

import matplotlib.pyplot as plt
plt.boxplot([df.loc[df.loc[:,'group']=='A','time1'],
                     df.loc[df.loc[:,'group']=='B','time1']],
                    labels=["A","B"])
plt.show()

plt.boxplot([df.loc[df.loc[:,'group']=='A','time2'],
                     df.loc[df.loc[:,'group']=='B','time2']],
                    labels=["A","B"])
plt.show()

plt.boxplot([df.loc[df.loc[:,'group']=='A','time3'],
                     df.loc[df.loc[:,'group']=='B','time3']],
                    labels=["A","B"])
plt.show()

(2) 结果解读

描述统计结果(图2)中，列出了各组观察变量的最小值和最大值，依据专业尚不能认为存在异常值；将group重新赋值(图3)后，绘制箱线图(图4-1~图4-3)发现，各时间点的分组箱线图中分别存在一个离群值，但根据专业尚不能判定为异常值。综上，本案例未发现需要处理的异常值，满足条件4。

2. 条件5判断(正态性检验)

(1) 软件操作

## 正态性检验 ##

from scipy import stats #导入stats模块
shapiro_test1 = stats.shapiro(df.loc[df.loc[:,'group']=='A','time1']) 
print(shapiro_test1)  #进行正态性检验
shapiro_test2 = stats.shapiro(df.loc[df.loc[:,'group']=='B','time1']) 
print(shapiro_test2)
shapiro_test3 = stats.shapiro(df.loc[df.loc[:,'group']=='A','time2']) 
print(shapiro_test3)
shapiro_test4 = stats.shapiro(df.loc[df.loc[:,'group']=='B','time2']) 
print(shapiro_test4)
shapiro_test5 = stats.shapiro(df.loc[df.loc[:,'group']=='A','time3']) 
print(shapiro_test5)
shapiro_test6 = stats.shapiro(df.loc[df.loc[:,'group']=='B','time3']) 
print(shapiro_test6)

(2) 结果解读

正态性检验结果(图5)显示A组3个时间点的P值分别为0.3597、0.7700、0.9191，B组3个时间点的P值分别为0.8127、0.8014、0.1820，均P＞0.1，提示各组数据均服从正态分布。此外，本案例也可以绘制Q-Q图，结果也提示各组数据均服从正态分布(读者可自行操作)。综上，本案例满足条件5。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现(链接)。

3. 条件6判断(方差齐性检验)

(1) 软件操作

leveneTestRes1 = stats.levene(df.loc[df.loc[:,'group']=='A','time1'],df.loc[df.loc[:,'group']=='B','time1'],center='mean') #进行方差齐性检验
print(leveneTestRes1)
leveneTestRes2 = stats.levene(df.loc[df.loc[:,'group']=='A','time2'],df.loc[df.loc[:,'group']=='B','time2'],center='mean')
print(leveneTestRes2)
leveneTestRes3 = stats.levene(df.loc[df.loc[:,'group']=='A','time3'],df.loc[df.loc[:,'group']=='B','time3'],center='mean')
print(leveneTestRes3)

(2) 结果解读

    Levene’s方差齐性检验结果(图6)显示，time1、time2、time3时，两组之间的方差齐性检验结果分别为F=1.3929、P=0.2533，F=1.7413、P=0.2035，F=0.0002、P=0.9876；提示每个时间点的两组之间都满足方差齐性。因此本案例满足条件6。关于方差齐性检验的更多内容请阅读医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现。

(三) 球形检验

1. 软件操作

## 宽数据转换为长数据 ##

df2 = df.melt(id_vars=['ID','group'])
df2

## 球形检验 ##

import pingouin
pingouin.sphericity(df2,dv='value',within='variable',subject='ID',method='mauchly')

2. 结果解读

莫奇来球形检验结果（SpherResults）见图8可知，W=0.884，P=0.330，表示数据满足球形假设。因此本案例可以直接采用非校正方法分析的结果。

(四) 交互作用判断

由于本案例有两个因素(一个为时间因素time，另一个为分组因素group)，因此需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。

## 重复测量方差分析 ##

pingouin.mixed_anova(data= df2,dv='value',within='variable',subject='ID',between='group')

分析结果见图9，time与group之间的交互作用结果为F_Interaction=0.106，P=8.995e-01，提示time与group之间的交互作用无统计学意义，本案例可直接以主效应分析结果进行判断，如果交互作用有统计学意义，则需要分析单独效应。

(五) 时间效应分析

time的检验结果(图9)可知，F_{time(variable)}=437.3，P＜0.001，表示体重变化具有时间变化趋势。

(六) 组间效应分析

    group的检验结果(图9)可知，F_group＜0.001，P=9.903e-01，表示两种饲料对家兔的增重效果差异无统计学意义。

(七) 事后检验(两两比较)

1. 软件操作

## 事后检验 ##

from statsmodels.stats.multicomp import pairwise_tukeyhsd
dc_sales_anova_post=pairwise_tukeyhsd(df2['value'],
                                      df2['variable'],
                                      alpha=0.05)
print(dc_sales_anova_post.summary())  #检验不同时间之间的差异

2. 结果解读

事后检验比较-时间 (TukeyHSD)结果(图10)中提供了各时间点两两比较的meandiff (均数差)、p.adj (调整后的P值)、lower (95%置信区间下限)、upper (95%置信区间上限)和reject(是否拒绝原假设)。可知，随着时间的延长，各时刻与time1时刻相比，均数差逐渐增大。time2比time1时均值增加了2.36 kg，至time3时增加了4.11 kg，time2到time3时增加了1.75 kg，差异均有统计学意义(reject=True)。

四、结论

本研究采用两因素重复测量方差分析比较A、B两种饲料对家兔的增重效果有无差别。通过专业知识判断，数据不存在异常值；通过Shapiro-Wilk检验，提示各组数据服从正态分布；通过Levene’s检验，提示每个时间点的两组之间都满足方差齐性；球形度检验提示，满足球形假设(W=0.884，P=0.330)；组别与time无交互作用(F_interaction=0.106，P=8.995e-01)，故进行主效应分析。

A、B两组在试验开始时、第1个月和第2个月时的体重分别为1.54±0.33、3.91±0.91、5.60±0.89 kg和1.51±0.27、3.86±0.58、5.67±0.91 kg。两因素重复测量方差分析结果显示，在试验开始时、第1个月和第2个月2组家兔的体重均呈上升趋势(F_time=437.3，P＜0.001)，但两种饲料对家兔的增重效果无差异(F_group＜0.001，P=9.903e-01)。进一步进行两两比较可知，随着时间的延长，第1个月和第2个月与试验开始时相比，体重增加均有统计学意义(reject=True)。综上可知，两种饲料对家兔的增重效果无差异。