當(dāng)重復(fù)測量時(shí)間序列數(shù)據(jù)不滿足來自正態(tài)分布總體或者協(xié)方差不符合球形性時(shí)帘腹，不能進(jìn)行重復(fù)測量方差分析蝌蹂。另外橡羞，當(dāng)結(jié)局變量屬于分類變量或等級變量時(shí)胖烛，也不能進(jìn)行重復(fù)測量方差分析。此時(shí)剃允，可以利用廣義估計(jì)方程進(jìn)行分析教翩。
1. 原始數(shù)據(jù)result.txt杆勇。
分析A、B兩種藥的減肥效果饱亿，每組選3個(gè)人蚜退，每個(gè)人服藥后第0、1彪笼、2钻注、3個(gè)月分別記錄下體重值。
內(nèi)容如下：
2.構(gòu)建廣義估計(jì)方程進(jìn)行分析配猫。
a.導(dǎo)入數(shù)據(jù)并進(jìn)行格式變換幅恋。
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.評價(jià)干預(yù)措施主效應(yīng)
geeglm1 <- geeglm(result ~ t0 + group + time, id = no, corstr = "independence", family = "gaussian", data = mydata)
summary(geeglm1)
anova(geeglm1)
其中，參數(shù)corstr設(shè)置作業(yè)相關(guān)矩陣類型泵肄，表示因變量的各次重復(fù)測量值兩兩之間相關(guān)性的大欣弧淑翼；參數(shù)family設(shè)置鏈接函數(shù)類型，需根據(jù)因變量Y的數(shù)據(jù)分布特征來選擇適合的模型品追，本例中因?yàn)橐蜃兞渴嵌窟B續(xù)數(shù)據(jù)玄括，所以選擇高斯分布，其它類型還有二元Logit分布诵盼、泊松分布等惠豺。
分析結(jié)果如下：
結(jié)果顯示，兩組體重值差異有統(tǒng)計(jì)學(xué)意義（P=2.6e-05）风宁，不同月份體重值有統(tǒng)計(jì)學(xué)差異（P=4.4e-16）。
c.評價(jià)干預(yù)措施的單獨(dú)效應(yīng)
geeglm2 <- geeglm(result ~ t0 + time + group : time, id = no, corstr = "independence", family = "gaussian", data = mydata)
summary(geeglm2)
分析結(jié)果如下：
結(jié)果顯示蛹疯，服用減肥藥第1戒财、2、3月兩組體重值均有統(tǒng)計(jì)學(xué)差異（P=0.00106捺弦、0.00028饮寞、3.8e-08）。
d.評價(jià)時(shí)間的單獨(dú)效應(yīng)
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)
分析結(jié)果如下：
結(jié)果顯示列吼，A幽崩、B兩種減肥藥均能使體重在1、2寞钥、3月逐月降低慌申。
