社會科學統計方法
迴歸係數假設檢定
迴歸係數假設檢定
蔡佳泓
2024
1 課程目標
2 迴歸假設
2.1 描述與推論
- \(E(Y|X)\)是我們要估計的期望值，在X固定的情況下，解釋Y的期望值。
- 最小平方法迴歸有一連串的假設必須滿足。在這些假設成立的前提下，進行區間估計、假設檢定等等。
- 母體迴歸模型可表示為：\(Y=\beta_{0}+\beta_{1}X+ \epsilon\)，其中Y是依變數(outcome)，X是自變數(regressor)，\(\beta_{0}\),\(\beta_{1}\)分別是截距以及斜率，\(\epsilon\)是誤差項，也是一個捕捉到所有會影響Y但不影響X的變數。
- 誤差(\(\epsilon\))也可說是Y與\(E(Y|X)\)之間的差距，也就是預測值(predicted value)與真實的變數值之間的差距，我們無法觀察到\(\epsilon\)。有人用\(u\)代表誤差，\(u\)來自residual，定義成\(Y-\hat{Y}\)。
- 因此，樣本估計的迴歸模型表示為：\(Y=\hat{\beta}_{0}+\hat{\beta}_{1}X+\hat{u}\)， \(\hat{\beta}_{0}\), \(\hat{\beta}_{1}\)分別是截距以及斜率。 \(\hat{u}\)代表殘差項，是對於誤差的估計，也就是X無法解釋的Y的部分。
2.2 最小平方法的假設
- 線性：參數是以線性的方式呈現(\(\color{rgb:red,4;yellow,3}{\mathrm {Linearity\hspace{.3em} in\hspace{.3em} parameters}}\))
- 隨機抽樣：資料以隨機的方式從母體抽出(\(\color{rgb:red,4;yellow,3}{\mathrm {Random \hspace{.3em}sampling}}\))
- 自變數X有變異量(\(\color{rgb:red,4;yellow,3}{\mathrm {Variation\hspace{.3em} in\hspace{.3em} X}}\))
- 殘差對X任何一值的條件機率平均值為0(\(\color{rgb:red,4;yellow,3}{\mathrm {Zero\hspace{.3em} conditional\hspace{.3em} mean}}\))
- 殘差對X任何一值的變異數相等(\(\color{rgb:red,4;yellow,3}{\mathrm {Homoskedasticity}}\))
- 殘差與X互相獨立而且呈常態分佈(\(\color{rgb:red,4;yellow,3}{\mathrm {Normality}}\))
- 殘差之間沒有自我相關，也就是\(\rm{Cov}(\epsilon_{i},\epsilon_{j})=0\quad i\neq j\)。
2.2.1 線性
2.2.2 隨機抽樣
- 會得到一個\(R^2=1\)的模型，無法估計標準誤。
- 無法預測新的資料。
- 資料太少使得參數估計會因為少數資料變動而不穩定。
2.2.3 自變數X有變異量
迴歸係數\(\hat{\beta}_{1}\)可用X的變異量與XY的共變量來表示，也就是寫成：
\[\hat{\beta_{1}}=\frac{\sum_{i=1}^n(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum_{i=1}^n(x_{i}-\bar{x})^2}\]
- 所以X必須有變異量，否則\(\hat{\beta_{1}}\)的分母等於0。
2.2.4 誤差對X任何一值的條件機率平均值為0
\(\epsilon\)代表所有影響Y，但是觀察不到的變數。如果有一個未觀察的變數，同時影響Y以及X，就違反這個假設。
例如薪資=\(\beta_{0}+\beta_{1}\)教育程度+\(\epsilon\)
在這個模型中，沒有觀察「能力」（\(q_{i}\)）而進入誤差項，\(\epsilon_{i}=\gamma q_{i}+v_{i}\)。如果\(Cov(q_{i},X_{i})\neq 0\), \(Cov(u_{i},X_{i})\neq 0\)。
如果\(\gamma>0\)，\(Cov(q_{i},X_{i})>0\)，迴歸係數估計會變大，也就是高估自變數與依變數的關係。
也被稱為內生性問題。
set.seed(02138)
X=rbinom(100, 20, prob=0.45)
Y=X+runif(100, 0, 0.5)
e=ifelse(X>=9, 1, 0)
OLS1=lm(Y ~ X)
#summary(OLS1)
OLS2=lm(Y+e ~ X)
#summary(OLS2)
stargazer(OLS1,OLS2,align = T,
title='不同依變數迴歸模型比較',
type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | ||
| Y | Y + e | |
| (1) | (2) | |
| X | 1.004*** | 1.177*** |
| (0.006) | (0.014) | |
| Constant | 0.213*** | -0.736*** |
| (0.057) | (0.127) | |
| Observations | 100 | 100 |
| R2 | 0.996 | 0.987 |
| Adjusted R2 | 0.996 | 0.986 |
| Residual Std. Error (df = 98) | 0.143 | 0.317 |
| F Statistic (df = 1; 98) | 25,878.000*** | 7,206.000*** |
| Note: | p<0.1; p<0.05; p<0.01 | |
- 根據Sandford Weisberg的Applied Linear Regression(p.36-37)，我們用
R的lm 函式所附的檢驗功能畫圖檢視預測值以及殘差項是否呈現特定趨勢，如果有的話，表示模型違反誤差的變異數相等之假設，但是我們也可以用來檢驗誤差的平均值是否為0的假設。圖2.1顯示，第一個依變數沒有其他變數作用的模型，並沒有違反誤差的平均值為0的假設：
plot(OLS1,1, yaxt='n')
axis(2, c(-0.6, -0.2, 0, 0.2, 0.6))
abline(h = 0, lty=2)
Figure 2.1: 條件平均值為0檢驗圖1
- 圖2.2則顯示，第二個依變數有其他變數作用的模型，可能違反誤差的平均值為0的假設，因為當X小於6或者大於12時，平均值也大於0或者小於0：
plot(OLS2,1, yaxt='n')
axis(2, c(-0.6, -0.2, 0, 0.2, 0.6))
abline(h = 0, lty=2)
Figure 2.2: 條件平均值為0檢驗圖2
2.2.5 誤差對X任何一值的變異數相等
- 以課本習題的資料為例，圖2.3顯示，變異數相等的假設似乎不成立：
X<-c(19,17,21,18,15,12,14,20)
Y<-c(1,3,1,1,2,3,2,1)
M1<- lm(Y ~ X); summary(M1)
##
## Call:
## lm(formula = Y ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.529 -0.335 -0.140 0.136 1.250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5000 1.2640 4.35 0.0048 **
## X -0.2206 0.0733 -3.01 0.0237 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.604 on 6 degrees of freedom
## Multiple R-squared: 0.602, Adjusted R-squared: 0.535
## F-statistic: 9.06 on 1 and 6 DF, p-value: 0.0237
car::spreadLevelPlot(M1)
Figure 2.3: 變異數相等假設檢驗圖1
##
## Suggested power transformation: 0.3847- 改用\(\tt{faithful}\)這筆資料，圖2.4顯示，變異數相等的假設似乎成立。
car::spreadLevelPlot(lm(eruptions~waiting, faithful))
Figure 2.4: 變異數相等假設檢驗圖2
##
## Suggested power transformation: 0.8039如果違反這個假設，標準誤有可能膨脹，也就是並不有效。
迴歸係數的條件變異數表示如下（後續說明如何得出）：
\[\boxed{\rm{Var}[\hat{\beta}_{1}|X]={\frac{\sigma_{u}^{2}}{\sum_{i=1}^{n}(x-\bar{x})^2}}=\frac{\sigma_{u}^{2}}{SST_{x}}}\]
\[\boxed{ \rm{Var}[\hat{\beta}_{0}|X] =\sigma_{u}^{2}(\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^{n}(x-\bar{x})^2}) =\frac{\sigma^2\sum_{i=1}^{n}x^2}{\sum_{i=1}^{n}(x-\bar{x})^2}}\]
\[\rm{where}\hspace{.4em}\rm{Var}[u|X]=\sigma_{u}^{2}\]
其中\(u\)代表殘差。
當\(\sigma_{u}^{2}\)小的時候，\(Var[\hat{\beta}_{1}|X]\)也會變小，也就是未知的變數對於係數的影響越小，離真實的迴歸線可能越近。
因為\(\sum_{i=1}^{n}(x-\bar{x})^2=(n-1)S_{x}^{2}\)，所以 n越大，\(Var[\hat{\beta}_{1}|X]\)越小，可能越近真實的迴歸線。
\(S_{x}^{2}\)越大，\(Var[\hat{\beta}_{1}|X]\)越小，離真實的迴歸線可能越近。
\(Var[\hat{\beta}_{1}|X]\)是\(Var[\beta_{1}|X]\)的估計之一，有可能因為n很大或者是\(Var[X]\)很大，使得\(Var[\hat{\beta}_{1}|X]\)變小
圖2.5似乎顯示殘差值在不同Y固定X的情況下，並不是隨機分佈。
m1 <- lm(ment ~ phd, data=pscl::bioChemists)
plot(fitted(m1), resid(m1))
abline(0,0, col='#FF11CC')
Figure 2.5: 殘差值與預測值散佈圖1
- \(\tt{olsrr}\)套件提供許多檢查迴歸假設的指令，例如我們用\(\tt{ols\_plot\_reside\_fit}\)這個函式畫預測值對殘差值的圖：
library(olsrr)
ols_plot_resid_fit(m1)
Figure 2.6: 用olsrr畫殘差值與預測值散佈圖1
set.seed(02138)
X<-rnorm(100, 0, 1); Y<-X+rnorm(100, 0, 1.1)
M1 <- lm(Y~X)
plot(fitted(M1), resid(M1), col='#22CDC0')
abline(0,0, col='#AA1100')
Figure 2.7: 殘差值與預測值散佈圖2
- 模擬另一筆資料，\(Y=X^2\)，殘差值與預測值散佈圖就不是隨機分佈，如圖2.8。
set.seed(02138)
X<-rnorm(100, 0, 1); Y<-X^2+runif(100, 0, 2)
M1 <- lm(Y~X)
plot(fitted(M1), resid(M1), col='#AE13C0')
abline(0,0, col='#AA1100')
Figure 2.8: 殘差值與預測值散佈圖4
2.2.6 誤差與X互相獨立而且呈常態分佈
\[Y|X \sim N(\beta_{0}+\beta_{1}X, \hspace{.4em}\sigma_{u}^{2})\]
- 也就是變異數同質性以及誤差項的平均值為0的假設。
image1<-here::here('Fig','TikZ_v3.png')
knitr::include_graphics(image1)
Figure 2.9: 誤差項的變異數圖
\[\hat{\beta}_{1}-\beta_{1}/\sqrt{\rm{Var[\hat{\beta}_{1}|X]}}=\frac{\hat{\beta}_{1}-\beta_{1}}{s.e._{\hat{\beta}}}\sim N(0,1)\]
- 或者：
\[\hat{\beta}\sim N(\beta_{1},\hspace{.3em}\rm{Var}[\hat{\beta_{1}}|X])\]
- 其中：
\[\rm{Var}[\hat{\beta_{1}}|X])=\sigma_{u}^2/\sum\limits_{i=1}^{n}(x_{i}-\bar{x})^2\]
3 中央極限定理(Central Limit Theorem)
#Two sampling
set.seed(02139)
datam1<-rep(NA,1000)
for (i in 1:1000){
datam1[i]<-mean(rnorm(100, 0, 1))
}
mean2<-function(x){
sum(x)/(length(x)+1)
}
datam2<-rep(NA,1000)
for (i in 1:1000){
datam2[i]<-mean2(rnorm(100, 0, 1))
}
#Two graphs
par(mfrow=c(1,2))
hist(datam1,xlab=expression(paste(hat(mu)[1])),
main=expression(paste(hat(mu)[1]==frac(1,n),sum(X[i],i=1,N))),prob=T,
breaks=seq(-1, 1,by=0.01))
hist(datam2,xlab=expression(paste(hat(mu)[2])),
main=expression(paste(hat(mu)[2]==frac(1,n+1),sum(X[i],i=1,N))),prob=T,
breaks=seq(-1, 1,by=0.01))
Figure 3.1: 常態分佈抽出的樣本平均值分佈
par(mfrow=c(2,2))
set.seed(02138)
res <-c()
s <- c()
simulations <- c(100, 1000)
samplesize <- c(100, 300)
for (i in 1:length(simulations)){
for (j in 1:length(samplesize)){
f=function(n = samplesize[j], mu=1, sigma = 1){
s = rnorm(samplesize[j], 1, 1)
(mean(s)-mu)/sqrt((sum(s-mean(s))^2)/sqrt(n))
}
#xvals = seq(-1, 1, 0.1)
hist(UsingR::simple.sim(simulations[i], f, samplesize[j], 1), probability=TRUE, breaks = 40, col='#FF9933', main = paste("N =", samplesize[j], ", Num. of Samples=", simulations[i]),
xlab=expression(mu))
# xvals = seq(-1, 1, .01)
#points(xvals, dnorm(xvals, 0, 1),type="l", lwd=2, col='#FF9933')
}
}
Figure 3.2: 不同樣本規模的樣本平均值分佈
# Set the parameters
lambda <- 0.5 # Rate parameter of the exponential distribution
sample_size <- 100 # Sample size
num_simulations <- 200 # Number of simulations
# Function to estimate lambda from a sample
estimate_lambda <- function(data) {
mean(data) # Estimator of lambda for exponential #distribution is the sample mean (1/lambda)
}
# Simulate data and estimate parameters
lambda_estimates <- replicate(num_simulations, {
# Generate a sample from exponential distribution
sample_data <- rexp(sample_size, rate = lambda)
# Estimate lambda from the sample
estimate_lambda(sample_data)
})
# Plot the distribution of estimated lambdas
hist(lambda_estimates, breaks = 30, main = "Distribution of 200 Estimated Lambdas",
xlab = "Estimated Lambda", col = "#00EEAA", border = "white")
Figure 3.3: 指數分佈抽出的200個樣本平均值分佈
# Set the parameters
lambda <- 0.5 # Rate parameter of the exponential distribution
sample_size <- 100 # Sample size
num_simulations <- 1000 # Number of simulations
# Function to estimate lambda from a sample
estimate_lambda <- function(data) {
mean(data) # Estimator of lambda for exponential #distribution is the sample mean (1/lambda)
}
# Simulate data and estimate parameters
lambda_estimates <- replicate(num_simulations, {
# Generate a sample from exponential distribution
sample_data <- rexp(sample_size, rate = lambda)
# Estimate lambda from the sample
estimate_lambda(sample_data)
})
# Plot the distribution of estimated lambdas
hist(lambda_estimates, breaks = 30, main = "Distribution of 1000 Estimated Lambdas",
xlab = "Estimated Lambda", col = "#112ee1", border = "white")
Figure 3.4: 指數分佈抽出的1000個樣本平均值分佈
\(V[Y]=V[\sum X_{i}]=\sum V[X_{i}]=np(1-p)\)
\(V[\frac{Y}{n}]=\frac{1}{n^2}V[Y]=\frac{1}{n^2}np(1-p)=\frac{p(1-p)}{n}\)
- 因此標準誤為\(\sqrt{\frac{p(1-p)}{n}}\)
set.seed(2019)
results =numeric(0)
k=100; p=0.25; mu=k*p # k is number of trials
for (i in 1:1000) {
S = rbinom(1, k, p)
results[i]=(S-mu)/sqrt(k*p*(1-p)) #
}
hist(results, probability = T, col='#3399FF', breaks=30, main='')
xvals = seq(-3,3,.01)
points(xvals,dnorm(xvals,0,1),type="l", lwd=2, col='#FF9933') 
Figure 3.5: 二元分佈抽出的樣本平均值分佈
3.1 迴歸係數的常態分佈
假設有迴歸模型如下：
\[Y=\hat{\beta}_{0} - \hat{\beta}_{1}X_{1}+u\]
\[\begin{align} Y=5 - 1x+u \tag{3.1} \end{align}\]
- 其中
\[u\sim N(0,4)\]
- 為了顯示\(\hat{\beta}_{0}, \hat{\beta}_{1}\)呈現常態分佈，我們先根據\(u\sim N(0,4)\)模擬200次，每次200個\(u_{i}\)。
- 我們也用均等分佈模擬X\(\sim Unif(1,10)\)。
- 將每一群模擬的\(u_{i}\)放入(3.1)得到\(Y_{i}\)，估計 \(\hat{\beta}_{0}\), \(\hat{\beta}_{1}\)。
set.seed(02138)
y<-rep(NA,40000)
x<-matrix(NA,4000,c(200,200))
for(j in 1:200){
x[,j]<-runif(200,1,10) }
ru<-matrix(NA,4000,c(200,200))
for(j in 1:200){
ru[,j]<-rnorm(200,0,4)}
y<-5-x+ru
lm.beta0<-rep(NA,200);lm.beta1<-rep(NA,200)
for(j in 1:200){ lm.beta0[j]<-lm(y[,j]~x[,j])$coefficients[1]
lm.beta1[j]<-lm(y[,j]~x[,j])$coefficients[2] }- 我們把模擬得到的\(\beta_{0}\)、\(\beta_{1}\)分別畫成長條圖，結果成常態分配。
dt <- data.frame(b1=lm.beta1, b0=lm.beta0)
b1<-ggplot2::ggplot(data=dt,aes(x=b1))+
geom_histogram(bins=40, fill='#FF6666')+
labs(x=expression(paste(beta,'1')))+
geom_density(alpha=.9)
b1
Figure 3.6: 迴歸係數長條圖1
b0<-ggplot2::ggplot(data=dt,aes(x=b0))+
geom_histogram(bins=40, fill='#FF1111')+
labs(x=expression(paste(beta,'0')))+
geom_density(alpha=.9)
b0
Figure 3.7: 迴歸係數長條圖2
4 假設檢定
步驟
- 設定虛無假設 \(\tt{H}_{0}\)
- 決定\(\alpha=\)Pr(否定\(\tt{H}_{0}\) | \(\tt{H}_{0}\))的大小。
- 根據樣本統計以及\(\tt{H}_{0}\)選擇\(t\)值。
- 假設\(\tt{H}_{0}\)為真，建立\(T\)的虛無分佈（常態分佈）。
- 從分佈表查拒斥值(critical values)，估計\(T\)觀察值有多難出現（被觀察到）。
- 如果從上述分佈觀察到\(T\)的機率，跟\(t\)值比拒斥值\(\alpha\)小的機率一樣少見，那麼可以拒斥\(\tt{H}_{0}\)。
- 否則的話，如果我們沒觀察到與\(\alpha\)一樣小的拒斥值，無法拒斥\(\tt {H}_{0}\)。
5 迴歸係數的檢定
\[ \hat{u}_{i}=y_{i}-\hat{y}_{i}=y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{i} \]
\[ 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 \]
- 這裏MSD是一種變異數的表示方式，只是分母用\(n\)而不是\(n-1\)。
\[ \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} \]
- 殘差的變異數的無偏估計表示為：
\[\begin{align*} \sigma^{2}_{u}=\frac{1}{n-2}\sum\limits_{i=1}^{n}\hat{u}_{i}^2 \tag{5.1} \end{align*}\]
5.1 迴歸係數與標準誤推導過程
5.1.1 迴歸係數
- 這是因為：
\[\begin{align*} \sum(x_{i}-\bar{x})\bar{y}&= (\sum x_{i}-n\bar{x})\bar{y}\\ &=\sum x_{i}\bar{y}-n\frac{\sum x_{i}}{n}\bar{y}\\ &=\sum x_{i}\bar{y}-\sum x_{i}\bar{y}\\ & = 0 \tag{5.2} \end{align*}\]
- 在(5.2)中，使用到總和(summation)的代數運算的公設之一：
\[\sum_{i=1}^{n} c=n\cdot c\]
- 所以：
\[\sum_{i=1}^{n} \bar{X}=n\cdot \bar{X}\]
\[\begin{align*} \hat{\beta}_{1}=\beta_{1}+\frac{\sum (x_{i}-\bar{x})u_{i}}{SST_{x}} \tag{5.3} \end{align*}\]
- 如果誤差項的平均數為0此一假設成立，那麼\(\hat{\beta}_{1}\)是\(\beta_{1}\)的無偏估計：
\[\begin{align*} E[\hat{\beta}_{1}|X]&=E[\beta_{1}|X]+E[\frac{\sum (x_{i}-\bar{x})u_{i}}{SST_{x}}|X]\\ &=\beta_{1}+\sum_{i=1}^{n}\frac{\sum (x_{i}-\bar{x})}{SST_{x}}E[u_{i}|X] \\ &=\beta_{1} \end{align*}\]
- 根據Law of Iterated Expectations，證明\(\hat{\beta}_{1}\)是\(\beta_{1}\)的無偏估計:
\[E[\hat{\beta}_{1}]=E[E[\hat{\beta}_{1}|X]]=\beta_{1}\]
- 在迴歸課程系列一提到可以求出\(\hat{\beta_{0}}\)：
\[ \hat{\beta_{0}} = \bar{y} - \hat{\beta_{1}}\bar{x} \]
- 連帶地可以求出\(\hat{y}\):
\[ \hat{y} = \hat{\beta_{0}} + \hat{\beta_{1}}x_{i} \]
5.1.2 迴歸係數標準誤(OLS estimator variance)
\[\hat{\beta}_{1}=\beta_{1}+\frac{\sum (x_{i}-\bar{x})u_{i}}{SST_{x}}\]
- 我們要求出\(\hat{\beta}_{1}\)的條件變異數，所以在上述式子左右各取條件變異數：
\[\begin{align*} V[\hat{\beta}_{1}|X]&=V[\beta_{1}|X]+V[\frac{\sum (x_{i}-\bar{x})u_{i}}{SST_{x}}|X]\\ &=0 + \sum \{\frac{(x_{i}-\bar{x})u_{i}}{SST_{x}}\}^2\hspace{.1em}\rm{Var}[u_{i}|X]\\ &=\frac{\sum (x-\bar{x})^2}{SST_{x}^2}\sigma_{u}^2\\ &=\frac{SST_{x}}{SST_{x}^2}\sigma_{u}^2=\frac{\sigma_{u}^2}{SST_{x}} \tag{5.4} \end{align*}\]
在上面的式子(5.4)的簡化過程中，使用到變異數的特性。因為\(\beta_{1}\)是常數，沒有變異數可言，所以Var[c+x] = Var[c]+Var[x]=0+Var[x]。而且Var\([c\times x]=c^2\)Var[x]。
根據式(5.4)，迴歸係數的條件變異數寫成：
\[\begin{align*} \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}} \tag{5.5} \end{align*}\]
\[\begin{align*} \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}} \tag{5.6} \end{align*}\]
- 根據(5.1)， \[ \sigma_{u}^{2}=\frac{1}{n-2}\sum\limits_{i=1}^{n}\hat{u}_{i}^2 \]
\(\blacksquare\) 上面的方程式(5.6)、(5.5)分別取開根號得到標準誤：\(\sqrt{\hat{V}[\hat{\beta}_{0}|X]}\)以及\(\sqrt{\hat{V}[\hat{\beta}_{1}|X]}\)
\[ u\sim N(0,\sigma_{u}^{2}) \]
\[ \mathrm{Y}|\mathrm{X} \sim N(\beta_{0}+\beta_{1}X, \sigma_{u}^{2}) \]
- 也就是變異數同質性以及誤差項的平均值為0的假設。
5.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}} \]
- 換句話說，\(\hat{\beta_{1}}\)的常態分配如下：
\[\hat{\beta_{1}}\sim N\Big(\beta_{1},\frac{\sigma^2}{\sum (x_{i}-\bar{x})^2}\Big)\]
- 另外，\(\hat{\beta_{0}}\)的常態分配如下：
\[\hat{\beta_{0}}\sim N\Big(\beta_{0},\frac{\sum x_{i}^2\sigma^2}{n\sum (x_{i}-\bar{x})^2}\Big)\]
5.3 迴歸係數的信賴區間
- 因為\(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\)
\[ t=\frac{\hat{\beta}_{1}-0}{\mathrm{SE}(\hat{\beta}_{1})} \]
- 圖5.1表示\(t\)分佈的雙尾拒斥域，也就是顯著水準：
par(bty='n', yaxt='n')
cord.x<-c(1.962, seq(1.962,3,0.01), 3)
cord.x2<-c(-3, seq(-3, -1.962,0.01), -1.962)
cord.y<-c(0, dt(seq(1.962, 3, 0.01), 1000), 0)
cord.z<- c(0, dt(seq(-3, -1.962, 0.01), 1000), 0)
curve(dt(x, 1000),xlim=c(-3,3),
main=expression(paste("Normal Density with"," ", t[alpha/2]=="0.025")) , ylab='', xlab='t value')
polygon(cord.x, cord.y, col="red3")
polygon(cord.x2, cord.z, col="red3" , add=T)
Figure 5.1: t分佈的雙尾拒斥域
par(bty='n', yaxt='n')
cord.x<-c(1.646, seq(1.646,3,0.01), 3)
cord.y<-c(0, dt(seq(1.646, 3, 0.01), 1000), 0)
curve(dt(x, 1000),xlim=c(-3,3),
main=expression(paste("Normal Density with"," ", t[alpha/2]=="0.05")) , ylab='', xlab='t value')
polygon(cord.x, cord.y, col="red3")
Figure 5.2: t分佈的右尾拒斥域
plot(function(x) dt(x, df = 300), -3, 3, ylim = c(0, .5),
main = "t - Density", yaxs = "i", col="white",ylab="Density",
xlab=expression(paste(X %~% tau[n])))
curve(dt(x, df = 3), -3,3, bty="l", col="blue", add=T, lwd=2)
curve(dt(x, df = 1), -3,3, bty="l", col="black", add=T, lwd=2)
curve(dt(x, df = 1000), -3,3, bty="l", col="red", add=T, lwd=2)
curve(dt(x, df = 15), -3,3, bty="l", col="green", add=T, lwd=2)
legend("topright", lty=c(1,1,1,1), lwd=c(2,2,2,2),
legend=c(expression(nu==1, nu==3, nu==15,nu==1000)),
col=c("black", "blue", "green","red"))
Figure 5.3: 不同自由度的t分佈
- 根據以上的\(t\)分佈，等我們計算標準誤之後，就可以檢定統計量是否落在雙尾或者單尾的拒斥域。
- 如果統計量\(T\)大於1.96，因為\(\texttt{pt(1.96, df = 1000)}\) \(\approx 0.975\)，那麼出現這個\(T\)值的機率應該小於0.05。當\(p<0.05\)，可以解釋為在95%的信心水準下，\(\hat{\beta}_{1}\neq 0\)，也就是得到拒斥\(\textit{H}_{0}: \beta_{1}=0\)的虛無假設的結論。
5.4 迴歸係數估計的自由度
5.5 \(t\)分佈的計算
R計算\(t\)分布的百分位的機率，也可以計算機率對應的百分位。首先計算特定機率的百分位如表5.1：
qtd <-data.frame(P=c('95%','97.5%','99%'),
t=c(qt(0.950, 1000),qt(0.975, 1000),qt(0.990, 1000)))
newqtd <- qtd %>% knitr::kable("html", caption='累積機率對應的t值') %>%
kableExtra::kable_styling(bootstrap_options = 'striped')
newqtd| P | t |
|---|---|
| 95% | 1.646 |
| 97.5% | 1.962 |
| 99% | 2.330 |
- 其次，計算機率對應的百分位如表5.2：
tt <- data.frame(p_168=pt(1.68, 1000), p_196=pt(1.96, 1000), p_300=pt(3.00, 1000))
colnames(tt)<-c("t=1.68","t=1.96","t=3.00")
newtt <- tt %>% kable("html", caption='t值對應的百分位') %>%
kable_styling(bootstrap_options = 'striped')
newtt| t=1.68 | t=1.96 | t=3.00 |
|---|---|---|
| 0.9534 | 0.9749 | 0.9986 |
6 例題
6.1 房價與房間數迴歸模型與檢定
- 房間數與房價有關嗎？表 6.1顯示房價（單位：千元）與房間數的資料。
TAB <- data.frame(price=c(300,250,400,550,317,389,425,289,389,559),
bedroom=c(3,3,4,5,4,3,6,3,4,5))
knitr::kable(TAB, format='simple', caption='房價與房間數資料')| price | bedroom |
|---|---|
| 300 | 3 |
| 250 | 3 |
| 400 | 4 |
| 550 | 5 |
| 317 | 4 |
| 389 | 3 |
| 425 | 6 |
| 289 | 3 |
| 389 | 4 |
| 559 | 5 |
- 估計OLS迴歸模型如表 6.2：
pricemodel<-lm(price ~ bedroom, data=TAB)
stargazer::stargazer(pricemodel, title='房價與房間數模型',
type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| price | |
| bedroom | 73.100** |
| (23.760) | |
| Constant | 94.400 |
| (97.980) | |
| Observations | 10 |
| R2 | 0.542 |
| Adjusted R2 | 0.485 |
| Residual Std. Error | 75.150 (df = 8) |
| F Statistic | 9.462** (df = 1; 8) |
| Note: | p<0.1; p<0.05; p<0.01 |
- 請問「每增加一間房間就可以多賣100,000元」的假設對嗎？
- 虛無假設是：\(\textit{H}_{0}:\beta_1\le 100\)。對立假設是：\(\textit{H}_{A}:\beta_1> 100\)。
- 假設m=100，計算T值等於(\(\hat{\beta}_{1}\)-m)/s.e.(\(\hat{\beta}_{1}\))：
T=(100-73.1)/23.76
cat("Test statistic:", T, "\n")Test statistic: 1.132
- 因為我們要驗證的單尾檢定是：\(\emph{H}_{A}:\beta_1> 100\)，所以計算p value時，\(\texttt{lower.tail=F}\)，也就是從分佈的「左邊」計算\(T\)值：
pvalue=pt(T, df = 8, lower.tail = T)
cat("p-value:", pvalue)p-value: 0.8548
- 結果顯示\(p>0.05\)，所以無法拒斥「每增加一間房間售價就可以增加100,000元」的虛無假設。
6.2 氣溫直減率檢定
- 海拔越高的地方，氣溫越低，兩者之間的反比關係稱為lapse rate（氣溫直減率），理論上lapse rate等於\(9.8^{\circ}\)C/Km，也就是每上升每一千公尺，空氣下降攝氏9.8度，或者轉換成\(5.34^{\circ}\)F/1000英呎。有一筆資料有高度與溫度兩個變數，請檢定\(5.34^{\circ}\)F/1000英呎這個理論。
TAB <- data.frame(elevation=c(600,1000,1250,1600,1800,2100,2500,2900),
temp=c(56,54,56,50,47,49,47,45))
knitr::kable(TAB, format='simple', caption='氣溫與高度資料')| elevation | temp |
|---|---|
| 600 | 56 |
| 1000 | 54 |
| 1250 | 56 |
| 1600 | 50 |
| 1800 | 47 |
| 2100 | 49 |
| 2500 | 47 |
| 2900 | 45 |
- 估計OLS迴歸模型如表6.4：
m1<-lm(temp ~ elevation, data=TAB)
stargazer::stargazer(m1, title='氣溫直減率檢定',
type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| temp | |
| elevation | -0.005*** |
| (0.001) | |
| Constant | 59.290*** |
| (1.717) | |
| Observations | 8 |
| R2 | 0.837 |
| Adjusted R2 | 0.810 |
| Residual Std. Error | 1.879 (df = 6) |
| F Statistic | 30.810*** (df = 1; 6) |
| Note: | p<0.1; p<0.05; p<0.01 |
- 表6.4顯示，迴歸係數為-0.00515，標準誤為0.000921。進行雙尾的\(t\)檢定，確認是否拒斥\(\hat{\beta}_{1}=-0.00534\)的虛無假設。因為T>0，我們要計算的\(p\)值應該是右尾累積的面積，所以設定\(\tt{lower.tail=F}\)，也就是\(P(T\leq a)\)。同時，因為是雙尾檢定，所以\(p\)要乘以2：
T=(-0.005115-(-0.00534))/0.000921
T[1] 0.2443
pvalue=pt(T, df =6, lower.tail = FALSE)
pvalue[1] 0.4076
pvalue*2[1] 0.8151
- 因為\(p>0.05\)，也就是說我們觀察到T值的機會高於5%，所以結論是我們無法拒斥這個虛無假設。
- 檢定的指令可參考Hypothesis Testing。
6.3 吳郭魚的產值與產量
現代統計學第二版 的習題14.10有表6.5的資料，請問產量（公噸）與產值（萬元）的關係為何？
dt <- data.frame(quantity=c(2537,2,26761,2500,111,1900,439,4,4485,
22778,1354,69,1492,4,2),
value=c(15332,33,122517,12306,569,12350,2525,56,24713,
117794,8578,451,11932,59,23))
knitr::kable(dt, format='simple', caption='吳郭魚的產值與產量資料')| quantity | value |
|---|---|
| 2537 | 15332 |
| 2 | 33 |
| 26761 | 122517 |
| 2500 | 12306 |
| 111 | 569 |
| 1900 | 12350 |
| 439 | 2525 |
| 4 | 56 |
| 4485 | 24713 |
| 22778 | 117794 |
| 1354 | 8578 |
| 69 | 451 |
| 1492 | 11932 |
| 4 | 59 |
| 2 | 23 |
- 估計OLS迴歸模型如表6.6：
m1<-lm(value ~ quantity, data = dt)
stargazer::stargazer(m1, title='吳郭魚的產值與產量模型',
type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| value | |
| quantity | 4.788*** |
| (0.103) | |
| Constant | 1,383.000 |
| (950.800) | |
| Observations | 15 |
| R2 | 0.994 |
| Adjusted R2 | 0.994 |
| Residual Std. Error | 3,259.000 (df = 13) |
| F Statistic | 2,156.000*** (df = 1; 13) |
| Note: | p<0.1; p<0.05; p<0.01 |
- 表6.6顯示，迴歸係數為4.7875，標準誤為0.1031。
- 去年吳郭魚每增加1公噸平均增加5.1萬元的產值，請問今年是否高於去年的幅度？
T=(4.78 - 5.1)/0.1031
T[1] -3.104
n=nrow(dt)
pvalue=pt(T, df = n-2, lower.tail = TRUE)
pvalue[1] 0.004193
- 因為\(p<0.05\)，顯示去年的5.1顯著的高於今年的4.78，或者4.78要高於5.1的情況不太可能發生，所以拒斥今年的產量影響產值的作用高於5.1的虛無假設。
6.4 大學教師性別與薪水
- 請問大學教師的薪水與性別有關嗎？使用
car裡面的Salaries資料：
library(carData)
fit1 <- with(Salaries, lm(salary ~ sex))
stargazer::stargazer(fit1, title='大學教師的薪水迴歸模型',
type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| salary | |
| sexMale | 14,088.000*** |
| (5,065.000) | |
| Constant | 101,002.000*** |
| (4,809.000) | |
| Observations | 397 |
| R2 | 0.019 |
| Adjusted R2 | 0.017 |
| Residual Std. Error | 30,035.000 (df = 395) |
| F Statistic | 7.738*** (df = 1; 395) |
| Note: | p<0.1; p<0.05; p<0.01 |
- 表6.7顯示：
- 當教師為男性，平均薪水為101003+14088=115090，當教師為女性，平均薪水為101002。
- 性別的迴歸係數\(t\)值為2.782，因為 \[ \frac{14088}{5065}=2.78 \]
- 當\(t\)=2.78，我們想知道對應到的機率值有多高？如果小於虛無假設的機率值，也就是我們觀察到虛無假設成立的機會非常小，也就是我們可以拒斥虛無假設。
h0=1-pt(1.96, 395)
h1=1-pt(2.78, 395)
h0; h1
## [1] 0.02535
## [1] 0.002848- 因此，我們可以拒斥兩者無關的虛無假設，並且得到性別會決定教師薪水的結論。圖6.1顯示性別與薪水的關係：
ggplot2::ggplot(data=Salaries,
ggplot2::aes(x=sex, y=salary))+
geom_point(col='#bb3311') +
stat_summary(geom = "line", fun = mean, group = 1) +
theme_bw()
Figure 6.1: 大學教師薪水與性別分佈圖
6.5 氣溫與壓力
- 在這個例題中，我們觀察殘差值與自變數的關係。首先我們讀取資料，然後估計OLS模型，並且畫圖6.2：
library(readxl)
file <- here::here('Data', 'pressure.xlsx')
pressure <- read_excel(file) #Upload the data
lmTemp = lm(Pressure~Temperature, data = pressure) #Create the linear regression
plot(pressure, pch = 16, col = "blue") #Plot the results
abline(lmTemp) #Add a regression line
Figure 6.2: 氣溫與壓力
- 然後觀察殘差與自變數的關係：
plot(lmTemp$residuals, pch = 16, col = "#3333EE")
- 加上一個平方的變數可以把殘差值轉變為隨機分佈：
lmTemp2 = lm(Pressure ~ Temperature + I(Temperature^2), data = pressure) #Create a linear regression with a quadratic coefficient
summary(lmTemp2) #Review the results
##
## Call:
## lm(formula = Pressure ~ Temperature + I(Temperature^2), data = pressure)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.605 -1.633 0.555 1.180 4.827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.75000 3.61559 9.33 3.4e-05 ***
## Temperature -1.73159 0.15100 -11.47 8.6e-06 ***
## I(Temperature^2) 0.05239 0.00134 39.16 1.8e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.07 on 7 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 0.999
## F-statistic: 7.86e+03 on 2 and 7 DF, p-value: 1.86e-12
plot(lmTemp2$residuals, pch = 16, col = "#333EEE")
- 這是因為溫度越高，壓力越大，但是到一定程度之後，兩者的線性關係 會變成非線性關係，所以要加進平方項，才能去掉殘差項裡面的作用。
7 檢定力(Power of test)
- \(\hat{\beta}\)，也就是\(X\)與\(Y\)的相關，在實驗設計則是實驗分組的差異
- \(\sigma\)，也就是\(Y\)的變異程度(variability)。變異程度越大，\(t\)值越小。
- \(1/\sqrt{n}\)，樣本數越大，\(t\)值也越大。這是因為樣本數越大，越能偵測到\(X\)的效果。
\(\blacksquare\)有關Type I跟Type II error的說明可見Jim Frost (https://statisticsbyjim.com/hypothesis-testing/types-errors-hypothesis-testing/)。
8 殘差的變異數異質性問題
8.1 什麼是殘差的變異數異質性問題？
set.seed(02138)
X<-rnorm(100, 0, 1); Y<-X+rnorm(100, 0, 1.1)
M1 <- lm(Y~X)
plot(fitted(M1), resid(M1), col='#22CDC0')
abline(0,0, col='#AA1100')Figure 8.1: 殘差值與預測值散佈圖2
已知\(\hat{\beta}_{1}=\frac{\sum (x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum (x_{i}-\bar{x})^2}=\frac{\sum (x_{i}-\bar{x})y_{i}}{SST_{x}}\)。
根據(5.3)，迴歸係數\(\hat{\beta}_{1}\)可從上面的程式改寫為：
\[\hat{\beta}_{1}=\beta_{1}+\frac{\sum (x_{i}-\bar{x})u_{i}}{\sum (x_{i}-\bar{x})^2}\]
假設： \[w_{i}\equiv \frac{\sum (x_{i}-\bar{x})}{\sum (x_{i}-\bar{x})^2}\]
根據式(5.4)， \[\begin{align*} V[\hat{\beta}_{1}|X]&=V[\beta_{1}|X]+V[\frac{\sum (x_{i}-\bar{x})u_{i}}{SST_{x}}|X]\\ &=0 + \sum \{\frac{(x_{i}-\bar{x})u_{i}}{SST_{x}}\}^2\hspace{.1em}\rm{Var}[u_{i}|X]\\ &=\frac{\sum (x-\bar{x})^2}{SST_{x}^2}\sigma_{u}^2\\ &=\frac{SST_{x}}{SST_{x}^2}\sigma_{u}^2=\frac{\sigma_{u}^2}{SST_{x}} \tag{5.4} \end{align*}\]
在推演過程中，我們假設：
\[\frac{\sum (x-\bar{x})^2}{SST_{x}^2}\sigma_{u}^2\]
- 裡面的\(\sigma_{u}^2\)對每一個X都相同。如果這個假設不成立，上面要改寫成：
\[\frac{\sum (x-\bar{x})^2}{SST_{x}^2}\sigma_{u}_{i}^2\]
- 所以得出：
\[\text{Var}(\hat{\beta_{1}})=\frac{\sigma_{u}_{i}^2}{SST_{x}}\]
- 變異數異質性(heteroskedasticity)會使得迴歸係數的標準誤不精確，也就是說迴歸係數的離散程度並不是最小，所以假設檢定的結果並不正確。
8.2 檢驗方式
- 首先用 Breusch-Pagan test確認有沒有自變數造成標準誤失真，假設 殘差是自變數的函數：
\[u_{i}^2=\alpha_{0}+\alpha_{1}z_{i,1}+\ldots +\alpha_{k}z_{i,k}+\nu_{i}\]
- 要檢驗的假設是：
\[\textit{H}_{0}:\alpha_{2}=\alpha_{3}=\ldots = \alpha_{k}=0\]
- 可以用\(\texttt{lmtest}\)套件檢驗如下：
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:data.table':
##
## yearmon, yearqtr
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Example data
set.seed(02139)
x1 <- rnorm(200)
x2 <- rnorm(200)
y <- 2*x1 + 3*x2 + rnorm(200)
# Fit linear regression model
model <- lm(y ~ x1 + x2)
# Perform Breusch-Pagan test
bptest(model)
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 8.3, df = 2, p-value = 0.02如果p<0.05，代表違反變異數同質性的假設。
也可用圖8.2檢測。
showtext::showtext_auto()
set.seed(02139)
x1 <- rnorm(200)
x2 <- rnorm(200)
y <- 2*x1 + 3*x2 + rnorm(200)
# Fit linear regression model
model <- lm(y ~ x1 + x2)
# visualization
plot(fitted(model), resid(model), col='#3322cc',
main = '檢測殘差變異數同質性假設散佈圖')
abline(0,0, col='#bb1100')
Figure 8.2: 殘差值與預測值散佈圖3
- 第二種檢測是the White test，做法是把自變數變成平方項，然後對殘差的平方項進行迴歸，再對該模型的殘差進行卡方檢定:
\[u_{i}^2=\alpha_{0}+\alpha_{1}z_{i,1}+\alpha_{2}z_{i,1}^2+\ldots +\alpha_{k}z_{i,k}+\alpha_{k+1}z_{i,k+1}^2+\nu_{i}\]
- 要檢驗的虛無假設是：
\[\textit{H}_{0}:\alpha_{2}=\alpha_{3}=\ldots = \alpha_{k}=0\]
data <- data.frame(y=y, x1=x1, x2=x2)
model2 <- lm(y ~ x1 + x2, data)
u2 <- model2$residuals^2
Ru2<- summary(lm(u2 ~ x1 + x2 + x1^2 + x2^2))$r.squared
LM <- nrow(data)*Ru2
p.value <- 1-pchisq(LM, 4)
cat("White Test's p value:", p.value, '\n')
## White Test's p value: 0.08185因為p>0.05，所以沒有違反殘差變異數同質性假設。
圖2.5似乎顯示殘差在不同Y固定X的情況下，並不是隨機分佈。我們用Breusch-Pagen test以及White test檢測：
mydata <- pscl::bioChemists
m1 <- lm(ment ~ phd + kid5, data = mydata)
u2 <- m1$residuals^2
Ru2<- summary(lm(u2 ~ mydata$phd + mydata$phd^2 +
mydata$kid5 + mydata$kid5^2))$r.squared
LM <- nrow(mydata)*Ru2
p.value <- 1-pchisq(LM, 4)
cat("White Test's p value:", p.value, '\n')
## White Test's p value: 0.03154
#alternative way
fit_Y <- m1$fitted.values
r2.u <- summary(lm(u2 ~ fit_Y + fit_Y^2))$r.squared
LM <- nrow(mydata)*r2.u
p.value <- 1-pchisq(LM, 2)
cat("White Test's p value:", p.value, '\n')
## White Test's p value: 0.03086
# skedastic package
library(skedastic)
white(m1)
## # A tibble: 1 × 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 21.3 0.000278 4 White's Test greater
#Breusch-Pagen test
bptest(m1)
##
## studentized Breusch-Pagan test
##
## data: m1
## BP = 11, df = 2, p-value = 0.005結果發現p<0.05，所以違反殘差變異數同質性假設。
還有一個Goldfeld-Quandt test，適用在資料可以分成兩組，檢驗兩組的變異數是否相等：
\[F=\frac{\hat{\sigma_{A}}^2/\sigma_{B}^{2}}{\hat{\sigma_{B}}^{2}/\sigma_{B}^{2}}\sim F_{N_{A}-k,N_{B}-k}\]
8.3 殘差異質性的標準誤
- 遇到殘差異質性時，我們得到的標準誤失真，heteroskedasticity-consistent standard errors (or robust errors)表示如下：
\[\text{Var}[\hat{\beta_{1}}]=\frac{n}{n-k}(\mathbb{X^{\prime}X})^{-1}\hat{\Sigma}(\mathbb{X^{\prime}X})^{-1}\]
- 其中
\[\hat{\Sigma}=\sum_{i=1}^{n}\mathbb{x_{i}}{x_{i}^{\prime}}u_{i}^{1}\]
- 用\(\texttt{sandwich}\)套件裡面的\(\texttt{coeftest}\)函數估計robust error。
8.3.1 實例一
- 用Germán Rodríguez所提供的20個拉丁美洲國家的出生率與家庭計劃、社會環境的資料為例計算穩健的標準誤。
library(haven)
dt <- read_dta("https://grodri.github.io/datasets/effort.dta")
dt$effortg = cut(dt$effort, breaks=c(min(dt$effort),5,15, max(dt$effort)), right=FALSE, include.lowest=TRUE, labels=c("Weak", "Moderate", "Strong"))
#OLS model
OLS1 <- lm(change ~ setting + effortg, data = dt)
library(lmtest)
bptest(OLS1)studentized Breusch-Pagan testdata: OLS1 BP = 3.6, df = 3, p-value = 0.3
library(sandwich)
robustmodel <- coeftest(OLS1, vcov = vcovHC(OLS1, type="HC1"))
stargazer::stargazer(OLS1, robustmodel, type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | ||
| change | ||
| OLS | coefficient | |
| test | ||
| (1) | (2) | |
| setting | 0.169 | 0.169*** |
| (0.106) | (0.045) | |
| effortgModerate | 4.144 | 4.144 |
| (3.191) | (3.191) | |
| effortgStrong | 19.450*** | 19.450*** |
| (3.729) | (2.567) | |
| Constant | -5.954 | -5.954** |
| (7.166) | (2.708) | |
| Observations | 20 | |
| R2 | 0.802 | |
| Adjusted R2 | 0.764 | |
| Residual Std. Error | 5.732 (df = 16) | |
| F Statistic | 21.550*** (df = 3; 16) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
- 標準誤變小，但是迴歸係數不變。
8.3.2 實例二
- 同樣地，用1985 Current Population Survey (CPS) 的資料為例計算穩健的標準誤，由性別、教育程度、工會會員等等變數預測薪資。
showtext::showtext_auto()
wages <- read_dta("https://grodri.github.io/datasets/wages.dta")
m1 <- lm(wages ~ education + workexp + unionmember + south + occupation + female,
data = wages)
# visualization
plot(fitted(m1), resid(m1), col='#3322cc',
main = '檢測殘差變異數同質性假設散佈圖')
abline(0,0, col='#bb1100')
Figure 8.3: 檢測殘差變異數同質性假設散佈圖
#OLS model
lm1 <- lm(wages ~ education + workexp + unionmember + south + occupation + female,
data = wages)
library(lmtest)
bptest(lm1)studentized Breusch-Pagan testdata: lm1 BP = 12, df = 6, p-value = 0.07
library(sandwich)
robustmodel <- coeftest(lm1, vcov = vcovHC(lm1, type="HC1"))
stargazer::stargazer(lm1, robustmodel, type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | ||
| wages | ||
| OLS | coefficient | |
| test | ||
| (1) | (2) | |
| education | 0.904*** | 0.904*** |
| (0.081) | (0.094) | |
| workexp | 0.104*** | 0.104*** |
| (0.017) | (0.019) | |
| unionmember | 1.444*** | 1.444*** |
| (0.523) | (0.516) | |
| south | -0.787* | -0.787* |
| (0.427) | (0.414) | |
| occupation | -0.062 | -0.062 |
| (0.125) | (0.151) | |
| female | -2.207*** | -2.207*** |
| (0.398) | (0.400) | |
| Constant | -3.362** | -3.362* |
| (1.466) | (1.750) | |
| Observations | 534 | |
| R2 | 0.270 | |
| Adjusted R2 | 0.261 | |
| Residual Std. Error | 4.416 (df = 527) | |
| F Statistic | 32.450*** (df = 6; 527) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
- 大部分的迴歸係數標準誤變小，但是迴歸係數不變。
9 作業
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(format='simple', caption='習題一資料') %>%
kable_styling(bootstrap_options = 'striped')
## Warning in kable_styling(., bootstrap_options = "striped"): Please specify
## format in kable. kableExtra can customize either HTML or LaTeX outputs. See
## https://haozhu233.github.io/kableExtra/ for details.| 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 |
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 |
\[\rm{y2000} = \beta_{0}+\beta_{1}\rm{y1970}\]
carData::Migration這筆資料，檢驗migrants跟distance之間的迴歸模型，有沒有符合殘差值與預測值無關的假設？
car::spreadLevelPlot函數確認預測值沒有特殊的分佈。
carData裡面的Salaries資料分析，回答大學教師的薪水與服務年資有關嗎？跟性別比起來，哪一個變數與薪水相關比較大？
10 更新日期
## 最後更新日期 05/17/2024