다변량 통계분석의 경우 두 변수 간의 상관관계에 아주 많이 의존합니다.  따라서 상관관계 두 연속형 변수 (continuous variable)의 상관관계를 분석하는 것이 매우 중요한데요, 상관관계 분석 기법으로는

 

통계량 분석

 

  (1) 공분산 (covariance)

  (2) 상관계수 (correlation coefficient)

 

 

그래프 분석 

 

  (1) 산점도 (scatter plot)

  (2) 산점도 행렬(scatter matrix plot)

  (3) 상관계수행렬(correlation coefficient plot)

 

등이 있습니다. 

 

그래프는 오른쪽 메뉴의 'R 그래프/시각화'에서 ggplot 을 이용해서 그래프 작성하는 방법을 이미 소개하였으므로, 이번 포스팅에서는 통계량 분석에 대해서만 알아보도록 하겠습니다.

 

 

R 분석 예제 데이터는 MASS 패키지 내 Cars93 데이터프레임의 고속도로연비(MPG.highway)와 무게(Weight) 연속형 변수를 사용하겠습니다.

 

> library(MASS)
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...

 

 

 

 

상식적으로 생각해보면 차 무게가 많이 나갈수록 고속도로연비가 안좋게 나올텐데요, 아래의 산점도를 보면 역(-)의 상관관계가 있음을 알 수 있습니다.

 

> # Scatter Plot of MPG.highway ~ Weight with linear regression line 
> plot(MPG.highway ~ Weight, data = Cars93, 
+      main = "Scatter Plot of MPG.highway ~ Weight", 
+      xlab = "Weight", ylab = "MPG.highway")
> abline(lm(MPG.highway ~ Weight, data = Cars93), col = "blue", lty = 2)
 

 

 

 

 

 

그럼, 먼저 공분산(covariance)을 알아보고, 그 다음으로 상관계수(correlation coefficient)를 알아보겠습니다.

 

 

(1) 공분산 (covariance) : cov()

 

 

공분산의 이해를 돕기 위해 x가 증가할 때 y도 역시 증가하는 경우를 가정해보면서 설명을 해보겠습니다.

 

- μ(x) 보다 큰 x값에 대응하는 y값도 또한 μ(y)보다 크므로 (X-μ(x)(Y-μ(y))는 양수가 됩니다.(+*+=+)

- μ(x) 보다 작은 x값에 대응하는 y값도 또한 μ(y)보다 작으므로 (X-μ(x)(Y-μ(y))는 양수가 됩니다.(-*-=+)

 

 

반대로, x가 증가할 때 y는 감소하는 경우를 가정해봅시다. 이 경우에는,

 

- μ(x) 보다 큰 x값에 대응하는 y값은 μ(y)보다 작으므로 (X-μ(x)(Y-μ(y))는 음수가 됩니다.(+*-=-)

- μ(x) 보다 작은 x값에 대응하는 y값은 μ(y)보다 크므로 (X-μ(x)(Y-μ(y))는 음수가 됩니다.(-*+=-)

 

즉, 공분산(covriance) 값이 '+'이면 x와 y 두 변수는 같은 방향으로 움직이며 (x가 증가하면 y도 증가, x가 감소하면 y도 감소), 공분산이 '-'이면 x와 y 두 변수는 서로 다른 방향으로 움직인다는 것을 알 수 있습니다.

 

 

R의 cov() 함수를 사용해서 고속도로연비(MPG.highway)와 무게(Weight)와의 공분산을 2가지 방법을 사용해서 구해보겠습니다.

 

 

> ##-- covariance 
> # way 1. cov(x, y)
> with(Cars93, cov(x=MPG.highway, 
+                  y=Weight, 
+                  use="complete.obs", 
+                  method=c("pearson")))
[1] -2549.655
> # way 2. cov(x, y=NULL)
> x <- Cars93[c("MPG.highway", "Weight")]
> with(Cars93, cov(x, 
+                  y=NULL, 
+                  use="complete.obs", 
+                  method=c("pearson")))
            MPG.highway     Weight
MPG.highway     28.4273  -2549.655
Weight       -2549.6546 347977.893

 

 

way 1. cov(x, y)를 사용(x, y는 벡터)하든지,  way 2. cov(x, y=NULL)을 사용(x는 데이터프레임)하든지 결과는 pearson covariance 가 -2549.6546 으로 음의 상관관계를 가지는 것으로 나왔습니다.

 

그런데 단위가 -2549.6546 이면 어느정도로 음의 상관관계라고 하는 것인지 감이 안잡힙니다. 이때 공분산을 표준화한 상관계수를 사용하면 두 변수간의 상관관계의 정도를 가늠할 수 있게 됩니다.

 

 

 

(2) 상관계수 (pearson correlation coefficient) : cor() 

 

 

두 변수간의 관련성을 나타내는 측도서 상관계수는 다음의 특징을 가집니다.

 

  1) 표준화변수들의 공분산은 상관계수가 된다.

  2) 상관계수는 -1 ≤ ρ ≤ 1 의 범위의 값을 가진다.

  3) 상관계수 ρ가 +1에 가까울수록 강한 양의 선형관계를 가지며,

      -1에 가까울수록 강한 음의 선형관계를 가진다.

  4) 상관계수 ρ = 0 이면 두 변수간의 선형관계는 없으며,

      0에 가까울수록 선형관계가 약해진다.

      (단, 비선형 관계를 가질 수는 있음. 그래프 분석 병행 필요)

  5) 위치 변환이나 척도 변환 후에도 상관계수는 변함이 없다. 

 

 

