1 課程目標

本週上課將介紹兩個變數的相關，包括兩個類別變數以及兩個連續變數。配合之前討論的視覺化，可以描述兩個變數之間的關係，例如：

with(anscombe, cor(x1, y1))
## [1] 0.8164205

2 類別變數

社會科學經常分析類別資料。類別資料之間的關聯性有各種分析方法，原則為減少實際觀察到的數值與理想上或者不知道另一個變數時之間誤差的比例，也就是 Proportional Reduction of Error (PRE)。PRE可表示為： \(PRE=\frac{E_{1}-E{2}}{E_{1}}\)

\(E_{1}\)表示不知道自變項(X)時、預測依變項(Y)的誤差。\(E_{2}\)表示知道X而預測Y時的誤差。

PRE越大，表示X與Y的關聯性越高。

根據資料的特性，使用不同的關聯性指標來表示兩個變數之間的相關程度：

  • 名目變數：Goodman-Kruskal’s \(\lambda\)(lambda)、 Goodman-Kruskal’s \(\tau\)(tau)

    • Goodman-Kruskal’s \(\lambda\)=\(\frac{S-R}{N-R}\)其中S代表每一列之中，數目最高的格子的總和，R代表總數最高的列的數值。N代表全部總數。最小是0，最大是1。

    Goodman-Kruskal’s \(\tau\)來自於不知道自變項(X)時、預測依變項(Y)的誤差，減去不知道依變項(Y)而預測自變項(X)時的誤差，再除以不知道自變項(X)時、預測依變項(Y)的誤差。

  • 順序變數：Goodman-Kruskal’s \(\gamma\)應用在分析對稱（也就是不分自變數與依變數）的關聯、Somers’D應用在不對稱的關係

    • 例如：

    學生 數學 英文
    Ariel 4 2
    Ben 3 3
    Carl 2 5
    Danny 1 4
    Elaine 5 1

    在上面的表格計算同序、不同序的數目，代入公式計算關聯係數。

  • 雙變數：卡方檢定：\(\mathcal{X}^2\)
  • 直方圖可以表示關聯性的視覺化。

    • 猜測

  • 我們經常要預測某個隨機變數的值，例如受訪者常用那一種語言。
  • 如果我們亂猜，也許會猜他們是說「閩南語」，因為大多數「本省閩南人」佔了七成。
  • 但是如果知道某些受訪者住在新竹縣，或許我們會修改答案為「客語」，因為許多客家人居住在當地。
  • 兩個變數之間的關聯性越強，表示我們可以用某一變數來正確預測另一個變數。

    • 卡方檢定

  • 卡方檢定計算隨機變數的分佈與理論上分佈之間的差距以及對應的機率值\(p\)。如果差距太大則拒斥觀察到的變數分佈等於理論分佈的假設，也就是隨機產生的。一般設定\(p<0.05\)。
  • 可分為適合度檢定(goodness of fit)以及獨立性檢定(test of independence)。

  • 為何要用\(\mathcal {X}^2\)分佈？卡方分佈適用於數個獨立變數的平方和，最初是研究變異數的學者發現。它是\(\Gamma\)分佈的一種。因此卡方分佈的最重要參數是自由度。

    卡方分佈

    自由度

  • 自由度的意義為在進行統計時有多少可以變化的參數。例如當我們已知\(n\)個數的平均數，只要知道\(n-1\)個數的數值時，最後一個數的大小會自動決定。
  • \(t\)分佈、\(\mathcal{X}^2\)分佈都會因為樣本大小不同而有不同的分佈。如果是進行適合度檢定，自由度是\(k-1\)，\(k\)是類別數。如果是進行雙樣本的獨立檢定，自由度是：\((n_{c}-1)\cdot (n_{r}-1)\)。因此我們當知道某幾個格子之後就可以得知其他的格子。例如一個2\(\times\) 2表格，只要知道其中一格就知道其他格子。

    • 適合度檢定

  • 適合度檢定比較觀察到的分佈與平均分配或者是已知的參數分佈。
  • 例如100位樣本中，52位是女性其他是男性，我們可以比較理論上男女各半的分佈，或者是比較男女比為47:53的假設，看看觀察到的分佈是否符合男女各半還是47:53這個假設。

    • 獨立性檢定

  • 獨立性檢定以列聯表決定兩個名目變數之間是否有統計上的相關。
  • 如果有兩個變數\(A\),\(B\)，當\(A\)從類別1變動到類別2，\(B\)也跟著改變，可以判斷這兩個變數不互相獨立。
  • 因為要用機率分佈進行假設檢定，所以列聯表內的期待值不能太小。
  • 樣本數如果太小，\(\mathcal {X}^2\)可能被高估，太大可能被低估。
  • 不表示兩個變數之間的關聯性大小，只呈現相關與否的證據。
  • \({H}_{0}:\)兩個變數獨立;\({H}_{1}:\)兩個變數不獨立
  • 期待值計算方式：

    • \(\hat{n_{ij}}=(RT\times CT)/ N\)

  • \(RT\)：包含該列的邊際的總和
  • \(CT\)：包含該列的邊際的總和
  • \(N\)：觀察值總數

    • 計算期待值 

    library(kableExtra)
