社會科學統計方法

雙母體假設檢定、變異數分析

1 課程目標

當我們想要了解兩個母體平均數有無差異時，可以檢驗母體平均數是否相等的假設，包括獨立兩母體平均數的檢定、成對兩母體檢定、兩母體比例差異的估計與檢定等等。如果我們想要檢驗三個或者以上的母體平均數是否相等，則使用變異數分析。

例如我們想知道不同種族的平均薪資是否有差異，可以先進行描述統計如下：

library(expss); library(ISLR); data("Wage")

tbl_caption<-cro_mean_sd_n(Wage$wage, Wage$race) %>%
  set_caption('不同種族的平均薪資、標準差')
tbl_caption

不同種族的平均薪資、標準差
	Wage$race
	1. White	2. Black	3. Asian	4. Other
Wage$wage
Mean	112.6	101.6	120.3	90
Std. dev.	41.7	37.2	46.4	29.2
Unw. valid N	2480	293	190	37

然後進行變異數分析：

m1 <- aov(wage ~ race, data=ISLR::Wage )
summary(m1)
##               Df  Sum Sq Mean Sq F value  Pr(>F)    
## race           3   63212   21071    12.2 5.9e-08 ***
## Residuals   2996 5158874    1722                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

結果顯示不同種族的平均薪資有統計上的顯著差異。我們也可以畫圖1.1顯示每個類別的薪資分佈以及相關的統計數字：

# plot
ggstatsplot::ggbetweenstats(
  data = ISLR::Wage,
  x = race,
  y = wage,
  xlab = "Race",
  ylab = "Wage",
  plot.type = 'box',
  type = "parametric",
  bf.message = F,
  p.adjust.method = 'none')

Figure 1.1: 不同種族的薪資

• 結果顯示白人與非裔之間薪資有顯著差異，而非裔與亞裔之間也有差異，亞裔與其他種族之間也有差異。

1.1 兩個獨立母體之樣本檢定

我們經常需要知道兩套樣本的差異是否達到統計上的顯著水準，例如兩種教學方法的成績、兩個餐廳的評價、兩種手術方法的患者恢復速度、兩個國家人民對於民主的滿意程度等等。

從兩個獨立母體抽樣，假設\(X_{i}\quad i=1,\ldots,n_{x}\)，\(Y_{i}\quad i=1,\ldots,n_{y}\)是獨立隨機樣本，虛無假設為：

\[H_{0}:\mu_{x}=\mu_{y}\]

• 對立假設則為：

\[H_{1}:\mu_{x}<\mu_{y},\mu_{x}>\mu_{y},\mu_{x}\neq \mu_{y}\]

• 樣本的平均數差異(\(\bar{X_{1}}-\bar{X_{2}}\))的抽樣分配趨近於常態分配。

\[T=\frac{(\bar{X}-\bar{Y})-E(\bar{X}-\bar{Y}|H_{0})}{SE(\bar{X}-\bar{Y}|H_{0})}\]

如果母體的變異數未知，只有樣本的變異數，我們可以用樣本變異數做為母體變異數的估計。

• 那麼如何計算兩套樣本的標準誤？假設兩個母體的變異數相等，\(\sigma=\sigma_{x}=\sigma_{y}\)，兩個母體的標準差公式為：

\[s_{pooled}^2=\frac{(n_{x}-1)s^{2}_{x}+(n_{y}-1)s^{2}_{y}}{n_{x}+n_{y}-2}\]

• 其中\(s^{2}_{x}=\frac{\sum(X-\bar{X})^2}{n-1}\)
• 如果\(n_{1}=n_{2}=n\)，

\[s_{pooled}^2=\frac{s^{2}_{x}+s^{2}_{y}}{2}\]

• 兩個母體的標準誤公式為：

\[SE(\mu_{x}-\mu_{y})=s_{p}\sqrt{\frac{1}{n_{x}}+\frac{1}{n_{y}}}\]

• 假設兩個母體的變異數不相等，兩個母體的標準誤公式為：

\[SE(\mu_{x}-\mu_{y})=\sqrt{\frac{s_{x}^{2}}{n_{x}}+\frac{s_{y}^{2}}{n_{y}}}\]

• 如果\(n_{1}=n_{2}=n\)，標準誤計算為：

\[SE(\mu_{x}-\mu_{y})=\sqrt{\frac{s_{x}^{2}+s_{y}^{2}}{n}}\]

