앤디 필드의 유쾌한 R 통계학 챕터 9 - 두 평균의 비교

이번 챕터에서는 t test에 대해 배운다.

유형은 여럿이지만 사용하는 함수는 t.test 하나고 사용도 어렵지 않다.

다만 paired 여부와 사용할 데이터를 여러모로 살펴보고 전처리하는데 주의해야 한다.

통계적 유의성과 별도로 효과크기도 측정하는 것이 좋다.

 

# 필요한 패키지 불러오기
library(ggplot2); library(pastecs); library(WRS);library(reshape)

# 데이터 불러오기
spiderLong <- read.delim('datafiles/spiderLong.dat',header = TRUE)
spiderWide <- read.delim('datafiles/spiderWide.dat',header = TRUE)

# 막대그래프와 오차 막대로 살펴보기
ggplot(spiderLong, aes(x=Group, y=Anxiety))+stat_summary(fun=mean, geom='bar')+stat_summary(fun.data = mean_cl_normal, geom= 'errorbar', position=position_dodge(), width=0.0)

# 각 개체의 평균과 총평균과의 차이 구하기
spiderWide$pMean <- (spiderWide$picture+spiderWide$real)/2
grandMean <- mean(c(spiderWide$picture, spiderWide$real))
spiderWide$adj <- grandMean - spiderWide$pMean

# 총평균과의 개체간 차이를 제거(모든 개체의 평균이 총평균과 같아진다)
spiderWide$picture_adj <- spiderWide$picture + spiderWide$adj
spiderWide$real_adj <- spiderWide$real + spiderWide$adj
spiderWide$pMean2 <- mean(c(spiderWide$picture_adj, spiderWide$real_adj))

# 위의 과정을 위한 함수 작성하기
rmMeanAdjust <- function(dataframe){
  varNames<-names(dataframe)
  pMean <- ((dataframe[,1]+dataframe[,2]))/2
  grandmean <- mean(c(dataframe[,1]+dataframe[,2]))
  adj <- grandmean-pMean
  varA_adj <- dataframe[,1]+adj
  varB_adj <- dataframe[,2]+adj
  output <- data.frame(varA_adj, varB_adj)
  names(output)<-c(paste(varNames[1], 'Adj', sep = '_'), paste(varNames[2], '_Adj', sep='_'))
  return(output)
}

# 조정된 값 시각화하기(오차 막대의 길이가 줄어든 것을 확인할 수 있다)
spider_adj <- spiderWide[,c('real_adj','picture_adj')]
spider_adj2 <- melt(spider_adj)
colnames(spider_adj2) = c('Group','Anxiety')
ggplot(spider_adj2, aes(x=Group, y=Anxiety))+stat_summary(fun=mean, geom='bar')+stat_summary(fun.data = mean_cl_normal, geom= 'errorbar', position=position_dodge(), width=0.0)

# 진행을 위해 데이터 덮어씌우기
Group <- gl(2, 12, labels=c('Picture', 'Real Spider'))
Anxiety <- c(30,35,45,40,50,35,55,25,30,45,40,50,40,35,50,55,65,55,50,35,30,50,60,39)
spiderLong<-data.frame(Group, Anxiety)

# 데이터 살펴보기
ggplot(spiderLong, aes(x=Group, y=Anxiety))+geom_boxplot()
by(spiderLong$Anxiety, spiderLong$Group, stat.desc, basic=FALSE, norm=TRUE)

# t검정 실행하기 - 독립표본 t검정
ind.t.test <- t.test(Anxiety~Group, data=spiderLong)
# 또는 ind.t.test <- t.test(spiderWide$real, spiderWide$picture)
ind.t.test

# 효과크기 확인하기
t <- ind.t.test$statistic[[1]]
df <- ind.t.test$parameter[[1]]
r <- sqrt(t^2/(t^2+df))
round(r,3) # 유의하지는 않지만 효과의 크기는 중간 정도이다

# t검정 실행하기 - 대응표본 t검정

dep.t.test <- t.test(spiderWide$real, spiderWide$picture, paired = TRUE)
dep.t.test

# 효과크기 확인하기
t <- dep.t.test$statistic[[1]]
df <- dep.t.test$parameter[[1]]
r <- sqrt(t^2/(t^2+df))
round(r,3) # 같은 데이터를 사용했지만 대응표본 t검정에서 효과의 크기가 더 크다

# 복습
penis <- read.delim('datafiles/penis.dat', header=TRUE)
ggplot(penis, aes(x=book, y=happy))+stat_summary(fun=mean, geom='bar')+stat_summary(fun.data = mean_cl_normal, geom= 'errorbar', position=position_dodge(), width=0.0)
penis_t <- t.test(happy~book, data=penis)
penis_t

p_t <- penis_t$statistic[[1]]
p_df <- penis_t$parameter[[1]]
p_r <- sqrt(p_t^2/(p_t^2+p_df))
round(p_r,3)

fnh <- read.delim('datafiles/field&hole.dat', header=TRUE)
fnh_melt <- melt(fnh)
ggplot(fnh_melt, aes(x=variable, y=value))+stat_summary(fun=mean, geom='bar')+stat_summary(fun.data = mean_cl_normal, geom= 'errorbar', position=position_dodge(), width=0.0)
fnh_t <- t.test(fnh$women, fnh$statbook, paired = TRUE)
fnh_t

f_t <- fnh_t$statistic[[1]]
f_df <- fnh_t$parameter[[1]]
f_r <- sqrt(f_t^2/(f_t^2+f_df))
round(f_r,3)