np_chA_R.Rmd

---
title: "기초적인 R 사용법"
subtitle: "비모수통계학 부록"
author: "정성규"
date: '마지막 업데이트: `r Sys.Date()`'
output:
  html_document: default
  pdf_document:
    latex_engine: xelatex
    keep_tex: true
mainfont: NanumMyeongjo
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
 

이 책은 기초적인 `R` 사용에 익숙한 독자를 가정하고 쓰였다. `R`을 처음 접하는 독자는 
`R`의 기본적인 사용법을 배우기 위한 가장 표준화된 참고도서인  
An Introduction to R (W. N. Venables, D. M. Smith
and the R Core Team 
https://cran.r-project.org/doc/manuals/r-release/R-intro.pdf )를 참고하도록 하자. `R`은 인터넷에서 참고자료를 많이 찾을 수 있다.  (한국어로 된 `R` 기초서도 여럿 있지만, 저자가 직접 살펴보지 않아 따로 추천하지 않는다.) 
이 절에서는 이 책에서 다루는 통계량의 계산, 재표집, 그리고 시각화를 위한 `R` 사용법을 설명한다. 

## 데이터 형태와 연산

`R`에 데이터를 입력하는 가장 간단한 방법은 벡터를 이용하는 방법이다. 예를 들어 제 3장 예 3.1의 회복기간 19, 22, 24,26의 데이터는 간단히 `c(19,22,24,26)`과 같이 입력할 수 있다. 이 때, `c()`는 함수로서 입력값들 (여기서는 `19,22,24,26`)을 결합("combine")하여 벡터의 형태로 만들어준다. 이렇게 생성된 벡터를 `x`라는 이름의 객체로  저장한다. 
```{r}
x = c(19,22,24,26)
x
```

행렬 형태의 데이터를 직접 생성하기 위해서는 `matrix()`함수를 이용한다. 24개의 값이 담겨있는 벡터 `x`를 이용하여 크기 $4\times 6$의 행렬을 `matrix()`를 이용하여 생성한다. 
```{r}
x = c(128, 128, 160, 142, 157, 179,
      122, 137, 177, 177, 160, 138,
      207, 188, 181, 164, 155, 175,
      120, 208, 199, 194, 177, 195)
matrix(x, byrow = TRUE, ncol=6)

```
함수 `matrix()`의 인자 `byrow = TRUE`를 지정하여, 벡터 `x`의 값들이 행별로 (즉 가로로 먼저) 채워지게 하였다. 인자 `byrow = FALSE`가 기본값이므로, `byrow` 지정을 생략하면 벡터의 값들이 열별로 (세로로 먼저) 채워진다. 
```{r}
X = matrix(x, nrow = 4)
X
```

통계분석을 위한 자료로서의 데이터는 데이터프레임의 형태로 저장된다. 예를 들어 `worldbank.csv`의 파일에 저장된 데이터를 읽으면 데이터프레임의 형태이다. 

```{r}
wbdata = read.csv("data/worldbank2020.csv")
str(wbdata)
```
`wbdata`에는 9개의 길이가 다른 벡터가 행렬로 저장되어 있으며, 각 벡터가 한 변수의 관측값들을 담고 있다. 
데이터프레임 `wbdata`의 세 번째 변수 `Score`만을 추출하려면 `wbdata$Score` 또는 `wbdata[,3]`등을 이용할 수 있다. `with()`를 이용하여 데이터프레임 내의 변수명을 직접 호출할 수 있다. 
```{r}
with(wbdata, lm(Score ~ GDP.per.capita))
```


`R`의 연산은 기본적으로 원소별로 실행된다. 예를 들어 두 벡터의 합과 곱은 원소별 합과 곱이다. 
```{r}
x = c(1,2,3,4)
y = c(5,6,7,8)
x + y
x * y
```

두 행렬의 곱(`*`) 역시 원소별 곱이다. 
```{r}
(X = matrix(1:6,nrow =3, ncol = 2))
X * X
```
두 행렬의 합 또는 차도 마찬가지로 원소별 곱이다. `X + X`와 `X-X`의 결과를 쉽게 상상할 수 있을 것이다. 
행렬 곱을 하기 위해서는 `%*%`을 이용한다. 
```{r}
(Y = matrix(1:6, nrow = 2, ncol = 3))
X %*% Y
```
행렬의 행과 열을 바꾸는 transpose는 함수 `t()`를 이용한다. 즉, $X'$는 `t(X)`로 구한다. 역행렬을 구하기 위해서는 함수 `solve()`를 이용하여, 만약 행렬 $X$의 역행렬이 존재한다면, $X^{-1}$는 `solve(X)`로 구한다. 이를 이용하여, 예를 들면 회귀분석에서의 햇 행렬 $X(X'X)^{-1}X'$를 계산할 수 있다. 
```{r}
(H = X %*% solve(t(X) %*% X) %*% t(X))
```
길이가 $3$인 벡터 `y`를 크기가 $3 \times 3$인 행렬 $H$와 곱한다면 어떻게 될까? 벡터 `y`를 크기가 $1 \times 3$인 행렬 (열벡터)로 보면 행렬곱 $yH$가 정의되고, 벡터 `y`를 크기가 $3 \times 1$인 행렬 (행벡터)로 보면 행렬곱 $Hy$가 정의된다. 행렬곱이 정의가 되도록 벡터 `y`를 행벡터 또는 열벡터로 자동으로 바꾸어준다. 
```{r}
y = c(0,1,2)
H %*% y # H times y 
y %*% H # y times H 
```

