1장 데이터 처리

KNOU/4 고급R 2017. 8. 30. 13:36
반응형



1. 자료 가져오기 

> setwd(choose.dir()) #대화형으로 폴더를 지정할 수 있음

> getwd()

[1] "C:/Users/amore/Documents/FOR ME/Data Scientist/KNOU/2017년 2학기/고급R활용"

> setwd("C:/Users/amore/Documents/FOR ME/Data Scientist/KNOU/2017년 2학기/고급R활용")

> getwd()

[1] "C:/Users/amore/Documents/FOR ME/Data Scientist/KNOU/2017년 2학기/고급R활용"

> savehistory(file='mywork.Rhistory')

> history()

> history(max.show=10)




2. 다양한 형태의 자료 가져오기

> insurance.data2 = read.table('insurance2.txt', header=T)

> insurance.data2[c(15:16),]

   id sex job religion edu amount salary

15 15   m  -9        3   4      6    110

16 16   f   1        3  -9      7     88


#-9를 NA로 설정

> insurance.data2 = read.table('insurance2.txt', header=T, na.strings = '-9')

> insurance.data2[c(15:16),]

   id sex job religion edu amount salary

15 15   m  NA        3   4      6    110

16 16   f   1        3  NA      7     88

> csv.data = read.csv('csv.txt',header=T) #구분자가 콤마

> tab.data = read.table('tab.txt', header=T,sep='\t') #구분자가 탭

> write.table(tab.data, file='C:/Rwork/tab.txt') #Rwork 폴더에 저장


#고정형식 텍스트 파일 읽기 

> fwf.data = read.fwf('insurance3.txt',

+                     widths=c(2,2,3,3,3,6,6),

+                     col.names=c('id','sex','job','religion','edu','amount','salary'))

> fwf.data[fwf.data$job==-9, 'job']=NA

> fwf.data[fwf.data$edu==-9, 'edu']=NA

> fwf.data[fwf.data$salary==-9, 'salary']=NA

> head(fwf.data)

  id sex job religion edu amount salary

1  1   m   1        1   3    7.0    110

2  2   m   2        1   4   12.0    135

3  3   f   2        3   5    8.5    127

4  4   f   3        3   5    5.0    150

5  5   m   1        3   3    4.5    113

6  6   m   2        1   2    3.5     95


> fwf2.data = read.fwf('insurance3.txt',

+                     widths=c(2,-2,-3,3,3,6,6),

+                     col.names=c('id','religion','edu','amount','salary'))

> head(fwf2.data)

  id religion edu amount salary

1  1        1   3    7.0    110

2  2        1   4   12.0    135

3  3        3   5    8.5    127

4  4        3   5    5.0    150

5  5        3   3    4.5    113

6  6        1   2    3.5     95





3. 통계 패키지 자료 가져오기

> library(foreign) #spss 파일

> ex1 = read.spss('ex1-1.sav',

+                 to.data.frame=T, use.value.labels = T)

Warning message:

In read.spss("ex1-1.sav", to.data.frame = T, use.value.labels = T) :

  ex1-1.sav: Unrecognized record type 7, subtype 18 encountered in system file

> ex1

  shock response count

1   yes      yes    25

2   yes       no    19

3    no      yes    31

4    no       no   141

> dim(ex1)

[1] 4 3

> mouse.data = ex1[rep(1:nrow(ex1),ex1$count),] #가중치

> head(mouse.data)

    shock response count

1     yes      yes    25

1.1   yes      yes    25

1.2   yes      yes    25

1.3   yes      yes    25

1.4   yes      yes    25

1.5   yes      yes    25

> dim(mouse.data)

[1] 216   3

> mouse.table=table(mouse.data$shock,mouse.data$response)

> mouse.table

     

      yes  no

  yes  25  19

  no   31 141

> summary(mouse.table)

Number of cases in table: 216 

Number of factors: 2 

Test for independence of all factors:

Chisq = 27.458, df = 1, p-value = 1.605e-07


 - 카이제곱 = 27.458 , p - value = 1.6 X 10마이너스 7승. 자극과 반응간의 유의한 관계가 있음.




4. RData 저장 및 가져오기

> save(ex1, file='c:/Rwork/ex1.Rdata')

> rm(ex1)

> load('C:/Rwork/ex1.Rdata')

> ex1

  shock response count

1   yes      yes    25

2   yes       no    19

3    no      yes    31

4    no       no   141

> load(file=file.choose()) #대화








참고문헌: 

[1] 고급 R 활용 (심송용, 이윤동, 이은경, 김성수 공저, KNOU PRESS)


반응형
Posted by 마르띤
,

세계 최고의 MBA는 무엇을 가르치는가

서평 2017. 8. 28. 13:49
반응형

URL: http://www.yes24.com/24/Goods/13315629?Acode=101


한 때 잠깐이나마 MBA를 꿈꾸던 시절이 있었다. 아니라고는 하지만, 아직도 많은 부분이 학벌로 평가되는 시기이고 MBA가 어느정도 보완 효과가 있다고 생각도 하였고, 좀 더 화려해 보이는 일들에 대한 욕심도 있고 하여 MBA를 잠깐 이나마 생각을 해 보았다. 학업에 대한 열망이라고 표현은 하지만 결국 MBA를 가는 사람은 결과적으로 더 나은 곳에서 더 나은 직장에 대한 열망으로 시작을 하였을테니 말이다.


일본 작가가 쓴 책이다. 늘 그렇지만 일본 작가는 참 좋은 기획을 한다. 일본 작가인 만큼 일본의 직장인들이 미국/유럽 등 유명 MBA를 수강하게 된 계기, 배운 것들, 현재 일하고 있는 직장에 대해 짧게 많은 케이스를 소개하고 있다. 그러면서 각 MBA 학교 특징, 커리큘럼, 학비 등의 정보를 짧게나마 소개도 하고 있다. 2014년에 출간된 책이라 2013년을 기준으로 쓰인 책이긴 하지만, 중간 중간 MBA 학교에서 제시된 토론 주제는 지금 생각해봐도 의미있는 내용들이 있어 참고해볼만 하다.


  • 하버드 비지니스 스쿨 : 리더십과 기업 윤리

 - 2003년 이라크 전쟁, 폭탄이 설치된 가옥이 있고 그 안에 이라크인 3명이 남아 있을 때, 부대장은 부하의 생명을 위험을 무릅쓰고 집 안으로 들어가야 하는가?


  • 스탠퍼드: 성장 기업의 매니지먼트

 - A와 B는 벤처기업의 공동차업자이다. 벤처캐피털이 주자 조건으로 이사를 A 한명이 맡아달라는 요청이 왔다. 내가 A라면 B에게 어떤 식으로 설명할 것이가?


  • 펜실베니아 와튼 스쿨: 영어 스피치로 호응을 얻는 기술 '스피치 라이팅'

"If the action changes, the reputation changes. The confidence changes if the reputation changes. The life changes if the confidence changes, This is my lesson at Wharton.


  • 노스웨스턴 대학 켈로그 스쿨 오브 매니지먼트: 변명하지 않기 위한 훈련

 - 항공업계에서 성공한 사례를 레스토랑 체인에 적용해보면 어떻게 될까? 스타벅스라는 말을 들으면 떠오르는 형용사는? 그 형용상를 지웠을 때 어떤 비지니스가 연상되는가?


  • 컬럼비아 비지니스 스쿨: 통합 마케팅 전략

 - 델타 항공은 제트블루 항공의 성공을 모델 삼아 저가항공 사업에 뛰어들어야 하는가?


  • 시카고대학 부스 비지니스 스쿨: 창업자 육성을 위한 컴피티션


  • 미시건 대학 로스 비지니스 스쿨: 글로벌 인재의 조건

 - 1994년 설립된 반야트리는 아시아에서 급성장을 이루었다. 어떻게 CSR를 다하면서 성장해 왔는가? 다른 호텔에 비해서 뭔가 특이한 점이 있는가?


  • 다트머스 대학 터크 비지니스 스쿨: 리버스 이노베이션

 - 중국 현지팀의 손으로 시장에 맞는 초음파 진단장치 개발 사례

 - 자기가 제어할 수 없는 외부 요인이 아닌 내부 요인에만 컨트롤 하기 - 인종 / 출신 국가 등


  • 듀크대학 퓨쿠아 비지니스 스쿨 : 숫자의 이면을 읽기

 - 실적 부진에 빠진 제조회사. 원료 조달, 마케팅, 제조 부분 책임자의 실적 평가


  • 런던 비지니스 스쿨: 역사와 전통을 팔기

 - 라이센싱 사업으로 여로 종류의 제품이 대량으로 판매되어 브랜드의 희석하가 시작되는 버버리가 취할 전략은? 


  • 인시아드: 심리적 요인을 찾아내서 자신을 바꾸자

 - 자신이 이제껏 바꾸고 싶어했지만 바뀌지 않는 습성이나 버릇은? 왜 바뀌지 않을까?

반응형

'서평' 카테고리의 다른 글

기업 위기 시스템으로 이겨라  (0) 2017.08.31
유라시아견문-제1권: 몽골 로드에서 할랄 스트리트까지  (0) 2017.08.07
구글신은 모든 것을 알고 있다 - 복잡계 네트워크, 정하웅 교수님  (0) 2017.07.23
불황을 넘어서  (0) 2014.07.13
월급전쟁  (0) 2014.07.13
Posted by 마르띤
,

유라시아견문-제1권: 몽골 로드에서 할랄 스트리트까지

서평 2017. 8. 7. 10:10
반응형

http://www.yes24.com/24/goods/30966027?scode=032&OzSrank=1



총 3권의 책인데, 내가 읽은 것은 그 중 제1권 <유라시아 견문> 중, 몽골의 올란바토르부터 인도네시아의 자카르타까지, 첫 1년의 기록이다.


유라시아의 사전적 의미는 뭘까? 라는 궁금증이 일어났다. Eurasia 유럽과 아시아를 하나로 묶어 부르는 이름이라 하고, 총 면적은 5492만 ㎢, 세계 전육지의 40%라는 광활한 영역을 가르킨다. 책 서문을 보면 서진에 서진을 거듭하여 터키까지 이른다는 표현이 있다. 그리고는 유럽과 아시아 사이의 공간적 장벽을 허물고, 전통과 근대 사이의 시간적 단층을 돌파해내고 싶다는 문장, 유라시아사를 쓰고 싶다는 작가의 말이 이 책의 존재를 정의내리고 있다. 


책이 두꺼워서 몇장 읽고 덮을 뻔 하였지만, 프롤로그에 저자만의 유라시아를 정의한 내용이 아니었다면 흥미를 쉽게 잃었을 것 같다.


