탐색적자료분석 5장

1. DISTRESS.DAT 는 severe idiopathic respiratory distress syndrome이 있는 신생아의 태어날 때 몸무게(Kg)이다. *표는 죽은 아이이다. 자료파일은 ASCII 형식이다.

가. 생존 여부에 관계없이 전체 몸무게의 자료로 상자그림과 줄기그림을 그리고 대칭화 재표현을 한 후의 그림들과 비교하여라.

데이터를 읽기위해 작업디렉토리 설정 edit함수도 이용가능함

setwd("C://Users/R/Desktop/4학년 1학기/탐자분/datafile")
library(readr)
read.any <- function(text, sep = "", ...) {
  
  encoding <- as.character(guess_encoding(text)[1,1])
  setting <- as.character(tools::file_ext(text))
  if(sep != "" | !(setting  %in% c("csv", "txt")) ) setting <- "custom"
  separate <- list(csv = ",", txt = "\n", custom = sep)
  result <- read.table(text, sep = separate[[setting]], fileEncoding = encoding, ...)
  return(result)
}

baby=read.any("DISTRESS.DAT",stringsAsFactors = F)
typeof(baby)

## [1] "list"

baby2=unlist(baby)
die=baby2[1:27]
live=as.numeric(baby2[28:50])
die=as.numeric(gsub("[*]", "", die))
all=c(die,live)

상자그림과 줄기그림을 그리기 앞서 DISTRESS.DAT파일을 불러 gsub함수를 이용하여 *을 제거하는 작업을 하였으며 그 후 죽은 신생아와 생존한 신생아를 나누어 각 변수(die,live)에 저장하였다.

boxplot(all,main="Weight of all babies",col=("Cyan"))

stem(all,1)

## 
##   The decimal point is at the |
## 
##   1 | 0111222233334
##   1 | 566677778888999
##   2 | 001223344
##   2 | 56666778
##   3 | 0024
##   3 | 6

boxplot(log(all),main="Weight of all babies (log)",col=("Cyan"))

stem(log(all),1)

## 
##   The decimal point is 1 digit(s) to the left of the |
## 
##    0 | 350267
##    2 | 0136674
##    4 | 14572444677
##    6 | 046601499
##    8 | 2289244469
##   10 | 04805
##   12 | 29

대칭화 전 자료의 상자그림과 줄기그림을 확인하면 작은 쪽에 자료들이 많이 몰려 있으며 시각적으로 왜도가 존재함을 확인할 수 있다. 수업시간에 배운 것과 같이 log변환을 먼저 적용하였고 이전 자료보다 위 수염과 아래 수염의 길이가 비슷해 졌고 상자 내부의 중간값 위치가 이전보다 더 가운데에 위치하게 되었다. 또한, 줄기그림을 통해 자료의 분포가 전보다 대칭을 이루고 있음을 볼 수 있다. 왜도를 계산하는 함수(R script 첨부)를 만들어 수치로 확인한 결과 왜도가 확연하게 감소하였다.

skewness = function(x) {
  lower=fivenum(x)[2]
  med=fivenum(x)[3]
  upper=fivenum(x)[4]
  value=((upper-med)-(med-lower))/((upper-med)+(med-lower))
  return(value)
}
skewness(all)

## [1] 0.1834862

skewness(log(all))

## [1] 0.0427223

나. 사망과 생존 집단으로 나누어 상자 그림을 그리고 산포(spread)가 같도록 하는 재표현을 찾아 상자 그림을 그리고 두 집단을 비교하여라.

위의 병이 있는 신생아 중 몸무게가 몇 Kg이면 사망할 가능성이 많다고 예측할 수 있겠는가?

HsprA <- fivenum(live)[4] - fivenum(live)[2]
HsprB <- fivenum(die)[4] - fivenum(die)[2]
Hspr <- c(HsprA, HsprB)
M <- c(median(live),  median(die))
plot(log(M), log(Hspr), main="Spread vs. Level plot") #aspect ratio

(RegrLine <- lm(log(Hspr) ~ log(M)))       # ( ) is needed to print coeff's

## 
## Call:
## lm(formula = log(Hspr) ~ log(M))
## 
## Coefficients:
## (Intercept)       log(M)  
##     -0.5157       0.6854

