이번 포스팅에서는 R Shiny를 사용하여 interactive하게 편하고 쉽게 두 연속형 변수 간 관계를 알아볼 수 있도록 

(1) 상관계수 (Correlation Coefficient)를 구하고, 

(2) 산점도 (Scatter Plot)을 그린 후, 

(3) 선형 회귀선 (Linear Regression Line) 을 겹쳐서 그려주는 

웹 애플리케이션을 만들어보겠습니다. 

 

 

R Shiny - Correlation Coefficient, Scatter Plot, Linear Regression Line

 

예제로 사용할 데이터는 MASS 패키지에 내장되어 있는 Cars93 데이터프레임으로서, Price, MPG.city, MPG.highway, EngineSize, Horsepower, RPM 의 6개 연속형 변수를 선별하였습니다. 

 

왼쪽 사이드바 패널에는 산점도의 X축, Y축에 해당하는 연속형 변수를 콤보박스로 선택할 수 있도록 하였으며, 상관계수와 선형회귀모형을 적합할 지 여부를 선택할 수 있는 체크박스도 추가하였습니다. 

 

오른쪽 본문의 메인 패널에는 두 연속형 변수에 대한 산점도 그래프에 선형 회귀선을 겹쳐서 그려서 보여주고, 상관계수와 선형회귀 적합 결과를 텍스트로 볼 수 있도록 화면을 구성하였습니다. 

 

library(shiny)

# Define UI for application that analyze correlation coefficients and regression

ui <- fluidPage(

   

   # Application title

   titlePanel("Correlation Coefficients and Regression"),

   

   # Select 2 Variables, X and Y 

   sidebarPanel(

     selectInput("var_x", 

                 label = "Select a Variable of X", 

                 choices = list("Price"="Price", 

                                "MPG.city"="MPG.city", 

                                "MPG.highway"="MPG.highway", 

                                "EngineSize"="EngineSize", 

                                "Horsepower"="Horsepower", 

                                "RPM"="RPM"), 

                 selected = "RPM"),

     

     selectInput("var_y", 

                 label = "Select a Variable of Y", 

                 choices = list("Price"="Price", 

                                "MPG.city"="MPG.city", 

                                "MPG.highway"="MPG.highway", 

                                "EngineSize"="EngineSize", 

                                "Horsepower"="Horsepower", 

                                "RPM"="RPM"), 

                 selected = "MPG.highway"), 

     

     checkboxInput(inputId = "corr_checked", 

                   label = strong("Correlation Coefficients"), 

                   value = TRUE), 

     

     checkboxInput(inputId = "reg_checked", 

                   label = strong("Simple Regression"), 

                   value = TRUE)

      ),

      # Show a plot of the generated distribution

      mainPanel(

        h4("Scatter Plot"), 

        plotOutput("scatterPlot"), 

        

        h4("Correlation Coefficent"), 

        verbatimTextOutput("corr_coef"), 

        

        h4("Simple Regression"), 

        verbatimTextOutput("reg_fit")

      )

   )

# Define server logic required to analyze correlation coefficients and regression

server <- function(input, output) {

  library(MASS)

  scatter <- Cars93[,c("Price", "MPG.city", "MPG.highway", "EngineSize", "Horsepower", "RPM")]

  

  # scatter plot

  output$scatterPlot <- renderPlot({

    var_name_x <- as.character(input$var_x)

    var_name_y <- as.character(input$var_y)

    

    plot(scatter[, input$var_x],

         scatter[, input$var_y],

         xlab = var_name_x,

         ylab = var_name_y,

         main = "Scatter Plot")

    

    # add linear regression line

    fit <- lm(scatter[, input$var_y] ~ scatter[, input$var_x])

    abline(fit)

    })

  

  # correlation coefficient

  output$corr_coef <- renderText({

    if(input$corr_checked){

      cor(scatter[, input$var_x], 

          scatter[, input$var_y])

      }

    })

  

  # simple regression

  output$reg_fit <- renderPrint({

    if(input$reg_checked){

      fit <- lm(scatter[, input$var_y] ~ scatter[, input$var_x])

      names(fit$coefficients) <- c("Intercept", input$var_x)

      summary(fit)$coefficients

    }

  })

}

# Run the application 

shinyApp(ui = ui, server = server)

 

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

 

728x90
반응형
Posted by Rfriend
,

