สหสัมพันธ์ (correlation)

ให้ \(X\) และ \(Y\) เป็นตัวแปรสุ่มดังนี้ \[H(x,y)=Pr(X\leq x,Y\leq y)=\int_{-\infty}^x\int_{-\infty}^y h(s,t)dsdt\] ถ้า \(X\) และ \(Y\) เป็นอิสระกัน(independent) จะได้ \[\begin{align*} H(x,y)&=Pr(X\leq x)\times Pr(Y\leq y)\\ &=\int_{-\infty}^x\int_{-\infty}^y h(s,t)dsdt\\ &=\int_{-\infty}^xf_x(s)ds\int_{-\infty}^y f_y(t)dt \end{align*}\] การวัดค่าสหสัมพันธ์เป็นกว่าดูว่าตัวแปรสุ่มสองตัวเป็นอย่างไร โดยสูตรทีี่ใช้ในการหาค่าสหสัมพันธ์ที่นิยม และจะพูดถึงในบทนี้คือ ค่าสหสัมพันธ์เพียร์สัน(Pearson correlation) และจะใช้อธิบายได้ดีถ้าตัวแปรสุ่ม \(X\) และ \(Y\) มีการแจกแจงร่วมแบบปกติ # สัมประสิทธิ์สหสัมพันธ์แบบเพียร์สัน มีนิยามและการคำนวณดังนี้ \[\rho_{xy}=\dfrac{cov(X,Y)}{\sigma_x\sigma_y}=\dfrac{E\left((X-\mu_x)(Y-\mu_y)\right)}{\sqrt{E(X-\mu_x)^2E(Y-\mu_y)^2}} \]

โดยที่

ด้วยการวัดด้วยค่าสัมประสิทธ์สหสัมพันธ์แบบเพียร์สัน ค่าที่ได้จะมีค่าดังนี้ \[-1\leq \rho_{xy}\leq 1\] เพื่อให้เห็นภาพ จะทำการสุ่มข้อมูลจากการแจกแจงปกติ 2 มิติ ที่มีค่าสหสัมพันธ์ที่แตกต่างกัน

# เรียกใช้ package MASS เมื่อการจำข้อมูลจากการแจกแจงปกติ 2 ตัวแปร
library(MASS)
# สหสัมพันธ์เท่ากับ -.90
Sigma <- matrix(c(1,-.9,-.9,1),2,2)
X1<-as.data.frame(mvrnorm(n = 1000, mu=rep(0, 2), Sigma))
colnames(X1)<-c("x","y")
ggplot()+geom_point(data=X1,aes(x=x,y=y))+
  labs(title="ค่าสหสัมพันธ์เท่ากับ -0.9")+
  theme(text = element_text(family = "TH Sarabun New"))

# สหสัมพันธ์เท่ากับ -.5
Sigma <- matrix(c(1,-.5,-.5,1),2,2)
X2<-as.data.frame(mvrnorm(n = 1000, mu=rep(0, 2), Sigma))
colnames(X2)<-c("x","y")
ggplot()+geom_point(data=X2,aes(x=x,y=y))+
  labs(title="ค่าสหสัมพันธ์เท่ากับ -0.5")+
  theme(text = element_text(family = "TH Sarabun New"))

# สหสัมพันธ์เท่ากับ 0
Sigma <- matrix(c(1,0,0,1),2,2)
X3<-as.data.frame(mvrnorm(n = 1000, mu=rep(0, 2), Sigma))
colnames(X3)<-c("x","y")
ggplot()+geom_point(data=X3,aes(x=x,y=y))+
  labs(title="ค่าสหสัมพันธ์เท่ากับ 0")+
  theme(text = element_text(family = "TH Sarabun New"))

# สหสัมพันธ์เท่ากับ 0.5
Sigma <- matrix(c(1,0.5,0.5,1),2,2)
X4<-as.data.frame(mvrnorm(n = 1000, mu=rep(0, 2), Sigma))
colnames(X4)<-c("x","y")
ggplot()+geom_point(data=X4,aes(x=x,y=y))+
  labs(title="ค่าสหสัมพันธ์เท่ากับ 0.5")+
  theme(text = element_text(family = "TH Sarabun New"))