abline(coef(RegrLine))

1 - coef(RegrLine)[2]

##    log(M) 
## 0.3145712

median(die)

## [1] 1.6

Spread-vs-level Plot을 통해 p값을 얻었고 루트변환과 log변환 중 더욱 적절한 방법을 적용하고자 하였다. 두 가지를 모두 시도해보기 전 H-spread 또한 산포의 정도를 나타내는 척도라 생각되어

#transformation X
(fivenum(live)[4]-fivenum(live)[2])-(fivenum(die)[4]-fivenum(die)[2])

## [1] 0.201

#log transformation
(fivenum(log(live))[4]-fivenum(log(live))[2])-(fivenum(log(die))[4]-fivenum(log(die))[2])

## [1] -0.04277151

#sqrt transformation
(fivenum(sqrt(live))[4]-fivenum(sqrt(live))[2])-(fivenum(sqrt(die))[4]-fivenum(sqrt(die))[2])

## [1] 0.02182665

#inverse transformation
(fivenum(-1/die)[4]-fivenum(-1/die)[2])-(fivenum(-1/live)[4]-fivenum(-1/live)[2])

## [1] 0.1047709

#inverse,sqrt transformation
(fivenum(-1/sqrt(die))[4]-fivenum(-1/sqrt(die))[2])-(fivenum(-1/sqrt(live))[4]-fivenum(-1/sqrt(live))[2])

## [1] 0.04321196

H-spread를 통해 산포를 확인하였고 두 집단의 산포 차이를 계산한 결과 루트변환이 두 집단의 산포를 가장 같도록 하는 재표현이라 판단되어 이를 토대로 루트변환한 상자그림을 그렸다.

boxplot(die,live,names=c("Die","Live"),main="Weight of babies (notch)",
        col=topo.colors(3,alpha=0.2), notch=T)

boxplot(sqrt(die),sqrt(live),names=c("Die","Live"),main="Weight of babies (notch, sqrt)",
        col=topo.colors(3,alpha=0.2), notch=T)
abline(1.36,0)
abline(median(die)-1.57*(fivenum(die)[4]-fivenum(die)[2])/length(die)^0.5,0)

원자료(왼쪽 그림)와 루트변환 자료(오른쪽 그림)에 대해 노치를 적용하여 그린 상자그림이다. 자료에서 노치를 적용하여 두 집단의 중간값의 차이가 유 의미한지 확인하고자 하였으나 판단이 애매하여 abline을 이용하였다. 선을 통해 확인한 결과 죽은 아이의 notch가 생존한 아이의 notch와 겹치지 않는 것을 확인할 수 있었다. 이는 죽은 신생아들의 중간 값과 생존한 신생아들의 중간 값이 95%의 신뢰도로 다르 다고 판단할 수 있을 것이다.

median(die)-1.57*(fivenum(die)[4]-fivenum(die)[2])/length(die)^0.5

## [1] 1.351031

median(die)+1.57*(fivenum(die)[4]-fivenum(die)[2])/length(die)^0.5

## [1] 1.848969

이 때 원자료의 죽은 신생아 notch구간을 구하면(1.3510,1.8489)이다. 위 신뢰구간을 바탕으로 강한 기준점을 잡으면 1.8489kg이하일 경우 생존하지 못할 가능성이 크다고 예상해 볼 수 있다.

2. R의 airquality 자료에서 Ozone 자료를 대칭화하기 위한 변환공식을 이용하여 변환 전과 변화 후를 비교하여라.

ozone=airquality$Ozone
summary(ozone)   

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1.00   18.00   31.50   42.13   63.25  168.00      37

stem(ozone)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##    0 | 1467778999
##    1 | 01112233334444666688889
##    2 | 0000111123333334478889
##    3 | 001222455667799
##    4 | 01444556789
##    5 | 0299
##    6 | 134456
##    7 | 13367889
##    8 | 024559
##    9 | 1677
##   10 | 8
##   11 | 058
##   12 | 2
##   13 | 5
##   14 | 
##   15 | 
##   16 | 8

boxplot(ozone, main="measurement of ozone")