dt <-data.frame(a=c("","" ,"X","" , "Total"),
                b=c("" ,"" , 1 ,    2,"" ),
                C=c("Y", 1, "f11", "f12", "R1"),
d=c("", 2, "f21"  , "f22" , "R2" ),
e=c("Total" , "" , "C1"  , "C2" , "N"  )
            )
colnames(dt)<-NULL

dt %>%
  kable("html") %>%
  kable_styling()
    Y Total
    1 2
    X 1 f11 f21 C1
    2 f12 f22 C2
    Total R1 R2 N
    根據上面表格計算期待值：
  • \(E11=\frac{R1\times C1}{N}\)
  • \(E12=\frac{R1\times C2}{N}\)
  • \(E21=\frac{R2\times C1}{N}\)
  • \(E22=\frac{R2\times C2}{N}\)

    • 計算卡方值公式

  • 觀察值以\(O_{ij}\)表示，期待值則以\(E_{ij}\)表示。
  • 兩者的差距越大，表示觀察值距離隨機分佈越遠，也就表示觀察值越不可能是兩者無關的隨機分佈。
  • \(\mathcal {X} ^2=\sum_{i}\sum_{j}\frac{(n_{ij}-\hat{n_{ij}})^2}{\hat{n_{ij}}}\)

    • 或者是

  • \(\mathcal{X}^2=\sum_{i}\sum_{j}\frac{(O_{ij}-E_{ij})^2}{E_{ij}}\)
  • 當交叉列表應用在\(2\times 2\)表格時，建議用Yate’s correction係數
  • \(\mathcal{X}^2=\sum_{i}\sum_{j}\frac{(|O_{ij}-E_{ij}|-0.5)^2}{E_{ij}}\)
  • 卡方分布的期望值等於自由度。而且自由度越接近30，越接近常態分佈。

    • 適合度檢定

    適合度檢定指的是變數的分佈是否符合特定分佈，例如單一分佈，也就是每一類別有一樣的觀察值數目。期待值計算方式為： \(E_{k}=\frac{N}{n}\)
  • \({H}_{0}\)為觀察值平均分佈在每一類別，如果拒斥這個假設，表示各類別之間的差距的確存在而非隨機分佈。
  • 自由度為\(k-1\)。因為觀察值的總和為\(N\)，只要知道\(k-1\)個類別的觀察值就知道最後一個類別的觀察值。
  • 另一方面，因為單一分佈的總和是固定的，所以在這一個限制下理論上的數值變化就喪失一個自由度，寫成\(p=1\)。\(k-p=k-1\)就是自由度。
  • 例如：抽出100位1990年之後的出生的受訪者，觀察到男性為52人、女性為48人，請問這是平均分佈嗎？
  • 計算為\(\mathcal{X}^2=\frac{(52-50)^2}{50}+\frac{(48-50)^2}{50}=0.16\)
  • 右尾卡方值參考表如下：
    • df
    1 2 3 4 5
    0.05 3.84 5.99 7.82 9.49 11.07
    0.01 6.63 9.21 11.34 13.28 15.09

    也可以用qchisq這個指令，設定p值與自由度之後，即可計算應該大於多少卡方值：

    qchisq(0.95, df=1)
    ## [1] 3.841459

    反過來，pchisq可以輸出一定的自由度下卡方值的機率，例如：

    1-pchisq(2.9, 1)
    ## [1] 0.08857955
  • 回到我們的例子，當卡方值為0.16，自由度為1時：可計算機率值等於0.68： 
    • 1-pchisq(0.16, 1)
    ## [1] 0.6891565
  • 查表可知機率值為0.68，因此無法拒斥資料來自母體的假設。
  • 我們也可以用卡方檢定測試觀察值是否符合特定的分佈。例如我們經常做加權檢定，想要知道是否樣本符合母體，虛無假設是樣本分佈等於母體分佈，以性別而言，自由度為1。假設母體的分佈是：
  • 男：\(\frac{9003150}{9003150+9184065}=0.495\)
  • 女：\(\frac{9184065}{9003150+9184065}=0.505\)
  • 假設我們觀察到1,067人，男女分佈如下，可計算對應的理論值：
    • 性別 個案數 理論值
    男 531 528.165
    女 539 538.835
    \(\mathcal{X}^2=\frac{(531-528.165)^2}{528.165}+\frac{(539-538.835)^2}{538.835}=0.015+0=0.015<3.84\)
  • 結論：我們無法拒斥\({H}_{0}\)：觀察值與母體有相同分佈在每一類別這個假設。

    • 如果我們想知道政黨支持與婚姻狀態是否互相獨立？以下表資料為例：

    婚姻 Total
    已婚 其他
    政黨支持 國民黨 141 77 218
    民進黨 91 72 163
    Total 232 149 381
  • 先寫語法轉出表格 
    • marri<-c(141,91,77,72)