## 패키지와 함수 

`R`의 패키지는 함수와 데이터셋 등의 객체를 모아 놓은 묶음이다. 
앞에서 소개한 `matrix()`, `t()`, `solve()` 등은 `R`의 `base` 패키지에 포함된 함수이다. 다른 함수 `read.csv()`는 `R`의 `utils` 패키지에 포함된 함수이다. 이 두 패키지는 `graphics`, `datasets` 등의 패키지와 함께 `R` 세션에 자동으로 포함된다. 만약 `R`에 포함되지 않은 패키지 내의 함수와 데이터셋을 이용하려면 `library(패키지 이름)`을 실행한다. 예를 들어 `tidyverse` 패키지 내의 함수를 이용하려면, 다음을 실행한다. 

```{r, warning=FALSE, message=FALSE}
library(tidyverse)
```

이 패키지를 처음 사용한다면, `install.packages("tidyverse")`와 같이 설치해 줄 수 있다. `tidyverse` 패키지는 사실 여러 패키지를 모아 놓은 패키지이다. 예를 들어 `ggplot2`, `dplyr` 패키지가 `tidyverse` 패키지에 포함되어 있다. 여러 패키지를 사용하다 보면 같은 이름의 함수 또는 데이터셋이 두 패키지에 모두 정의되어 있는 경우가 있다. 이 때, `패키지이름::함수이름()` 또는 `패키지이름::데이터셋이름`과 같이 어떤 패키지의 어떤 객체를 의미하는지를 확실하게 지정해 줄 수 있다. (이 때 패키지가 설치되어 있다면 `library()`로 로드하는 과정을 생략하고 직접 그 패키지 내의 객체를 불러올 수 있다.)

아래의 첫번째 함수 `select`는 패키지 `dplyr` (패키지 `dplyr`은 `library(tidyverse)`로 이미 `R` 환경에 로드되어 있다) 내의 함수이다.
```{r}
head( select(wbdata, Score) ) 
```

`select(wbdata, Score)`로 데이터셋 `wbdata` 내의 변수 `Score`만을 선택한다. 만약 `MASS` 패키지 내의 `select` 함수를 이용하려면  `MASS::select()`를 사용한다. 우리는 아직 `MASS` 패키지를 로드하지 않았지만, 다음과 같이 `::`를 통해 그 패키지 내의 함수를 이용할 수 있다.

```{r}
MASS::select(MASS::lm.ridge(Score ~ Generosity + Social.support, 
                            data = wbdata, lambda = seq(0,1,0.001)))
```

### 파이프 오퍼레이터 

