깜신 김종엽입니다.
이번 영상에서는 세 그룹 이상의 평균을 비교하고 난 뒤, 평균이 서로 동일하지 않다고 확인된 경우
서로 어떤 그룹 사이의 평균 차이가 있는지를 확인하는 방법을 알아보려고 합니다.
영어로는 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탄 (1) | 2017.07.25 |
| 여러 그룹의 평균 비교를 위한 다중비교-R스튜디오 따라하기, 깜신의 통계 왕초보 탈출 13탄 (4) | 2017.07.14 |
| 세 그룹의 등분산 검정과 평균 비교, R스튜디오 따라하기 -깜신의 통계 왕초보 탈출 12탄 (1) | 2017.07.11 |
| 정규 분포를 검정하는 방법 -깜신의 통계 왕초보 탈출 11탄 (0) | 2017.07.03 |
| 세 그룹 이상에서 평균을 비교하는 방법 -깜신의 통계 왕초보 탈출 10탄 (0) | 2017.06.26 |
안녕하세요.
미생물 생태학을 공부하는 학생인데 이 블로그에서 정말 도움을 많이 받았습니다.
기본적인 통계방법을 어떻게 적용하는지 잘 알게되었습니다.
궁금한 걸 댓글로 문의드리는게 맞는지 모르겠지만,
tukey나 mctp function을 사용해서 그룹간의 차이를 알게된 후, p값을 기반으로
boxplot 위에 문자로 차이를 나타내는 경우를 논문에서 많이 봤습니다.
(예를들어 차이없는건 같은 문자, 차이나는건 다른 문자로)
저도 그런 방식으로 표현하고 싶은데 구글링을 해도 쉽지가 않네요ㅜㅜ
최종적으로 plot을 통해 표현하는 방법까지 강의해주실 여력이 있으실까요?
감사합니다!
영상강의에서 welch's anova의 사후검정이랑 kruskal wallis H test의 사후검정이랑 뒤바뀐거같습니다
한번확인해주시면 감사하겠습니다
> install.packages('userfriendlyscience')
Installing package into ‘C:/Users/hp/Documents/R/win-library/4.0’
(as ‘lib’ is unspecified)
Warning in install.packages :
package ‘userfriendlyscience’ is not available (for R version 4.0.0)
이렇게 되는 데 해결 방법이 있을 까요..?
https://cran.r-project.org/src/contrib/Archive/userfriendlyscience/
저는 digest 문제는 아니었고, package ‘userfriendlyscience’ is not available (for R version 4.0.0)
라고 나오는 오류를 위의 링크의 tar.gz 파일 다운로드 받아서 수동 설치 하는 방법으로 해결하였습니다.
오랜시간 끙끙 댔기에 필요하신 분을 위해 남겨놓아요 ^^