새 천년 초원길과 바닷길의 복원은 판갈이의 출발이다. 백 년간 끊어지고 막혔던 동서의 혈로를 다시 뚫어 물류와 문류를 재가동하는 유라시아의 재활운동이다. 국격Border이 통로Gateway가 되고, 지리는 재발견되고 지도는 다시 그려진다. 작금의 모순과 균열을 미-중간의 패권 경쟁으로 오독해서는 안된다. 서구를 배타하지도, 근대롤 폄하하지도 않는다. 사물을 제자리에 돌려놓을 분이다. 유럽을 유라시아의 서단으로 지방화하는 것이다. 그리하여 이름을 바르게 불러주는 것이다. 유럽과 아시아, 근대와 전근대의 분단체제를 허물고 유라시아적 맥락으로 동서고금을 재인식하는 것이다. 유러브이 자만도, 아시아의 불만도 해소하는 대동세계의 방편이다. 한반도 동남단, 경주의 석굴암은 서역과 페르시아로 이어졌던 누천년 유라시아 연결망을 무묵히 증언하고 있다.


저자의 여정은 태국에서 '조선 지식인 유길준이 쓴 서유견문'을 읽는 것으로 시작한다. 태국 치앙라이에서 전혀 관계가 없을 것 같았던 흥화興華라는 간판을 발견하면서 동아시아의 냉전사로 이야기를 풀어나간다. 인민해방군이 남진하여 윈남성 성도인 쿤밍을 장악하였다. 1949년 10월 중화인민공화국이 설립되었지만, 국민당 잔군은 여전히 서남부 내륙 윈남성에는 그 세력들이 존재하였던 것이다. 항복을 거부한 일부는 남하하여 버마(미얀마)로 숨었고, 스스로를 반공구국군이라 개칭하였다.  미국은 반공구국군을 내륙으로 침투시켰다. 1953년 한반도 휴접협정에도 불구하고 전선의 총성은 멈추지 않자, 반공구국군 버마와 태국 산악지역에 살고 있는 소수민족도 끌여들였다. 하지만 유엔총회에서 버마의 영토주권을 침해하고 있는 중화민국의 잔군을 철수시키라는 요구도 있었고, 탈식민의 물결을 거스를 수 없던 장제스 역시 철군을 단행하였지만, 버마에서 쫓겨난 이들은 낯선 대만이 아닌 태국의 치앙라이의 마에살롱 마을을 선택하였다. 생계를 위해 태국과 버마 사이의 아편 무역 중계를 하였다. 해발 1800미터의 고산지대라 양귀비를 재배하기 적격이었다. 지금 이곳은 차밭으로 유명하다. 우롱차 차밭을 기반으로 지금은 커피, 약초 재배도 시작하였고 어느덧 호텔, 리조트, 레스토랑도 길게 늘어서고 있다. 시간이 지나면서 태국의 화교/화인 신문도 공급되고 있다. 또다른 중화망의 한 연결 고리이다.


태국에서 전혀 상관없어 보이던 중화민국의 연결고리를 찾아내어 네트워크 중화제곡, inter-localism이라는 키워드를 뽑아낸다. 태국의 영자신문 <방콕 포스트>에서 중국의 해양 실크로드 구상과 아세안의 미래를 토론하는 학술회의가 예고되어 있다며, 그 다음은 아세안의 이야기로 넘어간다. 모든 이야기가 이렇게 이어지니 단편적인 유라시아사가 아니라 정말 서문에서 밝힌 그대로 "근대와 전근대의 분단체제를 허물고 유라시아적 맥락으로 동서고금을 재인식하는 것"이라는 작가의 의도가 이해되었다. 


그 이후에도 

- 중국의 춘절이 크리스마스 같은 글로벌 축제가 될 수 있을까?

- 아시아 인프라 투자은행의 탄생

- 인도태평양의 균형자 인도네시아

- 유라시아의 축 몽골

- 아시아의 하늘을 잇는 저가항공사

- 실학자들의 나라 싱가포르

- 이슬람 경제의 메카 말레이시아

- 할랄산업

- 필리핀의 슬픈 민주주이, 근대화와 속국 민주주의

- 대동

- 시안의 미래는 장안이다, 서역의 출발점

- 서유기, 구도와 득도의 길

- 신장위구르의 중국화와 세계화

- 일대일로의 사상: 후안강과의 대화

- 동서고금의 교차로 : 카슈가르 

- 윈남에서 이슬람적 중국을 만나다


등 재미있는 주제로 한가득하다. 이중 몽골의 공산당 이야기, 새롭게 해석한 중국의 시안이야기, 서유기, 신장 위구르에 관한 이야기는 내겐 큰 흥미와 관심으로 이어져서, 관련 서적도 읽어 보고 싶게 만들고 있다.  반면 일대일오의 사상, 후안강과의 대화 내용은 다소 의문이 생겨 작가님에게 메일을 드렸고 작가님도 지금 다른 집필일로 답변을 당장 받진 못하였지만 앞으로의 작가님과 대화 그리고 앞으로 출판될 2권, 3권이 더 기대되기도 한다.


반응형

'서평' 카테고리의 다른 글

기업 위기 시스템으로 이겨라  (0) 2017.08.31
세계 최고의 MBA는 무엇을 가르치는가  (0) 2017.08.28
구글신은 모든 것을 알고 있다 - 복잡계 네트워크, 정하웅 교수님  (0) 2017.07.23
불황을 넘어서  (0) 2014.07.13
월급전쟁  (0) 2014.07.13
Posted by 마르띤
,

구글신은 모든 것을 알고 있다 - 복잡계 네트워크, 정하웅 교수님

서평 2017. 7. 23. 21:15
반응형

카이스트 명강 시리즈 중 하나. 카이스트 정하웅 교수님의 책. 예전에 한 번 읽긴 하였는데 최근 데이터 분석 공부를 하면서 신경망 그리고 텍스트마이닝의 네트워크 분석의 매력에 크게 빠져 있는 상태라 그런지 이 책을 다시 읽어보니 너무 재미있었다. 


우연히 티비를 틀었을 때 EBS의 <과학다큐 비욘드>에서 네트워크 과학을 소개하는 프로를 본적이 있었다. 물리학자 윤혜진 교수님이 진행이었는데, 이 책에는 카이스트 박사과정 윤혜진 학생으로 소개되어 다소 색다르기도 하였다. EBS 과학다큐 비욘드 내용을 잠시 소개하자면 아래와 같다.


“우리를 바라보는 새로운 시선, 네트워크 사이언스”

 

세상 만물을 물리로 풀어내는 과학자들이 나타났다. 도시가 어떤 모습으로 성장을 하는지, 부자들의 특징이 무엇인지를 과학으로 풀어내는 물리학자들. 이들은 어떤 시선으로 세상에 관심을 가지고 질문을 던지는 걸까? 이들이 바라보는 세상은 어떤 구조로 되어 있을까?

 

물리학자 윤혜진은 1870년대 이후 새로운 아이디어로 나온 특허는 드물다는 사실을 알아냈다. 기존의 아이디어가 조합되어 나온 ‘조합특허’가 대부분이라는 것이다. 세상을 바꾸는 혁신은 완전히 새로운 아이디어로 만들어지는 게 아니라는 거다.

 

세계 최대 복잡계 연구소인 산타페 연구소의 제프리 웨스트 교수는 도시가 ‘15%의 법칙’으로 움직인다는 것을 밝혀냈다. 도시가 두 배 성장할 때 어떤 것은 15% 더 성장하고, 어떤 것은 15% 덜 성장한다는 것이다. 세계의 어떤 도시라도 예외는 없다.

 

하버드대 ‘인터넷과 사회를 위한 버크만 클라인 센터’ 요하이 벵클러 교수는 지난해 미국 대선 보도를 받아들이는 이용자들의 이용패턴을 네트워크 지도로 그렸다. 미국 대선을 강타했던 ‘가짜뉴스’ 논란이 가짜뉴스 그 자체가 아니라 이용자들에게 달린 문제라는 것을 분석해냈다.

 

이들의 공통점은 모두 네트워크를 통해 세상을 바라본다는 것. 갈수록 복잡해지는 세상에서 네트워크 사이언스는 나무가 아닌 숲을 볼 수 있게 하는 길잡이가 되어 준다. 초연결사회, 연결을 거부할 수 없다면 내가 어디에, 어떻게 연결되어 있는지를 아는 것만으로도 나의 역할이 무엇인지, 내가 누구인지를 알 수 있다.

 

이외에도 네트워크를 통해 알아보는 부자의 특성, 가짜뉴스의 특징과 판별법 등 일상에 영향을 끼치는 과학에 대한 새로운 정보와 지식은 앞으로 살아갈 세상에서 우리가 나아가야할 방향을 제시해줄 것이다.

 


지금은 유튜브에서 너무 손쉽게 찾을 수 있어, 유튜브 링크도 찾아보았다.



복잡계 네트워크를 가장 쉽게 표현하는 것은 바로 고속도로와 항공망. 고속도로의 경우 한 도시마다 일반적으로 3-4개의 도로가 연결된 반면, 항공 네트워크의 경우 뉴욕, 시카고, LA의 경우 무수히 많은 항공망이 연결되어 있다. 


이를 다시 그래프로 표현하면 멱함수 분포 곡선으로 표현할 수 있는데, x축의 우측으로 가면 갈 수록 연결 횟수가 많은 소수의 점을 나타낸다. 이를 이용한 수많은 자연/사회 현상을 재미있게 소개한다. 


 - 국회의원 네트워크: 구글 검색을 통해 얻은 검색 결과를 바탕으로 허브 조사


의 그래프는 시각화도 잘 되어 있어 한참을 쳐다 보았다.


문득 드는 아이디어로는 화장품도 요즘은 멀티 브랜드, 멀티 채널인데 소비자가 2개 이상의 브랜드를 함께 쓴다고 가정하였을 때, 허브 브랜드가 존재할 것 같았다. 더 나아가 허브 제품도 있을 것이고. 허브 제품 구매자와 비구매자 간 비교 분석을 하거나, 허브 제품 비구매자에게 허브 제품 구매 유도하는 CRM을 하였을 때 어떤 반응을 보인지 분석하는 것도 의미 있어 보인다. 


이 책덕분에 복잡계 네트워크에 크게 관심을 가지게 되었다. 내가 만약 대학원을 가 석사 공부를 하게 된다면 이것을 전공으로 하고 싶을 정도로 말이다.


 - 책 링크: 구글 신은 모든 것을 알고 있다

 - 정하웅 교수님 홈페이지 링크 : http://stat.kaist.ac.kr/

 - 윤혜진 교수님(노스웨스턴대 경영학교 교수님) 홈페이지: http://hyoun.me/

 

참고문헌

- 알버트 라즐로 바라바시, "링크", 2002년

- 마크 뷰캐넌, "넥서스" , 2003년

- 던컨 와츠, "스몰월드", 2004년

- 마크 뷰캐넌, "세상은 생각보다 단순하다", 2004년

- 강병남, "복잡계 네트워크 과학: 21세기의 정보과학", 2010년