• 自由度：\(n_{1}+n_{2}-2\)

1.2 假設檢定步驟

1. 建立虛無假設
2. 假設雙母體變異數是否相等
3. 計算樣本標準差
4. 計算標準誤
5. 計算\(T\)值
6. 計算\(T\)的機率，自由度：\(n_{1}+n_{2}-2\)

1.2.1 雙母體平均數相等假設範例1

• 假設有兩種油漆，我們想知道這兩種油漆的效果是不是相似。兩種油漆的資料如下，

sample1<-c(19.7475, 20.5562, 16.0734, 18.1730, 19.8387, 
           14.5291, 19.4998, 18.8374, 16.6873, 14.7627,
           18.3869, 17.9287, 21.6973, 19.3439, 15.7374, 
           18.3563, 22.0878, 22.5745, 18.0030, 18.6004)

sample2<-c(19.4715, 18.8613, 17.4827, 21.6029, 23.0386, 
           20.6289, 19.9357, 19.8812, 16.6012, 21.7741, 
           21.9265, 18.4178, 20.4401, 20.1119, 22.9955,
           21.1385, 19.4969, 20.4448, 19.6675, 17.0984)

• 建立虛無假設： \[H_{0}: \mu_{x}=\mu_{y}\] li
• 先用箱形圖顯示兩套樣本的分佈，可以看出第一套樣本比較集中，第二套比較分散，如圖1.2：

dt=data.frame(sample1, sample2)
dt.m <- melt(dt) %>%
ggplot(aes(x=variable, y=value)) +
  geom_boxplot()  +
   labs(x="Groups", y='') +
  stat_summary(fun.y=mean, geom="point", shape=16, size=2, col='red') 
dt.m

Figure 1.2: 兩套樣本的盒型圖

• 也可以畫機率密度圖如圖1.3：

with(dt, plot(density(sample1), main='', col='#22ee00', ylim=c(0,0.25)))
with(dt, lines(density(sample2), lty=2, col='#aa11cc'))

Figure 1.3: 雙樣本機率密度圖

• 假設兩個母體的變異數相同。

• 根據以上公式，我們首先計算標準差以及標準誤：

n1<-length(sample1); n2<-length(sample2)

var1<-var(sample1); var2<-var(sample2)

sp=(var1+var2)/2
se <-sqrt(sp)*(sqrt(n1^-1+n2^-1)); se
## [1] 0.6422

• 然後計算\(T\)值：

T=(mean(sample1)-mean(sample2))/se
cat('T =',T, '\n')
## T = -2.304

• 接下來找出\(T\)的發生機率，自由度是\(n_{1}+n_{2}-2\)。因為\(T<0\)，所以拒斥域是落在左邊，如圖1.4，計算\(p\)值時不需要被1減。

Figure 1.4: 左尾累積機率圖

pvalue=pt(T, n1+n2-2)*2
cat('p value =', pvalue)
## p value = 0.02677

• 因為 \(p<0.05\)，因此我們不接受兩個母體的平均數相等的假設。
• 以\(\tt{t.test()}\)檢定假設：兩筆資料的平均數相等，或者是相減等於0：

 t.test(sample1, sample2, var.equal = T)
## 
##  Two Sample t-test
## 
## data:  sample1 and sample2
## t = -2.3, df = 38, p-value = 0.03
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.7798 -0.1797
## sample estimates:
## mean of x mean of y 
##     18.57     20.05

• 由以上的結果可以看出，\(T\)=-2.304，\(p<0.05\)，因此我們可以拒斥兩筆資料的平均數相等的假設。95%的區間是(-2.779, -0.179)，不包含0。結論是兩種油漆的效果不相同。

• 我們改假設變異數不等，得到類似結果。

 t.test(sample1, sample2, var.equal = F)
## 
##  Welch Two Sample t-test
## 
## data:  sample1 and sample2
## t = -2.3, df = 36, p-value = 0.03
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.7820 -0.1775
## sample estimates:
## mean of x mean of y 
##     18.57     20.05

1.2.2 單尾檢定

• 除了檢驗兩個樣本有相同的平均值，也可以檢驗一個大於另一個平均值。以上面資料為例，我們可以檢驗第一組的油漆效果是否比第二組來得好：

\[H_{0}:\mu_{x}>\mu_{y}\]

• \(t\)檢定的結果如下：

 t.test(sample1, sample2, var.equal = T, alternative = 'less')
