깜신 김종엽입니다.


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

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


영어로는 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)



신고