1 課程目標

本週上課將介紹迴歸的假設檢定。因為資料分析的樣本來自於母體，所以樣本的迴歸係數不同於母體的迴歸係數。例如我們想知道母體的平均數\(\mu\)，我們假設\(\hat{\mu}=\bar{y}=\frac{1}{n}\Sigma_{i=1}^{n}y_{i}\)。同理，我們想要系統化\(\hat{\beta}_{0}\)與\(\hat{\beta}_{1}\)對於\(\beta_{0}\)跟\(\beta_{1}\)的估計。

2 假設檢定

  • 假設檢定是要檢驗有關母體參數的假設。
  • 因為沒有辦法完全確定檢定結果無誤，所以必須假設願意接受多少判斷錯誤的機率。
  • \(\alpha\)是錯誤地拒絕假設的機率，\(\beta\)是錯誤地不拒絕虛無假設的機率。

    • 步驟

    1. 設定虛無假設 \(\tt{H}_{0}\)
    2. 決定\(\alpha=\)Pr(否定\(\tt{H}_{0}\) | \(\tt{H}_{0}\))的大小。
    3. 根據樣本統計以及\(\tt{H}_{0}\)選擇\(t\)值。
    4. 假設\(\tt{H}_{0}\)為真，建立\(T\)的虛無分佈（常態分佈）。
    5. 從分佈表查拒斥值(critical values)，估計\(T\)觀察值有多難出現（被觀察到）。
    6. 如果從上述分佈觀察到\(T\)的機率，跟\(t\)值比拒斥值\(\alpha\)小的機率一樣少見，那麼可以拒斥\(\tt{H}_{0}\)。
    7. 否則的話，如果我們沒觀察到與\(\alpha\)一樣小的拒斥值，無法拒斥\(\tt {H}_{0}\)。

    3 迴歸係數的檢定

  • 我們抽出一套樣本得到的迴歸係數可能高估母體的係數，但是下一套可能低估，不斷地抽樣之後，所得到的迴歸係數的平均值理論上等於母體的係數。
  • 但是我們只有一套樣本，所以只能假設估計的係數是無偏估計。
  • 如何估計OLS樣本統計的精確程度？也就是OLS樣本統計的變異數？利用殘差項：

    • \[ \hat{u}_{i}=y_{i}-\hat{y}_{i}=y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{i} \]

  • 估計\(Var[u]=\sigma^{2}_{u}\)，假設殘差項的分佈可以用Mean Squared Deviation (MSD)表示：

    • \[ MSD(\hat{u})\equiv\frac{1}{n}\sum\limits_{i=1}^{n}(\hat{u}_{i}-\bar{\hat{u}})^2=\frac{1}{n}\sum\limits_{i=1}^{n}\hat{u}_{i}^2 \]

  • 誤差項的變異數的無偏估計為：

    • \[ \begin{eqnarray} E[MSD(\hat{u})] & = & \frac{n-2}{n}\sigma^{2}_{u}\\ \sigma^{2}_{u} & = &\frac{n}{n-2}MSD(\hat{u}) \\ & =& \frac{n}{n-2}\cdot\frac{1}{n}\sum\limits_{i=1}^{n}\hat{u}_{i}^2 \\ & =& \frac{1}{n-2}\sum\limits_{i=1}^{n}\hat{u}_{i}^2 \end{eqnarray} \]

  • 因為估計\(\hat{u}_{i}\)需要\(\hat{\beta}_{0}\), \(\hat{\beta}_{1}\) 所以自由度為\(n-2\)

    • \(\frac{1}{n-2}\sum\limits_{i=1}^{n}\hat{u}_{i}^2\)代入： \[ \begin{eqnarray} $\hat{V}[\hat{\beta}_{1}|X] & = &\mathcal{\frac{\sigma_{u}^{2}}{\sum_{i=1}^{n}(x-\bar{x})^2}} \\ & = & \frac{\sigma_{u}^{2}}{SST_{x}} \end{eqnarray} \] \[ \begin{eqnarray} \hat{V}[\hat{\beta}_{0}|X] & = &\sigma_{u}^{2}\mathcal{\{\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^{n}(x-\bar{x})^2}}\}\\ & = &\mathcal{\frac{\sigma_{u}^2\sum_{i=1}^{n}x^2}{n\sum_{i=1}^{n}(x-\bar{x})^2}} \end{eqnarray} \] where \[ \sigma_{u}^{2}=\frac{1}{n-2}\sum\limits_{i=1}^{n}\hat{u}_{i}^2 \]

    \(\blacksquare\) 取開根號：\(\sqrt{\hat{V}[\hat{\beta}_{0}|X]}\)以及\(\sqrt{\hat{V}[\hat{\beta}_{1}|X]}\)
  • 最小平方法迴歸必須假設誤差項與自變項沒有相關，而且假設平均數為0的常態分佈：

    • \[ u\sim N(0,\sigma_{u}^{2}) \]

  • 以上假設意義為：

    • \[ \mathrm{Y}|\mathrm{X} \sim N(\beta_{0}+\beta_{1}X, \sigma_{u}^{2}) \]

    也就是變異數同質性以及誤差項的平均值為0的假設。

    4 標準誤

  • \(\beta_{0}\)與\(\beta_{1}\)的標準誤分別是：
    • \[ \mathrm{SE}(\hat{\beta_{0}})^2=\sigma^2[\frac{1}{n}+\frac{\bar{x}^2}{\Sigma_{i=1}^n(x-\bar{x})^2}] \] \[\mathrm{SE}(\hat{\beta_{1}})^2=\frac{\sigma^2}{\Sigma_{i=1}^n(x-\bar{x})^2} \]
  • 我們無法計算母體的\(\sigma^2\)，只能用\(\sqrt{RSS/(n-2)}\)估計，也就是residual standard error。如果我們估計\(\hat{\beta}_{1}\)，可代入計算\(\hat{Y}_{i}\)，以及計算\(\mathrm{SE}(\hat{\beta_{1}})\)為：

    • \[ \sigma^2=\frac{\Sigma(y-\hat{y})^2}{n-2} \] \[ \mathrm{SE}(\hat{\beta_{1}})=\sqrt{\frac{\sigma^2}{\sum(x_{i}-\bar{x})^2}} \]

    5 信賴區間

    \(95\%\)的信賴區間相當於兩個標準誤，係數加減兩個標準誤構成一個\(95\%\)的信賴區間：
    \[ \mathcal{[\hat{\beta_{1}}-2\cdot \mathrm{SE}(\hat{\beta_{1}}), \hspace{.3em}\hat{\beta_{1}}+2\cdot \mathrm{SE}(\hat{\beta_{1}})]} \]

    假設檢定

    \(\it{H}_{0}\): X與Y沒有關係，也就是\(\it{H}_{0}\): \(\beta_{1}=0\)

    \(\it{H}_{a}\): X與Y有某種關係，也就是\(\it{H}_{a}\): \(\beta_{1} \neq 0\)

  • 為了測試\(\hat{\beta_{1}}\)是否等於0，我們必須考慮\(\hat{\beta_{1}}\)的變異數，而我們用\(t\)分佈來檢驗：

    • \[ t=\frac{\hat{\beta}_{1}-0}{\mathrm{SE}(\hat{\beta}_{1})} \]

    5.1 \(t\)分佈的計算

  • R可以計算\(t\)分布的百分位的機率，也可以計算機率對應的百分位。首先計算特定機率的百分位： 
    • qt(0.975, 1000)
    ## [1] 1.962339
    qt(0.950, 1000)
    ## [1] 1.646379
    qt(0.900, 1000)
    ## [1] 1.282399
    qt(0.5, 1000)
    ## [1] 0
    qt(0.1, 1000)
    ## [1] -1.282399
    qt(0.05, 1000)
    ## [1] -1.646379
    qt(0.025, 1000)
    ## [1] -1.962339
  • 其次，計算機率對應的百分位： 
    • pt(1.68, 1000)
    ## [1] 0.9533652
    pt(1.96, 1000)
    ## [1] 0.9748634
    pt(2.00, 1000)
    ## [1] 0.9771148

    例題一

    請問大學教師的薪水與性別有關嗎？使用car裡面的Salaries資料：

    library(car)
    ## 
