[14-03-11]

1. 난수 생성 및 분포 함수

난수를 발생시키는 다양한 함수들은 다음과 같다. 자세한 목록은 http://www.stat.umn.edu/geyer/old/5101/rlook.html#dist를 참고하면 된다.
- 이항분포 : rbinom
- F 분포 : rf
- 기하분포 : rgeom
- 초기하분포 : rhyper
- 음이항분포 : rnbinom
- 정규분포 : rnorm
- 포아송 분포 : rpois
- t 분포 : rt
- 균등분포 : runif
위의 함수들은 “함수명(발생시킬 난수의 수, 파라미터1, 파라미터2, …)의 형식으로 호출된다.

rnorm(100, 0, 100)

##   [1]  -52.387  -75.956  -27.478   -9.157  -33.190  -10.232   40.785
##   [8]  137.488  114.543  -90.887 -105.818  115.696  209.621   70.454
##  [15]  169.191  -75.091   45.046  -50.854   70.437    8.152 -151.465
##  [22]  -72.537  -20.029  218.095  -12.332  160.019   58.769 -156.832
##  [29]  -73.633  -80.291  -38.332 -116.135  192.312   36.556   40.692
##  [36]  138.607   -1.217  220.009   87.936   31.348   64.837 -115.356
##  [43]   41.335  -79.014  -40.331  107.846  -93.921   23.212  119.622
##  [50]    3.926   53.602   33.005   46.322   -1.193  113.137  -52.094
##  [57] -173.893  -50.228  136.731  -88.384  -69.888  -28.890   19.346
##  [64]   -8.523  -43.269   30.387  172.145 -127.256 -155.719   71.371
##  [71]    6.641    5.661 -104.431  216.897   53.902   97.753 -151.889
##  [78]  137.263   13.545  -82.675  -83.899   32.189 -134.999   73.689
##  [85]  -29.760 -272.610  111.325   -9.863  -27.814   -6.417   59.362
##  [92]  -77.064  129.069    1.973   61.748   59.556  114.664  -78.234
##  [99] -134.718   74.186

정규분포에서 난수를 발생시킨후, 그 난수로 밀도그림(bin에서 경계가 끊기지 않도록 부드럽게 그린 확률밀도그래프)을 그려보자. 정규분포를 잘 따르고 있음을 알 수 있다.

plot(density(rnorm(1e+05, 0, 100)))

plot of chunk unnamed-chunk-2

난수를 생성하는 함수 외에, 누적분포함수, 분위수를 구하는 함수, 확률 밀도를 구하는 함수는 다음 그림과 같다. < 그림 >

dnorm(-1, 0, 1)

## [1] 0.242

dbinom(3, 10, 0.5)

## [1] 0.1172

dpois(3, 1)

## [1] 0.06131

pnorm(0)

## [1] 0.5

qnorm(0.5)

## [1] 0

2. 기초통계량

2.1 평균, 표본 분산, 표본 표준편차

평균, 표본 분산, 표본 표준편차는 각각 mean(), var(), sd()로 구한다.

mean(1:5)

## [1] 3

var(1:5)

## [1] 2.5

sum((1:5 - mean(1:5))^2)/(5 - 1)

## [1] 2.5

sd(1:5)

## [1] 1.581

sqrt(var(1:5))

## [1] 1.581

2.2 다섯 수치 요약

fivenum()으로 구한다.
summary()는 fivenum()과 유사하지만 다섯 수치 요약에 더해 평균까지 계산해 준다.

fivenum(1:11)

## [1]  1.0  3.5  6.0  8.5 11.0

summary(1:11)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     3.5     6.0     6.0     8.5    11.0

fivenum(1:4)

## [1] 1.0 1.5 2.5 3.5 4.0

summary(1:4)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    1.75    2.50    2.50    3.25    4.00

홀수일 경우 결과가 같지만, 짝수일 경우에는 다른 결과를 출력한다.
이유는 summary()가 1:4에서 제 1분위수를 구할 때, 1 + (4-1) X (¼)로 구하지만, fivenum()에서는 25%백분률의 위치가 1.75로 1과 2사이이므로 두개의 평균인 1.5로 잡기 때문이다.

다섯수치 요약을 별도의 함수로 하나씩 계산해보자.