## 
##  Two Sample t-test
## 
## data:  sample1 and sample2
## t = -2.3, df = 38, p-value = 0.01
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##    -Inf -0.397
## sample estimates:
## mean of x mean of y 
##     18.57     20.05

• 可以看出\(p< 0.05\)。我們拒斥第一筆資料平均數大於第二筆資料的平均數的假設，結論是第一組油漆的效果比較差。

• 雙樣本的檢定需要比較特別的標準誤以及自由度，檢定過程與單樣本的\(t\)檢定相同。

1.2.3 雙母體平均數相等假設範例2

在UsingR::babyboom資料中，有嬰兒的性別以及出生重量的變數，試問男嬰比女嬰重嗎？虛無假設是兩者相等，對立假設是男嬰比女嬰重。

\[H_{0}: \mu_{boy}=\mu_{girl}\]

\[H_{1}: \mu_{boy}>\mu_{girl}\]

• 我們假設男嬰、女嬰來自兩個不同母體，比較樣本的平均數以及標準誤，求出\(T\)值，進行檢定：

library(UsingR);data(babyboom)
male <-  dplyr::filter(babyboom, gender=='boy')

female <- dplyr::filter(babyboom, gender=='girl')

t.test(male$wt, female$wt, alternative='greater')
## 
##  Welch Two Sample t-test
## 
## data:  male$wt and female$wt
## t = 1.4, df = 28, p-value = 0.08
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -48 Inf
## sample estimates:
## mean of x mean of y 
##      3375      3132

• 因為p>0.05，所以結論是我們無法拒斥男嬰女嬰一樣重的假設。

1.3 雙樣本的比例檢定

• 如果兩個樣本都是來自二項分布，那麼我們可以檢驗比例是否相同，也就是檢驗以下的虛無假設：

\[H_{0}:p_{x}-p_{y}=0\]

• 如果把0換成任意數，等於是檢驗母體比例的差異是否達到某個程度。寫成：

\[H_{0}:p_{x}-p_{y}=D_{0}\]

• 以上的雙尾檢定可以改為單尾檢定。

標準誤公式為：

\[\sqrt{\frac{p_{1}(1-p_{1})}{n_{1}}+\frac{p_{2}(1-p_{2})}{n_{2}}}\]

• 計算T值或者Z值，然後加以檢定：

\[T=\frac{p_{1}-p_{2}}{SE}\]

1.3.1 雙樣本比例檢定範例

• 有一輛捷運，第一節車廂有50名乘客，第二節車廂有61人，第一節車廂有20位男性，第二節則有13位男性，請問兩節車廂有相同的男性比例嗎？用\(\tt{prop.test()}\)計算如下：

x1=20; n1=50; x2=13; n2=61
prop.test(c(x1,x2), c(n1, n2), correct = FALSE)
## 
##  2-sample test for equality of proportions without continuity
##  correction
## 
## data:  c(x1, x2) out of c(n1, n2)
## X-squared = 4.6, df = 1, p-value = 0.03
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.01659 0.35718
## sample estimates:
## prop 1 prop 2 
## 0.4000 0.2131

• 結果顯示，兩節車廂各有0.4以及0.21比例的男性。因為\(p<0.05\)，因此以\(\alpha=0.05\)的標準，我們拒斥兩節車廂有一樣男性比例的假設。

1.4 成對樣本檢定

前面的雙樣本檢定假定兩套樣本互相獨立，但是有不同的平均值與標準差。如果是想要確定兩個變數來自於同一套樣本，兩次測量的平均值相等，則是進行成對樣本檢定。

• 例如病患服用藥物前後的病況、選手接受訓練前後的表現、受訪者第一次受訪跟第二次受訪的回答等等，適合進行成對樣本檢定。

• 用R的\(\tt{t.test()}\)函數，可設定\(\tt{paired=T}\)，進行成對樣本檢定，虛無假設為兩次測量的平均值相等。

• 例如比較民進黨候選人在台北市第2選區在2016年的立委選舉與2019年的立委補選的各里的平均得票率是否相同：

filetaipei<-here::here('data','taipei201619.csv')
taipei<-read.csv(filetaipei, header=T, fileEncoding = 'BIG5')
library(dplyr)
taipei<-taipei %>% 
  mutate(DPP2016p=100*DPP2016/total2016, DPP2019p=100*DPP2019/total2019)