혼합모델 군집화(mixture model clustering)의 혼합 모델에는 어떠한 분포도 가능하지만 일반적으로 정규분포(Gaussian normal distribution)을 가장 많이 사용하며, 다음번 포스팅에서는 가우시안 혼합 모델(Gaussian misture model, 줄여서 GMM) 군집화를 다루겠습니다. .

 

혼합모델 군집화(mixture model clustering)에 들어가기에 앞서, 이번 포스팅에서는 정규분포(Gaussian normal distribution)에 대해서 먼저 몸풀기 차원에서 짚어보고 넘아가겠습니다.  준비운동 안하고 혼합모델 군집화에 바로 투입됐다가는 자칫 머리에 쥐가 나거나 OTL 하게 되는 불상사가 예상되기 때문입니다. ^^; 

 

이번 포스팅은 쉬어 가는 코너, 한번에 다 몰라도 좌절하지 않기 코너, 생각날 때 찾아보고 복기하는 코너로 여기면 좋겠습니다.

 

기계학습, 다변량 통계분석 책을 보다 보면 수식, 표기법 때문에 많이 어려워했던 기억이 있을 겁니다.  가장 기본이 되고 많이 사용하는 개념, 수식, 표기법을 중심으로 소개해보겠습니다.

 

 

1) 모집단(Population) vs. 표본집단(Sample)

 

알고자 하는 '전체 객체, 성분, 관측치를 포함하는 집단'을 모집단이라고 하며, 보통은 시간, 비용 상의 이슈로 인해서 모수(parameter)가 알려지지 않은(unknown) 경우가 많습니다.

 

모집단의 특성을 알기 위해 시간과 비용을 감당할 수 있는 범위 내에서 모집단으로 부터 모집단을 잘 대표할 수 있도록 표본(sample)을 추출해서 분석을 하게 됩니다.

 

 

모집단과 표본집단을 알았으니, 각 집단의 기대값과 공분산, 상관행렬, 상관계수는 어떻게 표기하는지 비교해서 살펴보겠습니다.

 

 

2) 모집단(Population)의 기대값과 공분산, 상관행렬, 상관계수

 

성분들이 확률변수로 이루어진 p x 1 확률벡터 X에 대한 평균 벡터(mean vector), 공분산 행렬(covariance matrix), 상관행렬(correlation matrix), 상관계수(correlation coefficient) 는 아래와 같습니다. 

 