- 알버트 라즐로 바라바시, "버스트: 인간의 행동 속에 숨겨진 법칙", 2010년


일반독자의 경우 링크나 넥서스, 과학 분야는 링크, 인문 사회 과학분야는 스몰월드, 대학 수준의 교재는 복잡계 네트워크 과학



<몇가지 용어들 정리>

1. 복잡계 네트워크: 점(vertex, node)과 연결선(edge, link)들로 이루어진 집합을 의미한다. 특히 연결선들이 들어오고 나오는 방향이 있는 경우 방향성 네트워크라고 하고, 각 연결선들에 가중치가 부여된 경우 가중치 네트워크라고 한다. '네트워크 이론'은 응용수학과 물리학 분야에서 다루는 이론으로, 수학의 그래프 이론에서 비롯하였다. 현재 전산학, 생물학, 경제학, 사회학 분야에 널리 적용된다 (출처: 위키피디아)

2. CSSPL : Complex Systems and Statistical Physics Lab 복잡계와 통계 물리학 연구실





반응형

'서평' 카테고리의 다른 글

세계 최고의 MBA는 무엇을 가르치는가  (0) 2017.08.28
유라시아견문-제1권: 몽골 로드에서 할랄 스트리트까지  (0) 2017.08.07
불황을 넘어서  (0) 2014.07.13
월급전쟁  (0) 2014.07.13
탐스 스토리  (0) 2014.07.13
Posted by 마르띤
,

Lesson 3: Explore One Variable

Python, R 분석과 프로그래밍/Data Analysis with R 2017. 7. 12. 10:23
반응형



#데이터 불러오기

getwd()[1] "C:/Users/amore/Documents/FOR ME/Data Scientist/Udacity/Data Analysis with R/eda-course-materials/lesson3"
list.files()
[1] "lesson3.Rmd"         "lesson3_student.rmd" "pseudo_facebook.tsv"
pf <- read.csv('pseudo_facebook.tsv', sep='\t')
names(pf)
 [1] "userid"                "age"                   "dob_day"              
 [4] "dob_year"              "dob_month"             "gender"               
 [7] "tenure"                "friend_count"          "friendships_initiated"
[10] "likes"                 "likes_received"        "mobile_likes"         
[13] "mobile_likes_received" "www_likes"             "www_likes_received"





#Histogram of User's Birthdays

library(ggplot2)qplot(x=dob_day, data = pf)

Here's some things that I noticed. On the first day of the month I see this huge bin of almost 8,000 people. This seems really unusual since I would expect most people to have the same number of birthday's across every day of the month.  



library(ggplot2)qplot(x=dob_day, data = pf) + 
  scale_x_continuous(breaks=1:31)





#Faceting

library(ggplot2)qplot(x=dob_day, data = pf) + 
  scale_x_continuous(breaks=1:31)+
  facet_wrap(~dob_month, ncol=4)

Now, you may have noticed some peaks in May or perhaps in October, but I think what's really interesting is this huge spike on January first. There's almost 4,000 users in this bin. Now, this could be because of the default settings that Facebook uses or perhaps users are choosing the first choice in the drop down menus. Another idea is that some users may want to protect their privacy and so they just go with January first by default. Whatever the case may be, I think it's important that we make our considerations in the context of our data. We want to look out for these types of anomalies. 



#Friend Count

ggplot(aes(x = friend_count),data = pf) +
   geom_histogram()





#Limiting the Axes

qplot(x = friend_count, data = pf, xlim = c(0,1000))



## alternative solution
qplot(x = friend_count, data = pf) + 
  scale_x_continuous(limits = c(0,1000))




#adjust the bin width

 qplot(x = friend_count, data = pf) + 
   scale_x_continuous(limits = c(0,1000), breaks = seq(0,1000,50))

## alternative solution
qplot(x = friend_count, data = pf) + 
  scale_x_continuous(limits = c(0,1000))



##splits up the data by gender 

qplot(x = friend_count, data = pf) + 
  scale_x_continuous(limits = c(0,1000), breaks = seq(0,1000,50)) + 
  facet_grid(gender~.)




#Omitting NA Obervations
qplot(x = friend_count, data = subset(pf, !is.na(gender)), 
      binwidth = 10) + 
  scale_x_continuous(limits = c(0,1000), breaks = seq(0,1000,50)) + 
  facet_grid(gender~.)

#equivalent ggplot syntax:

ggplot(aes(x = friend_count), data = pf) +
  geom_histogram() +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, 50)) +
  facet_wrap(~gender)




ggplot(aes(x = friend_count), data = subset(pf, !is.na(gender))) +
  geom_histogram() +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, 50)) +
  facet_wrap(~gender)



##statistics by gender

> table(pf$gender)female   male 
 40254  58574 
> by(pf$friend_count, pf$gender, summary)
pf$gender: female
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      37      96     242     244    4923 
--------------------------------------------------------------------- 
pf$gender: male
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      27      74     165     182    4917



### Tenure

qplot(x=tenure, data = pf, binwidth = 30,
      color = I('black'), fill = I('#099DD9'))



#Equivalent ggplot syntax:
ggplot(aes(x = tenure), data = pf) +
  geom_histogram(binwidth = 30, color = 'black', fill = '#099DD9')


## How would you create a histogram of tenure by year?


#create a histogram of tenure measured in years rather than in days





#It looks like the bulk of our users had less than two and a half years on Facebook.qplot(x = tenure / 365, data = pf, binwidth = 0.25, 
      color = I('black'), fill = I('#099DD9')) +
  scale_x_continuous(breaks = c(1,7,1), limits = c(0,7))
#Equivalent ggplot syntax:
ggplot(aes(x = tenure/365), data = pf) +
  geom_histogram(binwidth = .25, color = 'black', fill = '#F79420')



## Labeling Plots


qplot(x = tenure / 365, data = pf, binwidth = 0.25, 
      xlab = 'Number of years using Facebook',      ylab = 'Number of users in sample',
      color = I('black'), fill = I('#099DD9')) +
  scale_x_continuous(breaks = c(1,7,1), limits = c(0,7))


#Equivalent ggplot syntax:

ggplot(aes(x = tenure / 365), data = pf) +
  geom_histogram(color = 'black', fill = '#F79420') +
  scale_x_continuous(breaks = seq(1, 7, 1), limits = c(0, 7)) +
  xlab('Number of years using Facebook') +
  ylab('Number of users in sample')




## User Ages

> range(pf$age)
[1]  13 113
qplot(x = age, data = pf,binwidth = 5,
      xlab = 'User ages',
      ylab = 'Number of users in sample',
      color = I('black'), fill = I('blue')) +
  scale_x_continuous(breaks = c(0,113,5),limits= c(13,113))


#Equivalent ggplot syntax:

ggplot(aes(x = age), data = pf) +
  geom_histogram(binwidth = 1, fill = '#5760AB') +
  scale_x_continuous(breaks = seq(0, 113, 5))




## Transforming Data

> summary(pf$friend_count)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    31.0    82.0   196.4   206.0  4923.0 
> summary(log10(pf$friend_count))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   -Inf   1.491   1.914    -Inf   2.314   3.692 
> summary(log10(pf$friend_count + 1))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.505   1.919   1.868   2.316   3.692 
> summary(sqrt(pf$friend_count))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   5.568   9.055  11.090  14.350  70.160 
#Learn about how to use scales and how to create multiple plots on one page.  You'll also need to install and load the package gridExtra.
#Create 1 column thwi the folling thres histograms, friend count / transformed using  log 10 and square root.
> library(gridExtra)
> p1 <- qplot(x = friend_count, data = pf)
> p2 <- qplot(x = log10(friend_count+1), data = pf)
> p3 <- qplot(x = sqrt(friend_count), data = pf)
> grid.arrange(p1,p2,p3, ncol = 1)
#Equivalent ggplot syntax: we get the same outputs as we had before
pp1 <- ggplot(aes(x=friend_count), data = pf) + geom_histogram()
pp2 <- pp1 + scale_x_log10()
pp3 <- pp1 + scale_x_sqrt()

grid.arrange(pp1,pp2,pp3, ncol = 1)




## Add a Scaling Layer


> logScale <- qplot(x = log10(friend_count), data = pf)
> countScale <- ggplot(aes(x = friend_count), data = pf) +
  geom_histogram() +
  scale_x_log10()
> grid.arrange(logScale, countScale, ncol=2)

-> At the two plots, we can see that the difference is really in the x axis labeling. 
   Using scale_x_log10 will label the axis in actual friend_count. Where as using the 
   log10 wrapper will label the x axis in the log units. In general it is easier to 
   think about actual counts, so that's why people prefer using the scale_x_log10 as 
   a layer.

qplot(x=friend_count, data = pf) +
    scale_x_log10()




## Frequency Polygons : who has more friends on average men or women?


#This allows us to see the shape and the peaks of our distribution in more detail.qplot(x = friend_count, data = subset(pf, !is.na(gender)),
      binwidth = 10) +
  scale_x_continuous(lim = c(0,1000), breaks = seq(0,1000,50)) +
  facet_wrap(~gender)


qplot(x = friend_count, data = subset(pf, !is.na(gender)),
      binwidth = 10, geom='freqpoly', color = gender) +
  scale_x_continuous(lim = c(0,1000), breaks = seq(0,1000,50)) 



#But again, this plot doesn't really answer our question who has more friends on 
average men or women. Let's change the y-axis to show proportions instead of raw 
counts. 


qplot(x = friend_count, y = ..count.. / sum(..count..),
      data = subset(pf, !is.na(gender)),
      xlab = 'Friend Count',
      ylab = 'Proportion of Users with that firned count',
      binwidth = 10, geom='freqpoly', color = gender) +
  scale_x_continuous(lim = c(0,1000), breaks = seq(0,1000,50))






## Likes on the Web

qplot(x=www_likes, data = subset(pf, !is.na(gender)),
      geom = 'freqpoly', color = gender) +
  scale_x_continuous()+
  scale_x_log10()



# what's the www_like count for males?The first question is asking how many www_likes 
  there are in the entire data set for males.

> by(pf$www_likes, pf$gender, sum)

pf$gender: female
[1] 3507665
--------------------------------------------------------------------- 
pf$gender: male




## Box Plots


qplot(x = gender, y = friend_count,      data = subset(pf, !is.na(gender)),
      geom='boxplot')





## Adjust the code to focus on users who have friend counts between 0 and 1000.

qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom='boxplot',ylim = c(0,1000))




#same way
qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom = 'boxplot') +
scale_y_continuous(limits=c(0,1000))




-> Notice how the top of the box is just below 250, so it might be around 230. But this value might not be accurate for all of our data. use the ylim parameter or the scale_y_continious layer, we actually remove data points from calculations. 

# So a better way to do this is tu use the cord Cartesian layer to set the y limits instead.qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom = 'boxplot') +
  coord_cartesian(ylim= c(0,1000))



