4-3. Rmarkdown ANOVA 와 interactionplot

2019. 2. 21. 11:18
R Notebook

one-way ANOVA(요인 = 범주1개) 분석해보기

# 데이터 준비
anova_data<-data.frame(
  group=c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3),
  score=c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8,47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9,46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2))

# 데이터 확인
head(anova_data)
str(anova_data)
## 'data.frame':    30 obs. of  2 variables:
##  $ group: num  1 1 1 1 1 1 1 1 1 1 ...
##  $ score: num  50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 ...
# 1. 범주형 데이터를 (숫자로 나뉘어도) 문자형 변수로 변경***
# - as.character ( 칼럼인덱싱 ) 
anova_data$group <- as.character(anova_data$group)
str(anova_data)
## 'data.frame':    30 obs. of  2 variables:
##  $ group: chr  "1" "1" "1" "1" ...
##  $ score: num  50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 ...
# 2. 집단별 종속변수(숫자형) 데이터의 boxplot 확인 및 tapply(종속변수인덱싱, 범주칼럼인덱싱, 집계함수)로 집단별 종속변수의 평균과 표준편차.
boxplot(score~group,data=anova_data)
tapply(anova_data$score,anova_data$group, mean)
##     1     2     3 
## 51.46 46.94 46.35
tapply(anova_data$score,anova_data$group, sd)
##         1         2         3 
## 0.6834553 1.0415800 0.7763876
# 3. 아노바분석 전, 3집단 등분산성 
# 3_1. bartlett.test ***
bartlett.test( score ~ as.factor(group), data= anova_data) #  0.4368
## 
##  Bartlett test of homogeneity of variances
## 
## data:  score by as.factor(group)
## Bartlett's K-squared = 1.6565, df = 2, p-value = 0.4368
# 3_2. userfriendlyscience패키지의 oneway()함수로 등분산성 테스트 및 바로 oneway 아노바 적용
# install.packages("userfriendlyscience")
library(userfriendlyscience)

#      등분산인지 확인
oneway(factor(anova_data$group),y=anova_data$score,fullDescribe=T,correction=T)
## ### Oneway Anova for y=score and x=group (groups: 1, 2, 3)
## 
## Omega squared: 95% CI = [.78; .92], point estimate = .88
## Eta Squared: 95% CI = [.8; .92], point estimate = .89
## 
##                                    SS Df    MS      F     p
## Between groups (error + effect) 156.3  2 78.15 108.81 <.001
## Within groups (error only)      19.39 27  0.72             
## 
## 
## ### Welch correction for nonhomogeneous variances:
## 
## F[2, 17.55] = 136.29, p < .001.
## 
## ### Brown-Forsythe correction for nonhomogeneous variances:
## 
## F[2, 23.76] = 108.81, p < .001.
#      one-way anova 시행
oneway.test(score~group,data=anova_data,var.equal = T)
## 
##  One-way analysis of means
## 
## data:  score and group
## F = 108.81, num df = 2, denom df = 27, p-value = 1.199e-13
# 4. oneway 아노바 분석 aov( 종속변수 ~ 요인, data= ) : 그룹간 변동, 그룹내변동(Residual) 만 계산
# - H0 : Ma = Mb = Mc ( 모두 평균이 동일 )
# - Sum of Squares = 변동 = 편차(X-m)의 제곱 -> 제곱을 해야 합이 0이 안되므로.
# - group = 요인1에 대한 그룹간 변동 = ( 전체평균 - 그룹의 평균의  ) 의 제곱
# - Residuals = error = 그룹내 변동 =( 각 그룹별 x - 각 그룹별 평균) 의 제곱
# - aov()결과는 p-value가 안나온다.
aov( score ~ group, data = anova_data) # Sum of Square = 변동, Deg. of Freedom -> 변동의 평균 구해야하므로 필요.
## Call:
##    aov(formula = score ~ group, data = anova_data)
## 
## Terms:
##                   group Residuals
## Sum of Squares  156.302    19.393
## Deg. of Freedom       2        27
## 
## Residual standard error: 0.8475018
## Estimated effects may be unbalanced
# 5. summaray(  aov())를 통한  ANOVA 테이블 도출 : 여기서 p-value가 나온다.
# - Df : 자유도 -> F분포를 결정 + 변동들의 평균구할때 분모로 나눠짐
# - Sum Sq : 변동들
# - Mean sq : 변동들의 평균
# - F value : F검정 통계량 = effect(그룹간 변동의 평균) / error(그룹내 변동의 평균)
summary( aov( score ~ group, data = anova_data)  ) # 1.2e-13
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## group        2 156.30   78.15   108.8 1.2e-13 ***
## Residuals   27  19.39    0.72                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 6. 사후분석1 : laercio 패키지의 LDuncan ( 3집단 비교시 어느것이 차이가 있는지 알아야함)
# - aov()의 결과를 넣는다. 
# install.packages("laercio")
library(laercio)