## Attaching package: 'car'
    ## The following object is masked from 'package:dplyr':
## 
##     recode
    fit1 <- with(Salaries, lm(salary ~ sex))
summary(fit1)
    ## 
## Call:
## lm(formula = salary ~ sex)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -57290 -23502  -6828  19710 116455 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   101002       4809  21.001  < 2e-16 ***
## sexMale        14088       5065   2.782  0.00567 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30030 on 395 degrees of freedom
## Multiple R-squared:  0.01921,    Adjusted R-squared:  0.01673 
## F-statistic: 7.738 on 1 and 395 DF,  p-value: 0.005667
    1. 當教師為男性，平均薪水為101003+14088=115090，當教師為女性，平均薪水為101002。
    2. 性別的迴歸係數\(t\)值為2.782，因為 \[ \frac{14088}{5065}=2.78 \]
    3. 當\(t\)=2.78，我們想知道對應到的機率值有多高？如果小於虛無假設的機率值，也就是我們觀察到虛無假設成立的機會非常小，也就是我們可以拒斥虛無假設。
    h0=1-pt(1.96, 395)
h1=1-pt(2.78, 395)
h0; h1
    ## [1] 0.02534923
    ## [1] 0.002847915
  • 因此，我們可以拒斥兩者無關的虛無假設，並且得到性別會決定教師薪水的結論。

    • 例題二

    請問氣溫與風力大小、月份有關嗎？我們分析airquality資料：

    fit2 <- with(airquality, lm(Temp ~ Wind+factor(Month)))