marripid<-matrix(marri,2)
marripid
    ##      [,1] [,2]
## [1,]  141   77
## [2,]   91   72
  • 計算邊際機率 
    • rt<-margin.table(marripid,1)
ct<-margin.table(marripid,2)
newtable<-cbind(marripid,rt)
rbind(newtable,ct)
    ##             rt
##    141  77 218
##     91  72 163
## ct 232 149 232
  • 計算期待值 
    • n=sum(marri)
e11<-rt[1]*ct[1]/n; e21<-rt[2]*ct[1]/n
e12<-rt[1]*ct[2]/n;e22<-rt[2]*ct[2]/n
 
expec <- matrix(c(e11,e21,e12,e22),2)
print(expec)
    ##           [,1]     [,2]
## [1,] 132.74541 85.25459
## [2,]  99.25459 63.74541
  • 最後計算卡方值 
    • dt<-c((marripid[1,1]-expec[1,1])^2/e11,
            (marripid[2,1]-expec[2,1])^2/e21,
             (marripid[1,2]-expec[1,2])^2/e12,
             (marripid[2,2]-expec[2,2])^2/e22)
dt
    ## [1] 0.5133007 0.6865003 0.7992333 1.0689132
    chisq<-sum(dt)
chisq
    ## [1] 3.067948

    如果用R的函數計算，可以用chisq.test這個函數，得到Pearson’s 卡方檢定：

    marri<-c(141,91,77,72)
marripid<-matrix(marri,2)
chisq.test(marripid)
    ## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  marripid
## X-squared = 2.7075, df = 1, p-value = 0.09988

    因為Yates(1934)建議需要考慮討論2\(\times2\)表格的時候，卡方分配逼近不連續的樣本分配的問題，所以重新計算卡方值為\(\sum\frac{(|O-E|-0.5)^2}{E}\)

    # YATES <- min(0.5, abs(x - E))
dt<-c((abs(marripid[1,1]-expec[1,1])-0.5)^2/e11,
      (abs(marripid[2,1]-expec[2,1])-0.5)^2/e21,
      (abs(marripid[1,2]-expec[1,2])-0.5)^2/e12,
      (abs(marripid[2,2]-expec[2,2])-0.5)^2/e22)
dt
    ## [1] 0.4530003 0.6058532 0.7053428 0.9433419
    chisq<-sum(dt)
chisq
    ## [1] 2.707538

    可畫圖表示：

    z=curve(dchisq(x, df=1), from=0, to=6)
abline(v = chisq, lty = 3, lwd=1.5, col='red2')

    標準化殘差

  • 觀察值與期待值之間的差距稱為殘差。殘差越大表示觀察到的數目遠高於應該有的期待值，也表示往這個方向變動。
  • 標準化殘差就是考慮每一格所在的行與列所佔的比例進行標準化。
    • \(z_{ij}=\frac{O_{ij}-E_{ij}}{\sqrt{E_{ij}(1-{\mathrm {row\hspace{.2em} prop.}})(1-{\mathrm {column\hspace{.2em} prop.}})}}\)
  • 標準化殘差越大表示那一格的觀察值超過當虛無假設成立時的期待值越多，也就表示兩個變數的關聯越高。
  • 有一筆觀察值為1067的資料，教育程度分佈如下表，依據母體資料求出期待值以及殘差：。
    • 教育程度 理論值 殘差
    小學及以下 192.66 -6.66
    國中 228.42 -81.42
    高中 216.56 94.43
    專科 202.64 -65.64
    大學及以上 226.68 59.31
  • 由上可知高中、大學及以上教育程度這兩群人在樣本中太多。

    • 卡方檢定注意事項

  • 每一格的期待值最好大於5。
  • 卡方值不表示關聯性的大小。
  • \(\mathcal{X}^2\)代表卡方分佈在右邊的機率，假設\({H}_{0}\)為真。

    • 3 相關係數

  • 連續變數之間的關係可先用散佈圖表示，例如： 
    • library(ggplot2)