a1 <- aov( score ~ group, data = anova_data)
LDuncan( a1, "group" ) # 1번 그룹만 차이가 있는 그룹으로 나왔다.
## 
##  DUNCAN TEST TO COMPARE MEANS 
##  
##  Confidence Level:  0.95 
##  Dependent Variable:  score
##  Variation Coefficient:  1.75648 % 
##  
## 
##  Independent Variable:  group 
##   Factors Means   
##   1       51.46 a 
##   2       46.94  b
##   3       46.35  b
# 7. 사후분석2 : 튜키 정직유의검정(TukeyHSD) 및 튜키검정에 의한 평균 차이에 대한 95% 신뢰구간 plot
#    마찬가지로, aov()결과가 들어간다.
#    반드시 칼럼이 문자형이어야한다!
#    자동으로 Holm-Bonferroni correction을 통해 조절된 p-value ( p * nC2 )
TukeyHSD(aov(score~as.character(group),data=anova_data)) # 1-2 ,1-3 만 차이.. 1만 다른 집단!
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = score ~ as.character(group), data = anova_data)
## 
## $`as.character(group)`
##      diff       lwr        upr     p adj
## 2-1 -4.52 -5.459735 -3.5802652 0.0000000
## 3-1 -5.11 -6.049735 -4.1702652 0.0000000
## 3-2 -0.59 -1.529735  0.3497348 0.2813795
plot(TukeyHSD(aov(score~as.character(group),data=anova_data)))

COPD데이터로 one-way 연습 (no quiz)

  • 3집단을 가진 범주 : SMK_STAT_TYPE
  • 종속변수 : HMG