#function
tran = function(x) {
  lower=fivenum(x)[2]
  med=fivenum(x)[3]
  upper=fivenum(x)[4]
  value=1-((2*med*((upper-med)+(lower-med)))/((upper-med)^2+(lower-med)^2))
  return(value)
}


tran(ozone)

## [1] 0.03378238

대칭화하기 위한 변환공식(재표현 프린트 참조)은 함수로 만들어 사용하였으며 함수 내부코드는 R script에 첨부하였다. 이를 통해 얻어진 값은 0.0338로 0에 가깝다 생각되어 log변환을 이용하였다.

log변환 전에는 아래수염의 길이와 윗 수염의 길이 차이가 많이 났으며 2개의 특이값이 존재하였다. 또한 중간값의 위치가 상자내부에서 아래에 치우쳐 있음을 확인할 수 있다. log변환 이후 어느정도 대칭화가 이루어졌다. 위, 아래 수염의 길이 차이가 줄어들었고 중간 값이 상자 내부의 중앙으로 이동하였다. 다만 특이값이 아래에 1개 존재하였는데 이는 값이 1인 원본 데이터가 존재하여 나타난 것으로 판단하였다. (log1=0)

3. 2016년의 삼성전자 주식 가격을 인터넷에서 찾아 줄기잎 전시와 상자그림을 그려 대칭성을 점검하고 필요하다면 대칭화하는 변환을 시행착오방법으로 찾아라.

###어떤 시행 착오로 찾아 갔는지 중간 과정을 전부 보여라. 변환 후의 줄기잎 전시와 상자그림을 그려 전, 후를 비교하여라. 변환 전, 후의 skewness 값을 비교하여라.

2016년의 삼성전자 주식 가격은 한국거래소에서 다운 받아 사용하였다. 원자료에 대한 줄기잎그림과 상자그림은 아래와 같다.

setwd("C://Users/R/Desktop/4학년 1학기/탐자분/datafile")

read.any <- function(text, sep = "", ...) {
  
  encoding <- as.character(guess_encoding(text)[1,1])
  setting <- as.character(tools::file_ext(text))
  if(sep != "" | !(setting  %in% c("csv", "txt")) ) setting <- "custom"
  separate <- list(csv = ",", txt = "\n", custom = sep)
  result <- read.csv(text, sep = separate[[setting]], fileEncoding = encoding, ...)
  return(result)
}
sam=read.any("samsung.csv",header=T)
attach(sam)
boxplot(sam$Price,col="pink",main="Samsung stock price")

줄기잎그림을 확인해보면 완전한 대칭은 아니지만 예상했던 것보다 자료들이 어느정도 대칭을 이루고 있다. 또한, 이는 상자그림에서도 확인할 수 있다. 특이값은 존재하지 않으며 중간값이 제1분위수보다 제3분위수쪽에 가까운 것을 볼 수 있다. 이는 왼쪽으로 약간 늘어진 분포일 것이다. 가장 적절한 대칭화 변환을 찾기 위해 수업시간에 배운 것처럼 4가지 흔한 방법들을 적용해보았다.

stem(sam$Price,2)

## 
##   The decimal point is 4 digit(s) to the right of the |
## 
##   112 | 60012788
##   114 | 566802466
##   116 | 233488112255589
##   118 | 1570247
##   120 | 585
##   122 | 035
##   124 | 568903356
##   126 | 01345678999013599
##   128 | 00112256880022445666999
##   130 | 000582
##   132 | 3
##   134 | 
##   136 | 517
##   138 | 06889
##   140 | 06693
##   142 | 156001
##   144 | 580
##   146 | 045669
##   148 | 19
##   150 | 0276778
##   152 | 77033599
##   154 | 01355837889
##   156 | 167778888991357
##   158 | 5679902367788
##   160 | 068244689
##   162 | 00157999
##   164 | 00003345990023
##   166 | 5577
##   168 | 071
##   170 | 68
##   172 | 7
##   174 | 68929
##   176 | 627
##   178 | 02803589
##   180 | 2592

stem(log(sam$Price),0.5)