# สหสัมพันธ์เท่ากับ 0.9
Sigma <- matrix(c(1,0.9,0.9,1),2,2)
X5<-as.data.frame(mvrnorm(n = 1000, mu=rep(0, 2), Sigma))
colnames(X5)<-c("x","y")
ggplot()+geom_point(data=X5,aes(x=x,y=y))+
  labs(title="ค่าสหสัมพันธ์เท่ากับ 0.9")+
  theme(text = element_text(family = "TH Sarabun New"))

ข้อสังเกตุ ยิ่งค่าสหสัมพันธ์มีค่าใกล้่ \(-1\) หรือ 1 มากกว่าเท่าไหร่จะลักษณะของข้อมูลจะมีลักษณะเหมือนเส้นตรงมากขึ้น

ทฤษฎีที่สำคัญ คือถ้า \(X\) และ\(Y\) เป็นอิสระกันจะได้ \(\rho_{xy}=0\) แต่ถ้า \(\rho_{xy}=0\) ไม่ได้แสดงว่า \(X\) และ \(Y\) เป็นอิสระกัน สหสัมพันธ์ระหว่างตัวเอง หรือ \(cor(X,X=1)\) และ \(\rho_{xy}=\rho_{yx}\) ถ้าเราเก็บข้อมูล \((x_i,y_i),i=1,2,\cdots,n\) มาจำนวน \(n\) ชุดจะสามารถประมาณค่าสัมประสิทธ์สหสัมพันธ์แบบเพียร์สันได้โดยใช้่สูตร \[r_{x y}=\frac{\sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)\left(y_{i}-\overline{y}\right)}{\sqrt{\sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)^{2}}\sqrt{\sum_{i=1}^{n}\left(y_{i}-\overline{y}\right)^{2}}} \] ค่า \(\rho_{xy}\) เป็นค่าสหสัมพันธ์ของตัวแปรสุ่ม \(X\) และ \(Y\) แต่ค่า \(r_{xy}\) เป็นค่าสหสัมพะนธ์ของของตัวอย่างจาก ตัวแปร \(X\) และ \(Y\) ซึ่งค่า \(r_{xy}\) เป็นตัวแประมาณ \(\rho_{xy}\) นั้นเอง ตัวอย่างการคำนวณในอาร์โดยใช้ data(women) ที่มีตัวแปร คือ ส่วนสูง(height) และน้ำหนัก (weight)

# โหลดกรอบข้อมูล women
data(women)
# แสดงข้อมูลภายในต่างๆ ของกรอบข้อมูล women
str(women)
## 'data.frame':    15 obs. of  2 variables:
##  $ height: num  58 59 60 61 62 63 64 65 66 67 ...
##  $ weight: num  115 117 120 123 126 129 132 135 139 142 ...

การคำนวณด้วยสูตร

cov(women$height,women$weight)/(sd(women$height)*sd(women$weight))
## [1] 0.9954948

การคำนวณด้วยคำสั่งในอาร์

cor(women$height,women$weight)
## [1] 0.9954948

ค่าสหสัมพันธ์์ที่ได้มีสูงเข้าใกล้ แสดง มีความสัมพันธ์กันเชิงเส้นตรงที่มีความชันเป็นบวกสูงมาก ในกรณีที่ข้อมูลมีขนาด

มิติที่มากกว่าสองเป็นต้นไป การวัดสหสัมพันธ์จะถูกแสดงออกในรูปของเมตริกซ์ กรณีของข้อมูล women

cor(women)
##           height    weight
## height 1.0000000 0.9954948
## weight 0.9954948 1.0000000

หมายความว่า \[\begin{bmatrix} \rho_{height,height}&\rho_{height,weight}\\ \rho_{weight,height}&\rho_{weight,weight} \end{bmatrix}=\begin{bmatrix}1&.9955\\.995&1 \end{bmatrix}\] กรณีที่มีมิติมากกว่า 2 เช่นข้อมูล data(rock) เป็นตัวอย่างข้อมูล 4 มิติ