x <- 1:11
c(min(x), quantile(x, 1/4), median(x), quantile(x, 3/4), max(x))

##       25%       75%      
##  1.0  3.5  6.0  8.5 11.0

사분위범위는 IQR()함수를 사용해 계산한다.

IQR(1:10)

## [1] 4.5

quantile(1:10, c(1/4, 3/4))

##  25%  75% 
## 3.25 7.75

7.75 - 3.25

## [1] 4.5

2.3 최빈값(mode)

table()을 사용해 분할표를 만들고, which.max()를 사용해 최대값이 사용된 색인을 찾는다.

x <- factor(c("a", "b", "c", "c", "c", "d", "d"))
names(which.max(table(x)))

## [1] "c"

3. 표본 추출

3.1 단순 임의 추출(Random Sampling)

sample(데이터, 표본 크기, 복원추출여부, 가중치)의 형식으로 표본을 추출한다.

sample(1:10, 5)  # 비복원 추출

## [1] 10  5  8  1  9

sample(1:10, 5, replace = TRUE)  # 복원 추출

## [1] 3 2 1 8 8

가중치를 주어 표본을 추출 해보자.

sample(1:10, 5, replace = TRUE, prob = 1:10)

## [1] 2 4 7 9 7

sample을 주어진 데이터를 비복원추출로 추출하면 데이터의 순서를 섞는 목적으로 사용할 수 도 있다.

sample(1:10)

##  [1]  3  8  5  7  4 10  9  6  1  2

3.2 층화 임의 추출(Stratified Random Sampling)

데이터가 어떤 변수에 따라 분류했을 때 그 층화변수간 차이가 많이 날때 유용하다.
sampling 패키지의 strata()함수를 사용한다.

library(sampling)
x <- strata(c("Species"), size = c(3, 3, 3), method = "srswor", data = iris)
x

##        Species ID_unit Prob Stratum
## 16      setosa      16 0.06       1
## 26      setosa      26 0.06       1
## 33      setosa      33 0.06       1
## 72  versicolor      72 0.06       2
## 80  versicolor      80 0.06       2
## 81  versicolor      81 0.06       2
## 106  virginica     106 0.06       3
## 121  virginica     121 0.06       3
## 145  virginica     145 0.06       3

getdata(iris, x)

##     Sepal.Length Sepal.Width Petal.Length Petal.Width    Species ID_unit
## 16           5.7         4.4          1.5         0.4     setosa      16
## 26           5.0         3.0          1.6         0.2     setosa      26
## 33           5.2         4.1          1.5         0.1     setosa      33
## 72           6.1         2.8          4.0         1.3 versicolor      72
## 80           5.7         2.6          3.5         1.0 versicolor      80
## 81           5.5         2.4          3.8         1.1 versicolor      81
## 106          7.6         3.0          6.6         2.1  virginica     106
## 121          6.9         3.2          5.7         2.3  virginica     121
## 145          6.7         3.3          5.7         2.5  virginica     145
##     Prob Stratum
## 16  0.06       1
## 26  0.06       1
## 33  0.06       1
## 72  0.06       2
## 80  0.06       2
## 81  0.06       2
## 106 0.06       3
## 121 0.06       3
## 145 0.06       3

Species가 층변수이며, Species별로 3개씩 표본을 뽑는다. 층내에서의 추출 방법은 Simple random sampling with replacement, 즉 단순임의 비복원추출이다.

strata함수는 층별로 다른 수의 표본을 추출할 수 있다는 점이다.

strata(c("Species"), size = c(3, 1, 1), method = "srswr", data = iris)  # 복원추출

##        Species ID_unit    Prob Stratum
## 1       setosa       1 0.05881       1
## 7       setosa       7 0.05881       1
## 50      setosa      50 0.05881       1
## 89  versicolor      89 0.02000       2
## 140  virginica     140 0.02000       3

다수의 층을 기준으로 데이터를 추출할 수도 있다.

iris$Species2 <- rep(1:2, 75)  # 150개를 만든다.
strata(c("Species", "Species2"), size = c(1, 1, 1, 1, 1, 1), method = "srswr", 
    data = iris)

##        Species Species2 ID_unit Prob Stratum
## 27      setosa        1      27 0.04       1
## 2       setosa        2       2 0.04       2
## 73  versicolor        1      73 0.04       3
## 64  versicolor        2      64 0.04       4
## 139  virginica        1     139 0.04       5
## 140  virginica        2     140 0.04       6