# 데이터 준비
data <- read.csv("week4/4주차_코드/COPD_Lung_Cancer.csv")
str(data)
## 'data.frame':    944 obs. of  20 variables:
##  $ PERSON_ID            : int  10000065 10000183 10000269 10000471 10000788 10001096 10001122 10001260 10001546 10001919 ...
##  $ SEX                  : int  1 1 1 1 2 2 1 1 1 1 ...
##  $ DTH                  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ CTRB_PT_TYPE_CD      : int  8 9 8 1 7 3 10 10 8 8 ...
##  $ BMI                  : num  25.7 22.8 24.5 21.6 25.6 ...
##  $ BP_LWST              : int  90 100 70 100 54 60 100 73 90 70 ...
##  $ BLDS                 : int  112 105 76 70 113 137 89 98 73 233 ...
##  $ TOT_CHOLE            : int  177 177 156 228 199 182 243 284 221 133 ...
##  $ HMG                  : num  12.2 15.7 13.1 14 12.1 12 17 12.1 15 8.7 ...
##  $ FMLY_CANCER_PATIEN_YN: int  1 1 1 NA 1 1 1 1 1 1 ...
##  $ SMK_STAT_TYPE        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ DRNK_HABIT           : int  1 5 1 1 1 1 1 1 1 1 ...
##  $ EXERCI_FREQ          : int  1 1 1 1 1 1 1 2 1 1 ...
##  $ lungC                : int  0 0 1 0 0 0 0 0 0 0 ...
##  $ copd                 : int  0 1 1 1 0 0 0 0 0 1 ...
##  $ Asthma               : int  0 1 0 0 1 0 1 1 0 1 ...
##  $ DM                   : int  1 1 1 1 1 0 1 1 1 1 ...
##  $ Tub                  : int  0 0 1 0 0 0 0 0 0 0 ...
##  $ before_op_score      : int  44 44 37 50 31 34 31 45 35 43 ...
##  $ after_op_score       : int  52 53 60 44 43 49 41 51 59 42 ...
# 데이터 확인
tapply(data$HMG, data$SMK_STAT_TYPE, mean, na.rm=T)
##        1        2        3 
## 13.00319 13.78158 13.34242
tapply(data$HMG, data$SMK_STAT_TYPE, sd, na.rm=T)
##        1        2        3 
## 1.349230 1.600059 1.393871
# 1. 범주형인데 int/num 타입으로 된 것의 인덱스만 *** 가져와서, for문을 통해 
#    반복작업인 as.character로 변형해서 덮어씌우기*************8
index <- c(2, 3, 4, 10:18)
for (i in index) {
  data[,i] <- as.character( data[,i] )
}
str(data)
## 'data.frame':    944 obs. of  20 variables:
##  $ PERSON_ID            : int  10000065 10000183 10000269 10000471 10000788 10001096 10001122 10001260 10001546 10001919 ...
##  $ SEX                  : chr  "1" "1" "1" "1" ...
##  $ DTH                  : chr  "0" "0" "0" "0" ...
##  $ CTRB_PT_TYPE_CD      : chr  "8" "9" "8" "1" ...
##  $ BMI                  : num  25.7 22.8 24.5 21.6 25.6 ...
##  $ BP_LWST              : int  90 100 70 100 54 60 100 73 90 70 ...
##  $ BLDS                 : int  112 105 76 70 113 137 89 98 73 233 ...
##  $ TOT_CHOLE            : int  177 177 156 228 199 182 243 284 221 133 ...
##  $ HMG                  : num  12.2 15.7 13.1 14 12.1 12 17 12.1 15 8.7 ...
##  $ FMLY_CANCER_PATIEN_YN: chr  "1" "1" "1" NA ...
##  $ SMK_STAT_TYPE        : chr  "1" "1" "1" "1" ...
##  $ DRNK_HABIT           : chr  "1" "5" "1" "1" ...
##  $ EXERCI_FREQ          : chr  "1" "1" "1" "1" ...
##  $ lungC                : chr  "0" "0" "1" "0" ...
##  $ copd                 : chr  "0" "1" "1" "1" ...
##  $ Asthma               : chr  "0" "1" "0" "0" ...
##  $ DM                   : chr  "1" "1" "1" "1" ...
##  $ Tub                  : chr  "0" "0" "1" "0" ...
##  $ before_op_score      : int  44 44 37 50 31 34 31 45 35 43 ...
##  $ after_op_score       : int  52 53 60 44 43 49 41 51 59 42 ...
# 2. 3집단 등분산성 검정
bartlett.test( HMG ~ as.factor(SMK_STAT_TYPE), data = data) # 0.3899
## 
##  Bartlett test of homogeneity of variances
## 
## data:  HMG by as.factor(SMK_STAT_TYPE)
## Bartlett's K-squared = 4.3473, df = 2, p-value = 0.1138
# 3. 3집단 정규성 검정
unique(data$SMK_STAT_TYPE) # "1" "3" "2" NA 
## [1] "1" "3" "2" NA
data1 <- data[ data$SMK_STAT_TYPE == "1", ]
data2 <- data[ data$SMK_STAT_TYPE == "2", ]
data3 <- data[ data$SMK_STAT_TYPE == "3", ]