-> Here we will set the y limts from 0 to a 1000, notice how the top of the box has moved slightly closer to 250 for females.




## Box Plots, Quartiles, and Friendships

qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom = 'boxplot') +
  coord_cartesian(ylim = c(0,250))


> by(pf$friend_count, pf$gender, summary)
pf$gender: female
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      37      96     242     244    4923 
---------------------------------------------------------------------------- 
pf$gender: male
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      27      74     165     182    4917 
-> The third quartile of the 75% mark is at 244 and that's all the way up 
here(그래프를 보며). This means that 75% of female users have friend count below 244.
Or another way to say this is that 25% of female user have more than 244 friends.


## On average, who initiated more friendships in our sample: men or women?

> names(pf)
 [1] "userid"                "age"                   "dob_day"              
 [4] "dob_year"              "dob_month"             "gender"               
 [7] "tenure"                "friend_count"          "friendships_initiated"
[10] "likes"                 "likes_received"        "mobile_likes"         
[13] "mobile_likes_received" "www_likes"             "www_likes_received"   
qplot(x = gender, y = friendships_initiated,
      data = subset(pf, !is.na(gender)), geom = 'boxplot') +
  coord_cartesian(ylim = c(0,150))

> by(pf$friendships_initiated, pf$gender, summary)
pf$gender: female
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    19.0    49.0   113.9   124.8  3654.0 
---------------------------------------------------------------------------- 
pf$gender: male
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    15.0    44.0   103.1   111.0  4144.0 

-> On average, who initiated more friendships in our sample: men or women?
Women.


## Getting Logical : What percent of check in using mobile?

> head(pf$mobile_likes)
[1] 0 0 0 0 0 0
> summary(pf$mobile_likes)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0     4.0   106.1    46.0 25110.0 
> summary(pf$mobile_likes > 0)
   Mode   FALSE    TRUE    NA's 
logical   35056   63947       0 
> mobile_check_in <- NA
> pf$mobile_check_in <- ifelse(pf$mobile_likes>0, 1, 0)
> pf$mobile_check_in <- factor(pf$mobile_check_in)
> summary(pf$mobile_check_in)
    0     1 
35056 63947 
> sum(pf$mobile_check_in == 1) / length(pf$mobile_check_in) 
[1] 0.6459097
#what percent of check in using mobile? 65%




#reference

1. https://cn.udacity.com/course/data-analysis-with-r--ud651

2. http://www.cookbook-r.com/Graphs/Facets_(ggplot2)

3. https://en.wikipedia.org/wiki/Web_colors

4. http://ggplot2.tidyverse.org/reference/theme.html

5. http://lightonphiri.org/blog/ggplot2-multiple-plots-in-one-graph-using-gridextra

6. http://ggplot2.tidyverse.org/reference/scale_continuous.html

7. https://en.wikipedia.org/wiki/Linear_regression#Assumptions

8. https://en.wikipedia.org/wiki/Normal_distribution

9. https://www.r-statistics.com/2013/05/log-transformations-for-skewed-and-wide-distributions-from-practical-data-science-with-r/


반응형

'Python, R 분석과 프로그래밍 > Data Analysis with R' 카테고리의 다른 글

Lesson2: R Basic  (0) 2017.07.11
Posted by 마르띤
,

Lesson2: R Basic

Python, R 분석과 프로그래밍/Data Analysis with R 2017. 7. 11. 10:14
반응형


# 원하는 데이터만 발라서 보기

> data(mtcars)

> mean(mtcars$mpg)

[1] 20.09062

> subset(mtcars, mtcars$mpg >= 30 | mtcars$hp < 60)

                mpg cyl disp  hp drat    wt  qsec vs am gear carb

Fiat 128       32.4   4 78.7  66 4.08 2.200 19.47  1  1    4    1

Honda Civic    30.4   4 75.7  52 4.93 1.615 18.52  1  1    4    2

Toyota Corolla 33.9   4 71.1  65 4.22 1.835 19.90  1  1    4    1

Lotus Europa   30.4   4 95.1 113 3.77 1.513 16.90  1  1    5    2

> mtcars[mtcars$mpg >= 30 | mtcars$hp < 60, ] #column 전체를 보기 위해서는 콤마 필수

                mpg cyl disp  hp drat    wt  qsec vs am gear carb

Fiat 128       32.4   4 78.7  66 4.08 2.200 19.47  1  1    4    1

Honda Civic    30.4   4 75.7  52 4.93 1.615 18.52  1  1    4    2

Toyota Corolla 33.9   4 71.1  65 4.22 1.835 19.90  1  1    4    1

Lotus Europa   30.4   4 95.1 113 3.77 1.513 16.90  1  1    5    2




#reddit data 보기, 범주형 변수 그래프 그리기

> reddit <- read.csv('reddit.csv')

> library(ggplot2)

> qplot(data = reddit, x = age.range)





#오더 순서 정하기 setting levles of ordered factors solution


> reddit$age.range = ordered(reddit$age.range,levels=c("Under 18","18-24", "25-34","35-44","45-54","55-64","65 or Above"))

> qplot(data = reddit, x = age.range)



#alternative solution

> reddit$age.range = factor(reddit$age.range, levels=c("Under 18","18-24", "25-34","35-44","45-54","55-64","65 or Above"),ordered=T)

> qplot(data = reddit, x = age.range)




# practice

> nlevels(reddit$income.range)

[1] 8

> levels(reddit$income.range)

[1] "$100,000 - $149,999" "$150,000 or more"    "$20,000 - $29,999"   "$30,000 - $39,999"  

[5] "$40,000 - $49,999"   "$50,000 - $69,999"   "$70,000 - $99,999"   "Under $20,000"      

> qplot(data = reddit, x = income.range)



#아래같은 방법도 가능

> reddit$income.range = ordered(reddit$income.range, levels=c("Under $20,000" , "$20,000 - $29,999"   , "$30,000 - $39,999", "$40,000 - $49,999","$50,000 - $69,999" , "$70,000 - $99,999”, ”$100,000 - $149,999" , "$150,000 or more"))

> qplot(data = reddit, x = income.range)




#다른 예제

> tShirts <- factor(c('medium', 'small', 'large', 'medium', 'large', 'large'), levels = c('medium','small','large'))

> tShirts

[1] medium small  large  medium large  large 

Levels: medium small large

> qplot(x = tShirts)



> tShirts <- ordered(tShirts, levels = c('small', 'medium', 'large'))

> tShirts

[1] medium small  large  medium large  large 

Levels: small < medium < large

> qplot(x = tShirts)








참고

https://cn.udacity.com/course/data-analysis-with-r--ud651

https://cn.udacity.com/course/data-wrangling-with-mongodb--ud032

http://vita.had.co.nz/papers/tidy-data.pdf

http://courses.had.co.nz.s3-website-us-east-1.amazonaws.com/12-rice-bdsi/slides/07-tidy-data.pdf

http://www.computerworld.com/article/2497143/business-intelligence/business-intelligence-beginner-s-guide-to-r-introduction.html

http://www.statmethods.net/index.html

https://www.r-bloggers.com/

http://www.cookbook-r.com/

http://blog.revolutionanalytics.com/2013/08/foodborne-chicago.html

http://blog.yhat.com/posts/roc-curves.html

https://github.com/corynissen/foodborne_classifier

반응형

'Python, R 분석과 프로그래밍 > Data Analysis with R' 카테고리의 다른 글

Lesson 3: Explore One Variable  (0) 2017.07.12
Posted by 마르띤
,

3장 인자분석 - 주성분인자법(principal factor method) 사례

KNOU/3 다변량분석 2017. 3. 27. 19:59
반응형

### 5주차, 제3장 인자분석 ###

 

아래 자료는 무료 검진프로그램인 PHI를 개발하여 프로그램 유효서을 모니터링 한 자료이다. 11개 검진 항목, 128명의 자료가 이으며 R을 이용하여 주성분인자법(principal factor method)를 이용한 인자분석을 해 보자.

 

1. 자료 가져오기 및 요약 통계

> setwd('C:/Rwork')

> med.data = read.table('medFactor.txt',header=T)

> head(med.data)

lung muscle liver skeleton kidneys heart step stamina stretch blow urine

1   20     16    52       10      24    23   19      20      23   29    67

2   24     16    52        7      27    16   16      15      31   33    59

3   19     21    57       18      22    23   16      19      42   40    61

4   24     21    62       12      31    25   17      17      36   36    77

5   29     18    62       14      26    27   15      20      33   29    88

6   18     19    51       15      29    23   19      20      50   37    54

> boxplot(med.data)

 

2. 초기 인자분석 실행하기

> library(psych)

> library(GPArotation) #인자회전

> med.factor = principal(med.data,rotate='none') #rotate에 ‘none’, ‘varimax’,’quanrtmax’, ‘oblimin’ 등 방법을 사용

> names(med.factor)

[1] "values"       "rotation"     "n.obs"        "communality"  "loadings"     "fit"        

[7] "fit.off"      "fn"           "Call"         "uniquenesses" "complexity"   "chi"        

[13] "EPVAL"        "R2"           "objective"    "residual"     "rms"          "factors"    

[19] "dof"          "null.dof"     "null.model"   "criteria"     "STATISTIC"    "PVAL"       

[25] "weights"      "r.scores"     "Structure"    "scores"     

> med.factor$values #고유근

[1] 3.3791814 1.4827707 1.2506302 0.9804771 0.7688022 0.7330511 0.6403994 0.6221934 0.5283718 0.3519301

[11] 0.2621928

> plot(med.factor$values,type='b')

 

>> med.factor$values #고유근을 통해 세 번째 인자까지 고유근이 1 이상인 것을 알 수 있다. 스크리 그림에서는 4번째 인자다음부터 그래프의 기울기가 완만해지는 것을 볼 수 있다. 따라서 유효한 인자의 수는 3-4개로 생각할 수 있다.

 

3. 인자 회전

(1) 직교회전(orthogonal rotation)의 Varimax

> med.varimax = principal(med.data,nfactors=3,rotate='varimax')

> med.varimax

 

>> 공통성은 많은 문항(변수) 중에서 서로연관이 있는 일부 문항들을 선택하기 위하여 인자분석을 실시하는 경우 문항선택의 기준으로 이용된다. 예를 들면, 문항 100개를 이용하여 인자분석을 한 결과 공통성이 0.3 이하인 것이 40개라면 이 40개 문항은 다른 문항들과 공통점이 별로 없는 것으로 판단할 수 있다.


>> 고유분산은 u2 = 1 – h2(공통성) 공식으로 바로 구할 수 있다.lung의 경우 0.53 = 1 – 0.47


>> ss loading은 각 인자에 의해 설명되는 분산의 양을 나타낸다.

   RC1 2.39 = 0.662 + 0.112 + 0.782 + … + 0.252 + (-0.07)2


