본문 바로가기

깜신의 통계 이야기

여러 그룹의 평균 비교를 위한 다중비교-R스튜디오 따라하기, 깜신의 통계 왕초보 탈출 13탄




깜신 김종엽입니다.


이번 영상에서는 세 그룹 이상의 평균을 비교하고 난 뒤, 평균이 서로 동일하지 않다고 확인된 경우

서로 어떤 그룹 사이의 평균 차이가 있는지를 확인하는 방법을 알아보려고 합니다.


영어로는 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 자체가 무료 소프트웨어이고, 개발자들이 각자 만든 걸 스스로 공유하는 플랫폼이다보니, 상용화 프로그램들과는 다르게 각자의 컴퓨터 환경이 다르면 설치가 제대로 되지 않는 패키지들이 더러 있습니다. 


제보해주신 '사운드키'님께 감사드리고요.

제가 전달 받은 내용을 그대로 첨부합니다.



13탄 강의에서 Kruskal Wallis H Test 에서 세 그룹의 평균이 다른 경우 사후검정 방법 의 Windows 사용자 부분에서
install.packages('userfriendlyscience')   패키지 설치시 하나가 에러가 났습니다.

Warning in install.packages :
  cannot open URL 'https://cran.rstudio.com/bin/windows/contrib/3.4/digest_0.6.13.zip': HTTP status was '404 Not Found'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
Warning in install.packages :
  download of package ‘digest’ failed


digest 라는 패키지가 설치되지 않아서 
posthocTGH(InsectSprays$spray, y = InsectSprays$count, method = 'games-howell')  명령어가 에러가 났습니다.
원인은 다운로드 주소가 바뀌었더군요.. 폴더명이 미묘(?)하게 바뀌었습니다. -.-;;



그래서 digest 라는 패키지를 수동설치한 정보를 제보드립니다.

### digest 패키지 수동설치
# 수동설치방법 : http://rfriend.tistory.com/177

# version check
version
install.packages("D:/r_temp/digest_0.6.14.zip", repos = NULL, type = "win.binary")

설치후 
library(userfriendlyscience)
posthocTGH(InsectSprays$spray, y = InsectSprays$count, method = 'games-howell')

결과값이 정상적으로 나왔습니다.


감사합니다. ^^


그럼 모두 즐거운 R라이프 즐기시길 바랄게요~ ^^