각 층마다 동일한 갯수의 표본을 추출하고자 한다면 doBy패키지의 sampleBy()함수를 사용할 수 있다.
형식식은 sampleBy(Formula, 추출할 표본의 비율, 복원 추출 여부(기본값 FALSE), 데이터)이다.

library(doBy)

## Loading required package: survival
## Loading required package: splines
## 
## Attaching package: 'survival'
## 
## The following objects are masked from 'package:sampling':
## 
##     cluster, strata
## 
## Loading required package: MASS

sampleBy(~Species, 0.1, data = iris)

##               Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
## setosa.9               4.4         2.9          1.4         0.2     setosa
## setosa.26              5.0         3.0          1.6         0.2     setosa
## setosa.36              5.0         3.2          1.2         0.2     setosa
## setosa.38              4.9         3.6          1.4         0.1     setosa
## setosa.44              5.0         3.5          1.6         0.6     setosa
## versicolor.75          6.4         2.9          4.3         1.3 versicolor
## versicolor.76          6.6         3.0          4.4         1.4 versicolor
## versicolor.80          5.7         2.6          3.5         1.0 versicolor
## versicolor.86          6.0         3.4          4.5         1.6 versicolor
## versicolor.91          5.5         2.6          4.4         1.2 versicolor
## virginica.107          4.9         2.5          4.5         1.7  virginica
## virginica.136          7.7         3.0          6.1         2.3  virginica
## virginica.137          6.3         3.4          5.6         2.4  virginica
## virginica.142          6.9         3.1          5.1         2.3  virginica
## virginica.149          6.2         3.4          5.4         2.3  virginica
##               Species2
## setosa.9             1
## setosa.26            2
## setosa.36            2
## setosa.38            2
## setosa.44            2
## versicolor.75        1
## versicolor.76        2
## versicolor.80        2
## versicolor.86        2
## versicolor.91        1
## virginica.107        1
## virginica.136        2
## virginica.137        1
## virginica.142        2
## virginica.149        1

150의 10%인 15개를 각각 5개씩 추출하여 표본을 만든다.

3.3 계통 추출(Systematic Sampling)

처음 시작점을 정하고 일정한 간격으로 추출하는 방법이다.
데이터가 순서 모집단일 경우 효율성이 높지만, 주기성이 있다면 편향된 표본을 얻게 된다.
sampleBy()에서 systematic=TRUE옵션을 준다.

x <- data.frame(x = 1:10)
sampleBy(~1, frac = 0.3, data = x, systematic = TRUE)

##   [,1] [,2] [,3]
## 1    1    4    7

첫번째 인자는 표본을 추출할 그룹을 지정하는 formula이다. 층화추출에서는 그룹을 뜻하는 표현을 적어야하지만, 여기서는 그룹의 구분이 없으므로 상수 1을 사용하였다.

x를 사용하면 안된다.

sampleBy(~x, frac = 0.3, data = x, systematic = TRUE)

##    [,1]
## 1     1
## 2     2
## 3     3
## 4     4
## 5     5
## 6     6
## 7     7
## 8     8
## 9     9
## 10   10

4. 분할표

명목형 또는 순서형 자료의 도수를 표 형태로 기록한 것이다.

4.1 분할표의 작성

파이그래프를 그릴때 사용했던 것처럼, table()을 이용해 분할표를 작성할 수 있다.

table(c("a", "b", "b", "b", "c", "c", "d"))

## 
## a b c d 
## 1 3 2 1

xtabs()는 분할표를 만들 때, Formula를 사용해 데이터를 지정할 수 있다.

d <- data.frame(x = c("1", "2", "2", "1"), y = c("A", "B", "A", "B"), num = c(3, 
    5, 8, 7))
d

##   x y num
## 1 1 A   3
## 2 2 B   5
## 3 2 A   8
## 4 1 B   7

(xt <- xtabs(num ~ x + y, data = d))

##    y
## x   A B
##   1 3 7
##   2 8 5

num은 도수칼럼이며, x와 y에 따른 분할표를 만들었다.
num과 같은 도수를 나타내는 칼럼이 없고, 각 관찰 결과가 서로 다른 행으로 표현되어 있다면 '~ 변수 + 변수' 형태로 작성한다.