>> Proportion Var는 각 인자가 설명하는 총 분산의 비율을 의미한다. RC1은 총 분산의 22%, RC2가 19%, RC3는 14%, 아래 Cumulative var는 누적값. 세 인자에 의해 설명되어지는 변동은 총 변동의 56%.



>> 회전된 요인에 대한 변수들의 요인 적재값을 보면 첫 번째 인자는 lung, liver, kidneys, heart가 높은 값을 가지며, step은 상대적으로 높은 적재값을 가진다. (적재값이 뭘 의미하는지 잘 모르겠다 ㅠ)


두 번째 인자는 stamina, stretch, blow, urine의 값이 높다. 마지막으로 세번째 인자는 muscle과 skeleton에서 높은 값을 가진다. 따라서 첫 번째 인자는 생물의학(biomedical), 두번째 인자는 인체기능(performance), 세번째 인자는 근육골계통력(Muscular – skeletal strength)으로 특정 지을 수 있다.


>> 인자분석은 주성분분석과 달리 인자들의 결합으로 원변수들을 표현한다.

 Lung = 0.66 x factor1 + 0.12 x facotr2 + 0.16 x factor3 + error

 


인자점수 추정방법은 회귀분석, Barlett 등이 있는데, 아래는 회귀분석을 이용한 인자점수이다.

> head(med.varimax$scores)

             RC1        RC2         RC3

[1,] -0.11970907 -0.2574867 -0.74473196

[2,]  0.05696634 -0.4784611 -1.53199247

[3,] -0.59153602  0.8957933  1.52598533

[4,]  1.20919164  0.3179760 -0.42022213

[5,]  0.82291043  0.1513448 -0.02988744

[6,] -0.08606120  1.1451913  0.76290104 

 


인자분석의 시각화 몇가지 사례


> plot(med.varimax)

> plot(med.varimax$loadings[,c(1:2)])

> text(med.varimax$loadings[,c(1:2)],colnames(med.data))


> biplot(med.varimax$scores,med.varimax$loadings)

> fa.diagram(med.varimax)


> omega(med.data)

 





3. 인자 회전

(2) 사곽회전(Oblique rotation) - Oblimin

> med.oblimin = principal(med.data,nfactors=3,rotate='oblimin',scores=T,method='regression')

> med.oblimin

 


>> Oblimin을 이용한 인자분석에 대한 결과는 Varimax와 마찬가지로 세 인자에 대해 묶여지는 변수는 같으며 다만 인자 적재값이 차이가 나는 것을 볼 수 있다. Oblimin 방법을 이용한 결과에 대한 각 변수의 인자 모형은 아래와 같다.

Lung = 0.66 x factor1 + 0.01 x factor2 + 0.1 x facotr3

 

>> 회귀분석을 이용한 인자점수 함수식은 아래와 같다

Factor1 = 0.66 x lung + 0.02 x muscle + …+ (-0.09) x urine

 

인자점수는 아래와 같다.

> head(med.oblimin$scores)

TC1        TC2         TC3

[1,] -0.2325231 -0.2639103 -0.77351112

[2,] -0.1722964 -0.4641659 -1.54539480

[3,] -0.2852637  0.8308965  1.50251465

[4,]  1.1976234  0.4291411 -0.22225725

[5,]  0.8296036  0.2260879  0.09586509

[6,]  0.1766403  1.1289891  0.83993963 

 

 

4. 행렬도

> biplot(med.varimax)

 

 

참고문헌: 

[1] 다변량분석(김성수, 김현중, 정성석, 이용구 공저)

[2] R 프로그램에 기반한 다변량분석 및 데이터마이닝(이재길)


반응형

'KNOU > 3 다변량분석' 카테고리의 다른 글

2장 주성분분석 - heptathlon data / beer data  (0) 2017.03.27
Posted by 마르띤
,

2장 주성분분석 - heptathlon data / beer data

KNOU/3 다변량분석 2017. 3. 27. 19:13
반응형

[예제1] - heptathlon data



1. 자료 준비

> library(HSAUR2)

> library(tools)

> data("heptathlon")

> head(heptathlon)

> heptathlon은 R의 HSAUR2 패키지에 들어있는 자료로서, 1988년 서울 올림픽 여성 7종 경기에 대한 결과이다. Hurdles(110m 허들), highjump(높이뛰기), shot(포환던지기), run200m(200미터 달리기), longjump(멀리뛰기), javelin(창던지기), run800m(800미터 달리기),score(점수)로 구성되어 있다.


2. 자료 변형

> heptathlon$hurdles <- max(heptathlon$hurdles) - heptathlon$hurdles

> heptathlon$run200m <- max(heptathlon$run200m) - heptathlon$run200m

> heptathlon$run800m <- max(heptathlon$run800m) - heptathlon$run800m

>  작은 값일수록 좋은 점수이기 때문에 최대값에서 빼줌


 

3. 주성분분석

> library(stats)

> hep.data <- heptathlon[,-8] #주성분점수에서 예외

> heptathlon.pca <- princomp(hep.data,cor=T,scores=T) #cor = T, 상관계수 행령, F는 공분산 행렬, scores=T 주성분점수

> heptathlon.pca #변수가 7개이니, 주성분도 7개

Call:

  princomp(x = hep.data, cor = T, scores = T)

 

Standard deviations:

  Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6    Comp.7

2.1119364 1.0928497 0.7218131 0.6761411 0.4952441 0.2701029 0.2213617

 

7  variables and  25 observations.

> summary(heptathlon.pca)

> 주성분의 표준편차, 분산비율을 알수 있다. 첫 번째 주성분이 63.72%, 두 번째 주 성분이 17.06%의 분산 비율을 나타내고 있으며, 2개의 주성분이 총 변이의 80.8%의 정보를 가지고 있다.

 

> eig.val = heptathlon.pca$sdev^2 #주성분의 표준편차를 제곱

> eig.val #주성분 분산

Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6     Comp.7

4.46027516 1.19432056 0.52101413 0.45716683 0.24526674 0.07295558 0.04900101

> 제1주성분의 분산이 4.46, 제2주성분이 1.19, 제3주성분이 0.52로 유효한 주성분의 수는 2개임을 알 수 있다.

 

> par(mfrow=c(1,2))

> screeplot(heptathlon.pca)

> screeplot(heptathlon.pca,type='line',pch=19,main='Scree Plot of heptathlon.pca')

 

> 스크리 그림은 주성분의 고유값 크기 순으로 그린 것으로, 고유값이 1보다 큰 주성분값이 2개임을 알 수 있다.

 

> heptathlon.pca$loadings[,1:2]

             Comp.1      Comp.2

hurdles  -0.4528710  0.15792058

highjump -0.3771992  0.24807386

shot     -0.3630725 -0.28940743

run200m  -0.4078950 -0.26038545

longjump -0.4562318  0.05587394

javelin  -0.0754090 -0.84169212

run800m  -0.3749594  0.22448984

> 상관계수행렬을 이용한 경우 제1, 제2주성분은 다음과 같다.

PC1 = -0.453 x hurldes – 0.377 x highjump - ․․․ - 0.374 x rum800m

PC2 = 0.157 x hurldes + 0.248 x highjump + ․․․ + 0.224 x rum800m

 

여기서 제1주성분은 jvelin(창던지기)를 제외한 모든 변수의 절대값이 큰 값을 가지는 것으로 볼 때, 전반적인 체력 정도를 나타내는 성분이라고 할 수 있고, 제2성분은 jvelin의 계수가 다른 변수에 비해 상대적으로 절대값이 큰 것을 볼 때, 창던지기와 밀접한 관련이 있는 성분으로 파악할 수 있다.

 


4. 주성분 점수 및 행렬도 (biplot)

> head(heptathlon.pca$scores[,1:2])

                         Comp.1      Comp.2

Joyner-Kersee (USA) -4.206435 -1.26802363

John (GDR)          -2.941619 -0.53452561

Behmer (GDR)        -2.704271 -0.69275901

Sablovskaite (URS)  -1.371052 -0.70655862

Choubenkova (URS)   -1.387050 -1.78931718

Schulz (GDR)        -1.065372  0.08104469

> biplot(heptathlon.pca,cex=0.7,col=c('red','blue'),main='Biplot')

 

> 첫번째 주성분을 x축, 두 번째 주성분을 y축에 위치. Highjump, run800m, hurldes, longjump가 서로 

가까운 곳에 위치하고 벡터의 방향(화살표) 또한 비슷하다. Ru200m, shot이 가깝게 위치하고 있고, javelin은 다른 변수들과 다른 방향에 위치하고 있다.





[예제2] - beer data


99명의 소비자가 맥주 구매에 중요하게 생각하는 요인에 대해 점수를 준 엑셀 자료이다. 독립변수는 7개, 케이스 수는 99개인 자료이다. 점수는 0~100으로 되어 있다.

상관관계가 아닌 공분산을 활용해보자.

 

1. 자료 읽기

> beer<-read.csv('beer.csv',header=T)

> head(beer)

   cost size alcohol reputat color aroma taste

1   10   15      20      85    40    30    50

2  100   70      50      30    75    60    80

3   65   30      35      80    80    60    90

4    0    0      20      30    80    90   100

5   10   25      10     100    50    40    60

6   25   35      30      40    45    30    65

 

 

2. 상관계수행렬 및 산점도 행렬 보기

> round(cor(beer),2)

cost  size alcohol reputat color aroma taste

cost     1.00  0.88    0.88   -0.17  0.32 -0.03  0.05

size     0.88  1.00    0.82   -0.06  0.01 -0.29 -0.31

alcohol  0.88  0.82    1.00   -0.36  0.40  0.10  0.06

reputat -0.17 -0.06   -0.36    1.00 -0.52 -0.52 -0.63

color    0.32  0.01    0.40   -0.52  1.00  0.82  0.80

aroma   -0.03 -0.29    0.10   -0.52  0.82  1.00  0.87

taste    0.05 -0.31    0.06   -0.63  0.80  0.87  1.00

> 맥주의 가격(cost),크기(size), 알코올(alcohol)함량의 상관계수가 높아 서로 양의 상관관계가 있을 것으로 보인다. 또한 색(color), 향(aroma), taste(맛)이 서로 상관계수가 높게 나타났다. 반면 평판(reputation)은 비교적 다른 변수들과 상관계수가 높지 않으며, 다른 변수들과 모두 음의 상관 계수를 나타냈다.

 

> plot(beer[1:3])

 

 


3. 주성분분석 실행하기

> library(stats)

> beer.pca = princomp(beer,cor=F, scores=T) #cor = F 공분산

> beer.pca

Call:

  princomp(x = beer, cor = F, scores = T)

 

Standard deviations:

  Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6    Comp.7

39.117489 34.783073 17.974805  8.456192  6.479111  4.381410  2.579257

 

7  variables and  99 observations.

> summary(beer.pca)

 

 

> screeplot(beer.pca,type='l',pch=19)

> beer.pca$loadings[,1:3] #주성분계수