# โหลดกรอบข้อมูล rock
data(rock)
# แสดงข้อมูลภายในต่างๆ ของกรอบข้อมูล rock
str(rock)
## 'data.frame':    48 obs. of  4 variables:
##  $ area : int  4990 7002 7558 7352 7943 7979 9333 8209 8393 6425 ...
##  $ peri : num  2792 3893 3931 3869 3949 ...
##  $ shape: num  0.0903 0.1486 0.1833 0.1171 0.1224 ...
##  $ perm : num  6.3 6.3 6.3 6.3 17.1 17.1 17.1 17.1 119 119 ...

ก็สามารถใช้สูตร \[r_{x y}=\frac{\sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)\left(y_{i}-\overline{y}\right)}{\sqrt{\sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)^{2}}\sqrt{\sum_{i=1}^{n}\left(y_{i}-\overline{y}\right)^{2}}} \] โดยที่ \(x\) หรือ \(y\) คิือ area peri shape หรือ perm ถ้าแสดงผลออกมาเป็นเมตริกซ์สหสัมพันธ์จะได้

cor(rock)
##             area       peri      shape       perm
## area   1.0000000  0.8225064 -0.1821611 -0.3966370
## peri   0.8225064  1.0000000 -0.4331255 -0.7387158
## shape -0.1821611 -0.4331255  1.0000000  0.5567208
## perm  -0.3966370 -0.7387158  0.5567208  1.0000000

ถ้าต้องการเฉพาะ \(\rho_{area,peri}\) ก็สามารถใช้คำสั่ง

cor(rock)[1,2]
## [1] 0.8225064
cor(rock)[2,1]
## [1] 0.8225064
cor(rock)["area","peri"]
## [1] 0.8225064
cor(rock)["peri","area"]
## [1] 0.8225064

ในกรณีที่ต้องนำเสนอเมตริกซ์สหสัมพันธ์ ให้อยู่ในกราฟก็สามารถทำได้โดยใช้ package “corrplot”

#เรียกใช้ชุดคำสั่ง corrplot
library(corrplot)
## corrplot 0.84 loaded

สามารถวาดกราฟโดยใช้่คำสั่ง

# วาดกราฟสหสัมพันธ์
corrplot(cor(rock))

ถ้าต้องการแสดงเพียงครึ่งเดียวให้ใช้คำสั่ง

# วาดกราฟสหสัมพันธ์แบบเฉียงขึ้น
corrplot(cor(rock),type = "upper")

หรือ

# วาดกราฟสหสัมพันธ์แบบเฉียงลง
corrplot(cor(rock),type = "lower")

ถ้าต้องการให้แสดงตัวเลขด้วย

# วาดกราฟสหสัมพันธ์แบบเฉียงขึ้น และแสดงค่าสหสัมพันธ์ในกราฟ
corrplot.mixed(cor(rock))

corrplot น้ียังสามารถแสดงกราฟในรูปแบบอื่นๆ ได้อีกมาก ให้ผู้สนใจ ใช้คำว่า corrplot ค้นหาเพิ่มเติมได้จาก google

การทดสอบสมมุติฐาน

ในหาค่าสหสัมพันธ์นั้น แค่เราเชื่อหรือมีข้อมูลที่ยืมยันได้ \(X\) และ \(Y\) มีการแจกแจงร่วมแบบปกติ ถ้า \(\rho_{xy}\neq 0\) เราจะกล่าวว่า \(X\) และ \(Y\) มีความสัมพันธ์กันแบบเชิงเส้น ในการทดสอบสมมุติฐานจะตั้งสมมุติฐานดังนี้ \[\begin{align*} H_0: \rho_{xy}=0\\ H_a: \rho_{xy}\neq 0 \end{align*}\] จะยอมรับสมมุติฐานระดับนัยยะ \(1-\alpha\) \[t^{-1}_{df=n-2}(\frac{\alpha}{2})\leq\dfrac{r_{xy}}{\sqrt{1-r_{xy}^2}}\sqrt{n-2}\leq t^{-1}_{df=n-2}(1-\frac{\alpha}{2})\] จากข้อมูล women ถ้าต้องทดสอบว่า \(\rho=0\) หรือไม่สามารถทำได้โดยถ้าพิจารณาค่า \(\alpha =5\%\)