with(taipei, t.test(DPP2016p, DPP2019p, paired = T))
## 
##  Paired t-test
## 
## data:  DPP2016p and DPP2019p
## t = 26, df = 196, p-value <2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  11.07 12.90
## sample estimates:
## mean of the differences 
##                   11.98

• 結果顯示兩者有很大的差異，拒斥兩次選舉的平均值相減等於0的假設。

• 為了確認這個檢定結果，我們可以畫圖顯示兩次選舉的平均值。先整理成為長資料：

filetaipei<-here::here('data','taipei201619.csv')
taipei<-read.csv(filetaipei, header=T, fileEncoding = 'BIG5')

taipei<-taipei %>% 
  mutate(DPP2016p=100*DPP2016/total2016,
         DPP2019p=100*DPP2019/total2019)
dat<- data.frame(taipei$DPP2016p, taipei$DPP2019p)
colnames(dat)<-c('2016','2019')
dat.new <- reshape2::melt(dat, variable.name = "election", value.name = "DPP") 
dat.new <- dat.new %>% 
           group_by(election) %>%
           mutate(grp.mean=mean(DPP)) %>%
           mutate(Year=as.factor(election))

• 然後使用ggplot2畫分佈圖：

ggplot(dat.new, aes(x=DPP, color=Year)) +
  geom_density(position="identity", alpha=.4) +
  geom_vline(data=dat.new, aes(xintercept=grp.mean, color=Year),
             linetype="dashed") +
  labs(x="DPP's Vote Share", y="")

Figure 1.5: 民進黨2016大選及2019補選得票率

• 圖 1.5 顯示發現民進黨候選人在台北市第二選區的兩次選舉的得票率，的確有相當不同的分佈以及集中趨勢。

2 變異數分析

2.1 變異數分析基礎

變異數分析的功能是檢驗一個類別變數與其它連續變數(continuous variable)之間的關係，可運用One-way Anova(單因子變異數分析)，分析類別之間是否存在平均值的差異。

Two-way Anova(雙因子變異數分析)則是有兩個以上的類別變數，分析類別之間是否存在連續變數的平均值的差異。

Anova分析建立在三個假設：

1. 類別變數與連續變數互相獨立
2. 每一個類別內的連續變數呈常態分佈
3. 每一個類別內的連續變數的變異數相同

變異數分析的目標是比較每一組之間的離散程度跟每一類別（或是每一組）內部的離散程度。兩者相比越小，表示組跟組之間的差異越小於每一組內部的差異，因此分組的意義也越小，類別變數與連續變數的關係也越小。反之則越大。

• 變異數分析的虛無假設是：

\[H_{0}: \mu_{a}=\mu_{b}=\mu_{c}\]

• 如果每一個分組的平均數相等，那麼變異數來自於因為分組造成的差異，所以透過變異數分析，可以檢驗是否每一個分組有相同的平均數。
• 首先我們計算資料的全部的變異量，就是所有觀察值與全部的平均數的差異的平方，也就是：

\[SST=\sum_{i=1}\sum_{k=1}(x_{ij}-\bar{\bar{x}})^2\]

• 其中

\[\bar{\bar{x}}=\frac{\sum_{i=1}\sum_{k=1}x_{ik}}{n}\]

• 其次是 組間變異(SSB, sum of square between groups)。SSB為各組的樣本數乘以組平均減總平均的平方的加總。

\[SSB=\sum n_{k}(\bar{x_{k}}-\bar{x})^2\]

• 組內變異也稱為刺激的變異(SSTr)，也有人寫成SSF(sum of squares due to factor)。換句話說，是「組」這個變數帶來的差異。

• 然後是組內變異(SSW, sum of square within groups)，為各組的變異數的加總。

\[\begin{align*} SSW & =\sum_{k}\sum (x_{ki}-\bar{x_{k}})^2 \\ & =\sum_{k}(n_{k}-1)s_{k}^2 \end{align*}\]

• 有人把SSW寫成SSE(error sum of squares)，也就是組間差異解釋全部變異不足的部分就是SSE。

• 全部的變異量等於這兩者相加：

\[SST=SSB+SSW\]

• 每一組之間的平均離散程度用MSB(mean square between)表示，每一類別（或是每一組）內部的離散程度用MSW(mean square within)表示。計算方式為：

\[MSB=MSF=SSB/df(B) \quad df(B)=k-1\]

\[MSW=MSE=SSW/df(W) \quad df(W)=n-1-(k-1)=n-k\]

