깜신의 작은 진료소

###############################################################

####### 세 그룹(이상)의 평균을 비교하고 사후검정(다중비교)을 하는 방법 ######

###############################################################

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

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

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