d2 <- data.frame(x = c("A", "A", "A", "B", "B"), result = c(3, 2, 4, 7, 6))
xtabs(~x, d2)

## x
## A B 
## 3 2

x값이 A인 경우의 수와 B인 경우의 수를 센다.

행과 열의 합을 계산해 표시할 수 있다.

xt

##    y
## x   A B
##   1 3 7
##   2 8 5

margin.table(xt, 1)  # 각 행에 있는 도수의 합

## x
##  1  2 
## 10 13

margin.table(xt, 2)  # 각 열에 있는 도수의 합

## y
##  A  B 
## 11 12

margin.table(xt)  # 전체 도수의 합

## [1] 23

prop.table()은 분할표로부터 각 셀의 비율을 계산한다.(propotion)

xt

##    y
## x   A B
##   1 3 7
##   2 8 5

prop.table(xt, 1)  # 각 행에 대한 셀의 비율

##    y
## x        A      B
##   1 0.3000 0.7000
##   2 0.6154 0.3846

prop.table(xt, 2)  # 각 열에 대한 셀의 비율

##    y
## x        A      B
##   1 0.2727 0.5833
##   2 0.7273 0.4167

prop.table(xt)  # 전체에 대한 셀의 비율

##    y
## x        A      B
##   1 0.1304 0.3043
##   2 0.3478 0.2174

4.2 독립성 검정(Independence Test)

분할표의 행에 나열된 속성과 열에 나열된 속성이 서로에게 아무 영향이 없는 독립인지 아닌지를 테스트 한다.

MASS패키지의 survey데이터를 사용해 학생들의 성별에 따른 운동량에 차이가 있는지 독립성 검정을 해보자.

library(MASS)
data(survey)
str(survey)

## 'data.frame':    237 obs. of  12 variables:
##  $ Sex   : Factor w/ 2 levels "Female","Male": 1 2 2 2 2 1 2 1 2 2 ...
##  $ Wr.Hnd: num  18.5 19.5 18 18.8 20 18 17.7 17 20 18.5 ...
##  $ NW.Hnd: num  18 20.5 13.3 18.9 20 17.7 17.7 17.3 19.5 18.5 ...
##  $ W.Hnd : Factor w/ 2 levels "Left","Right": 2 1 2 2 2 2 2 2 2 2 ...
##  $ Fold  : Factor w/ 3 levels "L on R","Neither",..: 3 3 1 3 2 1 1 3 3 3 ...
##  $ Pulse : int  92 104 87 NA 35 64 83 74 72 90 ...
##  $ Clap  : Factor w/ 3 levels "Left","Neither",..: 1 1 2 2 3 3 3 3 3 3 ...
##  $ Exer  : Factor w/ 3 levels "Freq","None",..: 3 2 2 2 3 3 1 1 3 3 ...
##  $ Smoke : Factor w/ 4 levels "Heavy","Never",..: 2 4 3 2 2 2 2 2 2 2 ...
##  $ Height: num  173 178 NA 160 165 ...
##  $ M.I   : Factor w/ 2 levels "Imperial","Metric": 2 1 NA 2 2 1 1 2 2 2 ...
##  $ Age   : num  18.2 17.6 16.9 20.3 23.7 ...

head(survey[c("Sex", "Exer")])

##      Sex Exer
## 1 Female Some
## 2   Male None
## 3   Male None
## 4   Male None
## 5   Male Some
## 6 Female Some

Exer변수는 요인변수로 "자주”, “가끔”, “전혀안함"의 3개의 수준이 있다.

Sex와 Exer에 대한 분할표를 작성해 보자.

xtabs(~Sex + Exer, data = survey)

##         Exer
## Sex      Freq None Some
##   Female   49   11   58
##   Male     65   13   40

독립성 검정은 chisq.test()를 사용해 검정한다.

chisq.test(xtabs(~Sex + Exer, data = survey))

## 
##  Pearson's Chi-squared test
## 
## data:  xtabs(~Sex + Exer, data = survey)
## X-squared = 5.718, df = 2, p-value = 0.05731

p값이 0.05731로 유의수준 5%에서 '성별과 운동을 독립이다'라는 귀무가설을 기각할 수 없다.