## 
##   The decimal point is 1 digit(s) to the left of the |
## 
##   139 | 34444444
##   139 | 555566666777777777788888899999
##   140 | 0001122344444444
##   140 | 5555555555556666666666677777777777777888888889
##   141 | 03344
##   141 | 5555566667777778999
##   142 | 00000112233333444444
##   142 | 55555555666666666677777777778888888888888999999
##   143 | 00000000111111111112222223333344
##   143 | 5667778889999
##   144 | 0000000111

stem(sqrt(sam$Price))

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   106 | 133346770111234558889
##   108 | 11223344456799123489
##   110 | 25676678899
##   112 | 012344556666677891111222244556677888888
##   114 | 0000002455
##   116 | 8135
##   118 | 222336679244666
##   120 | 2348001127
##   122 | 05681222667889
##   124 | 111123346888991222222233345699
##   126 | 0111223444457800012233356
##   128 | 00011112223445556045569
##   130 | 0614
##   132 | 1224691345789
##   134 | 0112456

stem(-1/(sam$Price))

## 
##   The decimal point is 8 digit(s) to the left of the |
## 
##   -88 | 855430
##   -86 | 99333108755100
##   -84 | 9664433111987420
##   -82 | 9850830
##   -80 | 8633110
##   -78 | 887643211099888776422111100
##   -76 | 88665544332222000999652
##   -74 | 0
##   -72 | 3965
##   -70 | 655541108421
##   -68 | 999210533221
##   -66 | 527640
##   -64 | 99955422100998776422211
##   -62 | 988888887776541109998876666532000
##   -60 | 98877755000000099886666551
##   -58 | 76653162
##   -56 | 93221964321
##   -54 | 9987665432

stem(-1/sqrt(sam$Price))

## 
##   The decimal point is 5 digit(s) to the left of the |
## 
##   -94 | 21100
##   -93 | 877544332100
##   -92 | 8777554444333110
##   -91 | 98765410
##   -90 | 7544
##   -89 | 665543332110
##   -88 | 9998888877664444443322110000
##   -87 | 999888777777543
##   -86 | 6
##   -85 | 6421
##   -84 | 666553321
##   -83 | 987666210
##   -82 | 86666520
##   -81 | 6652222
##   -80 | 99888766665554211110
##   -79 | 9999999988877644433332211111
##   -78 | 998777666654411111110000
##   -77 | 998888532220
##   -76 | 9631
##   -75 | 766542100
##   -74 | 98776665433

par(mfrow=c(1,1))
hist(sam$Price,100,main="X")

hist(log(sam$Price),100,main="log")

hist((sqrt(sam$Price)),100,main="sqrt")

hist(-1/(sam$Price),100,main="-1/X")

hist(-1/sqrt(sam$Price),100,main="-1/sqrt(X)")

skewness = function(x) {
  lower=fivenum(x)[2]
  med=fivenum(x)[3]
  upper=fivenum(x)[4]
  value=((upper-med)-(med-lower))/((upper-med)+(med-lower))
  return(value)
}
skewness(Price)

## [1] -0.1333333

skewness(log(Price))

## [1] -0.1893818

skewness(sqrt(Price))

## [1] -0.1614364

skewness(-1/(Price))

## [1] -0.2446197

skewness(-1/sqrt(Price))

## [1] -0.2171244

줄기잎그림은 Plot형태로 확인할 수 없어 히스토그램을 그려 원자료와 여려 변환후의 분포를 확인하였다. 먼저 위에서 언급한 바와 같이 원자료가 어느정도 대칭을 이루고 있어 여러 변환들을 적용하여도 다른 문제들에 비해 눈에 띄게 대칭이 나타나고 있는 것을 확인하기는 어렵다. 그래서 대칭의 정도를 확인하는 왜도를 통해 어느 변환이 가장 대칭을 이루게 하는지 찾아보려고 하였다.

x=c()
for (i in seq(-5,5,0.1)) {
  if(i>0){
    
    x = c(x,skewness(Price^i))
  }
  
  else if(i<0){
    x = c(x,skewness(-(Price^i)))
  }
  
  else{
    x = c(x,skewness(log(Price)))
  }
}
y=cbind(seq(-5,5,0.1),x)
plot(y,xlab="p",ylab="skewness")
abline(0,0)

