R, Python 분석과 프로그래밍의 친구 (by R Friend)

> 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 
> # 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

> ##-- 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")
 

 구분

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)

 Pearson correlation coefficient 지정,

 가장 일반적으로 사용

 "kendall"

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

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

 "spearman"

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

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

> ##------------
> ## 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

# generate 50 random numbers from standard normal distribution X~N(0, 1)

set.seed(1004)

x <- rnorm(50, mean=0, sd=1)

df <- data.frame(matrix(x, nrow=5, ncol=10, byrow=TRUE))

df

# X1         X2         X3          X4          X5         X6        X7        X8         X9        X10

# 1 -0.6078408  0.7676292 -0.1642718 -0.02877676  0.01355787 -0.7157169 -1.477976 0.8280146 -0.8951765  1.8841561

# 2 -1.8952567  1.1259354  0.1494180 -0.21511242  2.63756778 -0.2678609  1.123417 0.4842308  0.6318069 -0.8180967

# 3 -0.7470136 -0.5132835  0.8709808  0.83062565 -1.03451768 -1.6140827  2.169269 0.6719930  0.4054343 -0.8316793

# 4  0.9159031  0.3905532 -0.6347263 -1.00658131  0.35622463  0.1083284 -1.712486 0.6846194 -0.6083049 -0.6192043

# 5 -0.7671233  0.3633518 -0.7325597  1.23058116  0.09165232  1.5699387  1.309473 0.9069780 -1.3729150 -0.1163513

# pearson's correlation coefficients in a long format

cor_long <- data.frame() # blank dataframe to save the corr results later

col_nm <- colnames(df)

n <- 1

for (i in 1:(ncol(df)-1)){

  for (j in (i+1):ncol(df)) {

    cor_long[n, 1] <- col_nm[i]

    cor_long[n, 2] <- col_nm[j]

    cor_long[n, 3] <- cor(df[,i], df[,j])

    n <- n+1

  }

}

# rename column names

colnames(cor_long) <- c("col_1", "col_2", "corr_coef")

cor_long

# col_1 col_2    corr_coef

# 1     X1    X2 -0.332162830

# 2     X1    X3 -0.448170165

# 3     X1    X4 -0.438930937

# 4     X1    X5 -0.497470094

# 5     X1    X6  0.100596534

# 6     X1    X7 -0.643972517

# 7     X1    X8  0.346978141

# 8     X1    X9 -0.465219120

# 9     X1   X10  0.051698753

# 10    X2    X3 -0.440364182

# 11    X2    X4 -0.454194089

# 12    X2    X5  0.847973001

# 13    X2    X6  0.342903949

# 14    X2    X7 -0.413224712

# 15    X2    X8 -0.235257160

# 16    X2    X9 -0.064519846

# 17    X2   X10  0.314659871

# 18    X3    X4  0.224054384

# 19    X3    X5 -0.167990551

# 20    X3    X6 -0.865225316

# 21    X3    X7  0.560975754

# 22    X3    X8 -0.500395969

# 23    X3    X9  0.805120716

# 24    X3   X10 -0.244220838

# 25    X4    X5 -0.419495665

# 26    X4    X6  0.170965784

# 27    X4    X7  0.755460153

# 28    X4    X8  0.482063832

# 29    X4    X9 -0.181887756

# 30    X4   X10 -0.001894032

# 31    X5    X6  0.239111079

# 32    X5    X7 -0.039280725

# 33    X5    X8 -0.632907538

# 34    X5    X9  0.356088708

# 35    X5   X10 -0.209073460

# 36    X6    X7 -0.074444629

# 37    X6    X8  0.457440745

# 38    X6    X9 -0.681812940

# 39    X6   X10 -0.023974644

# 40    X7    X8 -0.189050606

# 41    X7    X9  0.442853344

# 42    X7   X10 -0.552145749

# 43    X8    X9 -0.913090461

# 44    X8   X10  0.582373314

# 45    X9   X10 -0.556365998

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