Comp.1      Comp.2     Comp.3

cost     0.7344197  0.31868378 -0.1902272

size     0.3942553  0.34017182  0.1274959

alcohol  0.2831873  0.06888179  0.0370581

reputat -0.3356901  0.49057105 -0.7861697

color    0.2657171 -0.34379889 -0.3862210

aroma    0.1398211 -0.48510439 -0.3693624

taste    0.1488354 -0.42871339 -0.2062207

> screeplot을 보면 4번째 주성분부터 완만해짐로 유효한 주성분을 3개로 판단하여, 3개의 주성분 계수를 보자. 제1주성분을 보면 cost, size, alcohol의 점수가 높다. 비용, 크기, 알코올이 서로 연관이 있고, color, aroma, taste의 주성분점수가 상대적으로 낮게 측정되어 있다. 색, 향, 맛이 서로 연관되어 있다.

PC1 = 0.734 x cost + 0.397 x size + … + 0.148 x taste

 


4. 주성분 분석 점수 및 행렬도

> head(beer.pca$scores[,1:3])

Comp.1    Comp.2     Comp.3

[1,] -41.435562  40.03360   4.340612

[2,]  91.264626  23.06214   7.798284

[3,]  31.574345  15.79053 -34.501270

[4,]  -9.770913 -59.53112  -0.351844

[5,] -39.816498  37.52890 -16.165598

[6,]  -1.035114  22.08072  34.760924

> biplot(beer.pca,col=c('red','blue'),main = 'Biplot of Beer.pca')

 

 


참고 문헌

[1] R을 이용한 다변량 분석, 김성수 김현중 정성석 이용구 공저

반응형

'KNOU > 3 다변량분석' 카테고리의 다른 글

3장 인자분석 - 주성분인자법(principal factor method) 사례  (0) 2017.03.27
Posted by 마르띤
,

제3장 예측데이터의 기초분석

KNOU/3 예측방법론 2017. 3. 22. 20:05
반응형

### 3,4주차, 제3장 예측데이터의 기초 분석 ###

 

3.1 시계열의 요약

> gdp <- read.csv("gdpq.csv", header = TRUE)

> head(gdp)

gdp   gdpsa

1 12807.5 14934.6

2 15732.5 15246.5

3 14621.6 15472.1

4 18689.5 16198.0

5 14512.5 16664.4

6 17709.4 17180.3

> gdp_o  <- ts(gdp[,1]/1000, start=1970, frequency=4)

> gdp_o

Qtr1     Qtr2     Qtr3     Qtr4

1970  12.8075  15.7325  14.6216  18.6895

1971  14.5125  17.7094  16.5814  19.5058

1972  15.2014  18.4940  17.5132  21.5480

1973  17.1063  20.9626  20.6994  24.7473

1974  19.7445  23.2128  22.2597  26.1339

1975  20.0790  24.6205  23.9055  29.4487

1976  22.5622  28.2495  27.2866  33.1492

1977  24.5630  30.6101  30.7895  38.4318

1978  28.1520  34.5775  33.4854  40.9892

1979  32.1360  38.2981  36.1002  42.1807

1980  32.0372  36.8905  36.8476  40.1279

1981  33.0530  38.8709  39.6106  45.1715

1982  36.2445  41.7473  43.0336  48.6735

1983  40.6951  46.9821  49.2677  53.4268

1984  46.3244  52.1396  53.7155  56.9611

1985  49.3720  55.4260  57.4681  62.4993

1986  54.3321  61.6800  65.9805  70.2836

1987  62.3064  70.6516  72.8342  77.4281

1988  72.3188  76.2441  80.4584  87.2238

1989  75.6243  81.9830  85.7357  94.2551

1990  82.9977  90.4118  94.2292 101.3471

1991  91.5677 100.2760 102.4124 110.5688

1992  98.9458 107.5298 107.0952 114.5932

1993 103.2117 113.8639 115.0267 123.1619

1994 112.7728 123.0717 124.0430 135.3119

1995 123.1052 134.6932 136.4655 145.1600

1996 132.1325 144.3954 146.0147 155.6439

1997 139.4080 154.6306 155.6320 161.8583

1998 134.4830 143.3159 144.5975 154.1904

1999 143.0859 159.4288 161.7758 174.1675

2000 161.1069 174.1233 176.7639 182.6340

2001 167.9311 181.3261 181.9274 191.0444

2002 179.0416 194.0302 194.2672 206.5295

2003 185.3180 197.5801 198.0618 214.5984

2004 194.9231 209.3292 207.5780 220.4750

2005 200.2166 216.3471 216.9831 231.6941

2006 212.4693 227.4156 227.9013 242.2628

2007 221.9311 239.4580 239.1251 256.0003

2008 234.1789 249.8902 246.9623 247.4675

2009 224.3595 244.7084 249.5196 263.0376

2010 243.8225 263.1901 260.7457 275.9080

2011 254.2480 272.3620 270.2561 285.2294

2012 261.4804 278.8836 274.4844 289.3663

2013 265.2912 285.3123 283.5950 300.6547

> gdp_sa <- ts(gdp[,2]/1000, start=1970, frequency=4)

> gdp_sa

Qtr1     Qtr2     Qtr3     Qtr4

1970  14.9346  15.2465  15.4721  16.1980

1971  16.6644  17.1803  17.3029  17.1614

1972  17.5141  17.9938  18.2651  18.9837

1973  19.7739  20.4964  21.4395  21.8057

1974  22.8617  22.6609  22.9105  22.9179

1975  23.5435  24.1886  24.7458  25.5757

1976  26.7061  27.7776  28.2670  28.4967

1977  28.8486  30.2254  31.9596  33.3608

1978  32.9248  33.9149  34.2955  36.0687

1979  37.5335  37.7262  36.7787  36.6767

1980  36.6760  36.2477  36.8647  36.1147

1981  37.2882  38.7201  39.7069  40.9907

1982  40.9807  41.7302  42.6719  44.3161

1983  45.5381  46.9978  48.6324  49.2035

1984  51.1110  52.0170  52.8702  53.1425

1985  54.4390  55.6003  56.4163  58.3099

1986  59.6636  62.0513  64.7921  65.7692

1987  67.5395  70.8755  71.8664  72.9389

1988  77.9135  76.8991  79.7696  81.6630

1989  81.4808  82.6901  85.2974  88.1298

1990  89.1127  90.7379  93.9483  95.1870

1991  98.3701 100.0831 102.2049 104.1669

1992 105.8722 106.9636 106.9913 108.3369

1993 110.1757 112.9698 115.0582 117.0605

1994 120.1321 122.1115 124.0878 128.8680

1995 131.3819 133.6788 136.2208 138.1425

1996 141.2176 143.2933 145.5552 148.1203

1997 149.0865 153.3882 154.8591 154.1950

1998 143.4587 141.9956 143.8028 147.3296

1999 151.7891 157.9815 161.6553 167.0321

2000 170.3475 172.5127 176.8317 174.9363

2001 177.4443 179.7859 182.2794 182.7193

2002 189.0423 192.4812 195.1160 197.2289

2003 195.8857 196.0955 199.1527 204.4243

2004 206.1241 207.7387 208.5458 209.8967

2005 211.6290 214.7745 218.1001 220.7373

2006 224.3165 225.6878 229.0608 230.9839

2007 234.2216 237.6344 240.2943 244.3642

2008 246.6202 247.5462 247.8380 236.4944

2009 236.6559 242.6501 250.8293 251.4898

2010 256.9193 260.5657 262.1134 264.0678

2011 267.5257 269.7017 271.9174 272.9508

2012 275.1821 276.0057 276.1202 276.9067

2013 279.2230 282.3449 285.3485 287.9793

> gdp_gr <- ts((gdp_sa-lag(gdp_sa,-1))/lag(gdp_sa,-1)*100,  start=c(1970,2), frequency=4)

> gdp_sa-lag(gdp_sa,-1)

Qtr1     Qtr2     Qtr3     Qtr4

1970            0.3119   0.2256   0.7259

1971   0.4664   0.5159   0.1226  -0.1415

1972   0.3527   0.4797   0.2713   0.7186

1973   0.7902   0.7225   0.9431   0.3662

1974   1.0560  -0.2008   0.2496   0.0074

1975   0.6256   0.6451   0.5572   0.8299

1976   1.1304   1.0715   0.4894   0.2297

1977   0.3519   1.3768   1.7342   1.4012

1978  -0.4360   0.9901   0.3806   1.7732

1979   1.4648   0.1927  -0.9475  -0.1020

1980  -0.0007  -0.4283   0.6170  -0.7500

1981   1.1735   1.4319   0.9868   1.2838

1982  -0.0100   0.7495   0.9417   1.6442

1983   1.2220   1.4597   1.6346   0.5711

1984   1.9075   0.9060   0.8532   0.2723

1985   1.2965   1.1613   0.8160   1.8936

1986   1.3537   2.3877   2.7408   0.9771

1987   1.7703   3.3360   0.9909   1.0725

1988   4.9746  -1.0144   2.8705   1.8934

1989  -0.1822   1.2093   2.6073   2.8324

1990   0.9829   1.6252   3.2104   1.2387

1991   3.1831   1.7130   2.1218   1.9620

1992   1.7053   1.0914   0.0277   1.3456

1993   1.8388   2.7941   2.0884   2.0023

1994   3.0716   1.9794   1.9763   4.7802

1995   2.5139   2.2969   2.5420   1.9217

1996   3.0751   2.0757   2.2619   2.5651

1997   0.9662   4.3017   1.4709  -0.6641

1998 -10.7363  -1.4631   1.8072   3.5268

1999   4.4595   6.1924   3.6738   5.3768

2000   3.3154   2.1652   4.3190  -1.8954

2001   2.5080   2.3416   2.4935   0.4399

2002   6.3230   3.4389   2.6348   2.1129

2003  -1.3432   0.2098   3.0572   5.2716

2004   1.6998   1.6146   0.8071   1.3509

2005   1.7323   3.1455   3.3256   2.6372

2006   3.5792   1.3713   3.3730   1.9231

2007   3.2377   3.4128   2.6599   4.0699

2008   2.2560   0.9260   0.2918 -11.3436

2009   0.1615   5.9942   8.1792   0.6605

2010   5.4295   3.6464   1.5477   1.9544

2011   3.4579   2.1760   2.2157   1.0334

2012   2.2313   0.8236   0.1145   0.7865

2013   2.3163   3.1219   3.0036   2.6308

> gdp_gr

Qtr1         Qtr2         Qtr3         Qtr4

1970               2.088438927  1.479683862  4.691670814

1971  2.879367823  3.095821032  0.713608028 -0.817781990

1972  2.055193632  2.738936057  1.507741555  3.934279035

1973  4.162518371  3.653806280  4.601295837  1.708062222

1974  4.842770468 -0.878324884  1.101456694  0.032299601

1975  2.729743999  2.740034404  2.303564489  3.353700426