• 有了MSB, MSW，可以計算\(F\)值做為檢定用。

\[F=\frac{MSB}{MSW}\]

• 使用\(F\)分佈而不是\(t\)分佈檢驗假設成立與否。英國統計學家兼生物學家羅納德·費雪（Ronald Aylmer Fisher）發明了\(F\)檢驗。

• F 分佈的參數是兩個自由度，第一個自由度是\(k-1\)，也就是分組數減去1，第二個自由度是\(n-k\)，也就是個案數減去分組數。可畫圖 ?? 如下：

curve(df(x, df1=1, df2=2), from=0, to=5, ylab='')

Figure 2.1: 三種F分佈

curve(df(x, df1=5, df2=6), from=0, to=5, col='red', add=T)

Figure 2.2: 三種F分佈

curve(df(x, df1=20, df2=30), from=0, to=5, col='blue2', add=T)
abline(v=1, lwd=1.5, lty=2, col='green2')
legend('topright', c('df=1;df=2','df=5;df=6',
                     'df=20;df=30'),
            col=c('black',"red", "blue2"),
            lty=c(1,1,1))

Figure 2.3: 三種F分佈

• 圖 ??顯示自由度越大、\(F\)分佈圖形越接近常態分佈。

• R的\(F\)分佈有以下指令：

• df(x, df1, df2)回傳x百分比的機率 (y軸的高度)

• pf(q, df1, df2)回傳q百分比的累積機率

• qf(p, df1, df2)回傳p值的\(F\)值

• rf(n, df1, df2)回傳從\(F\)分佈抽樣的n個樣本

假設有10組資料，每一組資料各有6個觀察值，\(n=60\)，\(k=10\)，\(F\)值為1.8，請問\(p=\)？

• 計算如下，用\(\tt{pf()}\)求出特定\(F\)值的\(p\)值，也可以用\(\tt{qf()}\)求出臨界值：

alpha=0.05
pf(1.8, 9, 50, lower.tail = T)
## [1] 0.9084
qf(1-alpha/2, 9, 50)
## [1] 2.381

curve( df(x, df1=9, df2=50),
       xlim = c(0, 3),
       ylab = "Density",
      # main = "F機率密度",
       col='#e2a1cc', lwd=2)

n=60; k=10; alpha=0.05
upperbound <- qf(1-alpha/2, k-1, n-k)

x<-seq(upperbound, 4, by=0.1)
y=df(x, k-1, n-k)
cord.x1<-c(upperbound, x, 4)
cord.y1<-c(0, y, 0)

polygon(cord.x1, cord.y1,col='#20aaee')
text(2.48, 0.14, "2.5%")

Figure 2.4: F檢定圖

2.1.1 聯合信賴區間檢定方法

如果變異數分析發現全部的平均數不相等，可運用Bonferroni信賴區間。計算方式為：

\[\bar{X_{i}}-\bar{X_{j}}\pm t_{n-k,\alpha/2m}\sqrt{MSE}\sqrt{\frac{1}{n_{i}}+\frac{1}{n_{j}}}\]

其中m=\(k\choose{2}\)，例如有三組平均數比較，我們就比較三次，m=3，因為\(3\choose{2}\)=\(\frac{3!}{2!(3-2)!}=3\)。

2.2 變異數分析實例1

以Orange這筆資料為例，我們想檢驗每一種類型的樹的樹圍是否相等。先整理資料，計算平均樹圍以及標準差：

Orange$tree<-factor(Orange$Tree, levels=c("1","2", "3", "4","5"))

mydata <- Orange %>%
          group_by (tree) %>%
          summarise(
         count = n(),
         avg = mean(circumference, na.rm = TRUE),
         sd = sd(circumference, na.rm = TRUE)
  ) 
kable(mydata) %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

tree	count	avg	sd
1	7	99.57	43.29
2	7	135.29	66.32
3	7	94.00	42.98
4	7	139.29	71.90
5	7	111.14	58.86

• 可以看出五種平均樹圍之間有一些差異。

• 用箱形圖 2.5 檢視五種類型的樹圍：

p=ggplot(Orange, aes(x=tree, y=circumference)) 
p+geom_boxplot(data=Orange, aes(color=tree)) +
  theme_bw()

Figure 2.5: 五種樹的樹圍

• 箱形圖顯示平均數可能不相等。再用密度圖觀察是否常態分佈？

