社會科學統計方法
相關分析與卡方檢定
蔡佳泓
4/15/2022
1 課程目標
ggstatsplot::ggbarstats(
data = ISLR::Carseats,
x = Urban,
y = US,
#sampling.plan = "jointMulti",
title = "Store Location and US/Non-US",
xlab = "US or Not",
legend.title = "City",
#ggtheme = hrbrthemes::theme_ipsum_pub(),
ggplot.component = list(scale_x_discrete(
guide = guide_axis(n.dodge = 2))),
palette = "Set2",
messages = FALSE
)
Figure 1.1: 店面位置與是否在美國的柱狀圖
在總共258個美國店面當中，有72%位在都市，28%在郊區。而在142間美國以外店面當中，有有68%是位在都市，32%不是。卡方(\(\chi^2\))檢定顯示，\(p>0.05\)，表示是否位在美國以及是否位在都市這兩個變數，在統計上沒有達到顯著相關。這裡注意x=US/non-US，色塊表示City是Yes/No分別佔的比例，例如右邊的綠色指的是位於都市的店面佔了美國的68%的店面。 如果只是列表，可以用
janitor這個套件列出交叉列表，得到相同的結果，括號內是家數：
tab.carseats <- ISLR::Carseats |> janitor::tabyl(Urban, US) |>
janitor::adorn_percentages("col") |>
janitor::adorn_pct_formatting(digits = 2) |>
janitor::adorn_ns()
tab.carseats
## Urban No Yes
## No 32.39% (46) 27.91% (72)
## Yes 67.61% (96) 72.09% (186)- 也可以用
flextable這個套件，列出個案數、每一個格子佔的全部比例以及行、列為百分之百的比例如表 1.1。注意解讀每一個數字：
ftab <- ISLR::Carseats |> flextable::proc_freq('Urban', 'US')
ftab <- flextable::set_caption(ftab, "店面位置與是否在美國的交叉列表")
ftabUrban | US | |||
|---|---|---|---|---|
No | Yes | Total | ||
No | Count | 46 (11.5%) | 72 (18.0%) | 118 (29.5%) |
Mar. pct (1) | 32.4% ; 39.0% | 27.9% ; 61.0% | ||
Yes | Count | 96 (24.0%) | 186 (46.5%) | 282 (70.5%) |
Mar. pct | 67.6% ; 34.0% | 72.1% ; 66.0% | ||
Total | Count | 142 (35.5%) | 258 (64.5%) | 400 (100.0%) |
(1) Columns and rows percentages | ||||
2 類別變數的相關
\[PRE=\frac{E_{1}-E{2}}{E_{1}}\]
假設某一個縣市的小學教師之中有6成是研究所畢業，\(m=0.6\)。我們有\(100\times (1-0.6)\%=40\%\)誤認某一位老師是研究所學歷的可能。\(E_{1}=0.4\)
假設經過調查後，發現40歲以下的老師之中，有80%擁有研究所學歷，所以\(E_{2}=0.2\)。
知道老師的年齡會降低犯錯的程度，計算PRE=\(\frac{0.4-0.2}{0.4}=50\%\)
- 名目變數： \(\lambda\)
- 順序變數： \(\Gamma\)、Kendall’s\(\tau\)-b、Somers’ D、Spearman rank coefficient
- 檢定變數是否獨立：\(X^2\)
2.1 名目變數
2.1.1 \(\lambda\)範例
rawd0<-matrix(c("","","新聞台","","",
"","三立","民視","TVBS","總和",
"frequency","38","32","25","95"),
ncol=3, nrow=5)
kable(rawd0, booktabs=T, caption="觀看新聞台的次數分配") %>%
kable_styling(full_width=F, bootstrap_options = "striped", font_size=14) | frequency | ||
| 三立 | 38 | |
| 新聞台 | 民視 | 32 |
| TVBS | 25 | |
| 總和 | 95 |
rawd<-matrix(c("","","新聞台","","",
"","三立","民視","TVBS","總數",
"軍公教","10","5","15","30",
"私部門職員","15","10","5","30",
"勞工","5","15","5","25",
"農林漁牧","8","2","0","10",
"總數","38","32","25","95"),
ncol=7, nrow=5)
kable(rawd, booktabs=T, caption="職業與新聞台之交叉表") %>%
kable_styling(full_width=F, bootstrap_options = "striped", font_size=14) %>%
add_header_above(c(""," ", "職業" = 5)) | 軍公教 | 私部門職員 | 勞工 | 農林漁牧 | 總數 | ||
| 三立 | 10 | 15 | 5 | 8 | 38 | |
| 新聞台 | 民視 | 5 | 10 | 15 | 2 | 32 |
| TVBS | 15 | 5 | 5 | 0 | 25 | |
| 總數 | 30 | 30 | 25 | 10 | 95 |
Ld<-here::here('data','lambda.txt')
dt<-read.table(Ld, sep=';', header=T)
kable(dt, booktabs=T, caption="四種職業看新聞台之正確與錯誤") %>%
kable_styling(full_width=F, bootstrap_options = "striped", font_size=14) | Rating.of..Department | Modal.Category | Right.Guesses | Wrong.Guesses |
|---|---|---|---|
| 軍公教人員 (N=30) | TVBS | 15 | 15 |
| 私部門職員 (N=30) | 三立 | 15 | 15 |
| 勞工 (N=25) | 民視 | 15 | 10 |
| 農林漁牧 (N=10) | 三立 | 8 | 2 |
| 總數(N=95) | 53 | 42 |
- 計算\(\lambda\)需要知道\(E_{1}\)與\(E_{2}\)。\(E_{1}\)是單一變數的邊際機率，也就是不考慮其他事件時的次數分佈時，我們會犯錯的機會，等於：
\[\begin{align*} E_{1}&=N_{\text{total}}-N_{\text{mode}\hspace{.4em}\text{of}\hspace{.4em}\text{dependent}\hspace{.4em}\text{variable}}\\ &=95-38\\ &=57 \end{align*}\]
- \(E_{2}\)則是知道另一個變數與原來變數的交集時，我們會犯錯的機會，等於：
\[\begin{align*} E_{2} & = \sum(N_{\text{category}}-N_{\text{mode}\hspace{.4em} \text{of}\hspace{.4em}\text{category}}) \\ & = 95-53 \\ & = 42 \end{align*}\]
- 根據PRE的公式，\(\lambda\)值求得如下：
\[\begin{align*} \frac{E_{1}-E_{2}}{E_{1}} & = \frac{57-42}{57} \nonumber \\ & = \frac{15}{57} \nonumber \\ & = 0.2631579 \end{align*}\]
\[\lambda_{symmetric}=\frac{\sum Y_{x, modal}+\sum X_{y, modal}}{2N-(Y_{modal}+X_{modal})}\]
lambdatab=cbind(c(10,5,15),c(15,10,5),c(5,15,5),
c(8,2,0))
lambdatab
## [,1] [,2] [,3] [,4]
## [1,] 10 15 5 8
## [2,] 5 10 15 2
## [3,] 15 5 5 0
DescTools::Lambda(lambdatab, direction = "row")
## [1] 0.2632- 以另外一個矩陣資料為例，因為我們不知道A或者B是依變數，所以計算對稱與不對稱的\(\lambda\)值：
M <- as.table(cbind(c(1768,946,115), c(807,1387,438), c(189,746,288), c(47,53,16)))
dimnames(M) <- list(paste("A", 1:3), paste("B", 1:4))
M
## B 1 B 2 B 3 B 4
## A 1 1768 807 189 47
## A 2 946 1387 746 53
## A 3 115 438 288 16
# direction default is "symmetric"
DescTools::Lambda(M)
## [1] 0.2076
DescTools::Lambda(M, conf.level=0.95)
## lambda lwr.ci upr.ci
## 0.2076 0.1872 0.2281
DescTools::Lambda(M, direction="row")
## [1] 0.2241
DescTools::Lambda(M, direction="column")
## [1] 0.19242.2 順序變數
\[ \Gamma= \frac{N_{s}-N_{d}}{N_{s}+N_{d}} \]
TT<-matrix(c(""," ","Y"," ",
"","1","2","3","1","f11","f21","f31",
"2","f12","f22","f32",
"3","f13","f23","f33"),
ncol=5, nrow=4)
knitr::kable(TT, 'html',booktabs=T, caption="3X3表格") %>%
kable_styling(c("striped", "bordered")) %>%
add_header_above(c(""," ", "X" = 3)) | 1 | 2 | 3 | ||
| 1 | f11 | f12 | f13 | |
| Y | 2 | f21 | f22 | f23 |
| 3 | f31 | f32 | f33 |
- Ns=f11(f22+f23+f32+f33)+f12(f23+f33)+f21(f32+f33)+f22(f33)
- Nd=f13(f22+f21+f32+f31)+f12(f21+f31)+f23(f32+f31)+f22(f31)
- Ty=f11(f12+f13)+f12(f13)+f21(f22+f23)+f22(f23)+f31(f32+f33)+f32(f33)
- Tx=f11(f21+f31)+f21(f23)+f12(f22+f32)+f22(f32)+f13(f23+f33)+f23(f33)
TT<-matrix(c("","","教育程度","", "",
"","國中及以下", "高中五專","大學以上","總和",
"變差", "35", "76", "22", "133",
"差不多"," 73", "134", "56", "263",
"變好","30","54","30","114",
"總和","138", "264","108","150"),
ncol=6, nrow=5)
knitr::kable(TT, 'html', caption='學歷與對經濟的看法') %>%
kable_styling(c("striped", "bordered")) %>%
add_header_above(c(""," ", "經濟變好或是不好" = 3,"")) | 變差 | 差不多 | 變好 | 總和 | ||
| 國中及以下 | 35 | 73 | 30 | 138 | |
| 教育程度 | 高中五專 | 76 | 134 | 54 | 264 |
| 大學以上 | 22 | 56 | 30 | 108 | |
| 總和 | 133 | 263 | 114 | 150 |
- 在上面的表格計算同序、不同序的數目，代入公式計算關聯係數。
ns<-35*(134+56+54+30)+73*(54+30)+76*(56+30)+134*30
nd<-30*(76+22+134+56)+54*(22+56)+73*(76+22)+134*22
gamma <- (ns-nd)/(ns+nd)
cat('Gamma=', gamma)
## Gamma= 0.06752- 用
R計算Goodman-Kruskal’s \(\Gamma\)：
Gtable<-cbind(c(35,76,22),c(73,134,56),
c(30,54,30))
DescTools::GoodmanKruskalGamma(Gtable, direction="column")
## [1] 0.067522.2.1 Kendall’s \(\tau-b\)
\[\tau-\hspace{-.3em}b=\frac{N_{s}-N_{d}}{\sqrt{N_{s}+N_{d}+T_{x}}\cdot\sqrt{N_{s}+N_{d}+T_{y}}}\]
DescTools::GoodmanKruskalTau(Gtable)
## [1] 0.003948- 手動計算\(\tau-b\)如下：
tx<-35*(76+22)+76*22+73*(134+56)+134*56+30*(54+30)+54*30
ty<-35*(73+30)+73*30+76*(134+54)+134*54+22*(56+30)+56*30
taub<-(ns-nd)/sqrt((ns+nd+tx)*(ns+nd+ty))
cat('Tau-b=', taub)
## Tau-b= 0.04156- 用\(\tt{KendallTauB()}\)驗算：
DescTools::KendallTauB(Gtable, direction="column")
## [1] 0.04156- 跟\(\texttt{DescTools::GoodmanKruskalTau}\)進行分析得到的結果不太相同：
DescTools::GoodmanKruskalTau(Gtable)
## [1] 0.003948- 我們可以用\(\texttt{janitor::tabyl}\)以及\(\texttt{DescTools}\)分析
R的資料。先用\(\texttt{janitor::tabyl}\)顯示交叉表：
library(DescTools); library(janitor); data("d.pizza")
d.pizza %>% janitor::tabyl(area, operator) %>%
adorn_totals() %>%
adorn_percentages('row') %>%
adorn_pct_formatting(digits = 2) %>%
adorn_ns()
## area Allanah Maria Rhonda NA_
## Brent 32.28% (153) 32.28% (153) 35.23% (167) 0.21% (1)
## Camden 35.76% (123) 31.40% (108) 31.69% (109) 1.16% (4)
## Westminster 23.36% (89) 32.02% (122) 43.83% (167) 0.79% (3)
## <NA> 20.00% (2) 50.00% (5) 30.00% (3) 0.00% (0)
## Total 30.36% (367) 32.09% (388) 36.89% (446) 0.66% (8)- 可以看到列百分比，而每一格的個案數在括號內，而且還有遺漏值。如果不想看到遺漏值，可以設定
show_na=FALSE ：
library(DescTools); library(janitor); data("d.pizza")
d.pizza %>% janitor::tabyl(area, operator, show_na=FALSE) %>%
na.omit() %>%
adorn_totals() %>%
adorn_percentages('row') %>%
adorn_pct_formatting(digits = 2) %>%
adorn_ns()
## area Allanah Maria Rhonda
## Brent 32.35% (153) 32.35% (153) 35.31% (167)
## Camden 36.18% (123) 31.76% (108) 32.06% (109)
## Westminster 23.54% (89) 32.28% (122) 44.18% (167)
## Total 30.65% (365) 32.16% (383) 37.20% (443)- \(\texttt{janitor::tabyl}\)並沒有統計值。用\(\texttt{DescTools::Desc}\)重做一次，可以顯示更多統計結果，包括後面學到的卡方檢定以及G-test，這兩個統計都是檢驗兩個變數之間是否獨立，還有兩個圖形，分別顯示列與行的條件機率，列的條件機率在左邊的圖，行的條件機率在右邊的圖：
library(DescTools); library(janitor); data("d.pizza")
DescTools::Desc(area ~ operator, data=d.pizza)
## ------------------------------------------------------------------------------
## area ~ operator (d.pizza)
##
## Summary:
## n: 1'191, rows: 3, columns: 3
##
## Pearson's Chi-squared test:
## X-squared = 18, df = 4, p-value = 0.001
## Log likelihood ratio (G-test) test of independence:
## G = 18, X-squared df = 4, p-value = 0.001
## Mantel-Haenszel Chi-squared:
## X-squared = 8.7, df = 1, p-value = 0.003
##
## Contingency Coeff. 0.122
## Cramer's V 0.087
## Kendall Tau-b 0.073
##
##
## area Brent Camden Westminster Sum
## operator
##
## Allanah freq 153 123 89 365
## perc 12.8% 10.3% 7.5% 30.6%
## p.row 41.9% 33.7% 24.4% .
## p.col 32.3% 36.2% 23.5% .
##
## Maria freq 153 108 122 383
## perc 12.8% 9.1% 10.2% 32.2%
## p.row 39.9% 28.2% 31.9% .
## p.col 32.3% 31.8% 32.3% .
##
## Rhonda freq 167 109 167 443
## perc 14.0% 9.2% 14.0% 37.2%
## p.row 37.7% 24.6% 37.7% .
## p.col 35.3% 32.1% 44.2% .
##
## Sum freq 473 340 378 1'191
## perc 39.7% 28.5% 31.7% 100.0%
## p.row . . . .
## p.col . . . .
## 
Figure 2.1: DescTools交叉列表圖
Kendall’s \(\tau-b\) 在樣本數大於10的時候，抽樣分佈會形成標準常態分配，也就是可以計算Z值，並檢驗假設。
Kendall’s \(\tau-b\) 與 Spearman rank coefficient一樣， 可以用來計算兩個順序變數之間的相關程度，-1代表負相關，+1表示正相關。
set.seed(0000)
n <- 200 # Number of samples
residence <- factor(sample(c("Urban", "Rural"), n, replace = TRUE))
residence <- factor(residence, levels= c("Urban", "Rural"))
color <- factor(sample(c("Red", "Blue", "Green", "Darkgreen"), n, replace = TRUE))
color <- factor(color, levels = c("Red", "Blue", "Green", "Darkgreen"))
residencecolortable <- table(residence, color)
vcd::mosaic(color ~ residence, data = data.frame(residence, color), shade = TRUE)
Figure 2.2: 居住地與喜歡的顏色相關圖1
set.seed(02138)
n <- 150 # Number of samples
age <- c(factor(sample(c("20-35 yrs", "36-60 yrs"), n, replace = TRUE)), factor(rep("20-35 yrs", 10)), factor(rep('36-60 yrs', 40)))
age <- factor(age, levels= c("20-35 yrs", "36-60 yrs"))
color2 <- c(factor(sample(c("Red", "Blue", "Green", "Darkgreen"), n, replace = TRUE)), factor(rep("Darkgreen", 50)))
color2 <- factor(color2, levels = c("Red", "Blue", "Green", "Darkgreen"))
agecolortable <- table(age, color2)
vcd::mosaic(color2 ~ age, data = data.frame(age, color2), shade = TRUE)
Figure 2.3: 性別與喜歡的顏色相關圖2
- 計算\(\tau-b\)：
cat('Taub of gender and color:', DescTools::KendallTauB(residencecolortable, direction="row"))
## Taub of gender and color: -0.0393
cat('Taub of age and color:', DescTools::KendallTauB(agecolortable, direction="row"))
## Taub of age and color: 0.11792.2.2 Somers’ D
Somers’ D 也是計算兩個順序變數的相關程度，\(-1\leq D\leq 1\)。+1表示兩個變數正相關，一個變數的順序上升，另一個變數的順序也上升。-1表示兩個變數負相關，一個變數的順序下降，另一個變數的順序也下降。0表示兩個變數的順序變動沒有特定的模式。
跟\(\tau-b\)不同的地方是，Somer’s D考慮兩個變數之間所有的配對，包括順序正確以及順序相反的的關係。\(\tau-b\)只計算一致與不一致的順序之間的差異，並且針對平手的部分作出調整。
\(-1\leq \rm{Somers'}\hspace{0.2em} \rm{D}\leq 1\)。-1就是所有的配對都相反，1代表所有的配對都一致。
以下模擬X, Y1, Y2三個變數，X, Y1配對一致，X, Y2則完全不一致，得到的Delta值為1以及-1。
X<-c(rep(1,10), rep(2, 10), rep(3,15), rep(4,5))
Y1<-c(rep(1,10), rep(3, 10), rep(4,15), rep(6,5))
D1<-DescTools::SomersDelta(X,Y1)
Y2<-c(rep(6,10), rep(4, 10), rep(3,15), rep(1,5))
D2<-DescTools::SomersDelta(X,Y2)\[\rm{Delta}=(N_{s}-N_{d})/(N_{s}+N_{d}+T_{y})\]
★試計算以下的表格2.6的Somers’ D:
tab <- as.table(rbind(c(26,26,23,18,9),c(6,7,9,14,23)))
knitr::kable(tab,format = 'simple',
booktabs=T,caption="交叉列表範例")| A | B | C | D | E | |
|---|---|---|---|---|---|
| A | 26 | 26 | 23 | 18 | 9 |
| B | 6 | 7 | 9 | 14 | 23 |
- 用\(\texttt{DescTools}\)套件中的\(\texttt{SomersDelta}\)函式計算不同方向的Delta：
# Somers' D C|R
DescTools::SomersDelta(tab, direction="column", conf.level=0.95)
## somers lwr.ci upr.ci
## 0.4427 0.2786 0.6068
# Somers' D R|C
DescTools::SomersDelta(tab, direction="row", conf.level=0.95)
## somers lwr.ci upr.ci
## 0.2569 0.1593 0.3546- Goodman-Kruskal’s \(\Gamma\)並沒有考慮平手的情況，也就是自變數或是依變數的不同順序類別有相等的值。因此，我們計算Kendall’s \(\tau-b\)以及Somer’s D，兩者結果相近：
etable<-cbind(c(35,76,22),c(73,134,56),
c(30,54,30))
DescTools::SomersDelta(etable, direction="column")
## [1] 0.041632.3 小結
with(mtcars, stats::cor.test(cyl, gear,method='spearman'))
##
## Spearman's rank correlation rho
##
## data: cyl and gear
## S = 8535, p-value = 8e-04
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.5643
with(mtcars, stats::cor.test(cyl, gear,method='kendall'))
##
## Kendall's rank correlation tau
##
## data: cyl and gear
## z = -3.2, p-value = 0.002
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.5125- 結果顯示，兩者相當接近，大約是-0.5，也就是說我們可以選擇其中之一衡量順序變數的相關程度。
3 卡方檢定
卡方檢定可分為適合度檢定(goodness of fit)以及獨立性檢定(test of independence)。前者檢定我們感興趣的變數是否來自特定的機率分佈的假設，後者檢定兩個類別變數是否獨立的假設。
3.1 卡方分佈
v1<-here::here('Fig','chi2.jpg')
knitr::include_graphics(v1)
Figure 3.1: 4種自由度的卡方分佈機率密度圖
3.2 自由度
3.3 適合度檢定
3.4 獨立性檢定
\[\hat{n_{ij}}=(RT\times CT)/ N\]
- \(RT\)：包含該列的邊際的總和
- \(CT\)：包含該列的邊際的總和
- \(N\)：觀察值總數
3.5 計算期待值
tbyt<-matrix(c("","","Y","",
"","1","2","Total",
"1","f11","f21","C1",
"2","f12","f22","C2",
"Total","R1","R2","N"),
ncol=5, nrow=4)
knitr::kable(tbyt, 'html', caption='2X2表格') %>%
kable_styling(c("striped", "bordered")) %>%
add_header_above(c(""," ", "X" = 3)) | 1 | 2 | Total | ||
| 1 | f11 | f12 | R1 | |
| Y | 2 | f21 | f22 | R2 |
| Total | C1 | C2 | N |
- 根據表3.1計算期待值：
- \(E11=\frac{R1\times C1}{N}\)
- \(E12=\frac{R1\times C2}{N}\)
- \(E21=\frac{R2\times C1}{N}\)
- \(E22=\frac{R2\times C2}{N}\)
tbythy<-matrix(c("","","Y","",
"","1","2","Total",
"1","25","25","50",
"2","25","25","50",
"Total","50","50","100"),
ncol=5, nrow=4)
knitr::kable(tbythy, 'html', caption='2X2表格假想1') %>%
kable_styling(c("striped", "bordered")) %>%
add_header_above(c(""," ", "X" = 3)) | 1 | 2 | Total | ||
| 1 | 25 | 25 | 50 | |
| Y | 2 | 25 | 25 | 50 |
| Total | 50 | 50 | 100 |
tbythy<-matrix(c("","","Y","",
"","1","2","Total",
"1","50","0","50",
"2","0","50","50",
"Total","50","50","100"),
ncol=5, nrow=4)
knitr::kable(tbythy, 'html', caption='2X2表格假想2') %>%
kable_styling(c("striped", "bordered")) %>%
add_header_above(c(""," ", "X" = 3)) | 1 | 2 | Total | ||
| 1 | 50 | 0 | 50 | |
| Y | 2 | 0 | 50 | 50 |
| Total | 50 | 50 | 100 |
3.6 計算卡方值公式
\[\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，越接近常態分佈。
3.7 適合度檢定
\[E_{k}=\frac{N}{n}\]
\(H_{0}\)為觀察值平均分佈在每一類別，如果拒斥這個假設，表示各類別之間的差距的確存在而非隨機分佈。
自由度為\(k-1\)。因為觀察值的總和為\(N\)，只要知道\(k-1\)個類別的觀察值就知道最後一個類別的觀察值。
另一方面，因為單一分佈的總和是固定的，所以在這一個限制下理論上的數值變化就喪失一個自由度，寫成\(p=1\)。\(k-p=k-1\)就是自由度。
3.7.1 範例
抽出100位1990年之後的出生的受訪者，觀察到男性為52人、女性為48人，請問這是平均分佈嗎？
計算為
\[\mathcal{X}^2=\frac{(52-50)^2}{50}+\frac{(48-50)^2}{50}=0.16\]
3.8 卡方值對應之臨界值
現代統計學的第703-704頁提供拒斥域的機率值(\(\alpha\))對應的卡方分配的臨界值表，左邊每一列是自由度，每一欄\(\mathcal{X^2}_{\alpha}\)則是拒斥域\(\alpha\)的機率值，\(P(\mathcal{X^2}>\mathcal{X_{a}^2})\alpha\)。表格的由左而右，代表\(\alpha\)越來越小，非拒斥域的機率值越來越大。在同樣的自由度下，拒斥域越小，卡方分配的臨界值也就越大。
右尾卡方值參考表 3.4:
rtchi<-matrix(c("","conf. level","",
"", "5%","1%",
"1","3.84","6.63",
"2","5.99","9.21",
"3","7.82","11.34",
"4","9.49","13.28","5","11.07","15.09"),
ncol=7, nrow=3)
knitr::kable(rtchi, 'html',caption='右尾卡方值表') %>%
kable_styling(c("striped", "bordered")) %>%
add_header_above(c("","","d.f." = 5))| 1 | 2 | 3 | 4 | 5 | ||
| conf. level | 5% | 3.84 | 5.99 | 7.82 | 9.49 | 11.07 |
| 1% | 6.63 | 9.21 | 11.34 | 13.28 | 15.09 |
- 也可以用\(\tt{qchisq()}\)這個函式，設定\(\alpha\)值與自由度之後，即可計算應該大於多少卡方值，要注意的是\(\tt{qchisq()}\)應該要填\(1-\alpha\)，而不是\(\alpha\)：
alpha=0.05
qchisq(1-alpha, df=3)
## [1] 7.815- 反過來，\(\tt{pchisq()}\) 可以輸出一定的自由度下卡方值的機率，例如：
1-pchisq(2.9, 1)
## [1] 0.08858- 回到我們的例子，當卡方值為0.16，自由度為1時：可計算機率值等於0.68：
1-pchisq(0.16, 1)
## [1] 0.6892查表可知機率值為0.68，因此無法拒斥資料來自母體的假設。
我們也可以用卡方檢定測試觀察值是否符合特定的分佈。例如我們經常做加權檢定，想要知道是否樣本符合母體，虛無假設是樣本分佈等於母體分佈，以性別而言，自由度為1。假設母體的分佈是：
男：\(\frac{9003150}{9003150+9184065}=0.495\)
女：\(\frac{9184065}{9003150+9184065}=0.505\)
假設我們觀察到1,067人，男女分佈如表 3.5，可計算對應的理論值：
tmp<-matrix(c("性別","N","母體值",
"男", "531", "528.165",
"女","539","538.835"),
ncol=3, nrow=3)
knitr::kable(tmp, 'html',caption='男女性別分布') %>%
kable_styling(c("striped", "bordered")) | 性別 | 男 | 女 |
| N | 531 | 539 |
| 母體值 | 528.165 | 538.835 |
計算卡方值： \[\begin{align*} \mathcal{X}^2 & =\frac{(531-528.165)^2}{528.165}+\frac{(539-538.835)^2}{538.835} \nonumber \\ & =0.015+0 \nonumber \\ & =0.015 \nonumber \\ & <3.84 \end{align*}\]
結論：我們無法拒斥\(H_{0}\)：觀察值與母體有相同分佈在每一類別這個假設。
-如果我們想知道政黨支持與婚姻狀態是否互相獨立？以表 3.6 的資料為例：
| 婚姻 | 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.75 85.25
## [2,] 99.25 63.75- 最後計算卡方值：
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.5133 0.6865 0.7992 1.0689
chisq<-sum(dt)
chisq
## [1] 3.068- 如果用
R的函數計算，可以用\(\tt{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.7, df = 1, p-value = 0.1- 因為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.4530 0.6059 0.7053 0.9433
chisq<-sum(dt)
chisq
## [1] 2.708- 可畫圖表示如圖 3.2：
z=curve(dchisq(x, df=1), from=0, to=6)
abline(v = chisq, lty = 3, lwd=1.5, col='red2')
Figure 3.2: 自由度為1的卡方分佈
3.9 標準化殘差
- 觀察值與期待值之間的差距稱為殘差。殘差越大表示觀察到的數目遠高於應該有的期待值，也表示往這個方向變動。
- 標準化殘差就是考慮每一格所在的行與列所佔的比例進行標準化。
\[z_{ij}=\frac{O_{ij}-E_{ij}}{\sqrt{E_{ij}(1-{\mathrm {row\hspace{.2em} prop.}})(1-{\mathrm {column\hspace{.2em} prop.}})}}\]
標準化殘差越大表示那一格的觀察值超過當虛無假設成立時的期待值越多，也就表示兩個變數的關聯越高。
有一筆觀察值為1067的資料，教育程度分佈如表 3.7，依據母體資料求出期待值以及殘差：
| 教育程度 | 理論值 | 殘差 |
| 小學及以下 | 192.66 | -6.66 |
| 國中 | 228.42 | -81.42 |
| 高中 | 216.56 | 94.43 |
| 專科 | 202.64 | -65.64 |
| 大學及以上 | 226.68 | 59.31 |
- 由表 3.7 可知高中、大學及以上教育程度這兩群人在樣本中相對於母體而言太多，其他類別則是人太少，因此可能需要加權以符合母體。
★ 卡方檢定注意事項
- 變數交叉列表的每一格的期待值最好大於5，也就是個案數分佈要平均，才不會用其他的格子推論極少個案數的類別之關係。
- 卡方值\(\mathcal{X}^2\)不表示關聯性的大小，只能用來檢驗兩個變數是否互相獨立的假設。卡方值越大不代表兩個變數關聯越大。
4 運用卡方值產生關聯係數
4.1 二分類變數
- \(\Phi\)是卡方值產生的關聯係數之一，介於0與1之間。如果兩個變數都是二分類變數，我們可用卡方值代入\(\Phi\)的公式如下：
\[\Phi=\sqrt{\frac{\mathcal{X}^2}{N}}\]
4.1.1 範例
- 在\(\textbf{UsingR::too.young}\) 這筆資料中，有男性與女性的約會年齡，我們歸類為二分類，以卡方檢定得到卡方值，並且同時用\(\textbf{DescTools::Phi()}\)計算\(\Phi\)，得到結果如下，兩者相當接近：
male.age<-car::recode(UsingR::too.young$Male,
"10:35=1; 36:100=2")
female.age<-car::recode(UsingR::too.young$Female,
"10:30=1; 31:100=2; NA=1")
chisq.test(male.age, female.age, correct = F) #chisquared=28.283
##
## Pearson's Chi-squared test
##
## data: male.age and female.age
## X-squared = 28, df = 1, p-value = 1e-07
DescTools::Phi(male.age, female.age)
## [1] 0.5946
phi=sqrt(28.283/80); phi
## [1] 0.59464.2 多分類變數
- Cramer’s \(V\)也是卡方值產生的關聯係數之一，介於0與1之間。如果兩個變數都是多分類變數，我們可用卡方值代入Cramer’s \(V\)的公式如下：
\[V=\sqrt{\frac{\mathcal{X}^2}{N*min_{r-1,c-1}}}\]
- Cramer’s V的公式中，分母是觀察值乘以列或者行的數目比較小的其中之一，再減去1。
- 同時用\(\textbf{DescTools::CramerV()}\)計算\(V\)，得到結果如下，與代入卡方值計算結果相當接近：
newdf<- UsingR::too.young %>% mutate(Male.age=case_when(Male >=10 & Male<=30 ~ 1, Male >=31 & Male<=50 ~ 2,
Male >= 51 ~ 3)) %>%
mutate(Female.age=case_when(Female>=10 & Female<=25 ~ 1,
Female>=26 & Female<=40 ~ 2,
Female>=41 & Female<=100 ~ 3,
Female=NA ~ 2))
n = nrow(newdf)
#expss package
library(expss)
newdf %>%
# cro_cpct(female.age, male.age) %>%
calc_cro_cpct(Female.age , Male.age)%>%
set_caption('性別與約會年齡交叉列表')| 性別與約會年齡交叉列表 | |||
| Male.age | |||
|---|---|---|---|
| 1 | 2 | 3 | |
| Female.age | |||
| 1 | 100 | 20 | |
| 2 | 75 | 50 | |
| 3 | 5 | 50 | |
| #Total cases | 49 | 20 | 10 |
#chisq.test
ct <- chisq.test(newdf$Male.age, newdf$Female.age, correct = F) #chisquared=80.148
DescTools::CramerV(newdf$Male.age, newdf$Female.age)
## [1] 0.7253
v=sqrt(ct$statistic/(n*2)); v
## X-squared
## 0.72085 相關分析
5.1 相關係數
- \(X\)與\(Y\)兩個變數之間的線性相關可以用方程式 (5.1) 來估計：
\[\begin{align*} r & =\frac{\sum (X-\bar{X})(Y-\bar{Y})}{\sqrt{\sum (X-\bar{X})}\sqrt{\sum (Y-\bar{Y})}} \\ & = \frac{S_{XY}}{S_{X}S_{Y}} \tag{5.1} \end{align*}\]
其中
\[S_{XY}=\frac{\sum (X-\bar{X})(Y-\bar{Y})}{n-1}\]
\[S_{X}=\sqrt{\frac{\sum (X-\bar{X})^2}{n-1}}\]
\[S_{Y}=\sqrt{\frac{\sum (Y-\bar{Y})^2}{n-1}}\]
- 也可以計算為：
\[r=\frac{n(\Sigma xy)-(\Sigma x)(\Sigma y)}{\sqrt{(n\Sigma x^2-(\Sigma x)^2)(n\Sigma y^2-(\Sigma y)^2)}}\]
\(r_{XY}\)=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.4844
ggplot(df, aes(x=X, y=Y)) +
geom_point(fill='brown', size=5, shape=1) +
geom_smooth(method="lm", se=F, col='blue')
- 由於我們不知道母體的相關程度，所以用樣本估計。我們可以用\(t\)檢定來檢定兩變數之間有無線性相關。\(t\)檢定統計量為：
\[\frac{T}{\sqrt{\frac{1-r^2}{n-2}}} \sim\hspace{.15em}t_{n-2}\]
5.1.1 範例
\(\textbf{trees}\)這筆資料中有Girth（周長）與Volume（體積）兩個變數，我們想知道兩者的相關程度。
檢查是否兩個變數之間有線性相關可用散佈圖，例如圖 5.1：
ggplot(trees, aes(x=Girth, y=Volume)) +
geom_point(size=3, col="darkgreen") +
labs(caption="Data: trees") +
theme_bw()
Figure 5.1: Girth與Volume散佈圖
- 可以加上迴歸線表現兩個變數的相關，例如圖 5.2：
ggplot(trees, aes(x=Girth, y=Volume)) +
geom_point(col='darkgreen', size=3) +
geom_smooth(method="lm", se=F, col='blue')
Figure 5.2: Girth與Volume散佈圖加上迴歸線
- 用\(\textbf{ggstatsplot::ggscatterstats}\)套件，同時顯示散佈圖、迴歸線、相關係數以及檢定結果。見圖 5.3。
ggstatsplot::ggscatterstats(
data = trees,
x = Girth,
y = Volume,
xlab = "",
ylab = "",
title = "",
messages = FALSE
)
Figure 5.3: 結合相關分析的Girth與Volume散佈圖
- 也可以用非線性迴歸線（Locally Weighted Scatterplot Smoothing, LOESS）表示兩者的相關如圖 5.4：
ggplot(trees, aes(x=Girth, y=Volume)) +
geom_point(col='darkgreen', size=3) +
geom_smooth(method="loess", se=F) +
theme_bw()
Figure 5.4: Girth與Volume散佈圖加上非線性線
- 在觀察完兩個變數的關係後，以公式計算相關程度：
n=nrow(trees)
x<-trees$Girth; y<-trees$Volume
x_sd<-sd(trees$Girth)
y_sd<-sd(trees$Volume)
xy_sd<-sum((x-mean(x))*(y-mean(y)))/(n-1)
r=xy_sd/(x_sd*y_sd)
cat("correlation=",r)
## correlation= 0.9671- 再以
R計算：
with(trees, cor(Girth, Volume))
## [1] 0.96715.2 相關係數與迴歸係數
如同圖 5.3 所顯示，線性相關與迴歸係數跟\(R^2\)之間有關係。
相關係數與迴歸係數的關係可寫成：
\[\begin{align} \hat{\beta}=r\frac{S_{Y}}{S_{X}} \tag{5.2} \end{align}\]
方程式 (5.2) 顯示，當依變數\(Y\)相對於自變數\(X\)的離散程度越大，而且相關程度越大時，迴歸係數也越大。
相關係數與判定係數\(R^2\)的關係為：
\[\begin{align} r^2=R^2=\frac{\sum(\hat{Y}-\bar{Y})^2}{\sum (Y-\bar{Y})^2} \tag{5.3} \end{align}\]
\(\hat{Y}\) 是迴歸模型對\(Y\)的預測值。
- 方程式 (5.3) 說明，兩個變數的相關係數等於自變數解釋依變數的變異量的比例。
5.3 相關與因果關係
兩個變數的相關係數高表示一個變數如果增加，另一個變數可能增加或者下降，但是不表示一個變數的增加或者減少「造成」另一個變數增加或減少。例如，太陽在春分以及秋分時直射地球赤道，太陽光越往北越斜，所以住在赤道的人會比住在北極的人覺得熱，這是因果關係。但是天氣熱、飲料銷售量提高，這是相關，因為消費行為可能跟文化、廣告、價格有關。
在社會科學研究中，我們需要理論來說明相關，不能只是想當然爾。例如，教育與成就有正相關，但是兩者之間為什麼有高度相關？是因為知識可以轉換成收入？還是教育過程帶來人脈，而人脈又帶來收入？這些問題需要許多研究。
6 作業
- 有一筆資料的交叉列表如表 6.1 顯示，贊成女性有墮胎權利的民眾似乎傾向投給民主黨。請問墮胎議題立場與投票兩者的相關程度為何？
vtable <- matrix(c("","","投票","",
"","民主黨","共和黨","總和",
"贊成", "46","41","87",
"反對", "39", "73", "112",
"總和","85","114","199"),
ncol=5, nrow=4)
kable(vtable, 'html', caption='墮胎議題與投票選擇') %>%
kable_styling(full_width=F, bootstrap_options = "striped", font_size=13) %>%
add_header_above(c("","","墮胎議題" = 3))| 贊成 | 反對 | 總和 | ||
| 民主黨 | 46 | 39 | 85 | |
| 投票 | 共和黨 | 41 | 73 | 114 |
| 總和 | 87 | 112 | 199 |
請用
tabyl 顯示MASS::Cars93 裡面的Type與Origin兩個變數的交叉列表，請顯示個案數以及正確的列或者行百分比。此外，用DescTools 裡面的Lambda 計算Origin解釋Type的比例（提示：如果手算的話，總數扣掉Origin的兩個眾數當作分母）。表6.2的資料來自現代統計學的習題 12.10，請分析抽煙與否跟壽命長短有沒有相關？
stable <- matrix(c("","","壽命","","","",
"","50以下","50-60","60-70", "70-80", "80以上",
"吸菸", "38","47","43","32","34",
"不吸菸", "30", "55", "51","37","33"),
ncol = 4, nrow = 6)
kableExtra::kable(stable, 'html', caption='抽煙與壽命') %>%
kableExtra::kable_styling(full_width=F, bootstrap_options = "striped", font_size=13) | 吸菸 | 不吸菸 | ||
| 50以下 | 38 | 30 | |
| 壽命 | 50-60 | 47 | 55 |
| 60-70 | 43 | 51 | |
| 70-80 | 32 | 37 | |
| 80以上 | 34 | 33 |
- 請分析表6.3的政黨認同與性別的關係。
| Democrat | Independent | Republican | |
|---|---|---|---|
| M | 762 | 327 | 468 |
| F | 484 | 239 | 477 |
請使用
PP0797B2.sav這筆資料，並先把統獨立場(tondu)的9設定為遺漏值，然後顯示年齡(age)與統獨立場(tondu)之間的交叉列表，再進行卡方檢定，請問這兩個變數互相獨立嗎？根據上一題的資料，請問年齡的分佈是否與母體一致？假設母體的表格如表 6.4：
Table 6.4: 年齡分佈
20至29歲
30至39歲
40至49歲
50至59歲
60至69歲
N
351
488
610
374
235
%
17.06
23.71
29.64
18.17
11.42
請問在表格 6.5 中，部門與是否提前退休的人數有無統計上相關？
Table 6.5: 部門與退休人數
財務
工程
業務
行政
退休
26
23
18
9
在職
6
7
14
23
請先畫散佈圖表現
UsingR這個套件中的nym.2002這筆資料，裡面的age以及time之間的關係，並且計算兩者的相關程度。請計算並且畫圖表示\(\textbf{trees}\)資料中三個變數的相關程度：
請針對\(\textbf{mtcars}\)這筆資料的wt以及mpg繪製散佈圖，並且計算相關係數：
7 更新講義時間
## 最後更新時間: 2024-04-10 19:15:13