1976  4.419820376  4.012191971  1.761851276  0.812608342

1977  1.234879828  4.772501959  5.737558477  4.384285160

1978 -1.306923095  3.007155700  1.122220617  5.170357627

1979  4.061138882  0.513408022 -2.511517195 -0.277334435

1980 -0.001908569 -1.167793653  1.702176966 -2.034466576

1981  3.249369370  3.840088822  2.548547137  3.233191209

1982 -0.024395778  1.828909706  2.256639077  3.853121141

1983  2.757462863  3.205447746  3.478035142  1.174320001

1984  3.876756735  1.772612549  1.640233001  0.515034935

1985  2.439666933  2.133213321  1.467617980  3.356476763

1986  2.321561176  4.001937530  4.416990458  1.508054223

1987  2.691685470  4.939331798  1.398085375  1.492352476

1988  6.820228986 -1.301956657  3.732813518  2.373585927

1989 -0.223112058  1.484153322  3.153098134  3.320617041

1990  1.115286770  1.823758005  3.538102601  1.318491128

1991  3.344049082  1.741382798  2.120038248  1.919673127

1992  1.637084333  1.030865515  0.025896660  1.257672353

1993  1.697297966  2.536040161  1.848635653  1.740249717

1994  2.623942320  1.647686172  1.618438886  3.852272343

1995  1.950755812  1.748262127  1.901573024  1.410724353

1996  2.226034711  1.469859281  1.578510649  1.762286748

1997  0.652307618  2.885371915  0.958939475 -0.428841444

1998 -6.962806836 -1.019875407  1.272715493  2.452525264

1999  3.026886654  4.079607824  2.325462159  3.326089525

2000  1.984887935  1.271048885  2.503583794 -1.071866639

2001  1.433664711  1.319625370  1.386927451  0.241332811

2002  3.460499247  1.819116674  1.368860959  1.082894278

2003 -0.681036096  0.107103275  1.559036286  2.647014075

2004  0.831505843  0.783314518  0.388516921  0.647771377

2005  0.825310736  1.486327488  1.548414733  1.209169551

2006  1.621474939  0.611323732  1.494542461  0.839558755

2007  1.401699426  1.457081670  1.119324475  1.693714749

2008  0.923212156  0.375476137  0.117876986 -4.577022087

2009  0.068289143  2.532875791  3.370779571  0.263326493

2010  2.158934478  1.419278349  0.593976874  0.745631471

2011  1.309474309  0.813379799  0.821537276  0.380041880

2012  0.817473332  0.299292723  0.041484650  0.284839718

2013  0.836491136  1.118066921  1.063805296  0.921960340

> mean(gdp_gr);median(gdp_gr);var(gdp_gr);sd(gdp_gr)

[1] 1.720042

[1] 1.618439

[1] 2.977614

[1] 1.725576

> (1+mean(gdp_gr)/100)^4-1 #연평균 7.05% 성장

[1] 0.07059725

> plot(gdp_gr, ylab="경제성장률(전 분기대비)", xlab="연도", col="steelblue", main="경제성장률 추이")

> abline(h=0, lty=2, col="gray")

3.2 시계열의 분포

> hist(gdp_gr, breaks = 12, col = "lightblue", border = "black", freq=FALSE, main="경제성장률의 확률분포", xlab="", xlim=c(-10, 10))

> lines(density(gdp_gr)) #꼬리부분에 큰 값을 가지고 있다. IMF



> shapiro.test(gdp_gr) #p-value = 4.285e-06, 경제성장률은 정규분포를 따른다는 귀무가설을 기각,

 

 

Shapiro-Wilk normality test

 

data:  gdp_gr

W = 0.94723, p-value = 4.285e-06

 

 

3.3 자기 상관

시계열은 시간의 전개에 따른 시간영역 정보와 시간에 따라 반복되는 주파수 영역 정보를 가지고 있다. 시간 영역 정보를 파악하려면 시계열의 과거, 현재, 미래를 연결하는 구조를 이해해야 한다. 시간에 따른 시계열의 연결 구조를 이해하는 데에는 표본 자기 상관계수, 표본 부분자기 상관계수 등이 이용된다.

> set.seed(1)

> nn=length(gdp_o)

> wn=ts(rnorm(nn), start=1970, frequency=4)

> wn

Qtr1         Qtr2         Qtr3         Qtr4

1970 -0.626453811  0.183643324 -0.835628612  1.595280802

1971  0.329507772 -0.820468384  0.487429052  0.738324705

1972  0.575781352 -0.305388387  1.511781168  0.389843236

1973 -0.621240581 -2.214699887  1.124930918 -0.044933609

1974 -0.016190263  0.943836211  0.821221195  0.593901321

1975  0.918977372  0.782136301  0.074564983 -1.989351696

1976  0.619825748 -0.056128740 -0.155795507 -1.470752384

1977 -0.478150055  0.417941560  1.358679552 -0.102787727

1978  0.387671612 -0.053805041 -1.377059557 -0.414994563

1979 -0.394289954 -0.059313397  1.100025372  0.763175748

1980 -0.164523596 -0.253361680  0.696963375  0.556663199

1981 -0.688755695 -0.707495157  0.364581962  0.768532925

1982 -0.112346212  0.881107726  0.398105880 -0.612026393

1983  0.341119691 -1.129363096  1.433023702  1.980399899

1984 -0.367221476 -1.044134626  0.569719627 -0.135054604

1985  2.401617761 -0.039240003  0.689739362  0.028002159

1986 -0.743273209  0.188792300 -1.804958629  1.465554862

1987  0.153253338  2.172611670  0.475509529 -0.709946431

1988  0.610726353 -0.934097632 -1.253633400  0.291446236

1989 -0.443291873  0.001105352  0.074341324 -0.589520946

1990 -0.568668733 -0.135178615  1.178086997 -1.523566800

1991  0.593946188  0.332950371  1.063099837 -0.304183924

1992  0.370018810  0.267098791 -0.542520031  1.207867806

1993  1.160402616  0.700213650  1.586833455  0.558486426

1994 -1.276592208 -0.573265414 -1.224612615 -0.473400636

1995 -0.620366677  0.042115873 -0.910921649  0.158028772

1996 -0.654584644  1.767287269  0.716707476  0.910174229

1997  0.384185358  1.682176081 -0.635736454 -0.461644730

1998  1.432282239 -0.650696353 -0.207380744 -0.392807929

1999 -0.319992869 -0.279113303  0.494188331 -0.177330482

2000 -0.505957462  1.343038825 -0.214579409 -0.179556530

2001 -0.100190741  0.712666307 -0.073564404 -0.037634171

2002 -0.681660479 -0.324270272  0.060160440 -0.588894486

2003  0.531496193 -1.518394082  0.306557861 -1.536449824

2004 -0.300976127 -0.528279904 -0.652094781 -0.056896778

2005 -1.914359426  1.176583312 -1.664972436 -0.463530401

2006 -1.115920105 -0.750819001  2.087166546  0.017395620

2007 -1.286300530 -1.640605534  0.450187101 -0.018559833

2008 -0.318068375 -0.929362147 -1.487460310 -1.075192297

2009  1.000028804 -0.621266695 -1.384426847  1.869290622

2010  0.425100377 -0.238647101  1.058483049  0.886422651

2011 -0.619243048  2.206102465 -0.255027030 -1.424494650

2012 -0.144399602  0.207538339  2.307978399  0.105802368

2013  0.456998805 -0.077152935 -0.334000842 -0.034726028

> par(mfrow=c(2,1))