G=ggplot(Orange, aes(x=circumference, fill=tree, alpha=0.2))
G+geom_density() +
   xlim(c(0, 300)) +
   theme_bw()

Figure 2.6: 箱型圖與機率密度

• 圖 2.6 顯示部分類別可能不是成常態分佈，而且離散程度可能不相等。

  • 最後用aov函數進行Anova分析：

res.aov <- aov(circumference ~ tree, data = Orange)
# Summary of the analysis
summary(res.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)
## tree         4  11841    2960    0.88   0.49
## Residuals   30 100525    3351

• 結果顯示\(p>0.05\)，無法拒斥五種平均樹圍相等的假設。

• 我們可以畫圖檢視變異數相等的假設：

plot(res.aov, 1, cex=1.5, lwd=1.5)

Figure 2.7: 變異數圖形

• 圖 2.7 顯示，每一類別的平均值跟觀察值與該平均值之間的殘差(residual)並沒有明顯的關聯。因此，變異數相等的假設成立。

2.2.1 試算組間差、組內差以及F

我們也可以計算SSB, SSW, MSB, MSW，然後計算F=\(\frac{MSB}{MSW}\)

• 首先我們計算這5組的個案數。這邊要用到\(\tt{apply()}\)這個函式，它可以計算每一欄或者每一列的個案數、平均數、標準差等等，然後存成向量。

  • 首先整理資料：

#organize data
tmp <- Orange 
tmpdf<-data.frame(A=tmp$circumference[tmp$Tree==1],
                  B=tmp$circumference[tmp$Tree==2],
                  C=tmp$circumference[tmp$Tree==3],
                  D=tmp$circumference[tmp$Tree==4],
                  E=tmp$circumference[tmp$Tree==5])
tmpdf
##     A   B   C   D   E
## 1  30  33  30  32  30
## 2  58  69  51  62  49
## 3  87 111  75 112  81
## 4 115 156 108 167 125
## 5 120 172 115 179 142
## 6 142 203 139 209 174
## 7 145 203 140 214 177

• 然後計算每一組個案數：

g.n<-apply(tmpdf, 2, function(x) length(x))
g.n
## A B C D E 
## 7 7 7 7 7
n=sum(g.n); n
## [1] 35
k=length(g.n); k
## [1] 5

• 計算每一組的平均數：

#分組平均
g.mean<-apply(tmpdf, 2, function(x) mean(na.omit(x)))
g.mean
##      A      B      C      D      E 
##  99.57 135.29  94.00 139.29 111.14

• 然後我們可以計算組間變異(SSB)：

mean<-sum(tmpdf)/n
print(mean)
## [1] 115.9
SSB=sum((g.mean - mean)^2*g.n)
as.numeric(SSB)
## [1] 11841

• 再計算組內變異(SSW)，：

g.s<-apply(tmpdf, 2, function(x) sd(x))
g.s2<-g.s^2
SSW=sum(g.s2*(g.n-1))
as.numeric(SSW)
## [1] 100525

• 最後計算MSB以及MSW，就可以得到\(F\)值。

MSB=SSB/(length(g.n)-1)
as.numeric(MSB)
## [1] 2960
MSW=SSW/(sum(g.n)-length(g.n))
as.numeric(MSW)
## [1] 3351
Fstats=MSB/MSW; as.numeric(Fstats)
## [1] 0.8834

• 有了\(F\)值，計算兩個自由度：k-1, n-k，找出對應的\(p\)值：

pf(Fstats, k-1, n-k)
## [1] 0.5143

• 因為\(p > 0.05\)，所以接受\(\mu_{a}=\mu_{c}=\ldots=\mu_{k}\)的虛無假設。

2.3 變異數分析實例2

用chickwts這筆資料，請問不同飼料(feed)所得到的雞的重量(weight)的平均數是否不同？

• 我們用\(\tt{oneway.test()}\)這個函式，很快地計算\(F\)值以及檢定結果：

library(UsingR)
oneway.test(weight ~ feed, data=chickwts, var.equal = F)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  weight and feed
## F = 20, num df = 5, denom df = 30, p-value = 1e-08

• 結果顯示，\(p<0.05\)，所以結論為不同飼料餵的雞的重量有不同的平均數。

• 我們也可以用\(\tt{aov()}\)函式再做一遍，得到的結果一致。

