醫學統計與R語言:雙因素重複測量方差分析(Two-way repeated measures ANOVA)

微信公眾號：醫學統計與R語言
如果你覺得對你有幫助，歡迎轉發

輸入1：

setwd("C:\\Users\\mooshaa\\Desktop")
rma <- read.csv("rma.csv",header=T)
rma$group <- factor(rma$group,levels = c(1,2),labels = c("treatment","control"))
rma[sample(nrow(rma),5,replace=F),]

結果1：

   id     group before middle after
56 56   control     58     63    68
19 19 treatment     58     60    70
47 47   control     56     50    69
28 28 treatment     51     52    70
59 59   control     57     64    71

輸入2：

library(tidyr)
longrma <-  gather(rma,key=time,value=score,-id,-group)
longrma$time <- factor(longrma$time,levels = c("before","middle","after"))
longrma$group <- factor(longrma$group,levels = c("control","treatment"))
longrma[sample(nrow(longrma),5,replace=F),]

結果2：

    id     group   time score
179 59   control  after    71
163 43   control  after    51
14  14 treatment before    44
101 41   control middle    41
153 33   control  after    60

輸入3：

funx <- function(x){
  n <- length(x)
  mean <- mean(x)
  sd <- sd(x)
  se <- sd/sqrt(n)
  l1 <- mean-1.96*se
  u1 <- mean+1.96*se
  result <-c(n=n,mean=mean,sd=sd,se=se,l1=l1,u1=u1)
  return(result)
}
aggregate(longrma$score,list(longrma$time,longrma$group),FUN=funx )

結果3：

  Group.1   Group.2   x.n x.mean  x.sd  x.se  x.l1  x.u1
1  before   control 30.00  54.50  8.36  1.53 51.51 57.49
2  middle   control 30.00  60.27 10.54  1.92 56.50 64.04
3   after   control 30.00  64.37  8.60  1.57 61.29 67.44
4  before treatment 30.00  54.83  8.88  1.62 51.66 58.01
5  middle treatment 30.00  61.23  8.06  1.47 58.35 64.12
6   after treatment 30.00  74.90  7.91  1.44 72.07 77.73

輸入4：

library(nlme)
result.lme <- lme(score~time*group,random=~1|id,data=longrma)
anova(result.lme )

結果4：

             numDF denDF F-value p-value
(Intercept)     1   116    4516  <.0001
time            2   116      86  <.0001
group           1    58       5  0.0359
time:group      2   116      12  <.0001

shapiro.test(result.lme residuals);qqnorm(result.lmeresiduals)
qqline(result.lme $residuals)
plot(fitted(result.lme),residuals(result.lme),xlab = "Predicted Values",ylab = "Residuals")

輸入5：

library(ggplot2)
library(afex)    
afexout1 <- aov_ez("id", "score",  between=c("group"),within = c("time"),type=3,
                   data = longrma)    ##  return="Anova" 輸出多變量分類結果
summary(afexout1)

其它的表達方式：
 aov_4(score ~ group + (time|id),data = longrma)
 aov_car(score ~ group + Error(id|time),data = longrma)

結果5：

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

            Sum Sq num Df Error SS den Df F value  Pr(>F)    
(Intercept) 684870      1     8797     58 4515.55 < 2e-16 ***
group          700      1     8797     58    4.62   0.036 *  
time          6798      2     4576    116   86.17 < 2e-16 ***
group:time     980      2     4576    116   12.42 1.3e-05 ***
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1


Mauchly Tests for Sphericity

           Test statistic p-value
time                0.949   0.225
group:time          0.949   0.225


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

           GG eps Pr(>F[GG])    
time        0.952    < 2e-16 ***
group:time  0.952    1.9e-05 ***
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1

           HF eps Pr(>F[HF])
time        0.983   2.66e-23
group:time  0.983   1.49e-05

輸入6：

p1 <- afex_plot( afexout1 , x = "time", trace = "group",error = "within",
                   mapping = c("linetype","shape", "fill"),
                   data_geom = ggplot2::geom_boxplot, 
                   data_arg = list(width = 0.4))+theme_bw() + theme(legend.position="bottom")
p2 <- afex_plot( afexout1 , x = "group", trace = "time",error = "within",
                   mapping = c("linetype","shape", "color"),
                   point_arg = list(size = 2), line_arg = list(size = 1),
                   error_arg = list(size = 1, width = 0.1, linetype = 1),
                   factor_levels = list(time = c("治療前", "治療中", "治療後"),
                   group = c("對照組", "幹預組")), legend_title = "時間" )+ 
                   labs(y = "得分", x = "分組")+theme_classic() + theme(legend.position="bottom")

library("cowplot")
plot_grid(p1,p2)  
ggsave("p1.png", device = "png", width = 25, height = 12, units = "cm", 
         dpi = 600) 

結果6：

輸入7：

nice( afexout1, correction = "GG")

結果7：

Anova Table (Type 3 tests)

Response: score
      Effect           df    MSE         F ges p.value
1      group        1, 58 151.67    4.62 * .05     .04
2       time 1.90, 110.38  41.46 86.17 *** .34  <.0001
3 group:time 1.90, 110.38  41.46 12.42 *** .07  <.0001
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『+』 0.1 『 』 1

Sphericity correction method: GG 

輸入8：

library(lsmeans) ###compute least-squares means(predicted marginal means)
ref.grid(afexout1)
l1 <- lsmeans(afexout1,"group")
l2 <- lsmeans(afexout1,"time")
l3 <- lsmeans(afexout1,~time|group)
contrast(l1,method="pairwise")

結果8：

  contrast          estimate   SE  df t.ratio p.value
 control - treatment    -3.94 1.84 58 -2.149  0.0359 

Results are averaged over the levels of: time 

輸9：

 contrast(l2,method="pairwise",adjust = "bonferroni")

結果9：

 contrast        estimate   SE  df t.ratio p.value
 before - middle    -6.08 1.15 116  -5.310 <.0001 
 before - after    -14.97 1.15 116 -13.050 <.0001 
 middle - after     -8.88 1.15 116  -7.750 <.0001 

Results are averaged over the levels of: group 
P value adjustment: bonferroni method for 3 tests 

輸入10：

 contrast(l3,method="pairwise",adjust = "bonferroni")

結果10：

group = control:
 contrast        estimate   SE  df t.ratio p.value
 before - middle    -5.77 1.62 116  -3.560 0.0016 
 before - after     -9.87 1.62 116  -6.080 <.0001 
 middle - after     -4.10 1.62 116  -2.530 0.0384 

group = treatment:
 contrast        estimate   SE  df t.ratio p.value
 before - middle    -6.40 1.62 116  -3.950 0.0004 
 before - after    -20.07 1.62 116 -12.370 <.0001 
 middle - after    -13.67 1.62 116  -8.430 <.0001 

P value adjustment: bonferroni method for 3 tests 

輸入11：

   contrast(l2,method="poly")

結果11：

 contrast  estimate   SE  df t.ratio p.value
 linear        15.0 1.15 116 13.050  <.0001 
 quadratic      2.8 1.99 116  1.410  0.1613 

Results are averaged over the levels of: group