One-Way independent sample & Repeated Measure ANOVA R Example
Data: WOW_data
1.獨立樣本單因子變異數分析（One-way ANOVA）
某位研究者想了解在三種不同的種族下(ethnic)，他們的閱讀成績(reading)是否會有顯著差異?
用R讀取CSV檔，並將此資料命名為aov
aov <-read.csv("D:/104/ML_R/WOW_data.csv",
header=TRUE, sep=",")可用head( )看前幾行的data，或是用view( )看整個data
head(aov[c(1:14)])## X ID ethnic gender senson school_3 school_2 school_1 eco_3 eco_2
## 1 3 1041002 Human 0 1 206 204 102 0 0
## 2 15 1041007 Human 0 1 101 101 101 1 1
## 3 16 1041008 Human 0 1 103 103 103 0 0
## 4 30 1041012 Human 1 1 104 104 104 1 1
## 5 64 1041020 Human 0 1 107 107 107 0 0
## 6 66 1041021 Human 0 1 101 101 101 1 0
## eco_1 reading_3 math_3 WA_S_3
## 1 0 492 498 4.2
## 2 1 494 508 5.0
## 3 0 523 513 4.1
## 4 1 476 493 3.1
## 5 0 515 508 2.4
## 6 0 450 498 4.6利用[c( )]來選取欲分析的變項，並將它取代原本的aov資料檔,[]代表資料內
c( )是將選取的變項集合成數列,(對照data裡的欄位number=2,ethnic=3,reading_3=12)
aov<-aov[c(2,3,12)]
head(aov)## ID ethnic reading_3
## 1 1041002 Human 492
## 2 1041007 Human 494
## 3 1041008 Human 523
## 4 1041012 Human 476
## 5 1041020 Human 515
## 6 1041021 Human 450Summary( ) 可大略看某些統計資訊，如平均數、中位數…等。
summary(aov)## ID ethnic reading_3
## Min. :1041001 Human :87 Min. :423.0
## 1st Qu.:1041049 Orc :42 1st Qu.:489.0
## Median :1041097 Undead :64 Median :503.0
## Mean :1041097 Mean :499.3
## 3rd Qu.:1041145 3rd Qu.:512.0
## Max. :1041193 Max. :543.0str( ) 可看變數的資料型態為何，或是用class( )看指定變數的資料型。
其中「aov$ethnic」，表示在aov中的「種族」變數，在R中，一般是使用 $ (錢號)
去連接二者，但注意順序不要搞錯!!
str(aov)## 'data.frame': 193 obs. of 3 variables:
## $ ID : int 1041002 1041007 1041008 1041012 1041020 1041021 1041023 1041025 1041026 1041027 ...
## $ ethnic : Factor w/ 3 levels "Human ","Orc ",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ reading_3: int 492 494 523 476 515 450 510 515 522 529 ...class(aov$ID)## [1] "integer"這時候你可能會發現，id不是名義變項嗎? 怎麼會是整數?
透過as.factor()來將變項轉成因子就可以了
最後再用class( )檢查一下
aov$ID <-as.factor(aov$ID)
class(aov$ID)## [1] "factor"變異數同質性檢定，用bartlett.test(依變數 ~ 自變數,data=該筆資料)
attach()是將資料的變項進行快取
attach(aov)
bartlett.test(reading_3 ~ ethnic, data=aov)##
## Bartlett test of homogeneity of variances
##
## data: reading_3 by ethnic
## Bartlett's K-squared = 3.3976, df = 2, p-value = 0.1829進行ANOVA分析，用aov(依變數 ~自變數) 進行分析，並將分析結果儲存、並命名為 aov1，再利用summary() 看分析的報表，
發現p值< .01，達到顯著水準，表示3個水準間的均值不完全相同，至少有2個水準間的均值在統計上顯著的不相等。
aov1 <- aov(reading_3 ~ ethnic)
summary(aov1)## Df Sum Sq Mean Sq F value Pr(>F)
## ethnic 2 10758 5379 17.25 1.3e-07 ***
## Residuals 190 59233 312
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1繪製平均數比較圖與盒狀圖，幫助我們進行分析。
(請先載入car、sandwich、RcmdrMisc套件)，再利用plotMeans( )繪圖
library(car)
library(sandwich)
library(RcmdrMisc)with(aov, plotMeans(reading_3, ethnic, error.bars="se"))
with(aov, boxplot (reading_3 ~ethnic, error.bars="se"))
title(xlab="ethnic", ylab="score",main = "Box plot")
進行事後比較TukeyHSD( )發現Orc -Human、Undead -Orc兩組達顯著水準，Undead –Human則未達顯著水準。
TukeyHSD(aov1)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = reading_3 ~ ethnic)
##
## $ethnic
## diff lwr upr p adj
## Orc -Human -19.481117 -27.317880 -11.6443534 0.0000001
## Undead -Human -5.919361 -12.787910 0.9491884 0.1064997
## Undead -Orc 13.561756 5.279199 21.8443131 0.0004412結論
由以上結果，變異數同質性檢定結果不顯著(p=.182)，故可直接進行ANOVA分析，
分析結果達到顯著水準(F(2,190)=17.25,p<.001)，顯示三種種族間閱讀成績均值不完全相同，
至少有兩組間的平均數在統計上是有顯著差異的。根據事後比較的結果顯示：對於學生的閱讀成就表現來說，
三種種族的均值差異均達顯著，Human > Undead > Orc。
2. 相依樣本（或重複量測）單因子變異數分析 (Repeated measure ANOVA)
以下為三種種族各年(reading_1、reading_2、reading_3)閱讀成績的分數，研究者想了解三個時間點的成績是否有達到顯著差異。在讀取資料前，要先了解資料的型態是否和手邊的統計分析工具吻合；一般而言，重複量數設計的分析有以下兩種格式，而為了符合R軟體的
需求，我們需要重新編排資料的形式。
用R讀取CSV檔，並將此資料命名為RMaov_raw，選取所要分析的變項(reading_1、reading_2、reading_3)並另存為RMaov
RMaov_raw<-read.csv("D:/104/ML_R/WOW_data.csv",
header=TRUE, sep=",")
RMaov<-RMaov_raw[c(12,17,22)]
head(RMaov)## reading_3 reading_2 reading_1
## 1 492 478 472
## 2 494 486 478
## 3 523 505 500
## 4 476 457 431
## 5 515 502 500
## 6 450 444 429
利用stack( )進行資料編排，並將它取代原本的RMaov，目的是將原始資料轉換成R可以進行RM-ANOVA分析的形式
RMaov<-stack(RMaov)在RMaov裡新增一個變項RM_id，利用rep(要重複的變項,重複次數)將ID變項加到RMaov裡
RMaov$RM_id = rep(RMaov_raw$ID, 3)利用colnames( )為三個變項命名，最後head一下看是否成功
colnames(RMaov) = c("RM_read_score", "RM_time", "RM_id")
head(RMaov)## RM_read_score RM_time RM_id
## 1 492 reading_3 1041002
## 2 494 reading_3 1041007
## 3 523 reading_3 1041008
## 4 476 reading_3 1041012
## 5 515 reading_3 1041020
## 6 450 reading_3 1041021載入ez套件，使用ezANOVA( ) 進行相依樣本的球型檢定(別忘了要先將id變成factor喔!!)
RMaov$RM_id <-as.factor(RMaov$RM_id)
library(ez)
ezANOVA(data=RMaov, dv=.(RM_read_score), wid=.(RM_id), within=.(RM_time), detailed=TRUE)## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 192 138159768.37 210647.29 125929.3452 1.911016e-272
## 2 RM_time 2 384 47012.17 18005.16 501.3193 8.602855e-108
## p<.05 ges
## 1 * 0.9983477
## 2 * 0.1705412
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 RM_time 0.9358506 0.001779184 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 2 RM_time 0.9397176 1.756696e-101 * 0.9486901 2.020522e-102
## p[HF]<.05
## 2 *(由上述報表發現，球型檢定p<.05，因此本資料違反球型假定。)
載入套件繪製平均數圖
library(car)
library(sandwich)
library(RcmdrMisc)
with(RMaov, plotMeans(RM_read_score, RM_time, error.bars="se"))
進行RM-ANOVA統計分析
使用aov( )
Summary() 球型違反時要加上multivariate=FALSE
RMaov.aov <- aov(RM_read_score ~ RM_time + Error(RM_id/RM_time), data = RMaov)
summary(RMaov.aov,multivariate=FALSE)##
## Error: RM_id
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 192 210647 1097
##
## Error: RM_id:RM_time
## Df Sum Sq Mean Sq F value Pr(>F)
## RM_time 2 47012 23506 501.3 <2e-16 ***
## Residuals 384 18005 47
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(由上報表顯示，受測者內因子達顯著效果 ( F(2,384)=501.3, p<.001)，顯示三個時間點的學生閱讀成績均值不全相同。)
事後比較:
載入套件multcomp、nlme、MASS，建立RM-ANOVA模型，之後利用glht計算Tukey的事後比較，可以發現在每個時間點的測驗成績間顯著不同，
而且reading_3 > reading_2 > reading_1。
library(MASS)
library(nlme)
library(multcomp)Lme.mod<-lme(RM_read_score ~ RM_time,random = ~1|RM_id/RM_time,data = RMaov)
summary(glht(Lme.mod, linfct=mcp(RM_time="Tukey")))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = RM_read_score ~ RM_time, data = RMaov, random = ~1 |
## RM_id/RM_time)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## reading_2 - reading_1 == 0 11.6010 0.6971 16.64 <2e-16 ***
## reading_3 - reading_1 == 0 22.0622 0.6971 31.65 <2e-16 ***
## reading_3 - reading_2 == 0 10.4611 0.6971 15.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)結論
由上報表顯示，球型檢定不成立。因此要使用調整過後的ANOVA分析結果，由報表可以得知，受測者的內因子達顯著效果
(F(2,384)=501.3, p<.001) 表示三個時間點的學生閱讀成績均值不同。事後比較檢定結果發現，reading_3最高，
其次是reading_2，最後才是reading_1。
date: “2016年12月24日,第一版”
author: “邱浩恩”