공분산 행렬(covariance matrix)을 'Σ' 로 표기한다는 점 유념하시기 바랍니다.  다변량 정규분포 확률분포를 표기할 때 'Σ'가 나옵니다.  (* 주의 : 모두 summation 하라는 뜻이 아님 -.-").

 

 

 

위의 모집단(population)에 대한 표기법을 아래의 표본집단(sample)에 대한 표기법과 비교를 해가면서 유심히 살펴보시고 기억해놓으면 좋겠습니다. 

(강의 듣다보면 혼동해서 잘못쓰는 강사나 교재도 가끔 있음 -_-;) 

 

 

 

(3) 표본집단의 기대값과 공분산, 상관행렬, 상관계수

 

모집단으로부터 n개의 확률표본 X1, X2, ..., Xn에 대해 각 n x p 확률벡터 Xi = (Xi1, Xi2, ..., Xip)' , i = 1, 2, ..., n 의 표본 평균 벡터(sample mean vector), 표본 공분산 행렬(sample covariance matrix), 표본 상관행렬(sample correlation matrix), 표본 상관계수(sample correlation coefficient)는 아래와 같습니다.

 

표본 상관행렬의 우상(rigth upper)과 좌하(left down)의 값은 상호 대칭으로 같습니다. (가령, r12r21은 같은 값을 가짐)

 

 

 


 

 

정규분포는 변수의 개수(1개, 2개, p개)에 따라서 '일변량 정규분포 (univariate normal distribution)', '이변량 정규분포 (bivariate normal distribution)', 다변량 정규분포 (multivariate normal distribution)'로 구분할 수 있습니다.

 

일변량/이변량/다변량 정규분포 확률밀도함수(probability density function)의 표기법과 수식, 그리고 그래프를 살펴보도록 하겠습니다.

 

 

(4) 일변량 정규분포 (univariate normal distribution) vs. 다변량 정규분포 (multivariate normal distribution)

 

일변량 정규분포의 확률변수 X에 대한 확률밀도함수(probability density function) 식은 기억하기가 좀 복잡하고 어려운데요, 자꾸 보다보면 어느샌가 익숙해지고, 더 자주 보다보면 자기도 모르게 기억될 날이 올겁니다. (기억을 했다가도... 인간은 위대한 망각의 동물인지라 두어달만 안보면 다시... ㅋㅋ)

 

 

위의 표기에서 일변량 정규분포 대비 다변량 정규분포의 확률밀도함수의 어디가 서로 다른지 유심히 살펴보시기 바랍니다.  특히 다변량 정규분포의 확률밀도함수 표기법이 통계를 전공하지 않았거나, 선형대수 기본을 모르면 매우 생소하고 어렵게만 느껴질 것 같습니다. (왜 저런거냐고 묻는 댓글 질문은 정중히 사양해염... -_-*)

 

 

 

(5) 이변량 정규분포 (bivariate normal distribution)

 

3차원 이상 넘어가면 사람 머리로 인식하기가 곤란하니, 시각화로 살펴보기에 그나마 가능한 2개의 변수를 가지는 2차원의 '이변량 정규분포(bivariate normal distribution)'의 확률밀도함수(probability density function)에 대해서 알아보겠습니다.  위의 다변량 정규분포의 확률밀도함수 식에서 p = 2 를 대입하면 됩니다.

 

[ 이변량 정규분포 확률밀도함수 - 식 1 ]

 

 

 

위의 이변량 정규분포의 확률밀도함수를 다시 정리해보면 아래와 같습니다.

 

[ 이변량 정규분포 확률밀도함수 - 식 2 ]

 

 

 

'이변량 정규분포 확률밀도함수   식 2'를 보면 이건 뭐... 외계어가 따로 없습니다.  도대체 이렇게 복잡하면서도 쓰임새가 어마무시 많은 요물단지 정규분포를 정리한 천재들이 누구인지 궁금하시지요?

 

위키피디아를 찾아보니, 위대한 수학자 가우스(carl Friedrich Gauss)가 정규분포를 발견했고(그래서 영어로 Gaussian Normal Distribution 이라고 함), 라플라스(Marquis de Laplace)가 그 유명하고 중요하며 시험문제로 반드시 나오는 '중심극한의 정리(Central Limit Theorem)'를 증명하였습니다. 피어슨(Karl Pearson)이 20세기 들어서 '정규분포(Normal Distribution)'이라고 명명하는 것을 확산시키는데 일조했다고 나오네요.  가우스, 라플라스, 피어슨... 역시 천재 거물들 맞군요.

(Normal distribution 이외의 분포는 그럼 비정상(Abnormal) 이라는 소리? -_-?)

 

 

 

 

(* source : https://en.wikipedia.org/wiki/Normal_distribution)

 

 

 

수식으로만 봐서는 저게 뭔가 싶으실텐데요, R의 base 패키지에 내장되어 있는 persp() 함수를 가지고 위의 '이변량 정규분포 - 식 2'를 사용하여 3D plot 으로 그려보겠습니다.


먼저 "mvtnorm" 패키지의 dmvnorm() 함수로 다변량 정규분포 함수를 정의해서 3차원 그래프를 그리는 방법은 아래와 같습니다. 



# Multi-variable Normal Distribution

install.packages("mvtnorm")

library(mvtnorm)


x <- seq(from = -5, to = 5, by = 0.2)

y <- x


gaussian_func <- function(x, y) {

  dmvnorm(cbind(x, y))

}


persp(x, y, outer(x, y, gaussian_func), theta = 30, phi = 10)


 




이번에는 위 포스팅에서 소개한 다변량 정규분포 공식을 이용해서 직접 gaussian_func() 를 정의해보고, 이 함수를 사용해서 persp() 의 3D plot에 manipulate 패키지를 접목해서 동적(interactive plotting)으로 상하(phi), 좌우(theta) 로 돌려가면서 살펴보겠습니다. 

 

아래 첫번째 그래프의 노란색으로 표시해 놓은 단추를 누르면 interactive plot 의 slider 가 나타납니다.  theta(좌~우), phi(상~하)를 조절해가면서 3-D 그래프를 돌려서 보시면 '아, 이변량 정규분포 확률밀도함수가 이런거구나. 느낌 오네~' 하실겁니다.

 

 

##---------------------------------------
## (Gaussian) Mixture Model Clustering
##---------------------------------------

 

# bivariate nomral probability density function
mu1 <- 0   # setting the expected value of x1
mu2 <- 0   # setting the expected value of x2

 

s11 <- 3   # setting the variance of x1
s12 <- 4  # setting the covariance between x1 and x2
s22 <- 3   # setting the variance of x2

 

rho12 <- s12/(sqrt(s11)*sqrt(s22)) # setting the correlation coefficient between x1 and x2
rho12

 

x1 <- seq(-5, 5, length = 50) # generating the vector series x1
x2 <- seq(-5, 5, length = 50) # generating the vector series x2


# multivariate normal density function
gaussian_func <- function(x1, x2){
  term1 <- 1/(2*pi*sqrt(s11*s22*(1-rho12^2)))
  term2 <- -1/(2*(1-rho12^2))
  term3 <- (x1 - mu1)^2/s11
  term4 <- (x2 - mu2)^2/s22
  term5 <- 2*rho12*((x1 - mu1)*(x2 - mu2))/(sqrt(s11)*sqrt(s22))
  term1*exp(term2*(term3 + term4 - term5))
}


# calculating the density values
z_score <- outer(x1, x2, gaussian_func)
head(z_score)

 

 

##------------
# 3-D normal distribution density plot : persp() of {base} package
# with interactive plotting with Manipulate package in RStudio

 

install.packages("manipulate")
library(manipulate)

 

manipulate(persp(x1, x2, z_score,
                 theta = theta_x, # theta gives the azimuthal viewing direction
                 phi = phi_x, # phi gives the colatitude viewing direction
                 main = "Two dimensional Normal Distribution"), 
           theta_x = slider(10, 90, initial = 35), 
           phi_x = slider(10, 90, initial = 10))

 


 

# click the button at yellow circle position

 

[ theta : 35,   phi : 10 ]

 


 

[ theta : 60,   phi : 10 ]

 


[ theta : 80,   phi : 10 ]


[ theta : 90,   phi : 10 ]


 

[ theta : 35,   phi : 30 ]


[ theta : 35,   phi: 50 ]


 

[ theta : 35,   phi: 70 ]


 

[ theta : 35,   phi: 90 ]

 

 

 

 

마지막 그림이 phi = 90 으로서, 꼭대기 상단에서 비행기 타고 내려다본 그림인데요, 이걸 2-D 의 등고선 그림(contour plot)으로 그려보면 아래와 같습니다.

 

 

# 2-D contour plot
contour(x1, x2, z_score, xlab = "x1", ylab = "x2",
        main = "Bivariate(2 dimensional) Normal Distribution : contour plot")

 

 

 

 

 

참고로, 일변량 정규분포 그래프는 http://rfriend.tistory.com/102 의 이전 포스팅을 참고하세요.

 

[Reference]

1. 김재희, 'R 다변량 통계분석', 교우사, 2011

2. Wikipedia :  https://en.wikipedia.org/wiki/Normal_distribution

 

 

이번 포스팅에서 살펴본 '다변량 정규분포'에 대한 기본기를 바탕으로 해서, 다음번 포스팅에서는 '가우시언 혼합 분포 군집화 (Gaussian Mixture Model Clustering)'에 대해서 알아보겠습니다.

 

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

 

728x90
반응형
Posted by Rfriend
,

다변량 통계분석의 경우 두 변수 간의 상관관계에 아주 많이 의존합니다.  따라서 상관관계 두 연속형 변수 (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)

 Pearson 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

 

 



(3) 여러개의 숫자형 변수 간 상관계수 결과를 세로로 긴 형태(Long Format)로 저장하기


여러개의 숫자형 변수로 이루어진 데이터 프레임에서 숫자형 변수의 짝(pairs) 별로 상관계수를 구해서 세로로 길쭉한 형태(long format)으로 결과를 저장하고 싶다면 아래의 for loop 문을 사용한 예제를 참고하세요. 


표준정규분포로 부터 10개 칼럼, 5개 관측치의 난수를 생성하여 예제 숫자형 데이터 프레임을 만들었습니다.


for loop 순환문으로 칼럼 갯수만큼 i, j 를 돌리면 두개의 변수 쌍(pairs)이 두 번씩 출현할텐데요, 두 변수의 순서가 바뀌어도 상관계수 값은 동일하므로 두 변수 쌍은 한번만 계산하도록 for loop 문을 작성했습니다. 그리고 자기 자신과의 상관계수는 항상 '1'이므로 제외하도록 했습니다.  



# 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




R로 상관계수 행렬 그림을 그리는 방법은 https://rfriend.tistory.com/82 를 참고하세요. 

 

두 변수 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 패키지를 이용한 상관계수 행렬 그림 그리기

 

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

 

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

 


728x90
반응형
Posted by Rfriend
,