앤디 필드의 유쾌한 R 통계학 챕터 5 - 자료 가정

챕터 5에서는 R을 활용하여 표본 자료의 정규성, 등분산성을 검증한다.

 

정규성 검정은 stat.desc() 함수와 shapiroTest() 함수를 활용하고

등분산 검정은 leveneTest() 함수를 활용한다.

 

이 때 단순히 함수만 사용하는 것이 아니라 히스토그램, Q-Q plot, 분산비 등을 통해 다시 한 번 자료를 살피는 것도 필요하다.

 

rm(list=ls())

# 이번 챕터에서는 R을 통해 분포의 정규성과 분산의 동질성(homogeneity)을 확인하는 방법을 배운다.

library(car);library(ggplot2);library(pastecs);library(psych);library(Rcmdr)


# 눈으로 정규성 확인하기

  # 데이터 불러오기

dlf <- read.delim('datafiles/DownloadFestival(No Outlier).dat', header = TRUE)

  # 밀도 히스토그램 만들기

hist.Day1 <- ggplot(dlf, aes(day1)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine score Day 1', y = 'Density');hist.Day1
hist.Day2 <- ggplot(dlf, aes(day2)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine score Day 1', y = 'Density');hist.Day2
hist.Day3 <- ggplot(dlf, aes(day3)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine score Day 1', y = 'Density');hist.Day3

  # 정규분포 곡선 표시하기

hist.Day1 + stat_function(fun = dnorm, args = list(mean = mean(dlf$day1, na.rm = TRUE), sd = sd(dlf$day1, na.rm = TRUE)), colour = 'black', size = 1)
hist.Day2 + stat_function(fun = dnorm, args = list(mean = mean(dlf$day2, na.rm = TRUE), sd = sd(dlf$day2, na.rm = TRUE)), colour = 'black', size = 1)
hist.Day3 + stat_function(fun = dnorm, args = list(mean = mean(dlf$day3, na.rm = TRUE), sd = sd(dlf$day3, na.rm = TRUE)), colour = 'black', size = 1)

  # Q-Q plot 그리기

qqplot.day1 <- ggplot(dlf, aes(sample=dlf$day1)) + stat_qq(); qqplot.day1
qqplot.day2 <- ggplot(dlf, aes(sample=dlf$day2)) + stat_qq(); qqplot.day2
qqplot.day3 <- ggplot(dlf, aes(sample=dlf$day3)) + stat_qq(); qqplot.day3

  # 숫자로 정규성 판단하기(cbind, 인덱싱을 통해 여러 열의 값을 동시에 얻을 수 있다)
    # skew.2SE, kurt.2SE는 1, 1.29, 1.65보다 클 때 각각 p<.05, p<.01, p<.001에서 유의하다

describe(cbind(dlf$day1, dlf$day2, dlf$day3))
stat.desc(dlf[, c('day1', 'day2', 'day3')], basic = FALSE, norm = TRUE)

  # 전반적인 복습

rexam <- read.delim('datafiles/rexam.dat', header = TRUE)
rexam$uni <- factor(rexam$uni, levels= c(0:1), labels= c('Duncetown Univ', 'Sussex Univ'))

hist.exam <- ggplot(rexam, aes(exam)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Exam score', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(rexam$exam, na.rm = TRUE), sd = sd(rexam$exam, na.rm = TRUE)), colour = 'black', size = 1);hist.exam

hist.computer <- ggplot(rexam, aes(computer)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Computer score', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(rexam$computer, na.rm = TRUE), sd = sd(rexam$computer, na.rm = TRUE)), colour = 'black', size = 1);hist.computer

hist.lecture <- ggplot(rexam, aes(lectures)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Lecture score', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(rexam$lectures, na.rm = TRUE), sd = sd(rexam$lectures, na.rm = TRUE)), colour = 'black', size = 1);hist.lecture

hist.numeracy <- ggplot(rexam, aes(numeracy)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Numeracy score', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(rexam$numeracy, na.rm = TRUE), sd = sd(rexam$numeracy, na.rm = TRUE)), colour = 'black', size = 1);hist.numeracy

round(stat.desc(rexam[,c('exam', 'computer', 'lectures', 'numeracy')], norm = TRUE, basic = FALSE), digits = 3)

  # by 함수를 다양하게 활용하여 factor별 통계량 얻기

by(data = rexam$exam, INDICES = rexam$uni, FUN = describe)
by(rexam$exam, rexam$uni, stat.desc, basic = FALSE, norm = TRUE)
by(rexam$exam, rexam$uni, stat.desc, basic = FALSE, norm = TRUE)
by(rexam[, c('exam', 'numeracy')], rexam$uni, stat.desc, basic = FALSE, norm = TRUE)

  # factor별 그래프를 그리기 위해서는 데이터프레임을 나눈 뒤 각각 그린다

dunceData <- subset(rexam, rexam$uni == 'Duncetown Univ')
sussexData <- subset(rexam, rexam$uni == 'Sussex Univ')

hist.nueracy.dunce <- ggplot(dunceData, aes(numeracy)) + theme(legend.position = 'None') + geom_histogram(aes(y = ..density..), fill = 'white', colour = 'black', binwidth = 1) + labs(x = 'Numeracy', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(dunceData$numeracy, na.rm = TRUE), sd = sd(dunceData$numeracy, na.rm = TRUE)), colour = 'blue', size = 1);hist.nueracy.dunce

hist.nueracy.sussex <- ggplot(sussexData, aes(numeracy)) + theme(legend.position = 'None') + geom_histogram(aes(y = ..density..), fill = 'white', colour = 'black', binwidth = 1) + labs(x = 'Numeracy', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(sussexData$numeracy, na.rm = TRUE), sd = sd(sussexData$numeracy, na.rm = TRUE)), colour = 'blue', size = 1);hist.nueracy.sussex

hist.exam.dunce <- ggplot(dunceData, aes(exam)) + theme(legend.position = 'None') + geom_histogram(aes(y = ..density..), fill = 'white', colour = 'black', binwidth = 1) + labs(x = 'Exam', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(dunceData$exam, na.rm = TRUE), sd = sd(dunceData$exam, na.rm = TRUE)), colour = 'blue', size = 1);hist.exam.dunce

hist.exam.sussex <- ggplot(sussexData, aes(exam)) + theme(legend.position = 'None') + geom_histogram(aes(y = ..density..), fill = 'white', colour = 'black', binwidth = 1) + labs(x = 'Exam', y = 'Density') + stat_function(fun = dnorm, args = list(mean = mean(sussexData$exam, na.rm = TRUE), sd = sd(sussexData$exam, na.rm = TRUE)), colour = 'blue', size = 1);hist.exam.sussex


# 분포의 정규성 검정하기
  # 샤피로-윌크 검정은 표본의 점수와 해당 표본과 평균과 표준편차가 같은 정규분포에서 뽑은 점수를 비교한다. 검정 결과가 유의하지 않다면 표본의 분포는 정규분포를 따른다

  # shapiro.test 함수를 사용하지 않아도 stat.desc함수를 통해 normtest.W, normtest.p값을 확인할 수 있다

shapiro.test(rexam$exam)
by(rexam$exam, rexam$uni, shapiro.test)
by(rexam$numeracy, rexam$uni, shapiro.test)

  # Q-Q plot으로 시각적으로 다시 확인하기(표본이 크면 샤피로-윌크 검정을 통과하기 쉽다)

ggplot(rexam, aes(sample = exam)) + stat_qq()
ggplot(rexam, aes(sample = numeracy)) + stat_qq()

# 분산의 동질성 검정하기(정규성 검정 때와 같은 이유로 표본이 크면 분산비를 확인할 필요가 있다)

leveneTest(rexam$exam, rexam$uni) #두 번째 변수는 반드시 factor여야 하는 점에 주의!
leveneTest(rexam$numeracy, rexam$uni)


# 자료 변환

# 음수가 있는 경우 최댓값+1에서 각 값을 빼는 방식으로 모든 값을 양수로 바꾼다
# 이 때 값의 순서가 바뀌는 점을 꼭 인지한다
# 양수로 바뀐 값(혹은 양수만 있는 값)은 로그, 제곱근 또는 역수를 취해 변환한다
# 역수 변환을 쓰는 경우 대소관계가 다시 한 번 바뀐다다


  # 로그 변환

dlf$logday1 <- log(dlf$day1 + 1)
dlf$logday2 <- log(dlf$day2 + 1)
dlf$logday3 <- log(dlf$day3 + 1)

hist.logDay1 <- ggplot(dlf, aes(logday1)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Log Day 1', y = 'Density');hist.logDay1
hist.logDay2 <- ggplot(dlf, aes(logday2)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Log Day 2', y = 'Density');hist.logDay2
hist.logDay3 <- ggplot(dlf, aes(logday3)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Log Day 3', y = 'Density');hist.logDay3

  # 제곱근 변환

dlf$sqrtday1 <- sqrt(dlf$day1 + 1)
dlf$sqrtday2 <- sqrt(dlf$day2 + 1)
dlf$sqrtday3 <- sqrt(dlf$day3 + 1)

hist.sqrtDay1 <- ggplot(dlf, aes(sqrtday1)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Sqrt Day 1', y = 'Density');hist.sqrtDay1
hist.sqrtDay2 <- ggplot(dlf, aes(sqrtday2)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Sqrt Day 2', y = 'Density');hist.sqrtDay2
hist.sqrtDay3 <- ggplot(dlf, aes(sqrtday3)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Sqrt Day 3', y = 'Density');hist.sqrtDay3

  # 역수 변환

dlf$recday1 <- 1/(dlf$day1 + 1)
dlf$recday2 <- 1/(dlf$day2 + 1)
dlf$recday3 <- 1/(dlf$day3 + 1)

hist.recDay1 <- ggplot(dlf, aes(recday1)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Sqrt Day 1', y = 'Density');hist.recDay1
hist.recDay2 <- ggplot(dlf, aes(recday2)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Sqrt Day 2', y = 'Density');hist.recDay2
hist.recDay3 <- ggplot(dlf, aes(recday3)) + theme(legend.position = 'none') + geom_histogram(aes(y=..density..), colour = 'black', fill = 'white') + labs(x= 'Hygeine Sqrt Day 3', y = 'Density');hist.recDay3