> plot(wn, main="백색잡음계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> acf(wn, main="ACF 상관도표", col="steelblue")


 

>> 표준정규분포를 따르는 백색잡음계열의 상관도표는 위와 같다. 이를 보면 시차 0을 제외하고는 모두 점선 안의 작은 값을 갖는다. 시차의 1,2,3은 1,2,3년을 의미하며, 4개의 시차가 모여 1년이 된다. 백색잡음계열은 예측할 수 없는 계열이다. 시차 0을 제외한 모든 자기상관계수값이 점선 내에 있다. 따라서 백색잡음의 각 시차의 자기상관계수값이 5% 유의수준에서 0과 다르다고 할 수 없다.

 

 

 

> nn=length(gdp_sa)

> sin = ts(sin(1:nn/nn*12*pi), start=1970, frequency=4)

> par(mfrow=c(2,1))

> plot(sin, main="sin 계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> acf(sin, main="ACF 상관도표", col="steelblue")


>> sin 커브를 가지는 시계열의 상관도표를 구해보면 위와 같다. Sin 시계열의 자기상관계수는 sin커브와 같이 일정한 주기를 가지고 움직이는 것으로 나타났다.

 

> plot(gdp_o, main="GDP 원계열", xlab="", ylab="", col="steelblue")

> lines(gdp_sa,col='red')

> acf(gdp_o, main="ACF 상관도표", col="steelblue")

>> GDP원계열의 상관도표를 그려보면 위와 같다. 자기상관계수값이 서서히 감소하는 것으로 나타났다. 동 시계열은 강한 추세가 있는 불안정시계열임을 나타낸다. 강한 추세로 인해 GDP의 계절변동을 상관도표로 파악하기는 어렵다.

 

> plot(diff(log(gdp_o)), main="GDP로그차분", xlab="", ylab="", col="steelblue")

> acf(diff(log(gdp_o)), main="ACF 상관도표", col="steelblue")

>> GDP 원계열을 로그차분한 계열의 상관도표는 위와 같다. 이를 보면 자기상관계수값이 주기적으로 움직이는 것으로 나타났다. 이는 동 시계열의 로그차분을 통해 추세를 제거했지만, 계절변동은 그대로 남아 있음을 의미한다.

 

> plot(diff(log(gdp_o),4), main="GDP로그4차 차분", xlab="", ylab="", col="steelblue")

> acf(diff(log(gdp_o),4), main="ACF 상관도표 ", col="steelblue")

>>GDP 원계열을 로그 4차 차분한 계열의 상관도표를 그려보면 위와 같다. 자기상관계수값이 3차 이후 급격히 하락하는 것으로 나타났다. 이는 동 시계열의 로그 4차 차분을 통해 추세 및 계절 변동이 제거되었음을 의미한다.

 

 

륭(Ljung) – 박스(Box) 검정을 통해, 시계열에 m차까지의 자기상관관계가 존재하지 않는다는 귀무가설을 검정하는 것이다.

 - 귀무가설: H0 : p(1) = p(2) = … = p(m) = 0 잔차는 랜덤

> Box.test(wn,lag=8, type="Ljung")

 

Box-Ljung test

 

data:  wn

X-squared = 3.16, df = 8, p-value = 0.9239

 

> Box.test(sin,lag=8, type="Ljung")

 

Box-Ljung test

 

data:  sin

X-squared = 579.9, df = 8, p-value < 2.2e-16

 

> Box.test(gdp_o,lag=8, type="Ljung")

 

Box-Ljung test

 

data:  gdp_o

X-squared = 1249.8, df = 8, p-value < 2.2e-16

 

> Box.test(diff(log(gdp_o)),lag=8, type="Ljung")

 

Box-Ljung test

 

data:  diff(log(gdp_o))

X-squared = 756.05, df = 8, p-value < 2.2e-16

 

> Box.test(diff(log(gdp_o),4),lag=8, type="Ljung")

 

Box-Ljung test

 

data:  diff(log(gdp_o), 4)

X-squared = 205.63, df = 8, p-value < 2.2e-16

 

계열명

X2 검정통계량 값

유의확률

백색잡음계열

3.16

0.9239

sin계열

579.9

< 2.2e-16

GDP원계열

1249.8

< 2.2e-16

GDP 로그차분계열

756.05

< 2.2e-16

GDP 로그 4차 차분계열

205.63

< 2.2e-16

>> 백색잡음계열의 경우 p-value가 0.05보다 크므로 시계열에 8차까지의 자기상관관계가 존재하지 않는다는 귀무가설을 기각하지 못하지만, 나머지 계열은 모두 작은 값이어서 8차까지의 자기상관관계가 존재하지 않는다는 귀무가설을 기각하낟.

 

 

3.4 부분 자기 상관

표본 부분 자기 상관계수(sample partial autocorrelation coefficietns). 앞뒤 시차 영향을 제외한 순수한 영역

> gdp <- read.csv("gdpq.csv", header = TRUE)

> gdp_o  <- ts(gdp[,1]/1000, start=1970, frequency=4)

> gdp_sa <- ts(gdp[,2]/1000, start=1970, frequency=4)

> gdp_gr <- ts((gdp_sa-lag(gdp_sa,-1))/lag(gdp_sa,-1)*100,  start=c(1970,2), frequency=4)

> par(mfrow=c(2,1))

> plot(wn, main="", xlab="백색잡음계열", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> pacf(wn, main="Partial ACF 부분상관도표", col="steelblue") #main="", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")


>> 모두 점선 안(1.96 / root N)의 작은 값을 가지고 있다.

 

> plot(sin, main="sin 계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> pacf(sin, main="Partial ACF 부분상관도표", col="steelblue")


>> sin 계열의 부분자기상관계수는 앞의 2기에서 유의하게 큰 값을 가지고 있다.

 

> plot(gdp_o, main="GDP 원계열", xlab="", ylab="", col="steelblue")

> lines(gdp_sa,col='red')

> pacf(gdp_o, main="Partial ACF 부분상관도표", col="steelblue")


>> 부분자기상관계수값은 1차에서 1에 가까운 큰 값을 가지고 있다. 1차의 큰 값은 추세가 있음을, 주기적으로 크게 나타난 부분자기상관계수는 계절변동이 있음을 의미한다.

 

 

> plot(diff(log(gdp_o)), main="GDP 로그차분", xlab="", ylab="", col="steelblue")

> pacf(diff(log(gdp_o)), main="Partial ACF 부분상관도표 ", col="steelblue")


>> 부분자기상관계수값이 1차의 값은 작고 주기적으로 움직이는 것으로 나타났다. 이는 동 시계열의 로그차분을 통해 추세를 제거했지만, 계절 변동은 그대로 남아있음을 의미한다.



 

> plot(diff(log(gdp_o),4), main="GDP 로그 4차 차분", xlab="", ylab="", col="steelblue")

> abline(h=0,lty=2,col='grey')

> pacf(diff(log(gdp_o),4), main="Partial ACF 부분상관도표 ", col="steelblue")

>> GDP 원계열을 로그 4차차분한 계열의 부분상관도표를 그려보면 그림과 같다. 이를 살펴보면 부분자기상관계수값은 3차 이후 작은값을 가지는 것으로 나타났다. 이는 동시계열의 로그 4차 차분을 통해 추세 및 계절변동이 제거되었음을 의미한다.

 

> sin12 = ts(sin(1:nn/nn*12*pi), start=1970, frequency=4)

> sin36 = ts(sin(1:nn/nn*36*pi), start=1970, frequency=4)

> plot(cbind(sin12, sin36), main="sin12와 sin36 추이", xlab="", ylab="", col="steelblue")

> spectrum(cbind(sin12, sin36), spans=c(3,3), main="sin12와 sin36 함수 스펙트럼", col="steelblue")

>> 위의 그래프를 보면 첫번째 sin함수는 주기가 상대적으로 긴 함수, 즉 저주파 함수이고, 두번째 sin 함수는 상대적으로 주기가 짧은 함수, 즉 고주파 함수이다. 이를 평활화된 주기도(스펙트럼)을 구해보면 아래 그래프와 같은데 실선이 sin12, 점선이 sin36이다. 해당주파수(주기의 역수)의 값이 크게 나타난다. 저주파 sin12자료의 스펙트럼은 낮은 주파수에서 최고점이 나타나나, 고주파 sin36 자료의 스펙트럼은 높은 주파수에서 최고점이 나타난다.

 

> plot(sin12+sin36, main="sin12+sin36", xlab="", ylab="", col="steelblue")

> spectrum(sin12+sin36, spans=c(3,3), main="sin 12 + sin36 합성 스펙트럼", col="steelblue")

 

>> 두 sin 계열을 합성한 결과 스펙트럼은 다음과 같다. 두 계열 합성계열의 스펙트럼은 두 sin 계열에 해당되는 주파수에서 큰 값을 가지고 있어서 동 계열에 두 주파수에서의 변동요인이 있음을 알 수 있다.

 

> par(mfrow=c(2,1))

> plot(wn, main="백색잡음계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> wn=ts(rnorm(nn), start=1970, frequency=4)

> spectrum(wn, spans=c(3,3), main="스펙트럼", col=1:2)

>> 표준정규분포로부터 154개의 난수를 생성하여 보면 아무 특징이 없는 백색잡음계열이 생성된다. 백색잡음계열의 평활화한 스펙트럼은 그림과 같이 아무 특징이 없는데 y축 값을 보면 매우 작은값을 가지고 수평적으로 움직이고 있다. 백색잡음계열의 평활화한 스펙트럼은 이론적으로는 수평선이다.  다양한 peak가 나타난다는 것은 잔차가 random, 주파수에 특정 주기가 없다는 뜻이다.

 

> par(mfrow=c(2,1))

> plot(gdp_o, main="GDP 원계열", xlab="", ylab="")

> lines(gdp_sa, col=2)

> spectrum(cbind(gdp_o, gdp_sa), spans=c(3,3), main="스펙트럼(원계열+계정조정계열)", col=1:2)

 

>> GDP 원계열과 계절조정계열의 시계열도표와 스펙트럼은 다음과 같다. 원계열의 스펙트럼(실선)을 보면 저주파수와 계절주파수에서 큰 값을 가지는 것으로 나타났다. 여기서 저주파수는 추세변동요인의 존재를, 계절주파수의 큰 값은 계절변동요인의 존재를 의미한다. 반면 GDP 계절조정계열의 스펙트럼(점선)을 구해보면 큰 값이 사라짐을 알 수 있다. 즉 계절변동요인이 제거되었음을 말한다.

 

> dlgdp1 = diff(log(gdp_o))

> dlgdp2 = diff(log(gdp_o),4)

> dlgdp  = cbind(dlgdp1,dlgdp2)

> plot(dlgdp1, main="로그1차차분과 4차차분계", xlab="", ylab="", col="steelblue")

> lines(dlgdp2, col=2)

> spectrum( na.omit(cbind(dlgdp1,dlgdp2)), spans=c(3,3), main="스펙트럼(로그1차차분과 4차차분계열)", col=c("steelblue", "red"))

 

>> 로그1차 차분한 GDP 원계열(보라색)의 스펙트럼을 구해보면 아래 그림과 같다(보라색). 이를 보면 원계열에서의 스텍트럼에서 계절주파수의 큰 값은 남고 저주파수의 추세변동요인과 관련된 큰 값은 사라짐을 알 수 있다. 이는 차분이 시계열에서 추세변동을 제거하는 역할을 했음을 의미한다. 한편 로그 4차 차분 계열의 스펙트럼(점선)을 구해보면 추세변동요인과 계졀변동 요인에 해당되는 주파스의 밀도함수값이 작게 나타난반면, 순환변동에 해당하는 값이 상대적으로 크게 나타나 로그 4차 차분계열에서 순환변동을 파악할 수 있을 것으로 보인다.


출처: 예측방법론(이긍희, 이한식, 장영재 공저)

반응형

'KNOU > 3 예측방법론' 카테고리의 다른 글

제2장 예측데이터 - 시계열  (1) 2017.03.17
Posted by 마르띤
,

제2장 예측데이터 - 시계열

KNOU/3 예측방법론 2017. 3. 17. 19:04
반응형

### 2주차, 제2장 예측데이터 - 시계열 ###

 

> setwd("c:/Rwork/")

> gdp<-read.csv('gdpq.csv',header=T)

> head(gdp)

gdp   gdpsa #계절 조정

1 12807.5 14934.6

2 15732.5 15246.5

3 14621.6 15472.1

4 18689.5 16198.0

5 14512.5 16664.4

6 17709.4 17180.3

> gdp <- ts(read.csv("gdpq.csv", header = TRUE), start=1970, frequency=4)

> gdp

gdp    gdpsa

1970 Q1  12807.5  14934.6

1970 Q2  15732.5  15246.5

1970 Q3  14621.6  15472.1

1970 Q4  18689.5  16198.0

(…중간 생략…)

2013 Q1 265291.2 279223.0

2013 Q2 285312.3 282344.9

2013 Q3 283595.0 285348.5

2013 Q4 300654.7 287979.3

> plot(gdp[,1]/1000, ylab="GDP(조원)", xlab="연도", col="steelblue")

> lines(gdp[,2]/1000, col="red")


 

> ww=c(.5,1,1,1,.5) #5분기 이동평균

> gdpm5 = filter(gdp, sides=2, ww/sum(ww))

> dlgdp1 = diff(log(gdp)) #1차 차분

> dlgdp4 = diff(log(gdp),4) #4차차분

> plot(gdp, main="GDP", ylab=" ", xlab="YEAR") #gdp 원계열


> plot(gdpm5, main="GDPM5", ylab=" ", xlab="YEAR") #gdp m5  이동평균


> plot(dlgdp1, main="diff(log(GDP))", ylab=" ", xlab= "YEAR") #1차 차분


> plot(dlgdp4, main="diff(log(GDP),4)", ylab=" ", xlab= "YEAR") #4차 차분

 

 

 

 출처: 예측방법론(이긍희, 이한식, 장영재 공저, Knou press)

반응형

'KNOU > 3 예측방법론' 카테고리의 다른 글

제3장 예측데이터의 기초분석  (0) 2017.03.22
Posted by 마르띤
,

티스토리툴바