탐색적자료분석 8장
JUNG KYU HYEONG
2020 7 10
1. FRIDAY.DAT 자료는 영국에서 1986년 금요일에 발생한 자동차 사고 사망자 수이다. 첫 열은 하루를 24시간으로 표시한 것이고, 두 번째 열은 사망자 수이다. 시간에 따른 사망자 수를 평활하고 시간과 사망자 수의 관계를 서술하여라.
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoosetwd("C://Users/R/Desktop/4학년 1학기/탐자분/datafile")
friday=read.table("FRIDAY.DAT")
colnames(friday)=c("hour","death")
boxplot(friday$death,main="Number of Death")
summary(friday$death)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 138.0 858.8 2488.5 2310.1 3052.2 4923.0FRIDAY.DAT의 자료는 시간대에 따른 사망자 수의 데이터로 첫번째 열에는 시간대, 두번째 열에는 사망자 수의 값으로 이루어져 있어 열 이름(hour, death)을 지정해주었다. 시간대는 0-1과 같은 형식으로 이루어져 있어 Plot에서는 0-1은 1, 1-2는 2와 같은 형식으로 표현되었다.
상자그림을 통해 사망자 수를 살펴보았다 상자그림을 통해 많은 정보를 얻을 수는 없지만 자료의 분포는 어느정도 확인해 볼 수 있다. 특이 값은 나타나고 있지 않지만 최솟값과 최댓값의 차이가 꽤 있어 보인다.
또한, 중간 값이 3분위수와 가까운 값에 있음을 확인하였다. 만약 시간대의 흐름에 따라 최소에서 최대로 가는 방향이라면 자료를 해석하는데 있어 큰 어려움이 없지만 흐름과 다른 패턴의 점들이 분포한다면 이를 유의하여 시계열 자료를 분석해야 한다.
시간에 따른 사망자 수를 평활하는 방법은 여러가지가 존재하지만 수업시간에 배운 것처럼 3RSSH, twice방법과 4253H, twice방법을 이용하였다.
이 두 방법을 이용하기 전 Hanning이 어떤 의미를 갖는지 궁금하여 3RSS twice의 평활방법부터 시작하였다. 여러 평활은 smooth()함수를 통해 결과 값을 쉽게 얻을 수 있다. kind와 twiceit 변수를 이용하면 수업시간에 배운 대부분의 평활법을 사용할 수 있다.
a3rssh2=function(x){
n=length(x)
x3rss=smooth(x,kind="3RSS")
x3rssh <- vector("numeric",n)
x3rssh2 <- vector("numeric",n)
for (i in 2:(n-1)) x3rssh[i] <- x3rss[i-1]/4 + x3rss[i]/2 + x3rss[i+1]/4
x3rssh[1] <- x3rss[1]; x3rssh[n] <- x3rss[n]
rough=x-x3rssh
x3rss2=smooth(rough,kind="3RSS")
for (i in 2:(n-1)) x3rssh2[i] <- x3rss2[i-1]/4 + x3rss2[i]/2 + x3rss2[i+1]/4
x3rssh2[1] <- x3rss2[1]; x3rssh2[n] <- x3rss2[n]
end=x3rssh+x3rssh2
return(end)
}
ts.friday=ts(friday$death)
smooth.f=smooth(ts.friday,kind="3RSS",twiceit = T)
plot(friday$death,type="l",ylab="Death",lty="dotted",xlab="Hour",ylim=c(0,5000),lwd=2);par(new=T)
plot(smooth.f,xlab="",ylab="",ylim=c(0,5000),col="green",lwd=2) #3RSS3RSS
legend(x = 2, y = 5000, c("Raw", "3RSS twice"),
col = c("black","green"),
lty=c(3,1),lwd=2)
3RSS, twice방식을 통해 얻은 평활은 위와 같다. 점선은 원자료를 나타내며 초록색선이 평활을 통해 얻은 결과이다. 전체적인 추세는 어느정도 확인할 수 있지만 표시한 큰 값의 두 점의 특징이 잘 원 자료에 비해 뚜렷하지 않고 훨씬 평평해진 것을 확인할 수 있다.
만약 차이를 오차라고 생각한다면 적절한 평활이라 판단을 내릴 수 있지만 실제 특성을 잡지 못한 것은 아닌지 의심해봐야 한다.
h2=a3rssh2(friday$death)
plot(friday$death,type="l",ylab="Death",lty="dotted",xlab="Hour",ylim=c(0,5000),lwd=2);par(new=T)
plot(h2,type="l",ylab="Death",xlab="Hour",ylim=c(0,5000),col="red",lwd=2)
legend(x = 2, y = 5000, c("Raw", "3RSSH,twice"),
col = c("black","red"),
lty=c(3,1),lwd=2)
다음으로 3RSSH, twice방식을 적용한 평활의 결과이다. 추세선은 비슷하게 나타나고 있지만 위에서 표시한 두 점의 특징이 좀 더 실리게 되었다. 위에서 언급한 첫번째 지점에서 평평하던 점이 약간 상승하였고 두번째 지점은 확연하게 원자료와 가까워진 것을 확인할 수 있다.
또한 Hanning을 적용하기 전에는 평평한 부분이 종종 존재하였는데 Hanning을 통해 평평한 부분이 전에 비해 사라졌다. Hanning은 해당 시점의 가중치를 더 반영하는 특징을 가지는데 이는 Plot을 통해서도 확인해볼 수 있었다.
library(sleekts)
h4253=sleek(friday$death)
plot(friday$death,type="l",ylab="Death",lty="dotted",xlab="Hour",ylim=c(0,5000),lwd=2);par(new=T)
plot(h4253,type="l",ylab="Death",xlab="Hour",ylim=c(0,5000),col="blue",lwd=2)
legend(x = 2, y = 5000, c("Raw", "4253H,twice"),
col = c("black","blue"),
lty=c(3,1),lwd=2)
마지막으로 4253H, twice방법을 통해 얻은 결과이다. 무엇보다 이전의 평활에 비해 곡선의 형태가 두드러지는 것을 확인할 수 있다. 이는 원 자료에는 존재하지 않는 값으로 smoothing하였기 때문에 더 부드러운 곡선으로 나타나고 있다고 생각하면 될 것이다. 곡선의 형태로 나타나고 있기 때문에 이전의 평활법에 비해 추세를 쉽게 확인할 수 있다. 그렇기 때문에 4253H, twice 평활법이 가장 적절한 방법이라고 생각된다.
하지만 다른 평활법에 비해 8-9시 사이의 값, 16-17시 사이의 값 그리고 23-0시의 값이 원자료에 비해 많이 감소되었다. 이는 주의해야할 점이다. 8-9시, 16-17시는 출퇴근 시간으로 교통량이 많기 때문에 이 값들은 충분의 의미가 있는 값일수도 있지만 평활을 하면서 특징이 감소하게 되었다. 이러한 특징은 조금 더 반영될 필요성이 있어 보인다. 23-0시 또한 모든 평활법에서 급격히 감소하고 있다. 이 시간대에 급격하게 사망자 수가 늘어난 것은 전체적인 흐름에서 꽤나 특이한 특징이기 때문에 영국의 운전문화와 연관지어 해석할 필요성이 있다. 검색을 통해 영국의 몇몇 주차공간은 밤 12시 이후 주차 요금을 부과한다는 기사를 보았는데 이러한 모습이 약간 영향을 미치지 않았을까? 예상해 볼 수 있다. 새벽시간대에 사망자 수가 적은 것은 너무나도 당연한 일이다.
2. R의 lynx 자료를 평활하고 분석하여라.
lynx## Time Series:
## Start = 1821
## End = 1934
## Frequency = 1
## [1] 269 321 585 871 1475 2821 3928 5943 4950 2577 523 98 184 279 409
## [16] 2285 2685 3409 1824 409 151 45 68 213 546 1033 2129 2536 957 361
## [31] 377 225 360 731 1638 2725 2871 2119 684 299 236 245 552 1623 3311
## [46] 6721 4254 687 255 473 358 784 1594 1676 2251 1426 756 299 201 229
## [61] 469 736 2042 2811 4431 2511 389 73 39 49 59 188 377 1292 4031
## [76] 3495 587 105 153 387 758 1307 3465 6991 6313 3794 1836 345 382 808
## [91] 1388 2713 3800 3091 2985 3790 674 81 80 108 229 399 1132 2432 3574
## [106] 2935 1537 529 485 662 1000 1590 2657 3396boxplot(lynx)
자료를 분석하기에 앞서 lynx 자료는 1821년부터 1934년까지 캐나다에서 덫에 걸리는 스라소니의 수를 나타난 자료이며 시계열 자료로 작성되어 있다.
plot(lynx)
상자그림을 통한 자료요약을 살펴보면 4개의 특이 값이 존재하고 있는 것을 확인할 수 있었다. 또한, 최솟값과 최댓값 차이가 크기 때문에 몇몇 큰 값이 보일 것이라 예상된다. 1번에서 언급한 바와 같이 시간의 흐름에 따른 특성은 확인할 수 없기 때문에 자료의 Plot을 통해 다시 확인하였다.
plot(lynx,lty=2,ylim=c(0,7000),xaxt="n",xlab="year");par(new=T)
plot(sleek(lynx),xaxt="n",col="red",type="l",ylim=c(0,7000),ann=F);par(new=T)
h4253=a3rssh2(lynx)
plot(h4253,xaxt="n",xlab="year",type="l",ylim=c(0,7000),col="blue",ann=F)
axis(1,at=1:114,labels=1821:1934)
legend(x = 0.5, y = 7000, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
위 Plot은 2가지 평활법을 적용하여 나타낸 결과이다. 먼저 원 자료에 비해 최댓값과 최솟값의 폭이 많이 줄어들었음을 확인할 수 있다. 작은 값에서는 크게 변화가 없지만 값이 큰 자료들이 많이 줄어들었다.
추세를 확인함에 있어 훨씬 더 좋은 모양으로 변하였지만 전 문제에서 언급한 바와 같이 원 자료와 많이 차이나는 부분들은 특징을 날린 것이 아닌지 유의해야한다. 또한, 원자료에서 1911년과 1916년 사이의 값이 쌍봉으로 나타났는데 평활 후 하나의 봉우리만 나타나게 되었다. 2가지 평활법의 차이를 살펴보면 4253H방법에서 큰 값들의 감소폭이 크게 나타나고 있다. 덫에 걸리는 스라소니의 수는 스라소니의 전체 수와 비례한다고 볼 수 있다. 이를 토대로 위의 데이터를 분석한다면 약 10년을 주기로 스라소니의 수가 크게 증가하였다가 감소하고 있다.
더 정확한 분석이 필요한 부분이지만 약 1904년 폭발적으로 수가 증가한 후 10년 주기로 개체 수가 지난 때보다 증가폭이 덜한 것을 확인할 수 있다. 이는 환경 파괴로 인해 스라소니의 개체 수가 줄어든 것이 아닌가? 라는 생각이 들었다. 지난 학기 ‘시계열 분석’에서 배운 내용을 토대로 위의 해석을 덧붙이자면 위 자료는 주기를 가진 자료이기 때문에 이를 조정한 뒤 평활을 한다면 더욱 보기 좋은 결과 값이 나올 것이다. 위 결과로는 추세성 또한 확인하기 어렵기 때문에 위 자료를 가공한다면 더 많은 데이터의 특징을 얻을 수 있을 것으로 짐작된다.
3. 7장 숙제에 포함된 jumping 자료를 분석하여라.
setwd("C://Users/R/Desktop/4학년 1학기/탐자분/datafile")
jump=read.table("jumping.DAT")
colnames(jump)=c("year","highjump",
"polevault", "longjump", "tripplejump")
attach(jump)Jumping파일을 불러온 뒤 각 열에 변수지정이 되어있지 않아 우선적으로 각 열에 해당 변수(year, high jump, pole vault, long jump, tripple jump)를 지정해주었다. 또한 1916년, 1940년, 1944년의 자료가 누락되어 있다.
1916년에는 1차 세계대전이 발발되어 연기가 되었다가 1920년 올림픽이 개최됨으로서 개최가 최소 되었으며 1940년, 1944년의 경우 2차 세계대전으로 인해 취소되었다고 한다. 이러한 이유로 3개년의 값은 건너뛰고 Plot을 그렸다.
library(sleekts)
plot(highjump,xaxt="n",xlab="year",type="l",ylim=c(1.8,2.4),lty=2,main="Record of high jump");par(new=T)
plot(sleek(highjump),xaxt="n",col="red",type="l",ylim=c(1.8,2.4),ann=F);par(new=T)
h4253=a3rssh2(jump$highjump)
plot(h4253,xaxt="n",xlab="year",type="l",ylim=c(1.8,2.4),col="blue",ann=F)
axis(1,at=1:20,labels=jump$year)
legend(x = 0.5, y = 2.4, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
plot(polevault,xaxt="n",xlab="year",type="l",lty=2,ylim=c(3.3,6.0),main="Record of pole valut");par(new=T)
plot(sleek(polevault),xaxt="n",col="red",type="l",ylim=c(3.3,6.0),ann=F);par(new=T)
p4253=a3rssh2(jump$polevault)
plot(p4253,xaxt="n",xlab="year",type="l",ylim=c(3.3,6.0),col="blue",ann=F)
legend(x = 0.5, y = 6, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
axis(1,at=1:20,labels=jump$year)
plot(longjump,xaxt="n",xlab="year",type="l",lty=2,ylim=c(7,9.5),main="Record of long jump");par(new=T)
plot(sleek(longjump),xaxt="n",col="red",type="l",ylim=c(7,9.5),ann=F);par(new=T)
l4253=a3rssh2(jump$longjump)
plot(l4253,xaxt="n",xlab="year",type="l",ylim=c(7,9.5),col="blue",ann=F)
legend(x = 0.5, y = 9.5, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
axis(1,at=1:20,labels=jump$year)
plot(tripplejump,xaxt="n",lty=2,xlab="year",type="l",ylim=c(13,18),
main="Record of triple jump");par(new=T)
plot(sleek(tripplejump),xaxt="n",col="red",type="l",ylim=c(13,18),ann=F);par(new=T)
t4253=a3rssh2(jump$tripplejump)
plot(t4253,xaxt="n",xlab="year",type="l",ylim=c(13,18),col="blue",ann=F)
legend(x = 0.5, y = 18, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
axis(1,at=1:20,labels=jump$year)
프린트에 주어진 것처럼 2가지 방법(4253H, twice/3RSSH, twice)을 이용하여 평활한 결과는 위와 같다. 3RSSH방식이 원 자료를 더 잘 나타낸다고 생각되어 이를 바탕으로 자료를 해석하고자 하였다. 모든 종목에서 시간이 흐름에 따라 기록이 좋아지고 있다고 볼 수 있다. 시간에 따른 기술발전과 스포츠 분석 시스템이 발달되면서 모든 종목에서의 기록 향상은 당연한 결과일 것이다. 개별 종목을 살펴보면 Long jump의 원 자료 값이 변동이 크다는 것을 알 수 있으며 pole valut의 원 자료 값이 평활한 직선에 가장 잘 맞다고 판단된다.
마지막으로 long jump의 1968년 기록이 굉장히 큰 값을 가지고 있는 이유가 궁금하여 찾아본 결과 당시 미국의 멀리뛰기 선수가 세운 기록으로 당시 기록보다 21.5츠 앞서 경이로운 기록이었다고 한다. 이 기록은 23년이 지나서야 깨지게 되었다. 이러한 값들은 평활에 있어 하나의 특징으로 살려야 할지 아니면 이 시점에만 발생한 값이라 여겨야 할지 좀더 고민해볼 필요성이 있다. 위 Plot을 통해 각 종목을 개별적으로 판단하는 것은 전혀 문제가 없었지만 문제에서 어느 종목의 기록이 점점 기록 갱신하기 어려워지는지 묻고 있다. 여기서 의문점이 생기게 되었다.
fivenum(highjump)## [1] 1.800 1.940 2.035 2.235 2.380fivenum(tripplejump)## [1] 14.330 15.065 16.110 17.320 17.610high jump의 경우 1.8과 2.38사이에 값들이 존재하고 있는 반면 triple jump의 경우 14.33과 17.61사이에 값들이 존재하고 있다. 그렇다면 만약 절대적으로 1만큼의 변화가 있었다면 high jump의 경우에는 엄청난 변화이지만 triple jump의 경우 high jump만큼의 큰 변화는 아닐 것이다. 이러한 이유로 4개 종목의 표준화를 진행하였다. 표준화는 평균을 기준으로 얼마나 떨어져 있는지를 나타내는 값으로 보통 단위가 다른 2개 이상의 데이터가 주어졌을 때 적용한다. 표준화를 통해 같은 기준으로 볼 수 있게 된다. 또한, 표준화를 진행하면 데이터를 평평하게 하여 데이터의 진폭을 줄일 수 있다는 특징이 있다. 다만 유의해야 할 점은 표준화 후 값을 통해 원 자료의 값처럼 직접적인 수치 해석은 어렵다.
보통의 표준화 식은 (개별 값-평균값)/표준편차 이다. 물론 jumping 데이터에 큰 특이 값은 존재하지 않지만 이 표준화 방법은 데이터가 정규분포가 아닌 경우 적절하지 않을 수 있다. 이것의 대안으로 Min-Max Scaler와 Robust Scaler방식을 찾을 수 있었는데 Robust scale방법이 IQR을 이용하기 때문에 특이 값에 가장 민감하지 않고 그동안 배운 것들을 바탕으로 한다면 가장 적절한 표준화 방법이라 생각되어 Robust scale방법을 이용하였다. 즉, 각 종목들의 중앙값이 0 IQR이 1이 되도록 변환한 것이다.
(개별 값-중간값)/(3분위수-1분위수)
#표준화작업
stand= function(x){
std=(x-median(x))/(fivenum(x)[4]-fivenum(x)[2])
return(std)
}
highjump=stand(jump$`highjump`)
polevault=stand(jump$polevault)
longjump=stand(jump$longjump)
tripplejump=stand(jump$`tripplejump`)library(sleekts)
plot(highjump,xaxt="n",xlab="year",type="l",ylim=c(-2,2),lty=2,main="Record of high jump");par(new=T)
plot(sleek(highjump),xaxt="n",col="red",type="l",ylim=c(-2,2),ann=F);par(new=T)
h4253=a3rssh2(highjump)
plot(h4253,xaxt="n",xlab="year",type="l",ylim=c(-2,2),col="blue",ann=F)
axis(1,at=1:20,labels=jump$year)
legend(x = 0.5, y = 2, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
plot(polevault,xaxt="n",xlab="year",type="l",lty=2,ylim=c(-2,2),main="Record of pole valut");par(new=T)
plot(sleek(polevault),xaxt="n",col="red",type="l",ylim=c(-2,2),ann=F);par(new=T)
p4253=a3rssh2(polevault)
plot(p4253,xaxt="n",xlab="year",type="l",ylim=c(-2,2),col="blue",ann=F)
legend(x = 0.5, y = 2, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
axis(1,at=1:20,labels=jump$year)
plot(longjump,xaxt="n",xlab="year",type="l",lty=2,ylim=c(-2,2),main="Record of long jump");par(new=T)
plot(sleek(longjump),xaxt="n",col="red",type="l",ylim=c(-2,2),ann=F);par(new=T)
l4253=a3rssh2(longjump)
plot(l4253,xaxt="n",xlab="year",type="l",ylim=c(-2,2),col="blue",ann=F)
legend(x = 0.5, y = 2, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
axis(1,at=1:20,labels=jump$year)
plot(tripplejump,xaxt="n",lty=2,xlab="year",type="l",ylim=c(-2,2),main="Record of triple jump");par(new=T)
plot(sleek(tripplejump),xaxt="n",col="red",type="l",ylim=c(-2,2),ann=F);par(new=T)
t4253=a3rssh2(tripplejump)
plot(t4253,xaxt="n",xlab="year",type="l",ylim=c(-2,2),col="blue",ann=F)
legend(x = 0.5, y = 2, c("Raw", "4253H,twice", "3RSSH,twice"),
col = c("black","red","blue"),
lty=c(2,1,1),lwd=2)
axis(1,at=1:20,labels=jump$year)
y축을 동일하게 설정하여 쉽게 비교할 수 있게 하였다. 이번에는 각 평활법에 따라 크게 차이가 없지만 4253H, twice방식이 원 자료 값에 좀 더 적절하게 평활되었고 위 문제에서 언급한 바와 같이 부드러운 형태로 나타나 이 방식으로 자료를 분석해보고자 하였다. 하지만 각각 Plot을 그렸기 때문에 추이를 비교하는 것이 어려워 하나의 Plot에 표시해보았다. 또한, 위 Plot에서는 y축을 -2~2로 설정하였는데 너무 완만한 듯하여 아래에서는 -1~1로 변경하였다.
plot(lwd=2,sleek(highjump),xaxt="n",xlab="year",type="l",ylim=c(-1,1),col="red",ylab="",main="Record of Sports(4253H twice)");par(new=T)
plot(lwd=2,sleek(polevault),xaxt="n",xlab="year",type="l",ylim=c(-1,1),col="green",ann=F);par(new=T)
plot(lwd=2,sleek(longjump),xaxt="n",xlab="year",type="l",ylim=c(-1,1),col="blue",ann=F);par(new=T)
plot(lwd=2,sleek(tripplejump),xaxt="n",xlab="year",type="l",ylim=c(-1,1),col="black",ann=F)
axis(1,at=1:20,labels=jump$year)
legend(x = 0.5, y = 1, c("high jump", "pole vault", "long jump", "triple jump"),
col = c("red","green","blue","black"),
lty=c(1,1,1,1),lwd=2)
위 자료의 기울기를 통해 시간에 따라 얼마나 크게 기록이 변화하였는지 확인할 수 있다. 동 시간대에서 기울기가 큰 종목이 완만한 종목보다 기록향상이 더 크다고 생각하면 된다. 큰 값을 가질수록 각 종목의 중간 값에 비해 더 큰 폭으로 증가한 것이다. 사실 대부분 종목의 기울기가 비슷한 추이를 가지고 있어 두드러지는 특징이 보이지는 않지만 구간을 나누어 확인하니 몇몇 특징을 확인할 수 있었다.
가장 먼저 1968년 이후 triple jump의 기울기가 다른 종목에 비해 매우 완만해지고 있다. 이는 triple jump의 기록 갱신이 점점 어려워지고 있다고 볼 수 있다. 다른 종목은 위에서 언급한 바와 같이 분석기술과 장비의 발전으로 기록향상이 되고 있지만 triple jump의 기록향상 추이는 약간 다른 듯하다. 다음으로 1956년 이후 pole vault의 기울기가 이전보다 가파르게 나타나고 있다. 이 기간동안 해당 종목에서 어떠한 향상요인이 생겨났을 것으로 예상된다.