当重复测量时间序列数据不满足来自正态分布总体或者协方差不符合球形性时,不能进行重复测量方差分析。另外,当结局变量属于分类变量或等级变量时,也不能进行重复测量方差分析。此时,可以利用广义估计方程进行分析。
1. 原始数据result.txt。
分析A、B两种药的减肥效果,每组选3个人,每个人服药后第0、1、2、3个月分别记录下体重值。
内容如下:
2.构建广义估计方程进行分析。
a.导入数据并进行格式变换。
library(geepack)
mydata0 <- read.table("result.txt", sep="\t", header=TRUE, colClasses=c("factor", "factor", "factor", "numeric"))
mydata <- mydata0[mydata0$time!=0,]
mydata$t0 <- rep(mydata0[mydata0$time==0,4],each=3)
mydata$time <- as.factor(as.vector(mydata$time)) # 重置factor level
rownames(mydata) <- NULL # 重新排序行索引
mydata格式如下:
b.评价干预措施主效应
geeglm1 <- geeglm(result ~ t0 + group + time, id = no, corstr = "independence", family = "gaussian", data = mydata)
summary(geeglm1)
anova(geeglm1)
其中,参数corstr设置作业相关矩阵类型,表示因变量的各次重复测量值两两之间相关性的大小;参数family设置链接函数类型,需根据因变量Y的数据分布特征来选择适合的模型,本例中因为因变量是定量连续数据,所以选择高斯分布,其它类型还有二元Logit分布、泊松分布等。
分析结果如下:
结果显示,两组体重值差异有统计学意义(P=2.6e-05),不同月份体重值有统计学差异(P=4.4e-16)。
c.评价干预措施的单独效应
geeglm2 <- geeglm(result ~ t0 + time + group : time, id = no, corstr = "independence", family = "gaussian", data = mydata)
summary(geeglm2)
分析结果如下:
结果显示,服用减肥药第1、2、3月两组体重值均有统计学差异(P=0.00106、0.00028、3.8e-08)。
d.评价时间的单独效应
geeglm3_A <- geeglm(result ~ time, id = no, corstr = "independence", family = "gaussian", data = mydata0, subset = (group == "A"))
summary(geeglm3_A)
geeglm3_B <- geeglm(result ~ time, id = no, corstr = "independence", family = "gaussian", data = mydata0, subset = (group == "B"))
summary(geeglm3_B)
分析结果如下:
结果显示,A、B两种减肥药均能使体重在1、2、3月逐月降低。