패키지 `dplyr`을 로드하면 특별한 오퍼레이터 `%>%`를 사용할 수 있다. 파이프 오퍼레이터 `%>%`를 이용하여 `R` 코드를 읽기 편하게 바꾸어 줄 수 있다. 예를 들어, 관측값 벡터 `x`에서 평균 (`mean()`)을 구하려면, 
```{r}
mean(x)
```
를 통상적으로 이용하지만, 읽는 순서 그대로
```{r}
x %>% mean
```
로 사용할 수 있다. 파이프 ` %>% ` 전의 값(벡터, 행렬, 데이터프레임)이 파이프 ` %>% ` 후의 함수 `mean()`의 첫번째 입력값이 되는 구조이다. 파이프 오퍼레이터를 이용하여 `head( select(wbdata, Score) )`를 다시 쓰면 다음과 같다. 
```{r}
wbdata %>% select(Score) %>% head()
```
마지막 함수의 다른 입력값이 없다면, 괄호를 생략해도 된다. 

### 내가 만들어 쓰는 함수
같은 작업을 반복적으로 하기 위해 사용자가 함수를 정의해서 사용할 필요가 종종 있다. 
예를 들어 관측값 `x`의 블록 `block`별로 순위변환을 한다면 다음과 같은 과정을 거쳐 계산할 수 있다.  
```{r}
x = (1:6)^2;block = rep(c("a","b"), each = 3)
(temp.data = data.frame(x,block))
out = temp.data %>% dplyr::group_by(block) %>% dplyr::mutate(ranks = rank(x)) 
out$ranks
```
블록별 순위변환을 위해, `dplyr`의 두 함수 `group_by()`와 `mutate()`를 이용했다. 순서대로 설명하자면, `temp.data`의 관측값들을 `block`별로 그룹(`group_by()`)을 짓고, `ranks`라는 새로운 변수의 값을 각 그룹별로 기존의 변수 `x`를 변환(`mutate()`)하여 만들었다. 이 과정을 여러 데이터에 대한 반복하기 위해 다음과 같이 `rankinblock`의 이름을 가진 함수를 정의한다.
```{r}
rankinblock = function(x,block){
  require(dplyr)
  out = data.frame(x,block) %>% group_by(block) %>% mutate(ranks = rank(x))
  return(out$ranks)
}
```
`require(패키지이름)`은 `library(패키지이름)`과 같은 역할을 하지만 함수 내부에서 쓰이도록 디자인되었다. `?require`로 볼 수 있는 상세한 설명을 참조하자. 

이미 정의가 된 함수는 다음과 같이 사용한다. 
```{r}
rankinblock(x,block)
```

## 비모수 통계에서 자주 쓰이는 표현과 함수

### 순열과 조합

길이가 $n$인 수열 `1:n` 또는 벡터 `v`의 모든 개체들을 다른 순서로 뒤섞는 연산을 순열 (permutation)이라고 한다. 이 책에서는 "순열" 대신 "뒤섞기"라는 용어를 사용했다. $n$개의 개체들을 뒤섞는 모든 순열의 수는 
$$
n! = n(n-1)\cdots 2 \cdot 1
$$
이며 `R`에서는 
```{r}
n = 10
factorial(n)
```
과 같이 계산할 수 있다. 길이가 $n$인 벡터 `v`의 원소 중에서 $r$개의 개체를 골라 배열하는 것 역시 순열의 일부로 볼 수 있으며, 이러한 뒤섞기의 갯수는 $P(n,r) = n! / (n-r)! = n(n-1) \cdots (n-r+1)$이다. `gtools` 패키지의 `permutations()`를 이용하여 모든 뒤섞기를 생성할 수 있다. 

다음 코드에서는 길이 $n = 4$인 벡터 `v = 1:4`에서 $r  = 2$개의 개체를 골라 순서대로 배열하는 모든 뒤섞기를 생성한다. 
```{r}
gtools::permutations(n = 4, r = 2, v = 1:4)
```

제 I부에서 주로 다루는 정확한 순열검정에서는 사실 순열이 아닌 조합을 주로 이용한다. 길이가 $n$인 벡터 `v`에서 순서에 상관없이 $r$개의 개체를 선택하는 것이 곧 조합이다. 그 경우의 수는 
$$ {n \choose r} = \frac{n!}{n!(n-r)!}$$
이다. (${n \choose r}$는 "$n$ choose $r$"이라고 읽는다.) `R`에서는 다음과 같이 계산할 수 있다.
```{r}
choose(4,2)
```
다음 코드에서는 길이 $n = 4$인 벡터 `v = 1:4`에서 $r  = 2$개의 개체를 순서에 상관없이 고르는 모든 뒤섞기를 생성한다. 