• 不過，大部份時候\(\tt{oneway.test()}\)與\(\tt{aov()}\)這兩個函式得到的結果不同，所以在操作時要小心。建議直接使用\(\tt{aov()}\)進行變異數分析。

library(UsingR)
m1<-aov(Speed ~ Expt, data=morley)
summary(m1)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Expt         1  72581   72581      13 0.00048 ***
## Residuals   98 545444    5566                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.4 變異數分析實例3

• mtcars這筆資料有汽缸數(cyl)這個變數，請問4, 6,8這三種汽缸的汽車有重量(wt)的差別嗎？首先我們用\(\tt{aov()}\)進行變異數分析。

library(UsingR)
m1<-aov(mpg ~ cyl, data=mtcars)
summary(m1)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## cyl          1    818     818    79.6 6.1e-10 ***
## Residuals   30    308      10                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 根據以上結果，我們已經知道MSE=10.3，\(m=3\)。但是我們要查出自由度為\(n-k\)，\(t\)分配右尾機率為\(\alpha/2m\)的臨界值。

alpha=0.05;mse=10.3;m=3;k=3;n=32
tstar=qt(1-alpha/(2*m), n-k)
tstar
## [1] 2.541

mean.cyl<-mtcars %>% group_by(cyl) %>%
           summarise(obs=n(), mmpg=mean(mpg))
mean.cyl4<-mean.cyl$mmpg[1]
mean.cyl6<-mean.cyl$mmpg[2]
mean.cyl8<-mean.cyl$mmpg[3]
obs.cyl4<-mean.cyl$obs[1]
obs.cyl6<-mean.cyl$obs[2]
obs.cyl8<-mean.cyl$obs[3]
#4, 6
cat("mu1_mu2=(",(mean.cyl4-mean.cyl6)-tstar*sqrt(mse)*
      sqrt((1/obs.cyl4)+(1/obs.cyl6)),",",
(mean.cyl4-mean.cyl6)+tstar*sqrt(mse)*sqrt((1/obs.cyl4)+(1/obs.cyl6)),")","\n")
## mu1_mu2=( 2.978 , 10.86 )
#4,8
cat("mu1_mu3=(",(mean.cyl4-mean.cyl8)-tstar*sqrt(mse)*sqrt(1/obs.cyl4+1/obs.cyl8),",",
(mean.cyl4-mean.cyl8)+tstar*sqrt(mse)*sqrt(1/obs.cyl4+1/obs.cyl8),")","\n")
## mu1_mu3=( 8.278 , 14.85 )
#6,8
cat("mu2_mu3=(",(mean.cyl6-mean.cyl8)-tstar*sqrt(mse)*sqrt(1/obs.cyl6+1/obs.cyl8),",",
(mean.cyl6-mean.cyl8)+tstar*sqrt(mse)*sqrt(1/obs.cyl6+1/obs.cyl8),")")
## mu2_mu3=( 0.868 , 8.418 )

• 我們用\(\tt{pairwise.t.test()}\)測試聯合區間的假設，得到三個分組之間都有顯著差異的結果：

attach(mtcars)
pairwise.t.test(mpg, cyl, p.adj='bonf')
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  mpg and cyl 
## 
##   4     6   
## 6 4e-04 -   
## 8 3e-09 0.01
## 
## P value adjustment method: bonferroni

• R另外提供Bartlett以及Levene兩種檢定變異數相等的函數。Bartlett檢驗的做法是：

bartlett.test(mpg ~ cyl, data = mtcars)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  mpg by cyl
## Bartlett's K-squared = 8.4, df = 2, p-value = 0.02

• Levene檢驗則是：

library(car)
leveneTest(mpg ~ as.factor(cyl), data = mtcars)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)   
## group  2    5.51 0.0094 **
##       29                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 兩個檢驗結果都是拒斥變異數相等的假設。

• 雖然我們無法拒斥每一組平均數相等的假設，但是我們仍然可以用TukeyHSD函數檢驗哪些組別之間可能有差異：

A1<-aov(mpg ~ as.factor(cyl), data=mtcars)
TukeyHSD(A1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mpg ~ as.factor(cyl), data = mtcars)
## 
## $`as.factor(cyl)`
##        diff     lwr     upr  p adj
## 6-4  -6.921 -10.769 -3.0722 0.0003
## 8-4 -11.564 -14.771 -8.3565 0.0000
## 8-6  -4.643  -8.328 -0.9581 0.0112

• 結果發現每一組之間都有顯著的差異。