ggplot(anscombe, aes(x=x2, y=y1)) +
    geom_point(size=3, col="darkgreen") +
    labs(subtitle="Data: anscombe") +
    theme_bw()
    Scatter Plot

    Scatter Plot

    可以加上迴歸線表現兩個變數的相關：

    ggplot(anscombe, aes(x=x2, y=y1)) +
    geom_point(col='darkgreen', size=3) +
    geom_smooth(method="lm", se=F, col='blue')
    Scatter Plot with Regression Line

    Scatter Plot with Regression Line

    也可以用非線性迴歸線（Locally Weighted Scatterplot Smoothing, LOESS）表示兩者的相關：

    ggplot(anscombe, aes(x=x2, y=y1)) +
   geom_point(col='darkgreen', size=3) +
   geom_smooth(method="loess", se=F) +
   theme_bw()
    Scatter Plot with LOESS Line

    Scatter Plot with LOESS Line

  • 兩個變數之間的相關係數，計算方式為：
    • \(\mathcal{r=\frac{n(\Sigma xy)-(\Sigma x)(\Sigma y)}{\sqrt{(n\Sigma x^2-(\Sigma x)^2)(n\Sigma y^2-(\Sigma y)^2)}}}\)
  • 1代表正相關關係最高，X增加一個單位，Y也增加一個單位，或者Y增加一個單位，X也增加一個單位
  • -1代表負相關關係最高，X減少一個單位，Y增加一個單位
  • 0代表兩者沒有相關

    • 如果Y是X加上一個單位：

    df <- tibble(X=rnorm(20,0,1), Y=X+1)
with(df, cor(X, Y))
    ## [1] 1
    ggplot(df, aes(x=X, y=Y)) +
    geom_point(fill='brown', size=5, shape=1) +
    geom_smooth(method="lm", se=F, col='blue')

    如果把Y改為\(X^2\)：

    df <- tibble(X=rnorm(20,0,1), Y=X^2)
with(df, cor(X, Y))
    ## [1] 0.07247545
    ggplot(df, aes(x=X, y=Y)) +
    geom_point(fill='brown', size=5, shape=1) +
    geom_smooth(method="lm", se=F, col='blue')

    回到剛剛的例子：

    with(anscombe, cor(x2, y2))
    ## [1] 0.8162365
  • 一般用Pearson Product Moment Correlation (PPMC)或者是correlation coefficient \(r\) 代表相關係數，但是correlation of determination \(r^2\)，代表解釋一個變數的程度，例如\(r^2=0.85\)，代表有\(85\%\)的變異數被解釋。

    • 4 作業

    1. 請針對mtcars這筆資料的wt以及mpg繪製散佈圖，並且計算相關係數：

    2. 請使用PP0797B2.sav這筆資料，然侯顯示年齡(age)與統獨立場(tondu)之間的交叉列表，再進行卡方檢定，請問這兩個變數互相獨立嗎？

    3.根據上述的資料，請問年齡的分佈是否與母體一致？假設母體的表格如下：
    Category 20至29歲 30至39歲 40至49歲 50至59歲 60至69歲
    N 351 488 610 374 235
    % 17.05 23.71 29.64 18.17 11.41

    4. 請問以下這個表格中，部門與是否提前退休的人數有無統計上相關？

    tab <- as.table(rbind(Yes=c(26,23,18,9), No=c(6,7,14,23)))
names(dimnames(tab)) <- c("Retired", "Department")
tab
    ##        Department
## Retired  A  B  C  D
##     Yes 26 23 18  9
##     No   6  7 14 23

    5. 請先畫散佈圖表現UsingR這個套件中的nym.2002這筆資料，裡面的age以及time之間的關係，並且計算兩者的相關程度。