```{r}
gtools::combinations(n = 4, r = 2, v = 1:4)
```

두 모집단을 비교할 때, 첫 번째 모집단에서 $n_1$개의 표본 `x1`을, 두 번째 모집단에서 $n_2$개의 표본 `x2`를 얻었다고 가정하자. 만약 두 모집단이 사실은 하나의 모집단이었다면, 혼합표본 `c(x1,x2)`에서 $n_1$개의 표본을 임의로 고른 `x1*`가 첫 번째 모집단에서의 표본이었을 수도 있다. 이러한 경우의 수를 모두 고려하려면 순열(permutation)이 아닌 조합(combination)의 경우의 수를 따져야 한다. 조합의 경우에도 혼합표본을 뒤섞는 것으로 그 경우를 따져 볼 수 있으므로, "뒤섞기"라는 용어를 사용한다. 

```{r}
set.seed(1)
x1 = rpois(n = 3, lambda = 5) # n1 = 3
x2 = rpois(n = 2, lambda = 5) # n2 = 2
(x.combined = c(x1,x2))       # n = 5
gtools::combinations(n = 5, r = 3, v = x.combined, set = F)
```
`x.combined`에 같은 값 `4`가 두 개 포함되어 있다. `gtools::combinations()`의 인자 `set = F`는 입력된 `x.combined`에서 같은 값을 지우지 못하게 강제하는 역할을 한다. 

참고로, `set.seed(1)`을 사용하여 포아송분포에서의 랜덤 표집 (`rpois()`)의 결과가 항상 같게 나오도록 강제하였다. 분석의 재현을 위해 랜덤 표집 전 항상 `set.seed(숫자)`를 이용하는 것이 좋다. 

## 랜덤뒤섞기와 붓스트랩재표집
모든 뒤섞기를 고려하는 것은 표본의 수 $n$이 작을 때에만 가능하다. 크다고 볼 수 없는 표본 수 $n = 30$에서 $n_1 = 15$개를 고르는 모든 뒤섞기의 갯수는 ${30 \choose {15}} = 155117520$ (1억 오천만 개)이다. 따라서, 미리 정해놓은 $M = 1999$ 또는 $M = 9999$에 대해, ${n \choose n_1}$개의 가능한 뒤섞기 중 $M$개 만을 랜덤하게 표집하여 모든 가능한 뒤섞기를 "근사"하는 것이 일반적인 방법이다. 랜덤표집을 위해 `sample`함수를 이용한다. 

다음 코드은 벡터 `x.combined`에서 3개의 원소를 단순비복원추출한다. 
```{r}
set.seed(1)
sample(x.combined, size = 3)
sample(x.combined, size = 3)
sample(x.combined, size = 5)
```

복원추출이 아닌 비복원추출임을 강조하기 위해 인자 `replace = FALSE`를 추가했다. `sample()`함수의 `replace`인자의 기본값이 `FALSE`이므로, 결과는 그 전과 같다. 
```{r}
set.seed(1)
sample(x.combined, size = 3, replace = FALSE) 
```

제 II부에서 다루는 붓스트랩 재표집은 길이가 $n$인 표본 `x`에서 $n$개의 개체를 단순복원추출한다. 
```{r}
sample(x.combined, size = 5, replace = TRUE)
```

### 반복

랜덤순열검정과 붓스트랩 재표집시에는 단순비복원추출 또는 단순복원추출을 여러 번 반복한다. 반복을 실행하는 가장 기본적인 방법은 제어문 `for`를 이용하는 것이다. 

다음 코드에서는 데이터 프레임 `wbdata`의 변수 `Score`에서 50개의 개체를 단순비복원추출하여 평균을 구하는 과정을 $M = 1999$번 반복한다. 
```{r}
M = 1999
set.seed(2)
rand.mean = vector(length = M)
for(i in 1:M){
  rand.mean[i] = sample(wbdata$Score, size = 50) %>% mean
}
head(rand.mean)
```