2.5 二因子變異數分析

如果有兩個因子可能影響某一個變數，可以進行二因子變異數分析。例如學生成績可能受到母親職業以及母親學歷的影響，不同職業與不同學歷，造成學生成績的差異。

我們分別針對個別因子進行變異數分析，例如針對職業做變異數分析，然後針對不同學歷分析。首先我們讀取資料以及觀察變數欄位。

kidiq <- here::here('data','kidiq.dta')
kid <- sjlabelled::read_stata(kidiq)

head(kid)
##   kid_score mom_hs mom_iq mom_work mom_age
## 1        65      1 121.12        4      27
## 2        98      1  89.36        4      25
## 3        85      1 115.44        4      27
## 4        83      1  99.45        3      25
## 5       115      1  92.75        4      27
## 6        98      0 107.90        1      18

• 接下來進行變異數分析：

A2 <- aov(kid_score ~ mom_hs + mom_work, data=kid)
summary(A2)
##              Df Sum Sq Mean Sq F value  Pr(>F)    
## mom_hs        1  10125   10125   25.65 6.1e-07 ***
## mom_work      1    144     144    0.37    0.55    
## Residuals   431 170117     395                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 由以上可知，母親的學歷造成學生成績的差異，但是職業並沒有學生成績的差異。

兩個因子之間可能有交互作用，我們在未來再討論。

3 作業

3.1 雙母體統計估計

1. 請回答現代統計學二版的第11.4題。

2. 請回答現代統計學二版的第11.6題。

3. 請回答現代統計學二版的第11.11題，但是三個職業的樣本數各為100。

4. 某一調查中，18位受訪者對2位候選人的表現打成績如下表，老師想知道這兩位候選人得到的成績是否相同，請問如何檢定？

Respondent	Candidate.1	Candidate.2
1	75	66
2	80	74
3	86	74
4	70	74
5	80	89
6	69	79
7	59	83
8	86	85
9	74	80
10	80	80
11	88	60
12	88	59
13	76	49
14	90	42
15	91	30
16	70	72
17	68	60
18	55	88

5. 請用Moodle上面的「立委-A05-2（台北市）.ods」資料以及「候選人在各投開票所得票數一覽表」，整理台北市第二選區國民黨候選人在2016立委選舉以及2019補選的資料，然後先畫圖比較國民黨候選人在台北市第2選區在2016年的立委選舉與2019年的立委補選的各里的分佈，最後檢驗平均得票率是否相同。

3.2 變異數分析

1. 假設有三個月的里程數如下表3.1，請問月份與里程數有關嗎？

may=c(2166,1568,2233,1882,2019)
sep=c(2279,2075,2131,2009,1793)
dec=c(2226,2154,2583,2010,2190)
dt<-cbind(may, sep, dec)
knitr::kable(dt, 'simple', caption='里程數資料')

Table 3.1: 里程數資料

may	sep	dec
2166	2279	2226
1568	2075	2154
2233	2131	2583
1882	2009	2010
2019	1793	2190

2. 請建立以下的模擬資料，然後檢驗各組的平均數是否相等。如果不相等，請找出哪兩組有不同的平均數。

type1 <- c(303, 393, 396, 399, 398)
type2 <- c(322, 326, 315, 318, 320, 320)
type3 <- c(125, 351, 226, 119, 280)
wear <- list(type1=type1,type2=type2,type3=type3)
wear.st <- stack(wear)

3. 請使用UsingR::MLBattend這筆資料，檢驗在美國聯盟(AL)裡面，不同分區(division)之間的平均觀眾人數相等的虛無假設，並且視覺化這三個分區的觀眾人數分布。

4. 假設統一集團為觀察7-11便利超商在大都會，城市及鄉下之銷售額有無差異。則自3個區域隨機抽出幾家超商，並記錄如下：

seven <- rbind(metro=c(10,13,14,16,17,NA),
               city=c(5,8,10,11,12,14),
               country=c(6,8,9,13,NA,NA))
kable(seven)%>%
  kable_styling(bootstrap_options = "striped", full_width = F)

metro	10	13	14	16	17	NA
city	5	8	10	11	12	14
country	6	8	9	13	NA	NA

• 1. 請問SSB, SSW, MSF, MSE, F分別為何？
  • 2. 請問區域是否影響銷售額？

5. 請回答現代統計學二版的第13.8題。

4 更新內容日期

## 最後更新日期 05/14/2021