많은 p들을 대입하는 것은 시간이 너무 많이 소요되어 for 문을 이용하여 p에 -5부터 5까지 0.1 간격으로 넣어 왜도를 계산한 뒤 이를 Plot으로 나타냈고 p가 3일 때 왜도가 0에 가까워지고 있음을 확인하였다. (앞서 만든 skewness함수를 통해 확인한 결과 세제곱 변환을 하였을 때 왜도가 0에 근접하였다.)

stem(sam$Price^3)

## 
##   The decimal point is 17 digit(s) to the right of the |
## 
##   14 | 34455777011123444777899
##   16 | 1111222345679902569
##   18 | 234334557788
##   20 | 0112233444455679900001123445566777888999
##   22 | 0002467
##   24 | 48
##   26 | 132334488
##   28 | 02790223
##   30 | 245144557
##   32 | 5089
##   34 | 28990668
##   36 | 0025556799157889
##   38 | 04555666666891289
##   40 | 011234777880469
##   42 | 0024455691
##   44 | 00011114445889912
##   46 | 20224
##   48 | 047
##   50 | 75
##   52 | 2458
##   54 | 416
##   56 | 1462468
##   58 | 125825

boxplot(sam$Price^3,main="Samsung stock price (p=3)",
        col=gray.colors(3,alpha=0.1))

세제곱 변환의 줄기잎그림과 상자그림은 위와 같다.

왜도가 0에 근접하였기 때문에 당연히 가장 대칭적인 분포를 가진 변환일 것이라 생각했지만 결과는 예상과 다르게 상자내부는 적절하게 대칭을 이루게 되었지만 수염의 길이차이가 큰 상자를 확인하였다. 그 이유는 다음과 같다.

skew=((Hu-M)-(M-Hl))/((Hu-M)+(M-Hl)) (M=median,Hu=제3분위수,Hl=제1분위수) 왜도의 계산식을 보면 3가지(median, 제1,3분위수)만 쓰이고 있다. 상자그림의 경우 중간값과 상자 위 아래(제1분위수, 제3분위수)까지의 거리 차이가 최소화될 때 작은 왜도 값을 갖게 되는 것이다. 세제곱 변환의 상자그림을 확인해보면 위 수염과 아래 수염의 길이는 많이 다르지만 상자 내부의 중간값이 거의 가운데에 위치하고 있다.

즉, H-spread 내에서 가장 대칭적인 분포를 가진다고 볼 수 있다. 또한 대칭화 변환 공식을 이용하여 값을 구하면 p=3.31로 3에 가까운 값을 구할 수 있었다. 이 또한 p를 구하는 과정에서 3가지(median, 제1,3분위수)만 쓰였기 때문이다. 하지만 전체적인 데이터를 고려한다면 가장 대칭적인 분포를 가진다고 단언하기는 힘들다. 시각적으로 위에서 흔히 쓰이는 4가지 변환이 오히려 전체적인 분포는 더 대칭적으로 보인다.

par(mfrow=c(2,3))
boxplot(sam$Price,main="X",col=rainbow(1,alpha=0.1))
boxplot(log(sam$Price),main="log",col=rainbow(2,alpha=0.1))
boxplot((sqrt(sam$Price)),main="sqrt",col=rainbow(3,alpha=0.1))
boxplot(-1/(sam$Price),main="-1/X",col=rainbow(4,alpha=0.1))
boxplot(-1/sqrt(sam$Price),main="-1/sqrt(X)",col=rainbow(5,alpha=0.1))

원 자료, 4가지 변환 그리고 세제곱 변환의 상자그림을 확인해보자 세제곱 변환 상자그림의 중간값이 다른 상자그림에 비해 가장 상자 내부 중앙에 위치하고 있는 것을 확인할 수 있다.

전체적인 수염의 길이까지 고려한다면 오히려 원자료가 가장 대칭적인 분포를 띄고 있는 듯하다. 만약 제1분위수와 제3분위수 사이의 데이터를 중요시 여긴다면 세제곱변환이 가장 적절 할 것이다.

그렇지만 모든 데이터를 중요시 여겨 사용하고자 한다면 이 자료는 원본 데이터를 사용하는 것이 가장 나아 보인다. 실제로 세제곱 변환을 제외하면 원본 데이터가 가장 낮은 왜도를 가진다.