을이용한통계프로그래밍 Jinseog Kim 1 2011 년 3 월 22 일 1 Assistant Professor, Department of Statistics and Information Science, Dongguk University, Gyeongju, Korea. E-mail:jinseog.kim@gmail.com, Homepage: hdatamining.dongguk.ac.kr
제 1 장 머리말 본강의노트는처음 R 을접하는학생들에게 R 이어떻게동작하는지, 그리고어떻게 사용하는지를설명하기위하여만들어졌다. R 에서제공하는여러가지기능을사용 하려면, 처음에몇가지의개념들을이해할필요가있는데, 이강의노트에서는학생들 이이해하기쉽도록그개념들에대하여설명할것이다. R 을통계분석및그래픽스을지원하기위하여 Ross Ihaka 과 Robert Gentleman 이 개발하였고, AT&T Bell Lab 에서만든 S 언어에서유래된것으로소프트웨어이면서 프로그래밍언어이다. S 언어에서유래된또다른소프트웨어로는 Insightful 에서상 업용으로만든만든 S-PLUS 1 가있다. R 과 S 는소프트웨어디자인측면에서몇가지 차별성이있는데, 이에대하여는 Ihaka & Gentleman (1996) 2 의논문혹은 R-FAQ 3 를 참조하기바란다. R 은 GNU General Public 라이센스 4 에따라자유롭게사용과배포가가능하며, 현 재 R 의개발과배포는다수의통계학자및개발자들로구성된 R Development Core Team 이주도하고있다. R 은 C 언어와 Fortran 으로이루어진 source 코드혹은 Win- 1 http://www.insightful.com/products/splus/default.asp 2 Ihaka R. and Gentleman R. 1996. R: a language for data analysis and graphics. Journal of Computational and Graphical Statistics 5: 299 314. 3 http://cran.r-project.org/doc/faq/r-faq.html 4 http://www.gnu.org/ 1
dows, Linux, 그리고 Macintosh용으로만들어진바이너리코드로배포가되고있고, Comprehensive R Archive Network (CRAN) 인터넷사이트인 http://cran.r-project. org/ 에서구할수있다. R은통계분석뿐만아니라그래프에의한자료의표현을위하여많은함수를제공하는데, 이들그래픽함수를통해서얻어진이미지는다양한형태 (jpg, png, bmp, ps, pdf, emf, pictex, xfig) 로변환하여사용할수있다. 통계분석에관련된함수는스크린에분석결과를보여주기도하고결과를오브젝트의형태혹은파일로저장할수있으며, 이와연결된분석에사용할수도있다. R의프로그래밍언어의기능은사용자가 loop같은기능을사용하여여러개의자료를연속적으로분석하거나, 서로다른분석을위해작성한프로그램들을모아서복잡한분석에이용할수있도록한다. S 언어들은거의대부분 R에서동작하므로 R사용자들은인터넷상에서구할수있는많은 S 코드 5 들을구하여 R에서직접수행할수도있다. R을처음접할경우 R을읽히는데다소어렵다고느껴질수도있지만 R의강점은다른통계소프트웨어와달리매우유연한소프트웨어로알려져있다. 한예로 SAS 나 SPSS와같은전통적인통계분석소프트웨어들이분석결과를화면에출력해주는것과달리 R은결과를오브젝트형태로도저장하여주므로화면에출력된결과없이이오브젝트를가지고다른분석을하는데유용하게활용할수있다. R 초보자들은이점때문에간혹당황하는경우가있지만이런특징은실제로매우유용하며, 사용자들이 R수행의결과물중에서관심이있는필요한부분만발췌하여볼수도있게한다. 예를들어, 20개의회귀분석모형을만들어회귀계수를비교해보고자한다면회귀분석의결과들중회귀계수만을선택하여출력하면된다. 다른전통적인분석소프트웨어에서는이러한결과만을출력하는것이복잡하거나어렵다. 본강의노트가지금은 R이지니는다양한기능들을모두설명하고있지는못하지만, 학생들이 R을학습하는데도움이되기를바라면서, 본강의노트는저자의강의가지속되는한계속업데이트될것이다. 5 http://stat.cmu.edu/s/ 2
제 2 장 R 의설치 그림 2.1: R homepage Download the R system from http://www.r-project.org/. Figure 1 의좌측에있는 CRAN 을클릭 Figure 2-(1) 의우측에다운로드를위한 Mirror 사이트중에서하나를선택 3
Figure 2-(2) 의 우측에서 OS 를 선택 (Windows, Linux) Figure 2-(3) 의 우측에서 (Base, contrib) 중 택일, 처음인 경우는 Base 를 선택함 Figure 2-(4) 의 우측에서 가장 최근버전을 다운로드 한다. 통상적으로 윈도우버 젼인 셋업프로그램은 R-버전-win32.exe 라는 명칭으로 되어 있다. 다운로드된 파일을 실행하면 Windows 에 R 이 설치된다. 그림 2.2: R installation 제1절 Rgui R system 의 설치가 끝나면 윈도우즈의 바탕화면에 Rgui 를 실행할 수 있는 아이콘 (Shortcut) 이 생기는데, 이를 클릭하면 R 이 구동이되면서 Figure 2.3 와 같은 윈도가 생긴다. 이를 R-gui 라고 부른다. 4
그림 2.3: Rgui 제 2 절 패키지다운로드및설치 R시스템에는기본적으로설치된 package들이있다. 이들기본패키지를확인하려면 R-gui 에서 search() 명령을실행하면된다. > search() [1] ".GlobalEnv" "package:stats" "package:graphics" [4] "package:grdevices" "package:utils" "package:datasets" [7] "package:methods" "Autoloads" "package:base" 이들패키지들이설치된 path를확인하는명령 ( 함수 ) 1 는아래와같다. > searchpaths() [1] ".GlobalEnv" 1 명령과함수를혼용하여사용하는데이후로는함수라고부르겠다 5
[2] "C:/PROGRA~1/R/R-25~1.1/library/stats" [3] "C:/PROGRA~1/R/R-25~1.1/library/graphics" [4] "C:/PROGRA~1/R/R-25~1.1/library/grDevices" [5] "C:/PROGRA~1/R/R-25~1.1/library/utils" [6] "C:/PROGRA~1/R/R-25~1.1/library/datasets" [7] "C:/PROGRA~1/R/R-25~1.1/library/methods" [8] "Autoloads" [9] "C:/PROGRA~1/R/R-25~1.1/library/base" R 시스템이 SAS, SPSS 혹은 Splus와같은상용 (commercial) 통계소프트웨어에비하여지니고있는장점중의하나는자료의분석에필요한모듈이매우다양하다는점이다. 이는 R이공개소프트웨어이면서사용자가자신의모듈을쉽게시스템에접목할수있도록하는유연한 API 및제작도구를제공하고있기때문이다. 특정사용자가자체제작한모듈을 R 홈페이지에업로드하면이는모든사용자에게공개되어이를필요로하는일반사용자가다운로드하여사용할수있다. 특정사용자가자체제작한모듈을패키지라고부르며, R-gui는이를다운로드및설치를도와주는도구를제공한다. R-gui의메뉴바의 Packages를선택하고풀다운메뉴에서 Install package를선택하면다운로드사이트 (Figure 2.4 좌측 ) 및다운로드할수있는패키지들이나열된리스트창 (Figure 2.4 우측 ) 이나온다. 이리스트창에서원하는패키지를다운로드할수있다. 이리스트창에는단지패키지이름만나열되어있기때문에패키지의용도는스스로알아봐야한다. 패키지의용도를알아내는방법으로는 R홈페이지에서좌측윈도의 Package를클릭하면 Package list와그용도에대한설명이간략하게나와있으므로참고하기바란다. 아래는품질관리의관리도를그리기위한패키지인 qcc를다운로드하여설치를완료하였을때 R-gui에나타나는메시지이다. 네트워크의사정으로인해서간혹다운로드혹은설치에실패하는경우가있으므로이메시지를반드시확인하기바란다. 6
그림 2.4: Install the R contributed package > utils:::menuinstallpkgs() Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/rwin/bin /windows/contrib/2.5 trying URL http://bibs.snu.ac.kr/r/bin/windows/contrib/2.5/qcc_1.2.zip Content type application/zip length 318673 bytes opened URL downloaded 311Kb package qcc successfully unpacked and MD5 sums checked The downloaded packages are in C:\Documents and Settings\star\Local Settings\Temp\RtmpzkpPsb\downloaded_packages 7
updating HTML package descriptions 8
제 3 장 기본유틸리티 R은사용자의편의를위해서통계분석이외에아래와다양한유틸리티를제공하고있다. 아래는많이이용되는유틸리티에대한설명이다. ls(): 작업공간에있는객체 (objects) 를보여준다. 여기서객체는상수, 데이터프레임, 벡터, 행렬, 함수등 R에서정의된객체들이다. > ls() [1] "BostonHousing" "BostonHousing2" "cvt1" [4] "pred" "resid" "RSS" [7] "spam.p" "spam.t" "t1" [10] "t1.1" rm(): 작업공간에있는객체중에서 () 안에지정한특정객체 (objects) 를삭제한 다. > rm(bostonhousing2, cvt1, t1.1, spam.p) > ls() [1] "BostonHousing" "pred" "resid" "RSS" [5] "spam.t" "t1" 9
help(),?: () 안에지정한객체 ( 함수, 데이터 ) 의설명을아래 Figure 3.1 처럼보 여준다. help() 대신? 를사용할수있다. 그림 3.1: Windows help library(): () 안에지정한 package를로드한다. scan(): () 안에지정한외부파일을벡터형태로로드한다. read.table(), read.csv(): () 안에지정한외부파일을데이터프레임형태로로드한다. write.table(), write.csv(): () 안에지정한객체를외부파일로저장한다. save(), save.image(): () 안에지정한객체를 R data형식으로저장한다. x <- runif(20) y <- list(a = 1, b = TRUE, c = "oops") save(x, y, file = "xy.rdata") save.image() load(), unlink(): () 안에지정한 R data 파일을작업환경으로로드하거나, 메 모리에로드했던 Rdata 를 unload 한다. 10
load("all.rdata",.globalenv) unlink("all.rdata") sink(): R 코드를실행할때나타나는화면의출력내용을파일로저장려고할 때사용한다. 이경우화면에는출력되지않는다. 파일로저장되던내용을다시 화면으로출력하게하려면 sink() 라고하면된다. zz <- file("all.rout", open="wt") sink(zz) unlink(zz) sink() getwd(): 현재의작업디렉토리를출력해준다. > getwd() [1] "C:/Documents and Settings/jskim/ 바탕화면 / 강의 /statistics seminar(sna)" setwd(dir): 작업디렉토리를 dir 로바꾼다. > setwd("c:/work") > getwd() [1] "c:/work" setwd(dir): 타이틀바에현재작업디렉토리를나타내기. setwindowtitle(paste(":", getwd())) setstatusbar("error occurs") 11
제 4 장 객체 (R objects) R objects 에는아래와같은종류들이있고, 기본적으로 R 은 OOP(Object Oriented Programming) 개념에의해설계되었기때문에 R에서사용되는대부분이 object로간주될수있다. atomic ( 상수 ) vector matrix list data.frame function operator... 위의 R 오브젝트의분류중에서 atomic, vector, matrix, data.frame을데이터오브젝트라고부른다. 12
제 1 절 Data object 의 storage mode R에서모든데이터오브젝트는 4가지종류의저장타입 (mode or storage mode) 있다.. 아래는데이터오브젝트중 atomic 오브젝트의 4가지 mode에대한예이다. 이러한 4 종류의 mode는 vector, matrix, array 등에공통적으로적용된다. 숫자형 (numeric) > value <- 605 > value [1] 605 문자형 (character) > string <- "Hello World" > string [1] "Hello World" 논리형 (logical) > 2 < 4 [1] TRUE 복소수형 (complex number) > cn <- 2 + 3i > cn [1] 2+3i R 에서는데이터의 mode 를확인하는함수로 mode() 가있다. 13
> mode(value) [1] "numeric" > mode(string) [1] "character" > mode(2<4) [1] "logical" 데이터오브젝트는 class라는속성을부여할수있는데, 만일데이터오브젝트의 class 속성을따로부여하지않을경우 class속성은데이터의 mode와동일하다. > string <- "Hello World" > mode(string) [1] "character" > class(string) [1] "character" > class(string)<-"myclass" ## string의 class를 "myclass" 로부여 > class(string) [1] "myclass" 제 2 절 vector 벡터는하나이상의원소 (atomic) 로이루어진자료구조이다. 또한벡터를구성하고있는원소 (atomic) 는그형태 (mode) 가동일해야한다. 즉, (1, 2, a, b ) 와같은벡터는있을수가없다. 아래는벡터를생성하는코드이다. x1<-c(1,2,3,4) # 함수 c() 의이용 x2<-1:3 x3<-c("a", "B", "C") y <- c(x1, 0, x2) # 1,2,3,4, 0, 1,2,3 14
벡터의길이 (length) 를나타내는함수는 length() 가있다. length(x1) [1] 4 length(x3) [1] 3 length(x3) [1] 4 벡터를생성하는함수로는아래와같은것들이있다. rep(), seq(): # 2를 2번반복 > rep(2, 2) [1] 2 2 # vector (1,2) 의원소를각각 2번반복 > rep(c(1,2), each=2) [1] 1 1 2 2 #0과 1을등구간인 11개숫자로이루어진벡터를만듬 seq(0, 1, length=11) [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 #1에서 9까지 2씩증가하는숫자로이루어진벡터를만듬 seq(1, 9, by = 2) [1] 1 3 5 7 9 numeric(), double(), integer(), character(): 속성이 numeric, double, integer, 혹은 character 인벡터를괄호안의수만큼할당함 > x<-integer(length = 10) 15
> x [1] 0 0 0 0 0 0 0 0 0 0 R 데이터오브젝트중 numeric vector 는다음과같은 class 속성을가질수있다. 이는자료의속성을연속형, 범주형변수로구분하는것과같은맥락으로이해하면 된다. numeric: 연속형 factor: 범주형 ordered: 순서있는범주형 아래의표는모두 numeric vector 인데도불구하고 class 속성이서로다른경우에대한 예시이다. 표 4.1: 벡터오브젝트의 mode 및 class R code mode(x) class(x) x < c(1:10) numeric numeric x< factor(1:10) numeric factor x< ordered(1:10) numeric ordered factor R 데이터오브젝트의 component 를접근하는방법은아래와같이인덱스와 component 이름을이용한다. object [ arg1,..., argn ] object [[ arg1,..., argn ]] object $ tag # for vector, matrix, array # for list # for data.frame or named list 연습문제 16
x <- seq(length=51, from=-5, by=.2) y <- rep(x, times=5) z <- rep(x, each=5) 제 3 절 array( 배열 ) 과 matrix( 행렬 ) > Z <- array(data_vector, dim_vector) 예를들어만약 1부터 24까지의수를 3차원배열에넣는방법은아래와같다. > h<-1:24 > Z <- array(h, dim=c(3,4,2)) > Z,, 1 [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12,, 2 [,1] [,2] [,3] [,4] [1,] 13 16 19 22 [2,] 14 17 20 23 [3,] 15 18 21 24 행렬은 2차원배열이므로위와같은방법으로만들수있다. > x <- array(1:20, dim=c(4,5)) [,1] [,2] [,3] [,4] [,5] 17
[1,] 1 5 9 13 17 [2,] 2 6 10 14 18 [3,] 3 7 11 15 19 [4,] 4 8 12 16 20 또는간단히 matrix() 함수를이용해서아래와같이만들수도있다. X1 <- matrix(1:20, nrow=4, ncol=5) 대각행렬은 diag() 함수를이용하면좀더편하게만들수도있다. X2 <- diag(1, 10) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 0 0 0 0 0 0 0 0 0 [2,] 0 1 0 0 0 0 0 0 0 0 [3,] 0 0 1 0 0 0 0 0 0 0 [4,] 0 0 0 1 0 0 0 0 0 0 [5,] 0 0 0 0 1 0 0 0 0 0 [6,] 0 0 0 0 0 1 0 0 0 0 [7,] 0 0 0 0 0 0 1 0 0 0 [8,] 0 0 0 0 0 0 0 1 0 0 [9,] 0 0 0 0 0 0 0 0 1 0 [10,] 0 0 0 0 0 0 0 0 0 1 X2 <- diag(10) X2 <- diag(1:10) X2 <- diag(c(1,3,5,7,9)) [,1] [,2] [,3] [,4] [,5] [1,] 1 0 0 0 0 [2,] 0 3 0 0 0 18
[3,] 0 0 5 0 0 [4,] 0 0 0 7 0 [5,] 0 0 0 0 9 3.0.1 행렬의연산 element by element product (A 와 B 가크기가같은행렬일때 ) A * B matrix product( 행렬곱셈 ; A 의열수와 B 의행수가같은경우 ) A %*% B matrix inversion( 역행렬 ; A 가정방행렬인경우 ) A=matrix(c(1,2,3,4), ncol=2) solve(a) [,1] [,2] [1,] -2 1.5 [2,] 1-0.5 제 4 절 리스트 (list) List 는그구성요소 (component) 이또다른오브젝트들로이루어진오브젝트이다. 숫 자벡터, 논리값, 행렬, 문자, 배열, 함수등모든 R 오브젝트가리스트의 component 가 될수있다. 19
4.0.2 리스트의구성과수정리스트는 list() 함수를이용하여만들수있다. list(name_1=object_1,..., name_m=object_m) 여기서 name_1... name_m 은 [ 이름인수 ] 로리스트의콤포넌트의이름을지정한다. 또한 object_1... 은콤포넌트를위한오브젝트이다. 리스트를만드는간단한예를살펴보자. Lst <- list(name="fred", wife="mary", no.children=3, child.ages=c(4,7,9)) $name [1] "fred" $wife [1] "mary" $no.children [1] 3 $child.ages [1] 4 7 9 만약리스트가 4개의구성요소에대한접근방법은아래와같다. > Lst[[1]] [1] "fred" > Lst[[4]] [1] 4 7 9 리스트의구성요소에이름이붙여져있는경우, 그이름을이용하여접근할수있다. > Lst[["name"]] [1] "fred" > Lst$name 20
참고로 Lst[1], Lst[2] 는리스트에서첫번째혹은두번째 component 만으로이루어진서브리스트 (sublist) 를표현한다. > Lst[1] $name [1] "fred" > Lst[3:4] $no.children [1] 3 $child.ages [1] 4 7 9 리스트의콤포넌트의개수는 length() 함수를이용하여구할수있다. > length(lst) => 4 4.0.3 리스트의결합여러개의리스트는 c() 함수를이용하여연결할수있다. 이를통하여각리스트의구성요소들이차례로결합된형태의새로운리스트가생성된다. list1<-list(a1=1, b1=1:3) list2<-list(a2=2, b2=5:10) list3<-list(a3=3, b3=100:90) list.all <- c(list1, list2, list3) > list.all $a1 [1] 1 $b1 [1] 1 2 3 21
$a2 [1] 2 $b2 [1] 5 6 7 8 9 10 $a3 [1] 3 $b3 [1] 100 99 98 97 96 95 94 93 92 91 90 제 5 절 데이터프레임 (data frame) 데이터프레임은 data.frame 으로분류되는특별한리스트로서, 리스트에아래의제약을주어만든다. 데이터프레임의구성요소는반드시벡터, 펙터 (factor), 행렬, 리스트또는다른데이터프레임이어야한다. ( 이경우모두데이터프레임이가능하다는말은아님 ) 행렬, 리스트그리고데이터프레임이가진각각의행, 구성요소또는변수는새로운데이터프레임의행, 구성요소또는변수가된다. 벡터 ( 숫자, 문자등 ) 는데이터프레임의하나의열이된다. 데이터프레임에포함된변수는그길이가모두동일해야한다. 5.0.4 데이터프레임만들기데이터프레임을만들때는 data.frame() 함수를이용한다. name <-c("kim","lee","park","oh") sex<-c( f, m, f, m ) 22
income<-c(100,102,300,204) d1 <- data.frame(name=name, gender=sex, incom=income) > d1 name gender incom 1 kim f 100 2 lee m 102 3 park f 300 4 Oh m 204 리스트나행렬을통째로데이터프레임으로바꾸기위해서는 as.data.frame() 함수를 사용할수있다. 5.0.5 데이터프레임관련함수 앞줄보기 head(data) 변수명출력 names(data) 데이터차원출력 nrow(data) # 행 ncol(data) # 열 dim(data) ## 행, 열의차원 (dimension) 23
제 5 장 파일에서데이터읽어오기 제 1 절 read.table() 함수 대용량의데이터는 R 콘솔에서키보드를통해읽어들이는것보다는외부파일을통 해서읽어들이는것이편리하다. read.table() 함수는외부파일을읽어데이터프레임을 만들어주는함수이다. 이때외부파일을다음의형식을만족하도록작성할것을권 장한다. 파일의첫번째줄은변수명을지정한다. 다음줄부터는행식별자 ( 행라벨, row labels) 를입력하고그이후에는관측값을 변수명에대응하는순서대로입력한다. 아래는위의형식에의하여작성된외부파일 (houses.data) 의예이다. Price Floor Area Rooms Age Cent.heat 01 52.00 111.0 830 5 6.2 no 02 54.75 128.0 710 5 7.5 no 03 57.50 101.0 1000 5 4.2 no 04 57.50 131.0 690 6 8.8 no 24
05 59.75 93.0 900 5 1.9 yes... 위의예제데이터를불러들여데이터프레임 (HousePrice) 으로만드는 R code는아래와 같다. HousePrice<-read.table("houses.data") 행라벨입력을생략하고싶은경우, 즉외부파일이다음과같이입력된경우어떻 게이를일어들일수있을까? Price Floor Area Rooms Age Cent.heat 52.00 111.0 830 5 6.2 no 54.75 128.0 710 5 7.5 no 57.50 101.0 1000 5 4.2 no 57.50 131.0 690 6 8.8 no 59.75 93.0 900 5 1.9 yes... 위의경우에는다음과같은코드를이용하면읽어들일수있다. HousePrice <- read.table("houses.data", header=true) 제 2 절 read.csv() 함수 원본데이터를 Excel 파일로편집하는경우다음과같은방식으로 R data.frame 으로 불러들일수있다. as a CSV file (Comma Separated Values). CSV(Comma Separated Values) 형식으로저장 : 파일 Save As: 아래의함수를이용한다. 25
my.table=read.csv(file.choose()) ## 대화상자이용 my.table=read.csv("c:/xfile.csv") ## 파일명직접입력 제 3 절 Excel 에서직접읽어오기 Excel 파일을직접읽어데이터프레임으로만들기위해서는 RODBC package를설치하여야한다. library(rodbc); # filename<-paste(path, filename, ".xls", sep=""); channel <- odbcconnectexcel(file.choose()); x<-sqlfetch(channel, "Sheet1"); odbcclose(channel); 제 4 절 scan() 함수 아래와같이텍스트데이터파일 ( input.dat ) 입력되어있다고하자. 52.00 54.75 57.50 57.50 59.75 111.0 128.0 101.0 131.0 93.0 이러한파일을읽기위해서 scan() 함수를이용한다. inp <- scan("input.dat") 제 5 절 edit() 함수 * 기존의데이터 (olddata) 를수정할때의예 : 26
newdata<-edit(olddata) 새로운데이터를편집할때 : xnew<-edit(data.frame()) 27
제 6 장 R 오퍼레이터및내장함수 (built-in functions) 제 1 절 R operators - :Minus, can be unary or binary + :Plus, can be unary or binary * :Multiplication, binary / :Division, binary %% :Modulus, binary %/% :Integer division, binary < :Less than, binary > :Greater than, binary == :Equal to, binary >= :Greater than or equal to, binary <= :Less than or equal to, binary! :Unary not 28
: :Sequence, binary (in model formulae: interaction) ^ :Exponentiation, binary & :And, binary, vectorized && :And, binary, not vectorized :Or, binary, vectorized :Or, binary, not vectorized <- :Left assignment, binary, <<- :global assignment -> :Right assignment, binary = :Left ssignment, binary $ :List subset, binary 제 2 절 R built-in functions 제곱근지수함수로그함수최대값최소값합평균절대값누적연산삼각함수올림, 반올림... sqrt exp log(5), log2(5), log10(5), log(5, base=3) max, pmax min, pmin sum mean abs cummax, cummin, cumprod, cumsum sin, cos, tan ceiling, round, trunc, floor > sqrt(a) [1] 1.000000 1.414214 1.732051 2.000000 29
> exp(a) [1] 2.718282 7.389056 20.085537 54.598150 > log(a) [1] 0.0000000 0.6931472 1.0986123 1.3862944 > c <- (a + sqrt(a))/(exp(2)+1) > c [1] 0.2384058 0.4069842 0.5640743 0.7152175 > x1 <- seq(-2, 4, by =.5) > x1 [1] -2.0-1.5-1.0-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 > floor(x1) [1] -2-2 -1-1 0 0 1 1 2 2 3 3 4 > trunc(x1) [1] -2-1 -1 0 0 0 1 1 2 2 3 3 4 > a <- c(1,-2,3,-4) > b <- c(-1,2,-3,4) > min(a,b) [1] -4 > pmin(a,b) [1] -1-2 -3-4 제 3 절 Other built-in functions print(): Prints a single R object 30
> a<-c(5,3,6,2,4) > print(a) [1] 5 3 6 2 4 cat(): Prints multiple objects, one after the other > cat("mean of a is ",mean(a), "variance of a is ", var(a),"\n") mean of a is 4 variance of a is 2.5 unique():gives the vector of distinct values x<-c(1,5,1,3,5,7,5) unique(x) [1] 1 5 3 7 diff(): Replace a vector by the vector of first differences > diff(x) [1] 4-4 2 2 2-2 sort(): Sort elements into order, but omitting NAs order(): x[order(x)] orders elements of x, with NAs last rev(): reverse the order of vector elements > print(x) [1] 1 5 1 3 5 7 5 > sort(x) [1] 1 1 3 5 5 5 7 31
> order(x) [1] 1 3 4 2 5 7 6 > rev(x) [1] 5 7 5 3 1 5 1 32
제 7 장 R 프로그래밍 제 1 절 조건문 조건문에해당되는표현으로는다음의 3 가지종류가있다. if ( cond ) expr if ( cond ) expr1 else expr2 if ( cond1 ) expr1 else if( cond2 ) expr2 else expr3 예 : iris 데이터에서 Sepal.Length의 median을구하고 Sepal.Length를 median보다크면 L, 작으면 S 가되도록하여라. data(iris) n=length(iris$sepal.length) Sepal.Length.Cat = character(n) 33
Med=median(iris$Sepal.Length) for(i in 1:n) if(iris$sepal.length<med) Sepal.Length.Cat[i] = "S" else Sepal.Length.Cat[i] = "L"; 제 2 절 순환문 순환문의표현으로는다음세가지표현을사용한다. while ( cond ) expr repeat expr for ( var in list ) expr 위에서 while, repeat, for 구문안에서순환문을끝내기위한구문으로 break와이후의문장을건너뛰고다음순환을하라는구문인 next를사용할수있다. 이는 C 언어에서 break, 그리고 continue 와동일한용도로사용된다. 아래의 R 프로그램은정규분포에서난수 100개를발생시키고이중양수인값들만을합하는프로그램이다. x<-rnorm(100) ## 표준정규분포에서 100개의난수발생 sum.positive <- 0 for(i in 1:100) { if(x[i] < 0) next; sum.positive <- sum.positive + x[i] } 제 3 절 함수의작성 함수작성방법 34
function_name <- function(arg_1, arg_2,...){ expression...; return(...) } ( 예제 )mile 을 km 로바꾸는프로그램. miles.to.km <- function(miles) miles*8/5 > miles.to.km(175) # Approximate distance to Sydney, in miles [1] 280 만일 100, 200 300 miles를 kilometer로바꾼다면 > miles.to.km(c(100,200,300)) [1] 160 320 480 제 4 절 R 에서 C code 부르기 * We will explore using.c and.call with 7 code examples: Using.C /* usec1.c */ /* Calling C with an integer vector using.c */ void usec(int *i) { } i[0] = 11; The C function should be of type void. The compiled code should not return anything except through its arguments. To compile the c code, type at the command prompt: 35
R CMD SHLIB usec1.c The compiled code file name is usec1.so In R: > dyn.load("usec1.so") > a <- 1:10 # integer vector > a [1] 1 2 3 4 5 6 7 8 9 10 > out <-.C("useC", b = as.integer(a)) > a [1] 1 2 3 4 5 6 7 8 9 10 > out$b [1] 11 2 3 4 5 6 7 8 9 10 You have to allocate memory to the vectors passed to.c in R by creating vectors of the right length. The first argument to.c is a character string of the C function name. The rest of the arguments are R objects to be passed to the C function. All arguments should be coerced to the correct R storage mode to prevent mismatching of types that can lead to errors..c returns a list object. The second.c argument is given the name b. This name is used for the respective component in the returned list object (but not passed to the compiled code). II. 36
/* usec2.c */ /* Calling C with different vector types using.c */ void usec(int *i, double *d, char **c, int *l) { i[0] = 11; d[0] = 2.333; c[1] = "g"; l[0] = 0; } To compile the c code, type at the command prompt: RCMD SHLIB usec2.c to get usec2.so To compile more than one c file: RCMD SHLIB file1.c file2.c file3.c to get file1.so In R: > dyn.load("usec2.so") > i <- 1:10 # integer vector > d <- seq(length=3, from=1, to=2) # real number vector > c <- c("a", "b", "c") # string vector > l <- c("true", "FALSE") # logical vector > i [1] 1 2 3 4 5 6 7 8 9 10 > d 37
[1] 1.0 1.5 2.0 > c [1] "a" "b" "c" > l [1] "TRUE" "FALSE" > > out <-.C("useC", i1 = as.integer(a), d1 = as.numeric(d), c1 = as.character(c), l1 = as.logical(l)) > out $i1 [1] 11 2 3 4 5 6 7 8 9 10 $d1 [1] 2.333 1.500 2.000 $c1 [1] "a" "g" "c" $l1 [1] FALSE FALSE 38
제 8 장 R 그래픽스 다른통계분석도구들에비해서 R 이가진장점중의하나가뛰어난그래픽기능이다. 그래프를위한 R 기본함수는아래와같다. 고수준함수 : plot(), barplot(), boxplot(), hist(), pie(), persp() 저수준함수 점그리기 : points() 선그리기 : lines(), abline(), arrows() 문자출력 : text() 도형 : rect(), ploygon() 좌표축 : axis() 격자표현 : grid() x<-rnorm(100, sd=2) y<-0.3 + 2*x + rnorm(100, sd=1) plot(x) 39
그림 8.1: plot() example 제 1 절 barplot() 과 pie plot pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12) names(pie.sales) <- c("blueberry", "Cherry", "Apple", "Boston Cream", "Other", "Vanilla Cream") barplot(pie.sales) pie(pie.sales) 제 2 절 par() 함수를이용한그래프나누기 x<-rnorm(100) par(mfrow=c(1,2)) hist(x) plot(x) 40
0.00 0.05 0.10 0.15 0.20 0.25 0.30 Apple Cherry Blueberry Vanilla Cream Other Boston Cream Blueberry Cherry Apple Boston Cream Other Vanilla Cream 그림 8.2: barplot() and pie() examples Histogram of x Frequency 0 5 10 15 20 25 x 2 1 0 1 2 2 1 0 1 2 x 0 20 40 60 80 100 Index 그림 8.3: par() 함수이용 41
제 9 장 R 을이용한기초통계분석 제 1 절 확률분포 (Probability distributions) R에서확률분포와관련된계산을위해서몇가지의함수를제공하고있다. 확률분포와관련된기능으로는누적분포함수 (P (X x)), 확률밀도함수, quantile 함수 (min x {P (X x) > q}), 그리고랜덤넘버의생성이있다. 아래의표 (Table 9.1) 는잘알려진확률분포에대하여 R함수의기본형과입력값들이다. 확률분포와관련된네가지기능을구현하기위해서는위표에나온분포에따른기본형앞에다음과같은접두사를붙이면된다. p : 누적분포함수 (P (X x)) d : 확률밀도함수 (density) q : quantile 함수 (min x {P (X x) > q}) r : 랜덤넘버의생성아래의예제프로그램은 X t(13) 을따를때, P ( X > 2.43) 을구하는것이다. > ## 2-tailed p-value for t distribution 42
표 9.1: 확률분포에따른 R 함수의기본형확률분포 R 함수기본형 additional arguments beta beta shape1, shape2, ncp binomial binom size, prob chi-squared chisq df, ncp exponential exp rate F f df1, df2, ncp gamma gamma shape, scale geometric geom prob hypergeometric hyper m, n, k logistic logis location, scale negative binomial nbinom size, prob normal norm mean, sd Poisson pois lambda Student s t t df, ncp uniform unif min, max > 2*pt(-2.43, df = 13) [1] 0.0303309 제 2 절 기초통계 ( 기술통계 ; Descriptive statistics) 일변량자료를분석할때, 맨처음자료의수, 평균, 분산등요약통계량 (summary statistics) 을계산하는것부터시작한다. 이렇게함으로써자료의특징을개략적으로짐작할수있고, 이를바탕으로더세부적인분석을가능하게하기때문이다. 아래는요약통계를구하는함수 summary() 를이용하여기초분석을하는 R코드이다. > attach(faithful) > summary(eruptions) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.600 2.163 4.000 3.488 4.454 5.100 43
> fivenum(eruptions) [1] 1.6000 2.1585 4.0000 4.4585 5.1000 2.1 줄기 - 잎그림과상자그림 > stem(eruptions) The decimal point is 1 digit(s) to the left of the 16 070355555588 18 000022233333335577777777888822335777888 20 00002223378800035778 22 0002335578023578 24 00228 26 23 28 080 30 7 32 2337 34 250077 36 0000823577 38 2333335582225577 40 0000003357788888002233555577778 42 03335555778800233333555577778 44 02222335557780000000023333357778888 46 0000233357700000023578 48 00000022335800333 50 0370 44
2.2 Histogram A stem-and-leaf plot is like a histogram, and R has a function hist to plot histograms. par(mfrow=c(1,2)) hist(log10(islands)) hist(log10(islands), nclass=15, prob=true) lines(density(log10(islands), bw="sj"), col="blue") rug(log10(islands)) # show the actual data points 그림 9.1: Histogram example The bandwidth bw was chosen by trial-and-error as the default gives too much smoothing. (Better automated methods of bandwidth choice are available, and in this example bw = SJ gives a good result.) 2.3 분포가정의진단 : 정규성검정 표본자료의경험분포함수 (Empirical CDF) 는 ecdf() 함수를이용해서구한다. 45
x<-rnorm(200) x2<-runif(200) par(mfrow=c(1,2)) plot(ecdf(x), main="edf of x~norm(0,1)") plot(ecdf(x2), main="edf of x2~uniform(0,1)") 그림 9.2: 경험분포함수 EDF of x~norm(0,1) EDF of x2~uniform(0,1) Fn(x) 0.0 0.2 0.4 0.6 0.8 1.0 Fn(x) 0.0 0.2 0.4 0.6 0.8 1.0 3 2 1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 x x Figure 9.3의좌측그림을보면경험분포가표준적인분포와크게차이가남을알수있다. 따라서자료의값이 3 minutes보다큰것만을골라다시경험분포함수를그리고그위에정규분포의 CDF를도시하였다 (Figure 9.3의우측그림 ). 거의정규분포와비슷하다고보여진다. long <- eruptions[eruptions > 3] plot(ecdf(long), do.points=false, verticals=true) x <- seq(3, 5.4, 0.01) lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3) 정규성을체크하기위한좀더세밀한시각적인방법은 Quantile-quantile (Q-Q) plot 이다. 46
# creating simulated data y <- 4*x+rnorm(100) par(mfrow=c(1,2)) qqnorm(y,main="q-q plot of N(0,1) v.s. Sample(y)") qqline(y) qqplot(qt(ppoints(250), df = 3), y, xlab = "t(3) Quantiles", main="q-q plot of Student-t v.s Sample(y)") qqline(y) 그림 9.3: Q-Q plot Q Q plot of N(0,1) v.s. Sample(y) Q Q plot of Student t v.s Sample(y) Sample Quantiles 10 5 0 5 10 y 10 5 0 5 10 3 2 1 0 1 2 3 5 0 5 Theoretical Quantiles t(3) Quantiles 정규성에대한검토는위와같이그림으로나타내는것이외에정규성검정을통해서도할수있다. R은일변량자료의정규성점정을위하여 Shapiro-Wilk test, Kolmogorov- Smirnov test 등을제공하고있다. > shapiro.test(long) Shapiro-Wilk normality test data: long W = 0.9793, p-value = 0.01052 47
> ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long))) One-sample Kolmogorov-Smirnov test data: long D = 0.0661, p-value = 0.4284 alternative hypothesis: two.sided 제 3 절 제 4 절 추정및신뢰구간구하기 two-sample tests 두모집단의평균이같은지를검정하는전통적인방법이이표본 t-test이다. 아래의자료는얼음혼합물의잠재열량 (cal/gm) 을측정한자료이다 (Rice;1995, p.490). Method A: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97 80.05 80.03 80.02 80.00 80.02 Method B: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97 먼저두표본의분포를시각적으로비교하기위하여상자그림을그려본다. A <- c( 79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04, 79.97, 80.05, 80.03, 80.02, 80.00, 80.02) B <- c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95, 79.97) par(mfrow=c(1,2)) boxplot(list(a=a, B=B)) boxplot(list(a=a, B=B), horizontal = TRUE) 48
그림 9.4: Box plot 79.94 79.98 80.02 A B A B 79.94 79.98 80.02 위의상자그림을보면첫번째그룹 (A) 에서의열량이두번째그룹 (B) 보다높다고보여진다. 이를두모집단간에평균차가있는지에대한가설검정을통해서확인해보자. R에서 t-test를위한함수는 t.test() 이다. > t.test(a, B) Welch Two Sample t-test data: A and B t = 3.2499, df = 12.027, p-value = 0.00694 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.01385526 0.07018320 sample estimates: mean of x mean of y 80.02077 79.97875 49
R에서의 t 검정은기본적으로두집단의분산이서로다르다는전제에서수행한다. 두모집단의분산이같은지혹은다른지에대한검정을먼저해야한다. 이경우의검정방법은 F 검정이며, 이를위한 R 함수로는 var.test() 이있다. > var.test(a, B) F test to compare two variances data: A and B F = 0.5837, num df = 12, denom df = 7, p-value = 0.3938 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.1251097 2.1052687 sample estimates: ratio of variances 0.5837405 F-검정결과두모집단의분산이다르다고할만한근거는없다고판단된다. 따라서분산이동일한경우에대한 t-test를다음과같이한다. > t.test(a, B, var.equal=true) Two Sample t-test data: A and B t = 3.4722, df = 19, p-value = 0.002551 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.01669058 0.06734788 50
sample estimates: mean of x mean of y 80.02077 79.97875 위의모든 t 검정은기본적으로두표본이정규분포를따른다는가정에서출발한다. 따라서정규성을만족하는지를체크하여야한다. 정규성에대한진단은앞소절 (2.3) 에서처럼시각적인방법혹은가설검정을통해서할수있다. > plot(ecdf(a), do.points=false, verticals=true, xlim=range(a, B)) > plot(ecdf(b), do.points=false, verticals=true, add=true) 다음은두분포가서로동일한지를검정하는 Kolmogorov-Smirnov test 이다. > ks.test(a, B) Two-sample Kolmogorov-Smirnov test data: A and B D = 0.5962, p-value = 0.05919 alternative hypothesis: two-sided Warning message: cannot compute correct p-values with ties in: ks.test(a, B) 제 5 절 대응비교 (paired t-test) 이변량자료에대한대응비교 (paired t-test) 역시 t.test() 함수를이용한다. 다만, R 함수 t.test() 에는아래와같은많은인수들이있으며, 이함수의인수중에 paired=true 로지정하기만하면된다. t.test(x, y = NULL, 51
alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95,...) 제 6 절 비율검정및카이제곱검정 smokers <- c( 83, 90, 129, 70 ) patients <- c( 86, 93, 136, 82 ) prop.test(smokers, patients) ## Effect of simulating p-values x <- matrix(c(12, 5, 7, 7), nc = 2) chisq.test(x)$p.value # 0.4233 chisq.test(x, simulate.p.value = TRUE, B = 10000)$p.value # around 0.29! 52