summary(fit2)
    ## 
## Call:
## lm(formula = Temp ~ Wind + factor(Month))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8929  -4.4212   0.0799   3.4794  17.8880 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     74.1885     2.0544  36.112  < 2e-16 ***
## Wind            -0.7434     0.1488  -4.996 1.64e-06 ***
## factor(Month)6  12.5436     1.5943   7.868 7.18e-13 ***
## factor(Month)7  16.3621     1.6183  10.110  < 2e-16 ***
## factor(Month)8  16.3163     1.6239  10.047  < 2e-16 ***
## factor(Month)9  10.2792     1.5959   6.441 1.59e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.175 on 147 degrees of freedom
## Multiple R-squared:  0.5884, Adjusted R-squared:  0.5744 
## F-statistic: 42.03 on 5 and 147 DF,  p-value: < 2.2e-16
    1. 風力每增加一個單位，溫度減少0.7434個單位。標準誤為0.1488。\(t\)值為-4.996，拒斥風力與氣溫沒有關係的虛無假設。
    2. 相對於五月，六月的溫度增加12.5436度。標準誤為1.5943。

    6 作業

    1. 請根據以下的公式，計算\(\beta_{1}\)以及標準誤。
    \[ \hat{\beta_{1}}=\frac{\sum (x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum(y_{i}-\bar{y})^2} \] \[ \sigma^2=\frac{\Sigma(y-\hat{y})^2}{n-2} \] \[ \mathrm{SE}(\hat{\beta_{1}})=\sqrt{\frac{\sigma^2}{\sum(x_{i}-\bar{x})^2}} \] 

    library(kableExtra)
DT<-data.frame(X=c(4, 3, 5, 2, 4, 2, 2, 3, 2, 2, 2, 3, 5, 1, 1),
               Y=c(5, 5, 5, 3, 4, 3, 3, 4, 4, 5, 4, 5, 3, 2, 1))
DT %>%
  kable("html") %>%
  kable_styling()
    X Y
    4 5
    3 5
    5 5
    2 3
    4 4
    2 3
    2 3
    3 4
    2 4
    2 5
    2 4
    3 5
    5 3
    1 2
    1 1

    2. 接續上一題，如果增加了以下五筆資料，請重新估計迴歸模型：

    DT2<-data.frame(X=c(2, 1, 3, 4, 5),
               Y=c(3, 4, 5, 4, 4))
DT2 %>%
  kable("html") %>%
  kable_styling()
    X Y
    2 3
    1 4
    3 5
    4 4
    5 4