2017. 7. 14. 06:30ㆍ깜신의 통계 이야기
깜신 김종엽입니다.
이번 영상에서는 세 그룹 이상의 평균을 비교하고 난 뒤, 평균이 서로 동일하지 않다고 확인된 경우
서로 어떤 그룹 사이의 평균 차이가 있는지를 확인하는 방법을 알아보려고 합니다.
영어로는 multiple comparison 이라고 하니, 번역하면 다중비교가 될 것 같고요.
Post hoc test라고 쓰고, 사후검정이라고 읽기도 합니다.
이번 강의 내용은 코드가 많아 따라하기 어려우실까봐
코드 전체 내용을 함께 올려드립니다.
기초 통계 범위내에서는 이번 강의가 어쩌면 가장 복잡할 수도 있겠습니다.
돌려 이야기하면, 갈수록 더 어려워지는게 아니라,
다음 강의부터는 다시 쉬워진다는 의미일 수도 있겠죠?! ^^;;
어렵게 시작하신 공부,
포기하지 마시고,
끝까지 달려보시길 바랄게요.
감사합니다.
###############################################################
####### 세 그룹(이상)의 평균을 비교하고 사후검정(다중비교)을 하는 방법 ######
###############################################################
install.packages("moonBook")
require(moonBook)
moonBook::densityplot(LDLC~Dx, data=acs)
############### 정규분포를 확인하는 첫 번째 방법 ##########
shapiro.test(acs$LDLC[acs$Dx=="NSTEMI"])
shapiro.test(acs$LDLC[acs$Dx=="STEMI"])
shapiro.test(acs$LDLC[acs$Dx=="Unstable Angina"])
############### 정규분포를 확인하는 두 번째 방법 ##########
out = aov(LDLC~Dx, data=acs)
shapiro.test(resid(out))
############### 세 그룹의 등분산을 확인하는 방법 ##########
bartlett.test(LDLC~Dx, data=acs)
############### 정규분포하면서 등분산일 때 세 그룹의 평균 비교 #####
boxplot(LDLC~Dx, data= acs)
out = aov(LDLC~Dx, data= acs)
summary(out)
############### 정규분포는 하나 등분산은 아닐 때 세 그룹의 평균 비교 #####
oneway.test(LDLC~Dx,data=acs,var.equal=FALSE)
############### 세 그룹의 평균을 비교하는 비모수 방법 #####
kruskal.test(LDLC~Dx, data=acs) #이렇게 하면 오류가 납니다!!
kruskal.test(LDLC~factor(Dx), data=acs) #이게 올바른 방법입니다.
############## 사후 검정을 위한 방법 (정규분포하고 등분산일 때) ####
TukeyHSD(out)
############## 사후 검정을 위한 방법 (정규분포하고 등분산은 하지 않을 때) ####
#1) for windows
install.packages('userfriendlyscience')
library(userfriendlyscience)
data(InsectSprays)
densityplot(count~spray, data=InsectSprays)
posthocTGH(InsectSprays$spray, y = InsectSprays$count, method = 'games-howell')
#2) for mac
games.howell <- function(grp, obs) {
#Create combinations
combs <- combn(unique(grp), 2)
# Statistics that will be used throughout the calculations:
# n = sample size of each group
# groups = number of groups in data
# Mean = means of each group sample
# std = variance of each group sample
n <- tapply(obs, grp, length)
groups <- length(tapply(obs, grp, length))
Mean <- tapply(obs, grp, mean)
std <- tapply(obs, grp, var)
statistics <- lapply(1:ncol(combs), function(x) {
mean.diff <- Mean[combs[2,x]] - Mean[combs[1,x]]
#t-values
t <- abs(Mean[combs[1,x]] - Mean[combs[2,x]]) / sqrt((std[combs[1,x]] / n[combs[1,x]]) + (std[combs[2,x]] / n[combs[2,x]]))
# Degrees of Freedom
df <- (std[combs[1,x]] / n[combs[1,x]] + std[combs[2,x]] / n[combs[2,x]])^2 / # Numerator Degrees of Freedom
((std[combs[1,x]] / n[combs[1,x]])^2 / (n[combs[1,x]] - 1) + # Part 1 of Denominator Degrees of Freedom
(std[combs[2,x]] / n[combs[2,x]])^2 / (n[combs[2,x]] - 1)) # Part 2 of Denominator Degrees of Freedom
#p-values
p <- ptukey(t * sqrt(2), groups, df, lower.tail = FALSE)
# Sigma standard error
se <- sqrt(0.5 * (std[combs[1,x]] / n[combs[1,x]] + std[combs[2,x]] / n[combs[2,x]]))
# Upper Confidence Limit
upper.conf <- lapply(1:ncol(combs), function(x) {
mean.diff + qtukey(p = 0.95, nmeans = groups, df = df) * se
})[[1]]
# Lower Confidence Limit
lower.conf <- lapply(1:ncol(combs), function(x) {
mean.diff - qtukey(p = 0.95, nmeans = groups, df = df) * se
})[[1]]
# Group Combinations
grp.comb <- paste(combs[1,x], ':', combs[2,x])
# Collect all statistics into list
stats <- list(grp.comb, mean.diff, se, t, df, p, upper.conf, lower.conf)
})
# Unlist statistics collected earlier
stats.unlisted <- lapply(statistics, function(x) {
unlist(x)
})
# Create dataframe from flattened list
results <- data.frame(matrix(unlist(stats.unlisted), nrow = length(stats.unlisted), byrow=TRUE))
# Select columns set as factors that should be numeric and change with as.numeric
results[c(2, 3:ncol(results))] <- round(as.numeric(as.matrix(results[c(2, 3:ncol(results))])), digits = 3)
# Rename data frame columns
colnames(results) <- c('groups', 'Mean Difference', 'Standard Error', 't', 'df', 'p', 'upper limit', 'lower limit')
return(results)
}
games.howell(InsectSprays$spray, InsectSprays$count)
############## 사후 검정을 위한 방법 (정규분포마저 하지 않을 때) ####
install.packages("nparcomp")
require(nparcomp)
result = mctp(LDLC~Dx, data=acs)
summary(result)
강의 내용 중 install.packages('userfriendlyscience')를 설치하는 부분에 에러가 발생하는 분들이 계신 모양입니다.
R 자체가 무료 소프트웨어이고, 개발자들이 각자 만든 걸 스스로 공유하는 플랫폼이다보니, 상용화 프로그램들과는 다르게 각자의 컴퓨터 환경이 다르면 설치가 제대로 되지 않는 패키지들이 더러 있습니다.
제보해주신 '사운드키'님께 감사드리고요.
제가 전달 받은 내용을 그대로 첨부합니다.
posthocTGH(InsectSprays$spray, y = InsectSprays$count, method = 'games-howell')
결과값이 정상적으로 나왔습니다.
감사합니다. ^^
그럼 모두 즐거운 R라이프 즐기시길 바랄게요~ ^^
'깜신의 통계 이야기' 카테고리의 다른 글
| 상대위험도와 오즈비 개념 설명 - 깜신의 통계 왕초보 탈출 15탄 (1) | 2017.07.31 |
|---|---|
| R스튜디오에서 논문용 테이블 작성 따라하기 - 깜신의 통계 왕초보 탈출 14탄 (2) | 2017.07.25 |
| 세 그룹의 등분산 검정과 평균 비교, R스튜디오 따라하기 -깜신의 통계 왕초보 탈출 12탄 (1) | 2017.07.11 |
| 정규 분포를 검정하는 방법 -깜신의 통계 왕초보 탈출 11탄 (0) | 2017.07.03 |
| 세 그룹 이상에서 평균을 비교하는 방법 -깜신의 통계 왕초보 탈출 10탄 (0) | 2017.06.26 |