str(women)
## 'data.frame':    15 obs. of  2 variables:
##  $ height: num  58 59 60 61 62 63 64 65 66 67 ...
##  $ weight: num  115 117 120 123 126 129 132 135 139 142 ...
# คำนวณค่าสหสัมพันธ์แล้วนำไปเก็บไว้ในตัวแปรชื่อ corre
corre<-cor(women)[1,2] 
# ค่าสถิติที่คำนวณได้ มีค่ามากกว่า คอร์ไทล ที่ .025 ของการแจกแจง t(13) หรือไม่
corre*sqrt(15-2)/sqrt(1-corre^2) > qt(.05/2,15-2)
## [1] TRUE
# ค่าสถิติที่คำนวณได้ มีค่าน้อยกว่า คอร์ไทล ที่ .975 ของการแจกแจง t(13) หรือไม่
corre*sqrt(15-2)/sqrt(1-corre^2) < qt(1-.05/2,15-2)
## [1] FALSE

จากทดสอบพบ ค่าสถิติของ \(\rho\) มากกว่า \(t^{-1}_{df=n-2}(1-\frac{\alpha}{2})\) จึงปฎิเสธสมมุติฐานหลักที่เชื่อว่า \(rho=0\) ดังนั้น ตัวแปร weight และ height มีความสัมพันธ์กันแบบเชิงเส้น ในการกรณีที่ต้องการใช้คำสั่งที่มีอยู่แล้วในอาร์สามารถใช้ package Hmisc ช่วยในการทดสอบสมมุติฐานได้ โดยใช้คำสั่งดังนี้ ถ้าพิจารณาที่ข้อมูล women เหมือนเดิม ข้อมูลที่จะใช้กับ package นี้ต้องเป็น ข้อมูลหรือตัวแปรแบบเมตริกซ์เท่านั้น และควรกำหนดตัวแปรให้ใหม่เมื่อต้องการ

คำนวณาค่าสหสัมพันธ์และ ค่าทางสถิติ

library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
women_cor<-rcorr(as.matrix(women))
str(women_cor)
## List of 3
##  $ r: num [1:2, 1:2] 1 0.995 0.995 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "height" "weight"
##   .. ..$ : chr [1:2] "height" "weight"
##  $ n: int [1:2, 1:2] 15 15 15 15
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "height" "weight"
##   .. ..$ : chr [1:2] "height" "weight"
##  $ P: num [1:2, 1:2] NA 1.09e-14 1.09e-14 NA
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "height" "weight"
##   .. ..$ : chr [1:2] "height" "weight"
##  - attr(*, "class")= chr "rcorr"

จะได้ women_cor$r เป็น เมตริกซ์สหสัมพันธ์ ส่วน women_cor$P เป็นเมตริกซ์ของค่า p-value ใช้สำหรับตัดสินยอมสมมุติฐานหลัก ของการคำนวณค่า แถวและหลักต่างๆ

women_cor$r
##           height    weight
## height 1.0000000 0.9954948
## weight 0.9954948 1.0000000
women_cor$P
##              height       weight
## height           NA 1.088019e-14
## weight 1.088019e-14           NA

จะเห็นได้ว่า ค่า p-value=\(1.08819e1-14<.05\) แสดงว่าปฎิเสธสมมุติฐานหลัก ที่ระดับ \(\alpha=5\%\) จากใช้ package Hmisc สามารถนำมาประยุกต์ใช้การวาดกราฟสหสัมพันธ์ด้วย package corrplot ดังนี้ ถ้า ถ้าค่าสหสัมพันธ์ใดมีการยอมรับสมมุติฐานหลัก สามารถวาดกราฟให้ไม่ต้องแสดงผลออกมาได้ ดังนี้ จากข้อมูล mtcars ถ้าวาดกราฟสหสัมพันธ์แบบปกติ จะได้