R의 cor() 함수를 사용해서 고속도로연비(MPG.highway)와 무게(Weight)의 pearson 상관계수(pearson correlation coefficient)를 구해보면 다음과 같습니다.

 

 

> ##-- correlation coefficient
> # way 1. cor(x, y)
> with(Cars93, cor(x=MPG.highway, 
+                  y=Weight, 
+                  use="complete.obs", 
+                  method=c("pearson")))
[1] -0.8106581
> 
> 
> x <- Cars93[c("MPG.highway", "Weight")]
> cor(x, 
+     y=NULL, 
+     use="complete.obs", 
+     method=c("pearson"))
            MPG.highway     Weight
MPG.highway   1.0000000 -0.8106581
Weight       -0.8106581  1.0000000
> 
> 
> # options
> # use = c("everythig", "all.obs", "complete.obs", "na.or.complete", pairwise.complete.obs")
> # method = c("pearson", "kendall", "spearman")
 

 

고속도로연비(MPG.highway)와 무게(Weight)의 pearson 상관계수는 -0.8106581 으로서 매우 강한 음의 상관관계를 가지는 것으로 나왔습니다.

 

 

cov()와 cor() 함수에 공통으로 use와 method option이 사용되었는데요, 이에 대해 추가 설명을 하겠습니다.

 

 구분

option 

내용 

 use =

(결측값

처리)