shapiro.test(data1$HMG)
## 
##  Shapiro-Wilk normality test
## 
## data:  data1$HMG
## W = 0.99681, p-value = 0.1639
shapiro.test(data2$HMG) # 0.002898 정규성 만족안하지만, 한다고 가정
## 
##  Shapiro-Wilk normality test
## 
## data:  data2$HMG
## W = 0.94628, p-value = 0.002898
shapiro.test(data3$HMG)
## 
##  Shapiro-Wilk normality test
## 
## data:  data3$HMG
## W = 0.99057, p-value = 0.7168
# 3. 1way anova
one_way <- aov(HMG ~ SMK_STAT_TYPE, data=data)
summary(one_way) # 3.97e-06 -> 평균차이가 명확히 차이남
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## SMK_STAT_TYPE   2   47.8  23.915   12.61 3.97e-06 ***
## Residuals     892 1691.3   1.896                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 49 observations deleted due to missingness
# 4. 사후검정
library(laercio)
LDuncan(one_way, "group") # a, b, c  3집단 모두 평균이 다른 것으로 나오
## 
##  DUNCAN TEST TO COMPARE MEANS 
##  
##  Confidence Level:  0.95 
##  Dependent Variable:  HMG
##  Variation Coefficient:  10.50583 % 
##  
## 
##  Independent Variable:  SMK_STAT_TYPE 
##   Factors Means               
##   2       13.7815789473684 a  
##   3       13.3424242424242  b 
##   1       13.0031944444444   c

2way ANOVA 와 interaction plot

  • 종속변수 : BMI
  • 요인 : DM 과 SMK_STAT_TYPE
# 각 요인별 boxplot
par(mfrow=c(1,2))
boxplot(BMI ~ DM,data=data, 
        boxwex = 0.7,
        main="BMI~DM", 
        col = c("yellow","orange", "green"))
boxplot(BMI ~ SMK_STAT_TYPE,data=data, 
        boxwex = 0.5,
        main="BMI~Smkoing", 
        col = c("blue", "coral"))

# 요인1(DM이력-3가지집단), 요인2(흡연이력-3가지집단)에 따른 BMI 평균의 차이
two_way <- aov( BMI ~ DM + SMK_STAT_TYPE ,data = data)
summary(two_way) # DM : 1.66e-05 *** , SMK_STAT_TYPE : 0.136
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## DM              1    199  198.96  18.748 1.66e-05 ***
## SMK_STAT_TYPE   2     42   21.19   1.997    0.136    
## Residuals     896   9509   10.61                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 44 observations deleted due to missingness
# 요인2개를 순서바꿔서
two_way2 <- aov( BMI ~ SMK_STAT_TYPE + DM  ,data = data)
summary(two_way2) # DM :  3.1e-05 ***, SMK_STAT_TYPE : 0.0746
##                Df Sum Sq Mean Sq F value  Pr(>F)    
## SMK_STAT_TYPE   2     55   27.63   2.604  0.0746 .  
## DM              1    186  186.08  17.534 3.1e-05 ***
## Residuals     896   9509   10.61                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 44 observations deleted due to missingness
# interaction plot
# - interaction.plot( x축-요인1  -  범례 - 요인2  -  종속변수)
# 요인1-종속변수에  요인2가 영향을 끼치는지 보기
par(mfrow=c(1,2))
interaction.plot(data$SMK_STAT_TYPE, data$DM, data$BMI)
interaction.plot(data$DM, data$SMK_STAT_TYPE, data$BMI)

저작자표시

'한의대 생활 > └ 통계에 대한 나의 정리' 카테고리의 다른 글

4-5. 상관관계(계수, 산점도), t-test, anova 간단 복습  (0) 2019.02.21
4-4. Rmarkdown 칼럼명변경 +range() + 사용자정의함수로 scaling  (0) 2019.02.21
4-2. Rmarkdown T-test  (0) 2019.02.20
4-2. 3집단 이상의 (숫자형)분석 ANOVA  (0) 2019.02.19
4-1. 2개 집단의 평균 비교 - t-test  (4) 2019.02.19

+ Recent posts

티스토리툴바