corrplot(cor(mtcars),type = "upper")

ถ้าใช้คำสั่ง rcorr จาก package Hmisc ช่วยจะได้

mtcars_cor<-rcorr(as.matrix(mtcars))

ถ้าสหสัมพันธ์ระหว่างตัวแปรใด มีค่า p-value น้อยกว่า .1 ไม่ต้องแสดงในกราฟสหสัมพันธ์

corrplot(mtcars_cor$r, type="upper", p.mat = mtcars_cor$P, 
         sig.level = 0.1, insig = "blank")

ส่วนการศึกษาเรื่องกราฟสหสัมพันธ์เพิ่มเติมให้เป็นหน้าที่ของผู้อ่านต้องเรียนรู้ด้วยตัวเองต่อไป

สัมประสิทธ์สหสัมพันธ์บางส่วน(partial correlation)

สัมประสิทธ์สหสัมพันธ์บางส่วน จะถูกใช้มากในการศึกษาเรื่องการวิเคราะห์การถดถอยขั้นสูงต่อไป ถ้าจะหาค่านี้ จะทำการคำนวณเมื่อตัวแปรสุ่มที่เราสนใจศึกษามีขนาดตั้งแต่ 3 มิติขึ้นไปและจะทำกาคำนวณหาค่าสหสัมพันธ์เป็น คู่ๆ ไปทั้งหมด แล้วจึงค่อยมาคำนวณด้วยสหสัมพันธ์บางส่วนทีหลัง ด้วยสูตรการคำนวณนี้ ถ้าต้องแปรสุ่มที่สุ่มเราสน ใจมีขนาดมิติ \(n\) \(X_1,X_2,\cdots,X_n)\) ถ้าคำนวณค่าสหสัมพันธ์บางของ \(r_{ij|k}\) ส่วนของ \(X_i\) กับ \(X_j\) โดยให้ \(X_k\) คือ \[r_{i j | k}=\frac{r_{i j}-r_{i k} r_{j k}}{\sqrt{1-r_{i k}^{2}} \sqrt{1-r_{j k}^{2}}}\] จากข้อมูล mtcars ถ้าพิจารณาเฉพาะ 4 ตัวแปรแรกคือ mpg cyl disp และ hp จะได้

cor(mtcars[,1:4])
##             mpg        cyl       disp         hp
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475
## disp -0.8475514  0.9020329  1.0000000  0.7909486
## hp   -0.7761684  0.8324475  0.7909486  1.0000000

ดังนั้น สมมัติต้องการค่า \(r_{mpg,cyl|hp}\) จะได้ว่า \[r_{mpg,cyl|hp}=\frac{r_{mpg,cyl}-r_{mpg,hp} r_{cyl,hp}}{\sqrt{1-r_{mpg,hp}^{2}} \sqrt{1-r_{cyl,hp}^{2}}}\]

คำนวณอาร์ด้วยได้ดังนี้

# หาสหสัมพันธ์ระหว่าง mpg กับ cyl
cor_mpg.cyl<-cor(mtcars[,1:4])[1,2]
# หาสหสัมพันธ์ระหว่าง mpg กับ hp
cor_mpg.hp <-cor(mtcars[,1:4])[1,4]
 # หาสหสัมพันธ์ระหว่าง cyl กับ hp
cor_cyl.hp <-cor(mtcars[,1:4])[2,4]
# # หาสหสัมพันธ์บางส่วนของ mpg กับ cyl ภายใต้ hp
(cor_mpg.cyl-cor_mpg.hp*cor_cyl.hp)/(sqrt(1-cor_mpg.hp^2)*sqrt(1-cor_cyl.hp^2))
## [1] -0.5897431

ในบทนี้จะแนะนำการสูตรที่ใช้ในการคำนวณเท่านั้น

บทอื่นๆ