4.3 피셔의 정확성 검정(Fisher's Exact Test)

sample size가 작을 때는 독립성 검정이 잘못된 결과를 나태낼 수 있다. 기대빈도 E가 전체의 20%이상인 경우가 이에 해당한다.

xtabs(~W.Hnd + Clap, data = survey)

##        Clap
## W.Hnd   Left Neither Right
##   Left     9       5     4
##   Right   29      45   143

chisq.test(xtabs(~W.Hnd + Clap, data = survey))

## Warning: Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  xtabs(~W.Hnd + Clap, data = survey)
## X-squared = 19.25, df = 2, p-value = 6.598e-05

fisher.test(xtabs(~W.Hnd + Clap, data = survey))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  xtabs(~W.Hnd + Clap, data = survey)
## p-value = 0.0001413
## alternative hypothesis: two.sided

4.4 맥니마 검정(McNemar Test)

대응표본에 대해, 전후의 차이가 있는지 없는지를 검정하는 테스트이다.
mcnemar.test()로 수행한다.

Performance <- matrix(c(794, 86, 150, 570), nrow = 2, dimnames = list(`1st Survey` = c("Approve", 
    "Disapprove"), `2st Survey` = c("Approve", "Disapprove")))
Performance

##             2st Survey
## 1st Survey   Approve Disapprove
##   Approve        794        150
##   Disapprove      86        570

mcnemar.test(Performance)

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  Performance
## McNemar's chi-squared = 16.82, df = 1, p-value = 4.115e-05

p-value가 0.05보다 작으므로 유의수준5%에서 첫번째 검사와 두번째 검사에서의 찬반성향에 차이가 없다는 귀무가설을 기각한다. 즉, 두 검사간 비율에 차이가 발생했다고 할 수 있다.

5. 적합도 검정

데이터가 어떤 특정한 분포를 따르는지 가정할 수 있는지를 검정하는 방법이다.

5.1 Chi Square Test

분할표를 만들어 사용할 수 있다. 다만 독립성 검정과는 달리 분포로부터 기대밗을 계산하는 것이 차이다.

survey데이터에서 글씨를 왼쪽으로 쓰는 사람과 오른쪽으로 쓰는 사람의 비율이 30% : 70%인지의 여부를 분석해보자.

table(survey$W.Hnd)

## 
##  Left Right 
##    18   218

chisq.test(table(survey$W.Hnd), p = c(0.3, 0.7))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(survey$W.Hnd)
## X-squared = 56.25, df = 1, p-value = 6.376e-14

p-value가 0.05보다 작으므로 30% : 70%이라는 귀무가설을 기각한다.

5.2 Shapiro-Wilk Test

표본이 정규분포로부터 추출된 것인지 테스트하기 위한 방법이다.
shapiro.test()함수를 사용한다.

shapiro.test(rnorm(1000))

## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(1000)
## W = 0.9976, p-value = 0.1498

p-value가 0.05보다 크므로, 데이터가 정규 분포를 따른다는 귀무가설을 기각할 수 없다.

5.3 Kolmogorov-Smirnov Test

비모수 검정으로 경험적 분포함수와 비교대상이 되는 분포의 누적분포함수간의 최대 거리를 통계량으로 사용한다.
ks.test()함수를 사용해 수행한다.
귀무가설은 '주어진 두 데이터가 동일한 분포로부터 추출된 표본이다'이다.

ks.test(rnorm(100), rnorm(100))  # 정규분포와 정규분포의 표본

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(100) and rnorm(100)
## D = 0.09, p-value = 0.8127
## alternative hypothesis: two-sided

ks.test(rnorm(100), runif(100))  # 정규분포와 균일분포의 표분

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(100) and runif(100)
## D = 0.47, p-value = 5.099e-10
## alternative hypothesis: two-sided

'주어진 데이터가 특정분포로부터 뽑은 표본이다'를 검정하는 방법이다.

ks.test(rnorm(1000), "pnorm", 0, 1)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(1000)
## D = 0.0211, p-value = 0.7631
## alternative hypothesis: two-sided

pnorm은 평균이 0, 분산이 1인 정규분포의 누적분포함수이며, rnorm(1000)의 경험적 분포함수와 비교를 한다. p-value가 0.05보다 크므로 정규분포로부터 뽑은 표본이다라는 귀무가설을 기각할 수 없다.

5.4 Q-Q plot

자료가 특정 분포를 따르는지 시각적으로 검토하기 위해 Q-Q plot을 사용한다.
Quantile-Quantile Plot의 약자로 비교하고자하는 분포의 분위수끼리 좌표 평면에 표시하여 그린다.

qqnorm()은 sample과 sample이 정규분포를 만족한다고 가정할 때의 그에 해당하는 Z값을 좌표평면에 플롯한다.
qqline()은 Q-Q Plot에서 데이터가 만족해야하는 직선관계를 그린다.

x <- rnorm(1000, mean = 10, sd = 1)
qqnorm(x)
qqline(x, lty = 2)

plot of chunk unnamed-chunk-32

다른 분포에 대해서는 qqplot()함수를 쓰거나 car::qqPlot()함수를 사용한다.

6. 상관계수

변수간 독립성을 확인하는 수치이다.

6.1 피어슨 상관계수(Pearson Correlation Coefficient)

cor()함수를 사용한다.

cor(iris$Sepal.Width, iris$Sepal.Length)

## [1] -0.1176

Species를 제외한 모든 열의 상관계수행렬을 구해보자.

cor(iris[, 1:4])

##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length       1.0000     -0.1176       0.8718      0.8179
## Sepal.Width       -0.1176      1.0000      -0.4284     -0.3661
## Petal.Length       0.8718     -0.4284       1.0000      0.9629
## Petal.Width        0.8179     -0.3661       0.9629      1.0000

symnum()함수를 사용해 특정 범위 값을 문자로 치환하여 보기 쉽게 표시할 수도 있다.

symnum(cor(iris[, 1:4]))

##              S.L S.W P.L P.W
## Sepal.Length 1              
## Sepal.Width      1          
## Petal.Length +   .   1      
## Petal.Width  +   .   B   1  
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1

corrgram패키지는 상관계수를 시각화하는데 유용한 패키지이다.

library(corrgram)

## Loading required package: seriation

corrgram(cor(iris[, 1:4]), type = "corr", upper.panel = panel.conf)

plot of chunk unnamed-chunk-36

파란색은 양의 상관계수를 뜻하고 빨간색은 음의 상관 계수를 뜻한다. 색의 짙기는 상관계수의 크기를 뜻한다.

6.2 스피어만 상관계수(Spearman's Rank Correlation Coefficient)

두 값의 순위를 사용해 상관계수를 구하는 비모수, 비선형적 방법이다.
이산형 데이터, 순서형 데이터에 적용이 가능하다.
Hmisc::rcorr()을 사용해 구할 수 있다. 데이터 3, 4, 5, 3, 2, 1, 7, 5가 있다고 가정했을때, 순위는 다음과 같이 확인할 수 있다.

x <- c(3, 4, 5, 3, 2, 1, 7, 5)
rank(x)

## [1] 3.5 5.0 6.5 3.5 2.0 1.0 8.0 6.5

sort(rank(x))  # 보기 편하기 위해 정렬한 것뿐이다.

## [1] 1.0 2.0 3.5 3.5 5.0 6.5 6.5 8.0

이렇게 순위를 구한다음 스피어만 상관계수의 공식에 대입하여 구할 수 있다.
함수를 사용해 구해보자.

library(Hmisc)

## Loading required package: grid
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## 
## The following object is masked from 'package:seriation':
## 
##     panel.lines
## 
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, round.POSIXt, trunc.POSIXt, units

m <- matrix(c(1:10, (1:10)^2), ncol = 2)
m

##       [,1] [,2]
##  [1,]    1    1
##  [2,]    2    4
##  [3,]    3    9
##  [4,]    4   16
##  [5,]    5   25
##  [6,]    6   36
##  [7,]    7   49
##  [8,]    8   64
##  [9,]    9   81
## [10,]   10  100

rcorr(m, type = "pearson")$r

##        [,1]   [,2]
## [1,] 1.0000 0.9746
## [2,] 0.9746 1.0000

rcorr(m, type = "spearman")$r

##      [,1] [,2]
## [1,]    1    1
## [2,]    1    1

피어슨 상관계수는 0.97이지만, 스피어만 상관계수는 정확히 1로 나타난다. Y=X^2의 관계가 명확하게 존재하지만 선형관계를 나타내는 피어슨 상관계수는 0.97로 정확히 1이 아님을 알 수 있다.

6.3 켄달의 순위 상관 계수(Kendal's Rank Correlation Coefficient)

(X,Y)형태의 순서쌍으로 데이터가 있을 때, x가 커질수록 y가 커지면 concordant, x가 크지만 y가 작아진다면 discordant로 보는 것이다.
Kendall::Kendall() 함수를 사용한다.

library(Kendall)
Kendall(c(1, 2, 3, 4, 5), c(1, 0, 3, 4, 5))

## tau = 0.8, 2-sided pvalue =0.086

6.4 상관 계수 검정(Correlation Test)

cor.test()를 사용해 상관 계수의 유의성을 판단할 수 있다.
귀무가설은 '상관 계수가 0이다'이다.

cor.test(c(1, 2, 3, 4, 5), c(1, 0, 3, 4, 5), method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  c(1, 2, 3, 4, 5) and c(1, 0, 3, 4, 5)
## t = 3.928, df = 3, p-value = 0.02937
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1698 0.9945
## sample estimates:
##   cor 
## 0.915

cor.test(c(1, 2, 3, 4, 5), c(1, 0, 3, 4, 5), method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  c(1, 2, 3, 4, 5) and c(1, 0, 3, 4, 5)
## S = 2, p-value = 0.08333
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
## 0.9

cor.test(c(1, 2, 3, 4, 5), c(1, 0, 3, 4, 5), method = "kendall")

## 
##  Kendall's rank correlation tau
## 
## data:  c(1, 2, 3, 4, 5) and c(1, 0, 3, 4, 5)
## T = 9, p-value = 0.08333
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau 
## 0.8

7. 추정 및 검정

7.1 일표본 평균

t.test()는 모평균과 모평균의 95% 신뢰구간을 추청함과 동시에 가설검증을 수행한다.

x <- rnorm(30)
t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 1.555, df = 29, p-value = 0.1307
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.05723  0.42059
## sample estimates:
## mean of x 
##    0.1817

이 가설검정의 귀무가설은 '모평균이 0이다'라는 귀무가설이다. 다른 정규분포의 모평균에 대한 귀무가설은 mu를 이용해 지정해 주면 수행할 수 있다.

x <- rnorm(30, mean = 10)
t.test(x, mu = 10)

## 
##  One Sample t-test
## 
## data:  x
## t = -0.4615, df = 29, p-value = 0.6479
## alternative hypothesis: true mean is not equal to 10
## 95 percent confidence interval:
##   9.483 10.327
## sample estimates:
## mean of x 
##     9.905

7.2 독립 이표본 평균

t.test()에 있는 sleep데이터를 사용하여 독립 이표본을 만들어보자.

sleep

##    extra group ID
## 1    0.7     1  1
## 2   -1.6     1  2
## 3   -0.2     1  3
## 4   -1.2     1  4
## 5   -0.1     1  5
## 6    3.4     1  6
## 7    3.7     1  7
## 8    0.8     1  8
## 9    0.0     1  9
## 10   2.0     1 10
## 11   1.9     2  1
## 12   0.8     2  2
## 13   1.1     2  3
## 14   0.1     2  4
## 15  -0.1     2  5
## 16   4.4     2  6
## 17   5.5     2  7
## 18   1.6     2  8
## 19   4.6     2  9
## 20   3.4     2 10

sleep2 <- sleep[, -3]

식별번호인 ID를 제외하면 짝지어진 표본이 아닌, 서로 독립인 그룹 1,2 각각의 반응값으로 볼 수 있다.

수면시간 증가량(extra)의 평균을 구해보자. tapply()나, summaryBy()를 사용하자.

tapply(sleep2$extra, sleep2$group, mean)

##    1    2 
## 0.75 2.33

summaryBy(extra ~ group, sleep2)

##   group extra.mean
## 1     1       0.75
## 2     2       2.33

독립이표본에 대해서는 등분산검정을 먼저 실시해야한다.

var.test(extra ~ group, sleep2)

## 
##  F test to compare two variances
## 
## data:  extra by group
## F = 0.7983, num df = 9, denom df = 9, p-value = 0.7427
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1983 3.2141
## sample estimates:
## ratio of variances 
##             0.7983

p-value는 0.7427로 0.05보다 크므로 등분산이라는 가설을 기각할 수 없다. 따라서 독립표본 t-test를 수행할 수 있다.

t-test()에는 paired, var.equal가 있다. 각각 짝지어진 표본, 등분산을 뜻한다.

t.test(extra ~ group, data = sleep2, paired = FALSE, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  extra by group
## t = -1.861, df = 18, p-value = 0.07919
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.3639  0.2039
## sample estimates:
## mean in group 1 mean in group 2 
##            0.75            2.33

7.3 짝지은 이표본 평균

t.test()의 paired=TRUE로 지정하여 짝지은 이표번 검정을 수행할 수 있다.
with()를 통해 데이터 안의(with) 필드에 바로 접근할 수 있는 기능을 사용해서 group==1과 group==2로 나누고, t.test의 인자로 넘긴다

with(sleep, t.test(extra[group == 1], extra[group == 2], paired = TRUE))

## 
##  Paired t-test
## 
## data:  extra[group == 1] and extra[group == 2]
## t = -4.062, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4599 -0.7001
## sample estimates:
## mean of the differences 
##                   -1.58

7.4 이표본 분산

두개의 확률변수가 있을 때, 모분산을 비교해보자.
var.test함수를 사용해 분산을 비교한다.

with(iris, var.test(Sepal.Width, Sepal.Length))

## 
##  F test to compare two variances
## 
## data:  Sepal.Width and Sepal.Length
## F = 0.2771, num df = 149, denom df = 149, p-value = 3.595e-14
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2007 0.3825
## sample estimates:
## ratio of variances 
##             0.2771

7.5 일표본 비율

베르누이 시행에서 성공할 확률에 대한 비율을 구하는 것이 목표이다.
prop.test()를 사용하여 검정한다.
prop.test(x, n, p = NULL, alternative = c("two.sided”, “less”, “greater”), conf.level = 0.95, correct = TRUE)의 형태이며, p의 NULL값은 0.5이다. 양측검정을 할 것인지, 단측검정을 할것인지, 신뢰수준은 어느정도인지 정할 수 있다.

기본적으로 베르누이 시행을 여러번 반복한 이항분포는 시행횟수가 커질수록 그 성공횟수에 대한 분포가 정규분포로 근사한다.
따라서 모든 검정이 정규분포를 근사한다는 가정하에서 이루어 진다.

prop.test(42, 100)

## 
##  1-sample proportions test with continuity correction
## 
## data:  42 out of 100, null probability 0.5
## X-squared = 2.25, df = 1, p-value = 0.1336
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.3233 0.5229
## sample estimates:
##    p 
## 0.42

prop.test(42, 100, p = 0.4)

## 
##  1-sample proportions test with continuity correction
## 
## data:  42 out of 100, null probability 0.4
## X-squared = 0.0938, df = 1, p-value = 0.7595
## alternative hypothesis: true p is not equal to 0.4
## 95 percent confidence interval:
##  0.3233 0.5229
## sample estimates:
##    p 
## 0.42

반드시 정규분포로 근사해야하는 것은 아니다. binom.test()를 사용하면 이항분포에서 신뢰구간을 계산한다.

binom.test(42, 100)

## 
##  Exact binomial test
## 
## data:  42 and 100
## number of successes = 42, number of trials = 100, p-value = 0.1332
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.3220 0.5229
## sample estimates:
## probability of success 
##                   0.42

binom.test(42, 100, p = 0.4)

## 
##  Exact binomial test
## 
## data:  42 and 100
## number of successes = 42, number of trials = 100, p-value = 0.6843
## alternative hypothesis: true probability of success is not equal to 0.4
## 95 percent confidence interval:
##  0.3220 0.5229
## sample estimates:
## probability of success 
##                   0.42

7.6 이표본 비율

독립인 두집단 X,Y가 각각 서로 다른 이항분포를 따른 다고 할 때, 성공확률이 같은 지에 대한 검정이다.
prop.test()를 사용하며, 성공횟수와 시행횟수를 칼럼형태로 받는다.

prop.test(c(45, 55), c(100, 90))

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  c(45, 55) out of c(100, 90)
## X-squared = 4.307, df = 1, p-value = 0.03796
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.31185 -0.01037
## sample estimates:
## prop 1 prop 2 
## 0.4500 0.6111