EDA HIGHLIGHT
- Ch.1 R Language
- Ch.2 Introduction to EDA
- Ch.3 Stem and Leaf
- Ch.4 Numerical Summary and Box Plot
- Ch.5 Data Re-Expression
- Ch.6 QQ-Plot
- Ch.7 2원 자료, 빈도 표의 탐색
- Ch.8 시계열 자료의 탐색
- Ch.9 로버스트 선형회귀
- Ch.10 이변량 자료 탐색
library(tidyverse)
options(repr.plot.width=16, repr.plot.height=8)
- 생략
자료분석은 대체로 두 가지 단계로 나뉜다.
-
탐색적 자료분석(데이터의 구조와 특징 파악)
-
확증적 자료분석(모형, 재현성 평가)
-
EDA는 데이터 특징과 내재하는 구조적 관계를 알아내기 위한 기법들을 통칭한다.
-
데이터를 특정한 모형에 적합시키기 보다는 데이터를 있는 그대로 보려는 데에 중점을 둔다.
EDA에서는 네 가지 주제가 때로는 홀로, 때로는 얽혀서 나타난다.
-
저항성(resistance)의 강조(예를 들어 평균보다는 일부자료의 파손에 저항적인 중위수가 바람직한 대표값 측도로서 선호된다.)
-
잔차(residual) 계산
-
자료변수의 재표현(re-expression)을 통한 다각적 시도
-
그래프를 통한 현시성(revelation)
-
자료분석은 탐색적 자료분석과 확증적 자료분석의 두 단계로 나눌 수 있다.
-
EDA는 자료의 구조 및 특징의 파악을 목적으로 한다. 이를 위하여 효과적이고 신뢰성 있는 데이터의 요약과 그래프적 기법이 사용된다.
-
EDA의 네 개 주제는 저항성, 잔차, 재표현, 현시성이다.
-
통계적 모형은 '진실'로서가 아니라 '주요 사례'로서 의미가 있다. 또한 모형과 데이터는 일방통행이 아닌 쌍방통행으로 이해되어야 한다. 데이터와 분석도 사이클을 이룬다.
데이터의 값을 십 단위인 줄기(stem)와 일 단위인 잎(leaf)으로 분리한다.
- 장점
-
사분위수나 중앙값을 찾기 쉽다.
-
분포의 전체적인 모양(봉우리 개수, 대칭분포, 치우친 방향)을 쉽게 알 수 있다.
-
이상값 유무를 파악할 수 있다.
- 단점
- 자료가 많은 경우에는 부적합하고 자료의 수가 50개 이하일 때 적합하다.
exam1 = read.table('dataset/EDA/exam1.txt', header = TRUE)
head(exam1)
stem(exam1$score)
- 이 줄기그림의 주요 특징은 이봉분포의 모습을 보인다는 점이다. 이것은 자료가 2개의 군집으로 되어 있음을 말한다.
exam1[exam1$hw == 0, ]$score %>% stem # 과제 미제출자
exam1[exam1$hw == 1, ]$score %>% stem # 과제 제출자
과제 제출 여부에 따라서 줄기-잎 그림을 두 개 그려본 결과 흥미로운 점을 확인했다.
과제 미제출자에 경우 여전히 이봉분포를 띄고 있으나 과제 제출자의 경우 단봉분포를 띄고 있다.
이것은 과제물 미제출 그룹이 서로 이질적인 어떤 두 집단(50점대를 중심으로 하는 집단과 20점대를 중심으로 하는 집단)의 혼합임을 말하여 준다.
아마 전자는 학업습관이 불성실하나 학업능력은 우수한 집단이고 후자는 학업습관과 학업능력 모두 좋지 않은 학생들의 집단이 아닐까 생각해본다.
Tip :R에서 줄기 수를 줄이거나 늘이려면 stem() 함수 내 scale 파라미터를 조정하면 된다.
stem(exam1$score, scale = 0.5) # 줄기 수를 절반으로 줄인다.
이 줄기 그림은 단봉분포의 형태를 취한다. 이것은 너무 단순하여 이 자료의 주요 특성을 잃은 것으로 볼 수 있다.
즉 2개의 봉우리를 구분하지 못하고 1개만 본 것이다.
stem(exam1$score, scale = 2) # 줄기 수를 두 배로 늘인다.
줄기 수를 늘였더니 더 많은 봉우리를 볼 수 있다.
일반적으로 줄기 수를 늘이면 늘일수록 많은 수의 봉우리를 보게 되고 그 반대로 줄기 수를 줄이면 줄일수록 적은 수의 봉우리를 보게 된다.
hist(exam1$score)
| 줄기 - 잎 | 히스토그램 | |
|---|---|---|
| 정보 손실 | 손실되지 않음 | 손실됨 |
| 줄기 수 변환 | 쉬움 | 어려움 |
| 구간의 폭 | 조정 불가능 | 조정 가능 |
-
줄기 그림은 히스토그램과 마찬가지로 자료 분포의 특성을 그래프화한 것이다. 줄기 그림은 히스토그램에 비하여 정보의 보전 면에서 우수하며 쉽게 구간(줄기) 수를 늘이거나 줄일 수 있다. 그러나 구간(줄기)의 선정시 제약이 따른다.
-
적절한 줄기 그림을 그리기 위하여 여러 개의 그림을 그려보고 비교해 보아야 한다. 계획된 시행착오가 필요하다.
-
줄기 그림에서는 다음과 같은 자료의 특성을 관찰할 수 있다.
- 군집의 수
- 집중도가 높은 구간
- 대칭성 여부
- 자료의 범위 및 산포
- 특이점의 존재여부
- 한 쪽 꼬리가 긴 분포에서 평균값은 쉽게 휘둘리지만 중앙값은 쉽게 휘둘리지 않다. 따라서 중앙값이 대표값으로서 적합하다.
$\begin{cases} X_{\frac{N+1}{2}} \qquad if N = odd\\ \\ \dfrac{X_{\frac{N}{2}} + X_{\frac{N}{2} + 1}}{2} \qquad if N = even \end{cases}$
중간값의 깊이 $d(M)$은
$d(M) = \dfrac{(N+1)}{2}$
아래 4분위수 : $H_L$ 중간값 : $M$ 위 4분위수 : $H_U$
다섯 숫자 요약 = (min, $H_L, M, H_U$, max) = (최솟값, 제 1사분위수, 중앙값, 제 3사분위수, 최댓값)
summary(exam1$score)
왜도 = $SKEW = \dfrac{(H_U - M) - (M - H_L)}{(H_U - M) + (M - H_L)}$
왜도 < 0 이면 왼쪽으로 기울어진 분포
왜도 > 0 이면 오른쪽으로 기울어진 분포
low = quantile(exam1$score, c(1/2, 1/4, 1/8, 1/16), type = 8) %>% as.numeric %>% round(2)# 중간값, 아래 4분위수, 아래 8분위수, 아래 16분위수
high = quantile(exam1$score, c(1/2, 3/4, 7/8, 15/16), type = 8) %>% as.numeric %>% round(2) # 중간값, 위 4분위수, 위 8분위수, 위 16분위수
values = cbind(low, high, (low+high)/2, high-low)
colnames(values) = c('low', 'high', 'mid', 'spr')
rownames(values) = c('M', 'H', 'E', 'D')
values
- 표준화 변환이란 통상적으로 한 자료묶음의
평균이 0, 표준편차가 1이 되도록 하는 선형변환을 말한다.
$x_1, x_2, \dots, x_n$을 자료 값이라고 할 때 이것의 표준화변환 $z_1, z_2, \dots, z_n$은 다음과 같이 정한다.
$z_i = \dfrac{x_i - \bar{x}}{s_x}, \, i = 1,2, \dots,n. \qquad(1)$
그런데 (1)은 로버스트하지 않은 $\bar{x}$와 $s_x$에 의존하므로 EDA의 관점에서는 믿고 사용하기 어렵다. [^1]
- 왜냐하면 표본평균과 표본표준편차는 극단적인 이상점에 의해 크게 변동될 수 있기 때문이다.
- 그러나 중앙값 또는 사분위수범위(IQR)은 비교적 로버스트하다.
즉, 평균 $\bar{x}$ 대신에 중앙값 $med_x$를, 표준편차 $s_x$ 대신에 사분위수범위 $IQR$을 보정한 $\tilde{\sigma_x} = \dfrac{IQR}{1.35}$을 쓰는 것이 좋을 것이다.
따라서 로버스트 표준화 변환은 다음과 같다.
$\bar{z_i} = \dfrac{x_i - med_x}{\tilde{\sigma_x}}\, i = 1,2, \dots,n. \qquad(1)$
표준화 변환을 사용하는 예시 상황은 다음과 같다.
par(mfrow = c(1,2))
X_group <- rnorm(100, 40, 10)
Y_group <- c(rnorm(90,40,10), rnorm(10,80,5))
z_X <- (X_group-mean(X_group))/sd(X_group)
z_Y <- (Y_group-mean(Y_group))/sd(Y_group)
hist(z_X, breaks = seq(-6, 6, 0.5), freq = F, ylim = c(0, 0.7), main = 'Standardized X')
hist(z_Y, breaks = seq(-6, 6, 0.5), freq = F, ylim = c(0, 0.7), main = 'Standardized Y')
par(mfrow = c(1,2))
robust_z_X <- (X_group-median(X_group))/(IQR(X_group)/1.35)
robust_z_Y <- (Y_group-median(Y_group))/(IQR(Y_group)/1.35)
hist(robust_z_X, breaks = seq(-6, 6, 0.5), freq = F, ylim = c(0, 0.7), main = 'Robust Standardized X')
hist(robust_z_Y, breaks = seq(-6, 6, 0.5), freq = F, ylim = c(0, 0.7), main = 'Robust Standardized Y')
- 선형변환 $ax+b\, (a > 0)$은 분포의 형태를 바꾸지 않는다. 그러나 비선형변환은 분포의 형태를 바꾼다.
- 변환의 사다리는 $x^p$ 꼴의 파워(power, 멱승)형 변환을 일컫는데 변환의 사다리를 내려가면 $(p < 1)$ 오른쪽 꼬리가 짧아진다. $p=0$에 해당하는 변환은 로그변환이다.
- 자료의 재표현은 분포의 대칭화를 위하여, 또는 자료묶음들의 산포를 균일화하기 위한 목적으로 실행된다.
- 자료의 재표현은 자료 해석을 풍부하게 한다.
par(mfrow = c(1,2))
x1 <- rnorm(40, 100, 15)
qqnorm(x1)
qqline(x1)
x2 <- rnorm(4000, 100, 15)
qqnorm(x2)
qqline(x2)
par(mfrow = c(1,2))
x1 <- c(rnorm(20, 70, 15), rnorm(20, 130, 15))
qqnorm(x1)
x2 <- c(rnorm(2000, 70, 15), rnorm(2000, 130, 15))
qqnorm(x2)
par(mfrow = c(1,2))
x1 <- rnorm(38, 100, 15)
outliers <- c(25, 175)
qqnorm(c(x1, outliers))
x2 <- rnorm(3800, 100, 15)
qqnorm(c(x2, outliers))
par(mfrow = c(1,2))
x1 <- runif(40, 80, 120)
qqnorm(x1)
x2 <- runif(4000, 80, 120)
qqnorm(x2)
par(mfrow = c(1,2))
x1 <- c(rexp(20,1), -rexp(20,1))
qqnorm(x1)
x2 <- c(rexp(2000,1), -rexp(2000,1))
qqnorm(x2)
par(mfrow = c(1,2))
x1 <- exp(rnorm(40, 5, 1))
qqnorm(x1)
x2 <- exp(rnorm(4000, 5, 1))
qqnorm(x2)
par(mfrow = c(1,2))
x1 <- 1500 - exp(rnorm(40,5,1))
qqnorm(x1)
x2 <- 1500 - exp(rnorm(40,5,1))
qqnorm(x2)
leukemia <- c(1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)
leukemia_quant = ((1:length(leukemia)) - 1/3) / (length(leukemia) + 1/3)
x <- -log(1-leukemia_quant)
y <- sort(leukemia)
plot(y ~ x)
library(lattice)
qqmath(leukemia)
평균과 같은 추정치는 특이값에 의해 쉽게 영향을 받으므로 EDA적 관점에서 보면 바람직하지 않다.
따라서 평균을 중간값으로 대체하고 다듬을 필요가 있다.
R 함수
{r}
medpolish(world_temp)7장 이원자료빈도표,
7.4절
8장 시계자료탐색 표 데이터주고 3할 구하기 등, 스필트해보기 등,
149페이지 표3에 대한거 스필트해보기 등,
147페이지
19.5 봉우리인지 골자기인지, 봉우리를 스필트해라,
8.5 자기 상관과 가중치? 여기서 그림이 의미하는 바, hf? 자기상관 교차상관,
9장
165 lms lts ,
167 lts,
10장 그림을 그려서 알수잇는게 뭐냐
182 밀감 추정? 데이터 주고 커널 밀도 추정해봐라 계산해서 그려봐라 계산기 필요! 도자기 플룻 한번 보기
11장 가변량 자료 여러가지 데이터를 주고 평행자료 같은거 주면 해석하기 그림을 그리고 해석하라 알수잇는게 뭐냐 정도
world_temp <- read.table('dataset/EDA/WorldTemperature_Mean.txt', header = TRUE)
world_temp
polish <- medpolish(world_temp)
polish
comparison <- (matrix(polish$row, ncol = 1) %*% matrix(polish$col, nrow = 1))/polish$overall
comparison
plot(polish$residual ~ comparison, xlim = c(-15, 15), ylim = c(-15, 15))
round(polish$residuals[order(polish$row),],1)
house <- read.table('dataset/EDA/household.txt', header = TRUE)
house
house_mp <- medpolish(house)
house_mp
house_mp$residual
matrix(house_mp$row, ncol = 1) %*% matrix(house_mp$col, nrow = 1) / house_mp$overall
comparison <- matrix(house_mp$row, ncol = 1) %*% matrix(house_mp$col, nrow = 1) / house_mp$overall
plot(house_mp$residuals ~ comparison)
가법적이지 않다. (그렇다고 승법적인 것은 아니다.)
이제 어떻게 할 것인가? 로그 변환을 취한 뒤 다시 살펴보자
house_log_mp <- medpolish(log(house))
house_log_mp
comparison <- matrix(house_log_mp$row, ncol = 1) %*% matrix(house_log_mp$col, nrow = 1) / house_log_mp$overall
plot(house_log_mp$residuals ~ comparison)
이제는 경향선을 찾기 어렵다.
로그 변환 자료에 대한 중간값 다듬기 결과를 살펴보자.
round(house_log_mp$residuals, 2)
house_log_mp
UCBAdmissions
Tab1 <- UCBAdmissions[1,,]
Tab1
Tab2 <- UCBAdmissions[2,,]
Tab2
Tab <- Tab1 + Tab2
Tab
barplot(Tab, legend = rownames(Tab))
Tab_col <- apply(Tab, 2, sum)
Tab_col
Tab_c <- Tab %*% diag(1/Tab_col) * 100
Tab_c
colnames(Tab_c) <- c('A', 'B', 'C', 'D', 'E', 'F')
rownames(Tab_c) <- c('Male', 'Female')
Tab_c
barplot(Tab_c, legend = rownames(Tab_c))
par(mfrow = c(1,2))
barplot(t(Tab), beside = T, legend = colnames(Tab))
Tab_row <- apply(Tab,1,sum)
Tab_r <- diag(1/Tab_row) %*% Tab * 100
colnames(Tab_r) <- colnames(Tab_c)
rownames(Tab_r) <- rownames(Tab_c)
barplot(t(Tab_r), beside = T, legend = colnames(Tab))
- 남학생은 A, B에 50% 이상 집중되어 있는 것과 달리 여학생은 C, D, E, F에 90% 이상이 집중되어 있다.
Tab_M <- UCBAdmissions[,1,]
addmargins(Tab_M)
Tab_F <- UCBAdmissions[,2,]
addmargins(Tab_F)
이제 남학생의 합격률과 여성의 합격률을 비교해보자.
par(mfrow = c(1,2))
Tab_M_col <- apply(Tab_M, 2, sum)
Tab_M_C <- Tab_M %*% diag(1/Tab_M_col) * 100
colnames(Tab_M_C) <- c('A', 'B', 'C', 'D', 'E', 'F')
barplot(Tab_M_C, legend = rownames(Tab_M_C), main = 'Male')
Tab_F_col <- apply(Tab_F, 2, sum)
Tab_F_C <- Tab_F %*% diag(1/Tab_F_col) * 100
colnames(Tab_F_C) <- colnames(Tab_M_C)
barplot(Tab_F_C, legend = rownames(Tab_F_C), main = 'Female')
C, D, E, F 학과에서 남학생과 여학생의 합격률은 비슷하고
A, B 학과에서는 여학생의 합격률이 더 높은 것으로 보인다.
par(mfrow = c(1,2))
mosaicplot(~Dept+Gender, data = UCBAdmissions, color = T)
mosaicplot(~Gender+Dept, data = UCBAdmissions, color = T)
par(mfrow = c(1,2))
t <- seq(1:60)
s <- sin(2*pi*t/12) + sin(2*pi*t/30)
a <- rnorm(60, 0, 1)
x <- s+a
plot(t, x, type = 'l', lty = 'dotted')
lines(t, s, col = 'red')
plot(t, smooth(x), type = 'l', lty = 'dotted')
lines(t, s, col = 'red')
이동평균(이동중앙값)에서 $K$가 작으면 $smooth_t$ 계열이 평탄하게 된다. 반대로 $K$를 작게하면 $smooth_t$ 계열에 잦은 변동이 생긴다.
이동평균은 특이점에 취약하므로 평균을 중간값으로 대치한 이동중간값이 선호된다.
t <- 1:12
x <- c(2.1, 9.8, 19.5, 22.5, 16.6, 16.1, 18.5, -3.4, 8.9, -25.2, -14.0, -0.4)
runmed(x, k=3)
par(mfrow = c(1,2))
plot(t, x, type = 'l')
plot(t, runmed(x, k = 3), type = 'l')
smooth(x, kind = c("3RSS"), twiceit = T)
Telemetric <- read.table("dataset/EDA/Telemetric.txt", header = T)
head(Telemetric)
par(mfrow = c(2,2))
plot(Telemetric$temperature, type = 'l', lty = 'dotted')
plot(smooth(Telemetric$temperature), type = 'l', lty = 'dotted')
plot(smooth(Telemetric$temperature, kind = c('3RS3R'), twiceit = TRUE), type = 'l', lty = 'dotted')
plot(smooth(Telemetric$temperature, twiceit = TRUE), type = 'l', lty = 'dotted')
export <- read.table('dataset/EDA/Export_1988.txt', header = TRUE)
series <- ts(export$Series/1000, start = c(1988,1), frequency = 12)
par(mfrow = c(1,2))
plot(series)
plot(log(series))
decomp_out <- decompose(log(series))
decomp_out$season
round(decomp_out$season[1:12], 4)
plot(decomp_out$trend)
decomp_out$trend
x_tilde <- log(series) - decomp_out$seasonal # x_tilde = x - S_5
adjusted <- exp(x_tilde) # 변환
adj_series <- ts(adjusted, start = c(1988,1), frequency = 12)
adj_series
plot(adj_series)
geyser$waiting
par(mfrow = c(2,2))
library(MASS)
plot(geyser$waiting)
acf(geyser$waiting)
plot(geyser$waiting)
acf(geyser$duration)
round(acf(geyser$waiting)$acf[1:5],3)
round(acf(geyser$duration)$acf[1:5],3)
waiting과 duration 사이의 관계는 어떨까? 이것은 다음과 같이 정의되는 교차상관함수(CCF)로 알 수 있다.
ccf(geyser$waiting, geyser$duration)
round(cbind(ccf(geyser$waiting, geyser$duration)$lag, ccf(geyser$waiting, geyser$duration)$acf), 3)
par(mfrow = c(1,2))
plot(geyser$duration, geyser$waiting)
waiting_1 <- c(geyser$waiting[2:299], NA)
plot(waiting_1 ~ geyser$duration)
library(foreign)
ee <- read.spss('dataset/EDA/EEstock2000.sav')
x <- ee$change
par(mfrow = c(1,2))
plot(x, type = 'l', ylim = c(-15, 15))
acf(x, ylim = c(-1, 1))
par(mfrow = c(1,2))
m <- mean(x)
sd <- sd(x)
hist(x, freq=F, nclass = 20, main = 'Daily Stock Change')
curve(dnorm(x, m, sd), add = T)
ms <- read.spss('dataset/EDA/MSstock2000.sav')
x <- ms$change
m <- mean(x)
sd <- sd(x)
hist(x, freq = F, nclass = 20, xlim = c(-15, 15), main = 'Daily Stock Change')
curve(dnorm(x, m, sd), add = T)
set.seed(1234567)
x <- seq(1,10)
y <- 2.5 + 0.5*x + rnorm(10,0,1)
y[10] <- -10
par(mfrow=c(1,2))
plot(y~x)
library(MASS)
m0 <- lm(y ~ x)
m1 <- rlm(y ~ x)
m2 <- lqs(y ~ x)
plot(y ~ x)
abline(m0$coef)
abline(m1$coef, lty = 'dotted', col = 'blue')
abline(m2$coef, lty = 'dotted', col = 'red')
par(mfrow = c(1,2))
library(MASS)
data(geyser)
hist(geyser$duration)
density(geyser$duration)
hist(geyser$duration, freq = F)
lines(density(geyser$duration), lty = 2)
ker <- function(n
geyser$duration %>% length()
1/299 * sum(
d <- density(rnorm(10000))
plot(d)
approx(d$x, d$y, xout = c(-2, 0, 2))
d$x
[^1] 로버스트(robust) 한 통계량은 이상치/에러값으로 부터 영향을 크게 받지 않는 (건장한) 통계량