"everything" (default)

 결측값이 있을 경우 NA로 계산 결과 제시

 "all.obs"

 결측값이 있을 경우 오류 발생

 ("Error in cor(x = a, y = b, use = "all.obs", method = c("pearson")) : cov/cor에 결측치들이 있습니다

" error message)

 "complete.obs"

 결측값이 있는 case는 모두 제거된 상태에서 상관계수 계산

 "pairwise.complete.obs"

 상관계수가 계산되는 변수들만을 대상으로 결측값이 있는

 case 제거한  상관계수 계산

 method =

(상관계수
통계량)

 "pearson" (default)

 Peraon correlation coefficient 지정,

 가장 일반적으로 사용

 "kendall"

 Kendall의 순위상관계수 혹은 Kendall의 τ (tau) 지정,

 비모수 상관계수 계산 (정규성 불충족 시)

 "spearman"

 Spearman의 순위상관계수 혹은 Spearman이 ρ(rho) 지정,

 비모수 상관계수 계산 (정규성 불충족 시)

* kendall τ , Spearman ρ 수식은 생략함

 

 

 

아래에는 결측값(NA, missing value)이 포함되어 있을 경우 공분산이나 상관계수 계산 시에 결측값을 처리하는 옵션인 'use' 의 세부 옵션별("everything", "all.obs", "complete.obs", "pairwise.complete.obs")로 어떻게 다른지를 간단한 예제를 들어보았습니다.

 

a 벡터의 5번째 원소에 NA가 들어있고, c 벡터의 1번째 원소에 NA가 들어있는 예제입니다.  use 옵션을 좀더 명확히 이해하는데 도움이 되기를 바랍니다.

 

 - use = "everything" : NA가 포함되어 있으므로 결과는 "NA"

 - use = "all.obs" : NA가 포함되어 있으므로 다음과 같은 error message => "cov/cor에 결측치들이 있습니다"

 - use = "complete.obs" : a, b, c 벡터에서 한개라도 NA가 있으면 그 순번째의 a, b, c 벡터의 원소(만약 데이터 프레임이라면 해당 row 전체, 즉 해당 observation case의 모든 변수)는 통째로 삭제하라는 뜻이므로 a, b, c 벡터 모두의 1번째와 5번째 원소 삭제한 후에 상관계수 계산

 - use = "pairwise.compelte.obs" : a와 b의 짝 (pairwise b/w a and b)에서는 5번째 원소 삭제, a와 c 의 짝(pairwise b/w a and d)에서는 1번째와 5번째 원소 삭제, b와 c의 짝(pairwise b/w b and c)에서는 1번째 원소 삭제한 후에 상관계수 계산

 

> ##------------
> ## cov(), cor()
> ## use : an optional character string giving a method for computing covariances 
> ##       in the presence of missing values
> ##------------
> 
> # dataset for correlation analysis with NAs
> a <- c(10, 9, 8, 6, NA)
> b <- c(9, 6, 4, 3, 2)
> c <- c(NA, 13, 18, 20, 25)
> 
> x <- data.frame(a, b, c)
> 
> 
> ##------------
> # use = "everything"  
> #  ; (default setting) NAs will propagate conceptually, 
> #    i.e., a resulting value will be NA whenever one of its contributing observations is NA
> cor(x, 
+     y = NULL, 
+     use = "everything", 
+     method = c("pearson"))
   a  b  c
a  1 NA NA
b NA  1 NA
c NA NA  1
> 
> 
> ##------------
> # use = "all.obs" 
> # ; the presence of missing observations will produce an error
> cor(x, 
+     y = NULL, 
+     use = "all.obs", 
+     method = c("pearson"))
Error in cor(x, y = NULL, use = "all.obs", method = c("pearson")) : 
  cov/cor에 결측치들이 있습니다
> 
> 
> ##------------
> # use = "complete.obs" 
> #  ; missing values are handled by casewise deletion 
> #    (and if there are no complete cases, that gives an error)
> cor(x, 
+     y = NULL, 
+     use = "complete.obs", 
+     method = c("pearson"))
           a          b          c
a  1.0000000  0.9285714 -0.9078413
b  0.9285714  1.0000000 -0.9986254
c -0.9078413 -0.9986254  1.0000000
> 
> 
> # double check
> a <- c(10, 9, 8, 6, NA)
> b <- c(9, 6, 4, 3, 2)
> c <- c(NA, 13, 18, 20, 25)
> 
>   # use = "complete.obs" : delete 1st and 5th elements of a, b, c vectors
> a_complete.obs <- c(9, 8, 6)
> b_complete.obs <- c(6, 4, 3)
> c_complete.obs <- c(13, 18, 20)
> x_complete.obs <- data.frame(a_complete.obs, b_complete.obs, c_complete.obs)
> 
> cor(x_complete.obs, 
+     y = NULL, 
+     use = "all.obs", # because it is calculated after deleting 1st and 5th elements
+     method = c("pearson"))
               a_complete.obs b_complete.obs c_complete.obs
a_complete.obs      1.0000000      0.9285714     -0.9078413
b_complete.obs      0.9285714      1.0000000     -0.9986254
c_complete.obs     -0.9078413     -0.9986254      1.0000000
> 
> 
> ##------------
> # use = "pairwise.complete.obs" 
> #  ; the correlation or covariance between each pair of 
> #    variables is computed using all complete pairs of observations on those variables
> #    only works with the "pearson" method
> cor(x, 
+     y = NULL, 
+     use = "pairwise.complete.obs", 
+     method = c("pearson"))
           a          b          c
a  1.0000000  0.9221389 -0.9078413
b  0.9221389  1.0000000 -0.9824719
c -0.9078413 -0.9824719  1.0000000
> 
> # double check 
> a <- c(10, 9, 8, 6, NA)
> b <- c(9, 6, 4, 3, 2)
> c <- c(NA, 13, 18, 20, 25)
>   
>   # use = "pairwise.complete.obs" : delete 5th element of a and b pair vectors
> a_ab_pairwise.complete.obs <- c(10, 9, 8, 6)
> b_ab_pairwise.complete.obs <- c(9, 6, 4, 3)
> 
> cor(x = a_ab_pairwise.complete.obs, 
+     y = b_ab_pairwise.complete.obs, 
+     use = "all.obs", # because it is calculated after deleting 5th element
+     method = c("pearson"))
[1] 0.9221389
> 
>   # use = "pairwise.complete.obs" : delete 1st and 5th elements of a and c pair vectors
> a_ac_pairwise.complete.obs <- c(9, 8, 6)
> c_ac_pairwise.complete.obs <- c(13, 18, 20)
> 
> cor(x = a_ac_pairwise.complete.obs, 
+     y = c_ac_pairwise.complete.obs, 
+     use = "all.obs", 
+     method = c("pearson"))
[1] -0.9078413
> 
>   # use = "pairwise.complete.obs" : delete 1st element of b and c pair vectors
> b_bc_pairwise.complete.obs <- c(6, 4, 3, 2)
> c_bc_pairwise.complete.obs <- c(13, 18, 20, 25)
> 
> cor(x = b_bc_pairwise.complete.obs, 
+     y = c_bc_pairwise.complete.obs, 
+     use = "all.obs", 
+     method = c("pearson"))
[1] -0.9824719

 

 

 

 

두 변수 x, y별 상관계수를 산점도로 나타낸 아래의 그림을 참조하시기 바랍니다.  상관계수가 '0'이면 '선형관계'는 없는 것이지만 '비선형관계'는 있을 수도 있음에 유의하시기 바랍니다.  두번째 그림에서 보면 상관계수가 '1'인 것은 x, y가 선으로 되어있는 관계이며, 기울기는 상관없음을 알 수 있습니다.  y 값이 '0'인 선은 상관계수 값이 없습니다.

 

 

[ Examples of correlation coefficients ]

* source : https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient

 

 

 

상관관계(correlation) 분석할 때 인과관계(causality)와 혼동하는 경우를 종종 봅니다.  상관관계(correlation)는 단순히 두 변수 x와 y간의 상관성만을 보는 것이며, x와 y가 상관계수가 높다고 해서 x가 y의 원인이라거나 y가 x의 원인이라고 해석하면 안됩니다.

 

인과관계(causality)에는 아래처럼 6가지 종류의 인과관계가 있으며, 시간의 흐름(x chronologically precedes y), 영향을 끼치는 방향(x causally precedes y)을 함게 고려하여야 하므로 사용, 해석에 주의를 요합니다.

 

[ 6 different kinds of causality ]

 

There are six different kinds of causality within a model.

 

1. A direct causal relationship is one in which a variable, X, is a direct cause of another variable, Y (i.e. it is the immediate determinant of Y within the context of the theoretical system).

 

2. An indirect causal relationship is one in which X exerts a causal impact on Y, but only through its impact on a third variable, Z.

 

3. A spurious relationship is one in which X and Y are related, but only because of a common cause, Z.  There is no formal causal link between X and Y.

 

4. A bi-directional or reciprocal causal relationship is one in which X has a causal influence on Y, which in turn, has a causal impact on X.

 

5. An unanalyzed relationship is one in which X and Y or related, but the source of the relationship is unspecified.

 

6. A moderated causal relationship is one in which the relationship between X ad Y is moderated by a third variable.  In other words, the nature of the relationship between X and Y varies, depending on the value of Z.

 

From Jaccard, J., & Turrisi, R. (2003). Interaction Effects in Multiple Regression. Thousand Oaks, CA: Sage. 

 

* source : http://depts.washington.edu/methods/causality.html

 

 

 

위의 영어 설명을 예를 들어서 도식화하면 아래와 같습니다.

 

 

[ 인과관계의 종류 및 예시 ]

 

 

 

 

산점도, 산점도 행렬, 상관계수 행렬 그림 그리기는 아래 포스팅을 참조하시기 바랍니다.

 

ggplot2의 geom_point() 함수를 이용한 한개 변수의 산점도 그리기 

R Graphics 패키지 pairs() 함수를 사용한 산점도 행렬 그리기

corrplot 패키지를 이용한 상관계수 행렬 그림 그리기

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

Posted by R Friend R_Friend

이전 포스팅에서 연속형 변수(continuous variable)의 중심화 경향(central tendancy), 퍼짐 정도(dispersion), 분포형태 및 대칭정도 (distribution and symmetry) 등의 R 함수에 대해 알아보았습니다.

 

R에서는 이들 개별 통계량에 대해 개별 R 함수를 제공함과 동시에, 참 편리하게도 한방에(!) 연속형 변수의 기술통계량을 볼 수 있게 해주는 함수들이 있는데요,

 

(1) 연속형 변수 요약통계 한번에 보기

  :  summary(), stat.desc(), describe()

 

(2) 연속형 변수 그룹별(요인별) 요약통계 비교하기

  : tapply(var, factor, summary), aggregate(), summaryBy(), describe.by()

 

중에서 이번에는 R을 이용한 (2) 연속형 변수 그룹별(요인별) 요약통계 비교하기를 소개하겠습니다.

 

 

 

[ 연속형 변수 요약통계 한번에 보기 package & function 요약 ]

 

 category

 package

function 

statistics 

 (1)

descriptive

statistics

summary

base

summary() 

 min, 1Q, median, mean, 3Q, max

pastecs

 stat.desc()

 nbr.val, nbr.null, nbr.na, min, max, range, sum,

 median, mean, var, std.dev, coef.var,

 skewness, kurtosis, normtest.W, normtest.p,

 SE.mean, CI.mean.p

psych 

describe() 

 n, mean, std.dev, median, trimmed, mad,

 min, max, range, skew, kurtosis, se

 (2)

descriptive

statistics

comparison

by group (factor)

 base

tapply(var, factor, summary) 

 summary: min, 1Q, median, mean, 3Q, max,

 (user defined) functions...

 base

by() 

 (user defined) functions... 

 stats

aggregate() 

 (user defined) functions...

 doBy

summaryBy() 

 (user defined) functions... 

 psych

describeBy() 

 n, mean, sd, median, trimmed, mad,

 min, max, range, skew, kurtosis, se

 

 

 

MASS package 내 Cars93 데이터프레임의 차종(Type)별 가격(Price), 고속도로연비(MPG.highway) 연속형 변수에 대해서 요약통계를 계산해보겠습니다.

 

 
> library(MASS)
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...

 

 

 

 

 각 package별 R 함수 사용예는 아래와 같습니다.

 

(1) base package : tapply(var, factor, summary)

 

 

> ##-- (1) base package : with(dataset, tapply(variable, factor, summary))

> > with(Cars93, tapply(Price, Type, summary)) $Compact Min. 1st Qu. Median Mean 3rd Qu. Max. 11 13 16 18 21 32 $Large Min. 1st Qu. Median Mean 3rd Qu. Max. 18 20 21 24 27 36 $Midsize Min. 1st Qu. Median Mean 3rd Qu. Max. 14 17 26 27 34 62 $Small Min. 1st Qu. Median Mean 3rd Qu. Max. 7.4 8.6 10.0 10.2 11.3 15.9 $Sporty Min. 1st Qu. Median Mean 3rd Qu. Max. 10 14 17 19 22 38 $Van Min. 1st Qu. Median Mean 3rd Qu. Max. 16 19 19 19 20 23

 
> with(Cars93, tapply(MPG.highway, Type, summary))
$Compact
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     26      28      30      30      31      36 

$Large
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     25      26      26      27      28      28 

$Midsize
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     22      25      26      27      29      31 

$Small
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     29      33      33      36      37      50 

$Sporty
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     24      25      28      29      31      36 

$Van
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     20      21      22      22      23      24

 

tapply함수를 사용해서 특정 변수에 요인(factor)별 summary()함수를 적용하였습니다.  가격(Price)과 고속도로연비(MPG.highway) 를 한꺼번에 못하고 두번에 나누어서 적용해야 하며, 제시된 포맷도 일괄 비교해서 보기에는 조금 불편한 점이 있습니다.

 

 

 

(2) base package : by()

 

 

> ##-- (2) base package : by()
> fun_summary <- function(x, ...) {
+   c(n=sum(!is.na(x)), mean=mean(x, ...), sd=sd(x, ...))
+ }
> 
> by(Cars93[c("Price", "MPG.highway")], # dataset$variable
+    Cars93$Type, # by Group(Factor)
+    function(x) sapply(x, fun_summary, na.rm=T)) # function
Cars93$Type: Compact
     Price MPG.highway
n     16.0        16.0
mean  18.2        29.9
sd     6.7         2.9
-------------------------------------------------------------------------- 
Cars93$Type: Large
     Price MPG.highway
n     11.0        11.0
mean  24.3        26.7
sd     6.3         1.3
-------------------------------------------------------------------------- 
Cars93$Type: Midsize
     Price MPG.highway
n       22        22.0
mean    27        26.7
sd      12         2.5
-------------------------------------------------------------------------- 
Cars93$Type: Small
     Price MPG.highway
n       21        21.0
mean    10        35.5
sd       2         5.6
-------------------------------------------------------------------------- 
Cars93$Type: Sporty
     Price MPG.highway
n       14        14.0
mean    19        28.8
sd       8         3.6
-------------------------------------------------------------------------- 
Cars93$Type: Van
     Price MPG.highway
n      9.0         9.0
mean  19.1        21.9
sd     1.9         1.5

 

 

두번째는 base package의 by()함수를 사용하는 방법인데요, 위의 예에서는 관측값 개수 n, 평균 mean, 표준편차 sd 를 그룹(요인)별로 일괄 계산하는 script가 되겠습니다. 

 

by() 함수는 두 개 이상의 기술통계량을 그룹(요인)별로 계산하고 싶다면 사용자 정의 함수(user defined function)를 미리 만들어놓아야 하며, by() 함수 내에 sapply()를 적용해야 하므로 초보자가 이용하기에는 어려운 방법이 되겠습니다.  너무 복잡해서 외우기도 사실상 어렵습니다. -_-;

 

사용자 정의함수에서 '...' 은 '결측값 제거 후 계산'과 같은 추가 옵션을 설정할 수 있도록 하는 명령문입니다.

 

 

 

(3) stats package : aggregate()

 

 

> ##-- (3) stats package : aggregate()
> library(stats)
> fun_summary_2 <- function(x, ...) {
+   c(n=sum(!is.na(x)), mean=mean(x, ...), sd=sd(x, ...))
+ }
> 
> aggregate(Cars93[c("Price", "MPG.highway")], # dataset$variable
+           by=list(Car_Type=Cars93$Type), # by Group(Factor)
+           fun_summary_2, na.rm=T) # function
  Car_Type Price.n Price.mean Price.sd MPG.highway.n MPG.highway.mean MPG.highway.sd
1  Compact    16.0       18.2      6.7          16.0             29.9            2.9
2    Large    11.0       24.3      6.3          11.0             26.7            1.3
3  Midsize    22.0       27.2     12.3          22.0             26.7            2.5
4    Small    21.0       10.2      2.0          21.0             35.5            5.6
5   Sporty    14.0       19.4      8.0          14.0             28.8            3.6
6      Van     9.0       19.1      1.9           9.0             21.9            1.5

 

 

세번째로, stats package의 aggregate() 함수는 그룹(요인) 자리에 list() 로 들어가며, by()함수처럼 두개 이상의 통계량을 보려면 사용자 정의 함수를 미리 만들어서 사용해야 하므로 초보자에게 역시 쉽지 않습니다. 

 

통계량들이 제시되는 포맷이 tapply(summary), by()와는 달리, 위에서 보는 것처럼 row에 그룹(요인)이 위치하고 column에 각 통계량들이 나열되어 있습니다.

 

 

 

(4) doBy package : summaryBy()

 

 

> ##-- (4) doBy package : summaryBy()
> install.packages("doBy")
Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/doBy_4.5-13.zip'
Content type 'application/zip' length 3431577 bytes (3.3 MB)
downloaded 3.3 MB

package ‘doBy’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpqgHvNR\downloaded_packages
> library(doBy)
필요한 패키지를 로딩중입니다: survival
> 
> 
> summaryBy(Price + MPG.highway ~ Type, data=Cars93, FUN=c(mean, sd))

     Type Price.mean MPG.highway.mean Price.sd MPG.highway.sd
1 Compact         18               30      6.7            2.9
2   Large         24               27      6.3            1.3
3 Midsize         27               27     12.3            2.5
4   Small         10               35      2.0            5.6
5  Sporty         19               29      8.0            3.6
6     Van         19               22      1.9            1.5

 

 

 네번째로 doBy package의 summaryBy() 함수는 기본 표현법이 summaryBy( var1 + var2 ~ factor, data, FUN=c(,,)) 으로서 초보자도 사용하기에 매우 쉽고 직관적으로 되어있습니다.  분석결과가 제시되는 포맷도 row에 그룹(요인)이 오고 column에 기초통계량이 와서 일목요연하게 그룹 간 비교하기가 편하게 되어있어서 제가 제일 선호하는 함수입니다.

 

 

 

(5) psych package : describe.by()

 

 

> ##-- (5) psych package : describeBy() > library(psych) > > describeBy(Cars93[c("Price", "MPG.highway")], Cars93$Type, mat=FALSE) # list group: Compact vars n mean sd median trimmed mad min max range skew kurtosis se Price 1 16 18 6.7 16 18 5.3 11 32 21 0.78 -0.80 1.67 MPG.highway 2 16 30 2.9 30 30 3.0 26 36 10 0.48 -0.82 0.74 -------------------------------------------------------------------------- group: Large vars n mean sd median trimmed mad min max range skew kurtosis se Price 1 11 24 6.3 21 24 3.7 18 36 18 0.82 -1.0 1.91 MPG.highway 2 11 27 1.3 26 27 1.5 25 28 3 -0.07 -1.9 0.38 -------------------------------------------------------------------------- group: Midsize vars n mean sd median trimmed mad min max range skew kurtosis se Price 1 22 27 12.3 26 26 12.7 14 62 48 1.06 0.65 2.61 MPG.highway 2 22 27 2.5 26 27 2.2 22 31 9 0.11 -0.96 0.54 -------------------------------------------------------------------------- group: Small vars n mean sd median trimmed mad min max range skew kurtosis se Price 1 21 10 1.9 10 10 2.1 7.4 16 8.5 1 1.08 0.43 MPG.highway 2 21 35 5.6 33 35 4.4 29.0 50 21.0 1 0.22 1.22 -------------------------------------------------------------------------- group: Sporty vars n mean sd median trimmed mad min max range skew kurtosis se Price 1 14 19 8.0 17 19 4.3 10 38 28 1.0 -0.1 2.13 MPG.highway 2 14 29 3.6 28 29 4.4 24 36 12 0.4 -1.1 0.97 -------------------------------------------------------------------------- group: Van vars n mean sd median trimmed mad min max range skew kurtosis se Price 1 9 19 1.9 19 19 0.89 16 23 6.4 0.17 -0.6 0.63 MPG.highway 2 9 22 1.4 22 22 1.48 20 24 4.0 -0.05 -1.7 0.48 > >

> describeBy(Cars93[c("Price", "MPG.highway")], Cars93$Type, mat=TRUE) #matrix
             item  group1 vars  n mean   sd median trimmed   mad  min max range   skew kurtosis
Price1          1 Compact    1 16   18  6.7     16      18  5.34 11.1  32  20.8  0.780    -0.80
Price2          2   Large    1 11   24  6.3     21      24  3.71 18.4  36  17.7  0.815    -1.03
Price3          3 Midsize    1 22   27 12.3     26      26 12.68 13.9  62  48.0  1.060     0.65
Price4          4   Small    1 21   10  2.0     10      10  2.08  7.4  16   8.5  1.007     1.08
Price5          5  Sporty    1 14   19  8.0     17      19  4.30 10.0  38  28.0  1.029    -0.10
Price6          6     Van    1  9   19  1.9     19      19  0.89 16.3  23   6.4  0.166    -0.60
MPG.highway1    7 Compact    2 16   30  2.9     30      30  2.97 26.0  36  10.0  0.483    -0.82
MPG.highway2    8   Large    2 11   27  1.3     26      27  1.48 25.0  28   3.0 -0.068    -1.89
MPG.highway3    9 Midsize    2 22   27  2.5     26      27  2.22 22.0  31   9.0  0.106    -0.96
MPG.highway4   10   Small    2 21   35  5.6     33      35  4.45 29.0  50  21.0  1.021     0.22
MPG.highway5   11  Sporty    2 14   29  3.6     28      29  4.45 24.0  36  12.0  0.399    -1.06
MPG.highway6   12     Van    2  9   22  1.5     22      22  1.48 20.0  24   4.0 -0.049    -1.72
               se
Price1       1.67
Price2       1.91
Price3       2.61
Price4       0.43
Price5       2.13
Price6       0.63
MPG.highway1 0.74
MPG.highway2 0.38
MPG.highway3 0.54
MPG.highway4 1.22
MPG.highway5 0.97
MPG.highway6 0.48

 

 

 마지막으로 psych package의 describeBy() 함수의 경우 분석가가 임으로 요약통계량을 지정할 수 없게 되어있으며, describe()함수에서 제시되었던 요약통계량들(n, mean, sd, median, trimmed mean, mad, min, max, range, skewness, kurtosis, se 등)이 자동으로 그룹별로 제시됩니다.

 

mat=FALSE 옵션을 설정하면 list 형식으로 결과가 제시되며, mat=TRUE 옵션을 지정하면 matrix 형식으로 결과가 나타납니다.

 

위의 package별 function별 제공하는 요약통계량과 제시 포맷을 참고해서 필요로 하는 방법을 찾아서 사용하기 바랍니다.

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

Posted by R Friend R_Friend

지난 포스팅에서 연속형 변수(continuous variable)의 중심화 경향(central tendancy), 퍼짐 정도(dispersion), 분포형태 및 대칭정도 (distribution and symmetry) 등의 R 함수에 대해 알아보았습니다.

 

R에서는 이들 개별 통계량에 대해 개별 R 함수를 제공함과 동시에, 참 편리하게도 한방에(!) 연속형 변수의 기술통계량을 볼 수 있게 해주는 함수들이 있는데요,

 

(1) 연속형 변수 요약통계 한번에 보기

  :  summary(), stat.desc(), describe()

 

(2) 연속형 변수 그룹별(요인별) 요약통계 비교하기

  : tapply(var, factor, summary), aggregate(), summaryBy(), describe.by()

 

로 나누어서 소개해보겠습니다. 

 

먼저, 한국 아침 드라마 초간단 요약 사진을 보면서 자료 정리 및 요약의 효용에 대해 직관적으로 이해를 한 후에, R을 이용한 (1) 연속형 변수 요약통계 한번에 보기를 시작하겠습니다.

 

 

[ 한국 아침 드라마 초간단 요약, ㅋㅋ~ ]

(* 사진출처 : 미상)

 

 

[ 연속형 변수 요약통계 한번에 보기 package & function 요약 ]

 

 category

 package

function 

statistics 

 (1)

descriptive

statistics

summary

base

summary() 

 min, 1Q, median, mean, 3Q, max

pastecs

 stat.desc()

 nbr.val, nbr.null, nbr.na, min, max, range, sum,

 median, mean, var, std.dev, coef.var,

 skewness, kurtosis, normtest.W, normtest.p,

 SE.mean, CI.mean.p

psych 

describe() 

 n, mean, std.dev, median, trimmed, mad,

 min, max, range, skew, kurtosis, se

 (2)

descriptive

statistics

comparison

by group (factor)

 base

tapply(var, factor, summary) 

 summary: min, 1Q, median, mean, 3Q, max,

 (user defined) functions...

 base

by() 

 (user defined) functions... 

 stats

aggregate() 

 (user defined) functions...

 doBy

summaryBy() 

 (user defined) functions... 

 psych

describeBy() 

 n, mean, sd, median, trimmed, mad,

 min, max, range, skew, kurtosis, se

 

 

MASS package 내 Cars93 데이터프레임의 가격(Price), 고속도로연비(MPG.highway) 연속형 변수에 대해서 요약통계를 계산해보겠습니다.

 

 
> library(MASS)
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...

 

 

 

 

 각 package별 R 함수 사용예는 아래와 같습니다.

 

(1) base package : summary()

 

 

> ##--- (1) base package : summary()
> 
> # summary of one variable
> summary(Cars93$Price)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      7      12      18      20      23      62 
> 
> # summary of muliple variables
> summary(Cars93[c("Price", "MPG.highway")])
     Price     MPG.highway
 Min.   : 7   Min.   :20  
 1st Qu.:12   1st Qu.:26  
 Median :18   Median :28  
 Mean   :20   Mean   :29  
 3rd Qu.:23   3rd Qu.:31  
 Max.   :62   Max.   :50

 

 

base package의 summary() 함수는 중심화 경향과 퍼짐정도에 대해서 quick 하게 볼 수 있는 통계량들을 제공합니다.  base package에 들어있으므로 사용하기에 편리한 반면, 약간 부족한 듯한 면이 있습니다.

  •  summary() : min, 1Q, median, mean, 3Q, max

 

 

(2) pastecs package : stat.desc()

 

 

> ##--- (2) pastecs package : stat.desc()
> install.packages("pastecs")
Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/pastecs_1.3-18.zip'
Content type 'application/zip' length 1636171 bytes (1.6 MB)
downloaded 1.6 MB

package ‘pastecs’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpILyJlg\downloaded_packages
> library(pastecs)
필요한 패키지를 로딩중입니다: boot

다음의 패키지를 부착합니다: ‘boot’

The following object is masked from ‘package:survival’:

    aml

 

> # summary statistics of "Price", "MPG.highway" by stat.desc() function

> round(
+   stat.desc(Cars93[c("Price", "MPG.highway")], # dataframe[c(variable name)]
+                 basic = TRUE, # nbr.val, nbr.null, nbr.na, min, max, range, sum
+                 desc = TRUE, # median, mean, var, std.dev, coef.var
+                 norm = TRUE, # skewness, kurtosis, normtest.W, normtest.p
+                 p = 0.90 # SE.mean, CI.mean.p
+             )
+   , 2)
              Price MPG.highway
nbr.val       93.00       93.00
nbr.null       0.00        0.00
nbr.na         0.00        0.00
min            7.40       20.00
max           61.90       50.00
range         54.50       30.00
sum         1814.40     2705.00
median        17.70       28.00
mean          19.51       29.09
SE.mean        1.00        0.55
CI.mean.0.9    1.66        0.92
var           93.30       28.43
std.dev        9.66        5.33
coef.var       0.50        0.18
skewness       1.48        1.19
skew.2SE       2.97        2.38
kurtosis       3.05        2.30
kurt.2SE       3.08        2.32
normtest.W     0.88        0.92
normtest.p     0.00        0.00

 

 

 

stat.desc() 함수의 옵션 별로 제공하는 통계량은 아래와 같습니다.  IQR, quantile 빼고 다 있구요, 화면에 제시해주는 포맷도 아주 깔끔해서 보기에 좋아서 제가 제일 선호하는 함수이기도 합니다.

  • basic = TRUE : 관측치 개수, null 개수, NA 개수, 최소값, 최대값, 범위, 합
  • desc = TRUE : 중앙값, 평균, 분산, 표준편차, 변이계수
  • norm = TRUE : 왜도, 첨도, 정규성 검정통계량, 정규성 검정 P-value
  • p = 0.90 :  신뢰계수 90% (유의수준 10%) 값 => 90% 신뢰구간은 평균 ± CI.mean.0.9 값
    (위의 예 Price의 90% 신뢰구간은 19.51 ± 1.66)

 

 

 

(3) psych package : describe()

 

 
> install.packages("psych")
> library(psych)
> 
> describe(Cars93[c("Price", "MPG.highway")], 
+          na.rm = TRUE, # missing value is not included 
+          interp = TRUE, # method of median calculation
+          skew = TRUE, # skewness, kurtosis
+          ranges = TRUE, # ragne
+          trim = 0.1 # trimmed mean
+          )
            vars  n mean  sd median trimmed mad  min max range skew kurtosis   se
Price          1 93   20 9.7     18      18 8.3  7.4  62    54  1.5      3.0 1.00
MPG.highway    2 93   29 5.3     28      29 4.4 20.0  50    30  1.2      2.3 0.55
 

 

 psych package의 describe() 함수는 summary() 보다는 많고 stat.desc() 보다는 적은 기술통계량을 보여줍니다. 

 

  • describe() : 관측값 개수, 평균, 표준편차, 중앙값, 절삭평균, 평균절대편차, 최소값, 최대값, 범위, 
                    왜도, 첨도, 표준오차

 

summary()와 stat.desc() 함수가 column에 variable을 row에 statistics를 보여주는 반면, describe() 함수는 반대로 row에 variable을, column에 statistics를 보여줍니다. 

 

원하는 제시 포맷을 선택하고, 원하는 기술통계량을 옵션에서 선택해서 연속형 자료를 일목요연하게 정리, 요약하면 되겠습니다.

 

많은 도움 되었기를 바랍니다.  

질문은 댓글로 남겨주세요. 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

Posted by R Friend R_Friend