함수 `replicate()`을 이용하여 더 직관적이면서도 짧게 표현할 수 있다. 
```{r}
set.seed(2)
rand.mean = replicate(M, sample(wbdata$Score, size = 50) %>% mean)
head(rand.mean)
```
`replicate()`은 대체로 랜덤표집을 반복할 때 이용한다. 

행렬의 각 행벡터의 값에 대한 (또는 각 열벡터의 값에 대한) 계산을 반복할 필요도 있다. 예를 들기 위해, 크기 $10 \times 3$ 행렬인 `allperms`를 이용한다. 
```{r}
allperms = gtools::combinations(n = 5, r = 3, v = x.combined, set = F)
```

행렬 allperms의 각 행벡터 (가로벡터)의 값들의 표준편차를 계산하고자 한다. `for`제어문을 이용하는 것이 가장 기본적이다. 

```{r}
sd.allperms = vector(length = choose(5,3))
for(i in 1:choose(5,3)){
  sd.allperms[i] = allperms[i,] %>% sd()
}
round(sd.allperms,2)
```

덜 직관적이지만 `apply()`함수를 이용하여 `for` 제어문보다 빠른 계산을 할 수 있다. 다음 코드를 읽어보면, `allperms` 행렬의 각 행(1이 행, 2가 열에 해당한다)에 대해 표준편차 `sd()`를 계산한 뒤, 소숫점 두 자리로 정리한다. 
```{r}
apply(allperms, 1, sd) %>% round(2)
```

행렬 `allperms`의 각 열벡터(세로벡터)에 대해 변동계수 (표준편차 / 평균)을 계산해 보자. 변동계수를 계산한 뒤 소숫점 두 자리로 정리하는 함수를 `apply()`의 세번째 인자로 직접 지정해준다. 
```{r}
apply(allperms, 2, 
      function(x){
        cv = sd(x)/mean(x)
        round(cv,2)
      })
```

 

## 데이터시각화

데이터와 그 분포를 이해하는 가장 좋은 방법은 적절한 데이터 시각화이다. 이 책의 대부분의 그림은 `tidyverse` 패키지에 포함되어 있는 `ggplot2` 패키지의 `ggplot()`를 이용하여 생성하였다. 여기서 몇 가지 예를 보여준다. 더 배우고 싶다면 https://r4ds.had.co.nz/ 위캠의 R for Data Science 를 참조하자. 

### 일변량 데이터의 시각화
한 변수의 여러 관측값으로부터 그 분포를 시각화하는 가장 대표적인 방법은 히스토그램이다. 

```{r}
x = rnorm(100, mean = 1, sd = 1)
data.frame(x = x) %>% 
  ggplot(aes(x = x)) + 
  geom_histogram(aes(y = ..density..), binwidth = .5, 
                 fill = "grey", color = "white") +
  geom_rug()
```

데이터 `x`의 분포가 모수분포인 $N(1,1)$를 따르는지를 살펴보려면, 히스토그램에 $N(1,1)$의 밀도함수를 겹쳐 그리거나 분위수그림 (qqplot)을 생성한다. 분위수그림에 신뢰구간을 더한 QQ envelope plot을 `car::qqPlot()`으로 생성하는 것을 추천한다. 

```{r}
data.frame(x = x) %>% 
  ggplot(aes(x = x)) + 
  geom_histogram(aes(y = ..density..), binwidth = .5, 
                 fill = "grey", color = "white") +
  geom_function(fun = dnorm, args = list(mean = 1, sd = 1)) 
car::qqPlot(x, dist = "norm", mean = 1, sd = 1)
```

그룹이 있는 일변량 데이터의 그룹별 분포를 비교할 때는 상자그림을 이용한다. 
```{r}
farm = read.csv("data/smartfarm.csv")
farm %>% ggplot(aes(x = block, y = x, color = treatment)) + 
  geom_boxplot()
```

이변량 데이터의 시각화를 위해 산점도를 이용한다. 
```{r}
wbdata %>% ggplot(aes(x = GDP.per.capita, y = Score)) + 
  geom_point()
```


```{r do-not-run, echo = FALSE, eval = FALSE}
knitr::purl("np_chA_R.Rmd")
```