社會科學統計方法

複迴歸

1 課程目標

在介紹完單迴歸模型後，接下來將介紹複迴歸模型，以及變數的轉換。

2 複迴歸模型優點：

可以統計更多資訊，做為描述之用

改善模型的預測能力，解釋單一變數無法解釋的變異量

控制可能的「混淆」(confounding)變數，以利因果推論

估計更複雜的模型(例如：\(Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}^2\))

估計有交互作用變項的模型(例如：\(Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{1}X_{2}\))

單迴歸的模型寫成: \(E[Y|X_{1}=x_{1}]=\beta_{0}+\beta_{1}x_{1}\)

雙變數的迴歸模型可寫成: \(E[Y|X_{1}=x_{1}, X_{2}=x_{2}]=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}\)

3 複迴歸的假設

複迴歸與單迴歸模型相同，都做以下假設（參考現代統計學二版第625-626頁）：

1. 依變數為隨機變數，自變數\(X_{k}\)，\(k=1,\ldots,k\)。
2. 誤差項u的條件平均值為0。
3. 誤差項u的條件變異數為\(\sigma^2\)。
4. 誤差項服從常態分配。
5. 誤差項之間互相獨立。
6. 誤差項與自變數之間互相獨立。
7. 自變數之間沒有完美的線性相關。
8. 樣本數n>k+1

4 複迴歸係數的意義

根據上述假設，有X,Z兩個自變數的複迴歸模型可寫成：

\[E[Y|X,Z]=\beta_{0}+\beta_{1}X+\beta_{2}Z\]

• 估計的迴歸係數\(\hat{\beta}_{1}\)可計算為：

\[\hat{\beta}_{1}=\frac{\sum (z-\bar{z})^2\sum (x-\bar{x})(y-\bar{y})-\sum (x-\bar{x})(z-\bar{z})\sum (z-\bar{z})(y-\bar{y})}{\sum (x-\bar{x})^2\sum(z-\bar{z})^2-(\sum (x-\bar{x})(z-\bar{z}))^2}\]

如果要說明第一個變數的迴歸係數是\(X\)除去與\(Z\)相關部份之後的作用，\(\hat{\beta}_{1}\)改寫為：

\[\frac{\sum_{i}^{n}\hat{r}_{xz,i}y_{i}}{\sum_{i}^{n}\hat{r}_{xz,i}^2}\]

• 其中\(\hat{r}_{xz,i}\)代表用\(Z\)預測\(X\)的殘差：

\[X=\lambda+\delta\dot Z+r_{xz}\]

• 其中，\(\delta\)代表\(X,Z\)的相關程度
• \(\hat{r}_{xz}\)代表\(X\)與\(Z\)沒有相關的部份。
• \(\hat{r}_{xz}\)是\(X\)被\(Z\)解釋過後的部份，也就是淨作用。
• 所以\(\hat{\beta}_{1}\)代表\(X\)除去與\(Z\)相關部份之後的作用。

4.1 複迴歸的樣本統計的特性

與簡單OLS迴歸相同，複迴歸的樣本統計有以下的特性：

• 殘差項的平均值等於0，也就是\(\bar{\hat{u}}=0\)

• proof: \[\begin{align*} \bar{\hat{u}}&=\bar{y}-\hat{\beta}_{0}-\hat{\beta}_{1}\bar{x}-\hat{\beta}_{2}\bar{z} \\ & = \bar{y}-(\bar{y}-\hat{\beta}_{1}\bar{x}-\hat{\beta}_{2}\bar{z})- \hat{\beta}_{1}\bar{x}-\hat{\beta}_{2}\bar{z}=0 \end{align*}\]

• 這個特性可導出\(\bar{y} = \bar{\hat{y}}\)

• proof: \[\begin{align*} y_{i} & =\hat{y}_{i}+\hat{u}_{i}\Rightarrow \bar{y}\\ &=n^{-1}\sum_{i=1}^{n}y_{i}=n^{-1}\sum_{i=1}^{n}(\hat{y}_{i}+\hat{u}_{i})\\ &=\bar{\hat{y}}-\bar{\hat{u}}=\bar{\hat{y}} \end{align*}\]

• \(Cov(x,\hat{u})=Cov(z,\hat{u})=0\)

• 可導出\(Cov(\hat{y},\hat{u})=0\)

• (\(\bar{x},\bar{y},\bar{z}\))一定會落在迴歸平面上。

4.2 複迴歸係數與矩陣

我們可以用矩陣表示迴歸模型為（參考Bremer以及史丹佛大學講義）：

\[\bf{y}=\bf{X}\beta+\bf{\epsilon}\]

• 其中

\[ \bf{y}= \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{bmatrix} \]

\[ X=\begin{bmatrix} 1 & X_{11} & X_{12} &\ldots & X_{1k} \\ 1 & X_{21} & X_{22} &\ldots & X_{2k} \\ \vdots & \vdots & \vdots &\ldots &\vdots \\ 1 & X_{n1} & X_{n2} &\ldots & X_{nk}\\ \end{bmatrix} \]

• 迴歸係數則是表示為：

\[ \bf{\beta}= (\beta_{0},\beta_{1},\beta_{2}\ldots,\beta_{k})=\begin{bmatrix} \beta_{0} \\ \beta_{1} \\ \beta_{2} \\ \vdots \\ \beta_{k} \end{bmatrix} \]

\[ \bf{\epsilon}= \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{bmatrix} \]

• 假設所有的殘差項等於0，模型可表示為：

\[\begin{align*} \bf{\hat{y}}&=\bf{X}\hat{\beta} \tag{4.1} \end{align*}\]

• 因為理論上自變數與殘差應該等於0，在矩陣上表示為\(\bf{X'\epsilon}=0\)，因此式(4.1)改為：

\[\begin{align*} \bf{X'}(\bf{y}-\bf{X}\hat{\beta}) & =0\\ \Longleftrightarrow \bf{X'y}-\bf{X'X}\hat{\beta}&=0\\ \Longleftrightarrow \bf{X'X}\hat{\beta}&=\bf{X'y} \tag{4.2} \end{align*}\]

• 因此，迴歸係數的估計式為(4.3)：

\[\begin{align*} \hat{\beta}&=(\bf{X'X})^{-1}\bf{X'y} \tag{4.3} \end{align*}\]

• 式(4.3)也可以從normal equation推導而來。normal equation 的矩陣為：

\[(\bf{X'X})\hat{\beta}=\bf{X'}y\]

• 因為

\[\begin{align*} \bf{y}&=\bf{X}\beta+\bf{\epsilon}\\ &=\bf{X}\hat{\beta}+u \tag{4.4} \end{align*}\]

• 在式(4.4)中，u為殘差。所以：

\[\begin{align*} (\bf{X'X})\hat{\beta}&=\bf{X'}(\bf{X}\beta+\bf{\epsilon})\\ &= (\bf{X'X})\hat{\beta}+\bf{X'}u \end{align*}\]

• 可以推導出:

\[\bf{X'}u=0\]

• 換句話說，自變數與殘差沒有相關。還可以推導出預測值\(\hat{y}\)與殘差沒有相關。

• \[\hat{y}u=0\]

• 接下來是推導\(\rm{Var}(\hat{\beta})\)。因為：

\[\rm{Var}(\hat{\beta})=E[\hat{\beta}^2]-E[\hat{\beta}]E[\hat{\beta}']\]

• 而且：

\[\hat{\beta}=(\bf{X'X})^{-1}\bf{X'y}\]

• 所以

\[\begin{align*} Var(\hat{\beta})&=[(\bf{X'X})^{-1}\bf{X'y})^2]-\beta^2 \end{align*}\]

• 經過整理後， \(Var(\hat{\beta})= \sigma^2\bf{(X'X)^{-1}}\)（詳細過程）

• 其中\(\sigma^2\)可用\(\hat{\sigma}^2\)來估計，也就是殘差項的平方和除以n-k：

\[\hat{\sigma}^2=\frac{\sum u^2}{n-k}\]

• 如果只有一個自變數，k=2，因為要加上一個常數項，也就是單迴歸模型講到的\(\sigma^2=\frac{\sum (y-\hat{y})^2}{n-2}\)。

我們以\(\tt{UsingR::SAT}\)這筆資料為例，用各州的平均教育花費以及生師比預測SAT的分數，估計OLS迴歸模型如表4.1：

m1<-lm(total ~ expend+ratio, data=UsingR::SAT)
stargazer::stargazer(m1,  title='SAT分數迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 4.1: SAT分數迴歸模型

	Dependent variable:
	total
expend	-22.310***
	(7.956)
ratio	-2.295
	(4.784)
Constant	1,136.000***
	(107.800)
Observations	50
R2	0.149
Adjusted R2	0.113
Residual Std. Error	70.480 (df = 47)
F Statistic	4.114** (df = 2; 47)
Note:	p<0.1; p<0.05; p<0.01

如果我們用矩陣的方式來處理，先建立自變數的矩陣：

dat<-UsingR::SAT
X<-as.matrix(cbind(1,dat$expend,dat$ratio))
head(X)
##      [,1]  [,2] [,3]
## [1,]    1 4.405 17.2
## [2,]    1 8.963 17.6
## [3,]    1 4.778 19.3
## [4,]    1 4.459 17.1
## [5,]    1 4.992 24.0
## [6,]    1 5.443 18.4

• 然後建立Y的矩陣：

Y<-as.matrix(dat$total)
head(Y)
##      [,1]
## [1,] 1029
## [2,]  934
## [3,]  944
## [4,] 1005
## [5,]  902
## [6,]  980

• 最後根據估計式(4.3)求出迴歸係數：

beta_hat<-solve(t(X)%*%X)%*%t(X)%*%Y
beta_hat
##          [,1]
## [1,] 1136.336
## [2,]  -22.308
## [3,]   -2.295

• 對照表4.1，可以看出係數都相同。

4.2.1 標準誤矩陣

我們應用\(\hat{\sigma}^2=\frac{\sum u^2}{n-k}\)來計算迴歸係數的標準誤。計算時我們可以從\(\tt{anova(model)}\)找出\(\frac{\sum u^2}{n-k}\)，也就是表4.2：

anatable=anova(m1)

knitr::kable(anatable, format='simple', caption='迴歸模型變異數分析表')

Table 4.2: 迴歸模型變異數分析表

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
expend	1	39722	39722	7.9974	0.0069
ratio	1	1143	1143	0.2301	0.6337
Residuals	47	233443	4967	NA	NA

• 或者是直接計算\(\sum(y-\hat{y})^2/n-k\)，也就是

sum((dat$total-fitted(m1))^2)/47
## [1] 4967

• 因為\(Var(\hat{\beta})= \sigma^2\bf{(X'X)^{-1}}\)，所以我們解\(\bf{X'X}\)得到此一變異數與共變數矩陣(variance-covariance)的對角線，乘上\(\hat{\sigma}^2\)如下：

var_betaHat <- anova(m1)[[3]][3]* solve(t(X) %*% X)

sqrt(diag(var_betaHat))
## [1] 107.803   7.956   4.784

• 對照表4.1，可以看出係數的標準誤都相同。

• 迴歸係數的變異數與共變數矩陣如下：

\[ E[(\hat{\beta}-\beta)(\hat{\beta}-\beta)']=\begin{bmatrix} \rm{var}(\hat{\beta}_{1}) & \rm{cov}(\hat{\beta}_{1},\hat{\beta}_{2}) &\ldots & \rm{cov}(\hat{\beta}_{1},\hat{\beta}_{k}) \\ \rm{cov}(\hat{\beta}_{1},\hat{\beta}_{2}) & \rm{var}(\hat{\beta}_{2}) &\ldots & \rm{cov}(\hat{\beta}_{2},\hat{\beta}_{k}) \\ \vdots & \vdots & \vdots &\vdots \\ \rm{cov}(\hat{\beta}_{k},\hat{\beta}_{1}) & \rm{cov}(\hat{\beta}_{k},\hat{\beta}_{2}) &\ldots & \rm{var}(\hat{\beta}_{k})\\ \end{bmatrix} \]

• 由上可知，\(\hat{\beta}\)的標準誤就是變異數-共變數矩陣的對角線元素的開根號。

\[\hat{\beta}\sim N(\beta, \sigma^2\bf{(X'X)^{-1}})\]

5 複迴歸模型

學生的成績與母親學歷以及母親的智商有關嗎？先分別估計單迴歸模型：

library(foreign)
file<-here::here('data','kidiq.dta')
iq<-read.dta(file)
f1 <- with(iq, lm(kid_score ~ mom_hs) )
stargazer::stargazer(f1,  title='母親學歷與學生成績',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.1: 母親學歷與學生成績

	Dependent variable:
	kid_score
mom_hs	11.770***
	(2.322)
Constant	77.550***
	(2.059)
Observations	434
R2	0.056
Adjusted R2	0.054
Residual Std. Error	19.850 (df = 432)
F Statistic	25.690*** (df = 1; 432)
Note:	p<0.1; p<0.05; p<0.01

• 根據表5.1，估計結果寫成：
  \[\rm{學生成績} = 77.548 + 11.77\cdot \rm{母親高中畢業與否}\]

f2 <- with(iq, lm(kid_score ~ mom_iq))
stargazer::stargazer(f2,  title='母親智力與學生成績',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.2: 母親智力與學生成績

	Dependent variable:
	kid_score
mom_iq	0.610***
	(0.059)
Constant	25.800***
	(5.917)
Observations	434
R2	0.201
Adjusted R2	0.199
Residual Std. Error	18.270 (df = 432)
F Statistic	108.600*** (df = 1; 432)
Note:	p<0.1; p<0.05; p<0.01

• 估計結果寫成：
  \[\rm{學生成績} = 25.799+ 0.609\cdot \rm{母親智力成績}\]

接下來估計雙變數模型如下：

\(\widehat{\mathrm {Kid's\hspace{.25em} Scores}}=\hat{\beta}_{0}+\hat{\beta}_{1}(\mathrm {Mom\hspace{.25em} has\hspace{.25em} HS\hspace{.25em} degree})+\hat{\beta}_{2}(\mathrm {Mom's\hspace{.25em} IQ})\)

• 當\(X_{1}=1\)（母親有高中學歷），

\[\begin{align*} \widehat{\mathrm {Kid's\hspace{.25em} Scores}} & = \hat{\beta}_{0}+\hat{\beta}_{1}\cdot 1+\hat{\beta}_{2}(\mathrm {Mom's\hspace{.25em} IQ}) \\ & = \hat{\beta}_{0}+\hat{\beta}_{1}+\hat{\beta}_{2}(\mathrm {Mom's\hspace{.25em} IQ}) \end{align*}\]

• 當\(X_{1}=0\)（母親沒有高中學歷），

\[\begin{align*} \widehat{\mathrm {Kid's\hspace{.25em} Scores}} & = \hat{\beta}_{0}+\hat{\beta}_{1}\cdot 0+\hat{\beta}_{2}(\mathrm {Mom's\hspace{.25em} IQ}) \\ & = \hat{\beta}_{0}+\hat{\beta}_{2}(\mathrm {Mom's\hspace{.25em} IQ}) \end{align*}\]

• 也就是說，這兩個模型的斜率一樣，但是當\(X_{1}=1\)時，截距多了 \(\hat{\beta}_{1}\)。

• 進行估計如下：

f3 <- with(iq, lm(kid_score ~ mom_hs + mom_iq))
stargazer::stargazer(f3,  title='母親智力及學歷與學生成績',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.3: 母親智力及學歷與學生成績

	Dependent variable:
	kid_score
mom_hs	5.950***
	(2.212)
mom_iq	0.564***
	(0.061)
Constant	25.730***
	(5.875)
Observations	434
R2	0.214
Adjusted R2	0.210
Residual Std. Error	18.140 (df = 431)
F Statistic	58.720*** (df = 2; 431)
Note:	p<0.1; p<0.05; p<0.01

• 估計得到的模型為： \[\hat{Y}=25.73+5.95\cdot \rm{母親高中畢業與否} +0.56\cdot \rm{母親智商}\]

• 迴歸係數分別代表的意義為：

\(\hat{\beta}_{0}\)代表母親沒有高中畢業的迴歸線的截距

\(\hat{\beta}_{1}\)代表母親有無高中畢業的兩條迴歸線的垂直距離

\(\hat{\beta}_{2}\)代表兩條線的斜率

• 圖5.1顯示學生的母親有高中畢業比起沒有高中畢業，平均成績高了一些。

ggplot(iq, aes(x=mom_iq, y=kid_score, group=mom_hs)) +
    geom_point(aes(color=factor(mom_hs)), shape=1, size=2)  +
    geom_abline( slope=0.56, intercept=25.73, col="indianred", size=1.2) +
    geom_abline(slope=0.56, intercept=25.73+5.95, col="blue", size=1.2) +
    scale_color_discrete(name = "Mother with High School Degree",
       labels=c("No High School Degree","High School Degree")) +
    theme_bw() +
    theme(legend.position = "top") +
    labs(x="Mother's IQ", y="Kid's Score")

Figure 5.1: 雙變數迴歸模型

以上結果顯示雙變數的迴歸模型的估計方式與解釋方式，自變數其中之一是虛擬變數。

5.1 可以比較不同迴歸模型的係數嗎？

• 比較不同迴歸模型的係數，被稱為group comparison problem。 Clogg, Petkova, and Haritou (1995) 討論如何比較不同線性迴歸模型的係數，見Clogg, Clifford C., Eva Petkova, and Adamantios Haritou. 1995. Statistical Methods for Comparing Regression Coefficients between Models. American Journal of Sociology 100,5:1261-93. 而Paul D. Allison (1995) 提出修正的方式。見Allison, Paul D. 1995. The Impact of Random Predictors on Comparisons of Coefficients Between Models: Comment on Clogg, Petkova, and Haritou. American Journal of Sociology 100,5:1294-1305. Kuha and Mills (2017)帶入因果推論的概念來討論二元迴歸模型的係數比較。見Kuha, Johni and Mills, Colin. (2017) On group comparisons with logistic regression models. Sociological Methods & Research.
• 根據Kuha and Mills (2017)，如果我們正確地測量依變數是連續變數，就可以比較不同模型的係數。但是如果自變數與群體的性質有關，也就是不同的母體有不同的自變數的分佈，那麼係數的比較可能會有問題。

6 兩個連續變數

接下來我們改用母親的智商與年齡這兩個連續變數預測依變數：學生成績，估計模型：

fit4 <- with(iq, lm(kid_score ~ mom_iq + mom_age))
stargazer::stargazer(fit4,  title='母親智力及年齡與學生成績',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 6.1: 母親智力及年齡與學生成績

	Dependent variable:
	kid_score
mom_iq	0.604***
	(0.059)
mom_age	0.388
	(0.326)
Constant	17.600*
	(9.084)
Observations	434
R2	0.204
Adjusted R2	0.200
Residual Std. Error	18.260 (df = 431)
F Statistic	55.080*** (df = 2; 431)
Note:	p<0.1; p<0.05; p<0.01

\[ \hat{Y}=17.6+0.6\cdot \rm{母親智商} +0.39\cdot \rm{母親年齡} \]

雙連續變數的迴歸係數分別代表的意義為：

\(\hat{\beta}_{0}=17.6\)代表母親智商與年齡等於0的截距（沒有意義）

\(\hat{\beta}_{1}=0.6\)代表對於同樣的兩個觀察值，母親年齡相同但是母親智商有一單位的差異預期帶來的學生分數的差異。

\(\hat{\beta}_{2}=0.39\)代表對於同樣的兩個觀察值，母親智商相同但是母親年齡有一單位的差異預期帶來的學生分數的差異。

• 因為對Ｘ進行微分可得\(\beta_{1}\)：

\[ \mathcal{\frac{\partial (y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2})}{\partial X_{1}}}=\beta_{1} \]

• 所以\(\hspace{.3em}\beta_{1}\hspace{.3em}\)代表當其他變數在同一個值或是水準時，對於\(\hspace{.3em}y\hspace{.3em}\)的淨作用。

7 交互作用

• 一個自變數除了與依變數相關，該自變數會透過另一個自變數影響依變數，我們稱之為交互作用。
• 例如有慢性病的病人，如果感染冠狀病毒，反應會特別比沒有慢性病的感染者嚴重。或者有公民責任感的人，如果又支持特定政黨，可能比其他的政黨支持者，更想去投票。

7.1 交互作用原理

假設\(Y\)因為第二個變數\(Z\)而有差異，而且自變數\(X\)也因為\(Z\)有差異，我們需要考慮交互作用：

\[\widehat{\mathrm {Kid's\hspace{.25em} Scores}}=\hat{\beta}_{0}+\hat{\beta}_{1}(\mathrm {Mom\hspace{.25em} has\hspace{.25em} HS\hspace{.25em} degree})+\hat{\beta}_{2}(\mathrm {Mom's\hspace{.25em} IQ}) + \hat{\beta}_{3}(\mathrm{Mom's\hspace{.25em} IQ}\cdot \mathrm{Mom's\hspace{.25em} HS}) \]

重新估計：

fit5 <- with(iq, lm(kid_score ~ mom_iq + mom_hs+mom_iq:mom_hs))
fit5.1<-with(iq, lm(kid_score[mom_hs==1] ~ mom_iq[mom_hs==1]))
fit5.2<-with(iq, lm(kid_score[mom_hs==0] ~ mom_iq[mom_hs==0]))

library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
stargazer(fit5.1, fit5.2, fit5, align=TRUE, title='三個模型比較',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 7.1: 三個模型比較

	Dependent variable:
	kid_score[mom_hs == 1]	kid_score[mom_hs == 0]	kid_score
	(1)	(2)	(3)
mom_iq[mom_hs == 1]	0.485***
	(0.065)
mom_iq[mom_hs == 0]		0.969***
		(0.157)
mom_iq			0.969***
			(0.148)
mom_hs			51.270***
			(15.340)
mom_iq:mom_hs			-0.484***
			(0.162)
Constant	39.790***	-11.480	-11.480
	(6.663)	(14.600)	(13.760)
Observations	341	93	434
R2	0.143	0.294	0.230
Adjusted R2	0.140	0.286	0.225
Residual Std. Error	17.660 (df = 339)	19.070 (df = 91)	17.970 (df = 430)
F Statistic	56.420*** (df = 1; 339)	37.870*** (df = 1; 91)	42.840*** (df = 3; 430)
Note:	p<0.1; p<0.05; p<0.01

請計算完整的模型(3)，當母親高中程度變數分別為1與0，係數是否分別等於模型(1)與(2)？

• 重新畫圖7.1：

ggplot(iq, aes(x=mom_iq, y=kid_score, group=mom_hs)) +
    geom_point(aes(color=factor(mom_hs)), shape=1, size=2)  +
    geom_abline( slope=0.969, intercept=-11.482, col="indianred", size=1.2) +
    geom_abline(slope=0.485, intercept=39.786, col="blue", size=1.2) +
    scale_color_discrete(name = "Mother with High School Degree",
   labels=c("No High School Degree","High School Degree")) +
    theme_bw() +
    theme(legend.position = "top") +
    labs(x="Mother's IQ", y="Kid's Score")

Figure 7.1: 雙變數迴歸模型

• 代入不同的值，會得到不同的交互作用與自變數的值。

7.2 用Anova測試交互作用是否顯著

• Anova的型2測試提供模型的殘差平方和會不會因為變數一個被拿掉而其他留下的時候，發生下降的現象。如果有顯著的不同，代表該變數有作用。我們用學生成績的資料為例：

options(contrasts = c("contr.treatment", "contr.poly"))

fit5 <- with(iq, lm(kid_score ~ mom_iq + mom_hs+mom_iq:mom_hs))

library(car)
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.2
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
Anova(fit5, type="II")
## Anova Table (Type II tests)
## 
## Response: kid_score
##               Sum Sq  Df F value Pr(>F)    
## mom_iq         28504   1   88.26 <2e-16 ***
## mom_hs          2380   1    7.37 0.0069 ** 
## mom_iq:mom_hs   2878   1    8.91 0.0030 ** 
## Residuals     138879 430                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 因為交互作用達到統計上的顯著水準，所以拿掉這個交互作用項會影響模型與資料的符合程度，因此可以說兩個群體之間有不同的斜率。

• 同樣的方法適用在glm模型。

• National Institute of Statistical Sciences 在這篇文章探討如何測試迴歸模型有二元變數時的斜率相等假設，請大家參考。 # 變數的淨作用

• 我們可以從淨作用的角度詮釋迴歸係數。

• 假設迴歸模型為：

\[Y=\beta_{0}+\beta_{1}X+\beta_{2}Z+\epsilon\]

• 估計的迴歸係數可寫成\[\hat{\beta}_{1}=\frac{\sum_{i}^{n}\hat{r}_{xz,i}y_{i}}{\sum_{i}^{n}\hat{r}_{xz,i}}\]

• 其中\(\hat{r}_{xz,i}\)代表用\(Z\)預測\(X\)的殘差：

\[X=\lambda+\delta\dot Z+r_{xz}\]

• 也就是說：

\(\delta\)代表\(X,Z\)的相關程度

\(\hat{r}_{xz}\)代表\(X\)與\(Z\)沒有相關的部份。

\(\hat{r}_{xz}\)是\(X\)被\(Z\)解釋過後的部份，也就是

所以\(\hat{\beta}_{1}\)代表\(X\)除去與\(Z\)相關部份之後的作用。

7.3 範例：工黨得票率的迴歸模型

以pscl套件的UKHouseOfCommons為例，v2是工黨的得票率，v3是自由民主黨的得票率，labinc代表該選區現任者為工黨議員。工黨得票率受到自由民主黨的得票率影響，因為兩個黨的意識形態很接近。此外，現任者如果是工黨，這次選舉的得票率應該比較好。 先估計完整的複迴歸模型如表7.2：

library(pscl); library(stargazer)
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
labormodel<-lm(v2 ~ v3 + labinc , data=UKHouseOfCommons)
stargazer::stargazer(labormodel,  title='工黨得票率迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 7.2: 工黨得票率迴歸模型

	Dependent variable:
	v2
v3	-0.998***
	(0.041)
labinc	0.206***
	(0.009)
Constant	0.495***
	(0.010)
Observations	521
R2	0.770
Adjusted R2	0.769
Residual Std. Error	0.085 (df = 518)
F Statistic	865.700*** (df = 2; 518)
Note:	p<0.1; p<0.05; p<0.01

• 表7.2顯示估計得到的模型可寫成：

\[\rm{v2}=0.495-0.998\cdot \rm{v3}+0.206\cdot \rm{labinc}\]

估計第二個自變數預測第一個自變數的模型，把殘差值存起來，預測原來的依變數。

library(pscl)
m2 <- lm(v3 ~  labinc , data=UKHouseOfCommons)
resd.m2 <- residuals(m2)
m3 <- lm(v2 ~ resd.m2, data=UKHouseOfCommons)

stargazer::stargazer(m2, m3, align=T,  title='用殘差值建立的工黨得票率迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 7.3: 用殘差值建立的工黨得票率迴歸模型

	Dependent variable:
	v3	v2
	(1)	(2)
labinc	-0.080***
	(0.009)
resd.m2		-0.998***
		(0.073)
Constant	0.211***	0.361***
	(0.005)	(0.007)
Observations	521	521
R2	0.131	0.263
Adjusted R2	0.129	0.261
Residual Std. Error (df = 519)	0.091	0.153
F Statistic (df = 1; 519)	78.250***	185.000***
Note:	p<0.1; p<0.05; p<0.01

• 表7.3顯示：

\[v3=0.361-0.998\cdot \mathrm{residual\hspace{.1cm}of\hspace{.1cm}model\hspace{.1cm}2}\]

得到的係數等於複迴歸模型的第一個迴歸係數。這說明第一個自變數的係數代表去掉第二個自變數作用的淨作用。

library(pscl)
m2 <- lm(v3 ~  labinc , data=UKHouseOfCommons)
resd.m2 <- UKHouseOfCommons$v3-fitted(m2)

beta1<-sum(resd.m2*UKHouseOfCommons$v2)/sum((resd.m2)^2)
beta1

[1] -0.9978

8 \(R^2\)與F檢定

8.1 \(R^2\)的意義

\(R^2\)代表模型解釋的總變異量佔依變數的總變異的比例，所以有人稱為goodness-of-fit，模型適合度，也有人稱為判定係數(coefficient of determination)，表示如下：

\[\begin{align*} R^2 & =\frac{\sum\limits_{i=1}^{n}(\hat{Y}-\bar{Y})^2}{\sum\limits_{i=1}^{n}(Y-\bar{Y})^2}\\ &=1-\frac{\sum\limits_{i=1}^{n}(Y-\hat{Y})^2}{\sum\limits_{i=1}^{n}(Y-\bar{Y})^2}\\ &= 1-\frac{\sum\limits_{i=1}^{n}e^2}{\sum\limits_{i=1}^{n}(Y-\bar{Y})^2}\\ &=\frac{SSR}{SST} \tag{8.1} \end{align*}\]

如果有太多自變數，可能高估\(R^2\)，因此有人建議adjusted \(R^2\)：

\[\rm{adj}\hspace{.2em}R^2=1-\frac{\sum\limits_{i=1}^{n}(Y-\hat{Y})^2}{\sum\limits_{i=1}^{n}(Y-\bar{Y})^2}\times\frac{n-1}{n-k-1}\]

簡單來說，就是1減掉殘差值平方和除以y的變異數，然後乘以n-1，最後除以n-k-1

8.2 預測值

• 如果我們的模型是\(Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}\)，我們有\(\hat{\beta}_{0}\), \(\hat{\beta}_{1}\), \(\hat{\beta}_{2}\)三個係數，而預測值\(\hat{Y}=\bf{X'}\beta\)，其中自變數表示為：

\[ X'=\begin{bmatrix} 1 & X_{11} & X_{21} \\ 1 & X_{12} & X_{22} \\ \vdots & \vdots &\vdots \\ 1 & X_{1n} & X_{2n} \\ \end{bmatrix} \]

• 迴歸係數則是表示為：

\[ \bf{\beta}= (\hat{\beta}_{0},\hat{\beta}_{1},\hat{\beta}_{2})=\begin{bmatrix} \hat{\beta}_{0} \\ \hat{\beta}_{1} \\ \hat{\beta}_{2} \end{bmatrix} \]

• 因為\(\bf{X'}\)是\(n\times 3\)矩陣，\(\bf{\beta}\)是\(3\times 1\)矩陣，所以相乘得到\(n\times 1\)矩陣。\(\hat{Y}=X'\beta=\)

\[ \begin{bmatrix} \hat{\beta}_{0}\times 1 +X_{11}\hat{\beta}_{1} + X_{21}\hat{\beta}_{2} \\ \hat{\beta}_{0}\times 1 + X_{12}\hat{\beta}_{1} + X_{22}\hat{\beta}_{2} \\ \vdots\\ \hat{\beta}_{0}\times 1 +X_{1n}\hat{\beta}_{1} + X_{2n}\hat{\beta}_{2} \end{bmatrix} \]

• 產生預測值\(\hat{Y}\)的方法：

set.seed(02139)
#data
X1 <- rnorm(10, 0, 1)+runif(10, 0, 1)
X2 <- rnorm(10, 0, 1)+runif(10, 0, 0.5)
Y <- rnorm(10, 0, 1)
#model
M0 <-lm(Y ~ X1 + X2)
#predicted value
as.vector(t(M0$coefficients[c(1, 2,3)])%*%t(cbind(1,X1,X2)))

[1] -0.4156 -0.4394 -0.5445 -0.8770 -0.1836 -0.6284 -0.2466 -0.5551 -1.0363 [10] -0.4813

• R可以直接產生預測值，也就是\(\tt{predict(model)}\)：

as.vector(predict(M0))

[1] -0.4156 -0.4394 -0.5445 -0.8770 -0.1836 -0.6284 -0.2466 -0.5551 -1.0363 [10] -0.4813

有了預測值，我們可以計算\(R^2\)。根據(8.1)以及表7.2，計算如下，或者用自訂函數計算：

#calculate R2 step-by-step
y = pscl::UKHouseOfCommons$v2
mean.y<-mean(y)
sum.res<-sum((y-labormodel$fitted.values)^2)
var.y <- sum((y - mean.y)^2)
r2 <- 1- (sum.res/var.y); r2
## [1] 0.7697

#function
r2_general <-function(preds, actual){ 
  return(1- sum((preds - actual) ^ 2)/sum((actual - mean(actual))^2))
}
preds<-labormodel$fitted.values
actual <- y
r2_general(preds,y)
## [1] 0.7697

8.3 F檢定

F檢定可檢定模型中的所有自變數對於依變數是否有聯合解釋的能力，虛無假設為「\(\beta_{1}=\beta_{2}=\ldots=\beta_{k}=0\)」。對立假設為至少有一個係數不為0。

F檢定為：

\[F=\frac{MSR}{MSE}\sim F_{k,n-k-1}\]

• 其中MSR=SSR/k，k為自變數的數目。SSR=\(\sum(\hat{Y}-\bar{Y})^2\)。
• MSE=SSE/(n-k-1)=(SST-SSR)/(n-k-1)。SST=\(\sum(Y-\bar{Y})^2\)。SSE=\(\sum(Y-\hat{Y})^2\)。

決策原則為如果\(F>F_{k,n-k-1}\)，則拒斥虛無假設。

計算F值可以借用\(R^2\)，因為：

\[R^2=1-\frac{1}{1+\frac{F\times k}{n-k-1}}\]

• 所以F值可以寫成：

\[F=\frac{R^2}{1-R^2}\times \frac{n-k-1}{k}\]

• 連同前面的\(R^2\)函數，整合成F值函數如下：

#calculate F step-by-step
#function
r2_general <-function(preds, actual){ 
  return(1- sum((preds - actual) ^ 2)/sum((actual - mean(actual))^2))
}
#model
labormodel<-lm(v2 ~ v3 + labinc , data=UKHouseOfCommons)
#R2
preds<-labormodel$fitted.values
actual <- pscl::UKHouseOfCommons$v2
r2=r2_general(preds, actual)
cat('R squared:', r2, '\n')
## R squared: 0.7697
n=nrow(UKHouseOfCommons);n
## [1] 521
k=2
r2new=r2/(1-r2)
#F2
F_stat<-function(r2new, n, k){
  return((r2new*(n-k-1)/k))
}
cat('F-statistics:', F_stat(r2new, n, k))
## F-statistics: 865.7

---

9 複迴歸例題

9.1 氣溫的影響因素

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

fit2 <- with(airquality, lm(Temp ~ Wind+factor(Month)))
stargazer::stargazer(fit2,  title='氣溫的迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 9.1: 氣溫的迴歸模型

	Dependent variable:
	Temp
Wind	-0.743***
	(0.149)
factor(Month)6	12.540***
	(1.594)
factor(Month)7	16.360***
	(1.618)
factor(Month)8	16.320***
	(1.624)
factor(Month)9	10.280***
	(1.596)
Constant	74.190***
	(2.054)
Observations	153
R2	0.588
Adjusted R2	0.574
Residual Std. Error	6.175 (df = 147)
F Statistic	42.030*** (df = 5; 147)
Note:	p<0.1; p<0.05; p<0.01

• 表9.1顯示：

1. 風力每增加一個單位，溫度減少0.7434個單位。標準誤為0.1488。\(t\)值為-4.996，拒斥風力與氣溫沒有關係的虛無假設。
   
2. 相對於五月，六月的溫度增加12.5436度。標準誤為1.5943。

• 計算\(R^2\)：

y=airquality$Temp
mean.y<-mean(y)
sum.res<-sum((y - fit2$fitted.values)^2)
var.y <- sum((y - mean.y)^2)
r2 <- 1- (sum.res/var.y); r2
## [1] 0.5884

#function
r2_general <-function(preds, actual){ 
  return(1- sum((preds - actual) ^ 2)/sum((actual -    mean(actual))^2))
}
preds<-fit2$fitted.values
actual <- y
r2_general(preds,y)
## [1] 0.5884

9.2 新生兒體重的迴歸模型

Using R for Introductory Statistics的例題10.5要建立一個新生兒體重迴歸模型如下：

\[\rm{wt} \sim gestation+age+ht+wt1+dage+dht+dwt\]

• 估計結果如表9.2：

library(UsingR)
## Warning: package 'UsingR' was built under R version 3.6.2
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 3.6.2
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: HistData
## Warning: package 'HistData' was built under R version 3.6.2
## Loading required package: Hmisc
## Warning: package 'Hmisc' was built under R version 3.6.2
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.6.2
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.6.2
## Loading required package: Formula
## Warning: package 'Formula' was built under R version 3.6.2
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:xtable':
## 
##     label, label<-
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## 
## Attaching package: 'UsingR'
## The following object is masked _by_ '.GlobalEnv':
## 
##     iq
## The following object is masked from 'package:survival':
## 
##     cancer
wt.lm1 <- with(babies, lm(wt ~ gestation+age+ht+wt1+dage+dht+dwt))
stargazer(wt.lm1, align=TRUE, title='新生兒體重迴歸模型1',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 9.2: 新生兒體重迴歸模型1

	Dependent variable:
	wt
gestation	0.013*
	(0.007)
age	0.056
	(0.102)
ht	0.511***
	(0.123)
wt1	-0.006
	(0.004)
dage	0.024
	(0.077)
dht	0.123
	(0.133)
dwt	-0.004
	(0.005)
Constant	73.040***
	(11.520)
Observations	1,236
R2	0.022
Adjusted R2	0.016
Residual Std. Error	18.090 (df = 1228)
F Statistic	3.873*** (df = 7; 1228)
Note:	p<0.1; p<0.05; p<0.01

• 估計模型後，發現自變數有99或者999拒答的情形，以及篩選gestation小於350的研究對象，我們可以用\(\tt{subset}\)這個指令，設定子資料的條件。 重新估計結果如表9.3：

library(UsingR)
wt.lm2 <- lm(wt ~ gestation+age+ht+wt1+dage+dht+dwt, data = babies,
    subset = (gestation<350 & age< 99 & ht< 99 & wt1< 999 & dage< 99 & dht< 99 & dwt < 999))
stargazer(wt.lm2, align=TRUE, title='新生兒體重迴歸模型2',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 9.3: 新生兒體重迴歸模型2

	Dependent variable:
	wt
gestation	0.462***
	(0.040)
age	0.138
	(0.188)
ht	1.216***
	(0.285)
wt1	0.029
	(0.034)
dage	0.059
	(0.165)
dht	-0.066
	(0.270)
dwt	0.078**
	(0.033)
Constant	-105.500***
	(23.350)
Observations	700
R2	0.215
Adjusted R2	0.207
Residual Std. Error	16.410 (df = 692)
F Statistic	27.090*** (df = 7; 692)
Note:	p<0.1; p<0.05; p<0.01

圖9.1顯示，誤差項的變異程度應該不會隨著自變數而變化。

tmp <- data.frame(fv=fitted(wt.lm2), rv=resid(wt.lm2))
ggplot2::ggplot(data=tmp, aes(x=fv, y=rv)) +
    geom_point() +
    geom_hline(yintercept = 0, col='#FFAA11') +
    labs(x='fitted value', y='residual value', caption='UsingR::babies') 

Figure 9.1: 新生兒體重迴歸模型2的檢驗

10 轉換變數

迴歸假設的其中之一是參數成線性組合，也就是說迴歸係數只有一次項，不能有二次項。例如\(Y=\beta_{0}+\beta_{1}^{2}X+u\)違反線性假設。

在John Fox所寫的An R and S-Plus Companion to Applied Regression介紹如何用次方項轉換自變數，讓自變數變成對稱或者常態分配，也就是\(x'=x^{p}\)，當p=2，\(x'=x^2\)，當p=1/2，\(x'=\sqrt{x}\)，當p=-1，\(x'=\frac{1}{x}\)，以此類推。Fox主張設定\(x^{0}=\rm{log}_{e}x\quad\forall\hspace{.3em} x>0\)。

另外一個常用的轉換是：

\[\rm{logit}(p)=\frac{p}{1-p}\]

• logit轉換0到1的變數為趨近\(-\infty\)到\(\infty\)，也就是

\[(0,1)\rightarrow(-\infty, \infty)\]

10.1 確認常態分佈

我們可以用\(\tt{car::qqPlot()}\)確認變數的每一個分位符合常態分配。也有人用\(\tt{qqline()}\)與\(\tt{qqnorm()}\)。例如我們模擬一筆資料，圖10.1顯示部分的觀察值並沒有落在常態分佈的範圍內：

set.seed(02139)
y<-rnorm(30, 0, 1)
car::qqPlot(y)

Figure 10.1: qqplot散佈圖

## [1] 16  1

我們改確認\(\tt{carData::Prestige\$income}\)這個變數是否符合常態分配，請見圖10.2。

car::qqPlot(carData::Prestige$income)

Figure 10.2: qqplot散佈圖1

## [1]  2 24

• 我們用\(\tt{bcPower()}\)這個指令轉換變數，\(\tt{bcPower()}\)指的是以下列的公式轉換變數：

\[\begin{align*} x'=x^{p}=\begin{cases} \frac{x^p-1}{p} & \rm{when}\hspace{.3em} p\neq 0\\ \rm{log}{x} & \rm{when}\hspace{.3em} p = 0 \end{cases} \tag{10.1} \end{align*}\]

圖??的第2張圖似乎比較接近常態分佈，也就p=1/3。

library(car)
par(mrow=c(2,2))
for (p in c(-1/2, 1/3, 1/2,  1)){
  qqPlot(bcPower(carData::Prestige$income, p),
  ylab=paste('transformed income, p=', p))
}

Figure 10.3: qqplot散佈圖2

Figure 10.4: qqplot散佈圖2

Figure 10.5: qqplot散佈圖2

Figure 10.6: qqplot散佈圖2

• 我們可以用以下的圖10.7來驗證(10.1)，例如p=1/2，產生的圖就像圖10.7的第3張圖。

p=1/2
y=carData::Prestige$income
qqPlot((y^p-1)/p)

Figure 10.7: qqplot散佈圖3

## [1]  2 24

10.2 轉換線性

我們可以運用類似的方法轉換雙變數，使其成線性的關係。跟轉換單一變數相同，我們用雙迴圈進行試誤，找到比較好的轉換，如圖??。第一張圖似乎是比較好的轉換，也就是(0, 1/3)

par(mfrow=c(2, 3))

for (p1 in c(0,  1)){
  for (p2 in c(1/3, 1/2, 2/3)){
  scatterplot(bcPower(carData::UN$ppgdp,p1), bcPower(carData::UN$infantMortality,p2),
  xlab=paste('transformed PP GDP, p=', p1),  ylab=paste('transformed infant mortality, p=', 
              sprintf("%.2f",p2)))
  }
}

Figure 10.8: 雙變數轉換後散佈圖

Figure 10.9: 雙變數轉換後散佈圖

Figure 10.10: 雙變數轉換後散佈圖

Figure 10.11: 雙變數轉換後散佈圖

Figure 10.12: 雙變數轉換後散佈圖

Figure 10.13: 雙變數轉換後散佈圖

10.3 實例：美軍受傷人數

首先我們以\(\tt{ggpubr::gghistogram}\)的直方圖顯示美軍在各個戰役受傷人數的人數，圖10.14顯示大多數的觀察值落在0附近，少數是超過1萬，平均值則是數千。

dt<-here::here('data','killwound.txt')
mili <- read.table(dt, sep=',', header = T)
#with(mili, hist(Wounded, breaks = 20,
#        xlim=c(0, 7e+05)))

ggpubr::gghistogram(mili,x='Wounded', bins = 40,
      fill='lightgray',add='mean', rug = TRUE)
## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.
## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.

Figure 10.14: 受傷人數直方圖

• 圖10.15則以點狀圖排列受傷人數的戰役，以第二次世界大戰最多。

library(lattice)
w_order<-order(mili$Wounded, decreasing=F)
mili_order<-mili[w_order,]
dotplot(mili_order$Location ~ mili_order$Wounded, cex=1.1)

Figure 10.15: 受傷人數點狀圖1

• 圖10.16則以\(\tt{lattice}\)套件的點狀圖排列受傷人數的戰役：

library(lattice)
dotchart(mili_order$Wounded,labels=mili_order$Location,pch=16,cex=1,
         xlab="Number of American Soliders Wounded",aspect=0.8)

Figure 10.16: 受傷人數點狀圖2

如果受傷人數變數不經轉換，跟死亡人數的散佈圖呈現線性關係如圖10.17，但是許多觀察值聚集在左下角：

dt<-here::here('data','killwound.txt')
mili <- read.table(dt, sep=',', header = T)
with(mili, plot(Wounded ~ Killed, main="",xlab="Number of American Soliders Killed",
 ylab="Number of American Soliders Wounded"),pch=16)
mod1<-with(mili,lm(Wounded ~ Killed))
abline(mod1, col="red", lwd=2)
legend("topleft",legend=expression(paste("Wounded=-0.002+1.21 Killed"," ",R^2==0.86)),
       lty=1,col="red")

Figure 10.17: 受傷與死亡人數散佈圖1

• 圖10.18顯示取log之後的受傷人數的點狀圖，可以看到取log之後的確把數字變成線性。但是因為受傷人數有0，所以不能用\(\tt{bcPower}\)這個指令。

log.wounded <- log(mili$Wounded)
log.killed<-bcPower(mili$Killed, 0)
newdt<-data.frame(log.wounded, log.killed,mili$Wounded,mili$Killed,mili$Location)
is.na(newdt)<-sapply(newdt,is.infinite)
#newdt

worder<-order(newdt$log.wounded,decreasing=F)
nnewdt<-newdt[worder,]
#jpeg("logwounded.jpg")
dotchart(nnewdt$log.wounded,labels=newdt$mili.Location,pch=16,cex=1,
         xlab="Number of American Soliders Wounded")

Figure 10.18: log轉換變數之點狀圖

• 圖10.19顯示兩個變數取log之後，得到的迴歸模型有更高的判定係數，\(\beta_{1}\)則下降，顯示取log是必要的。

with(newdt, plot(log.wounded ~ log.killed, main="",xlab="log(Number of American Soliders Killed)",
                 ylab="log(Number of American Soliders Wounded(",pch=16))
mod3<-with(newdt,lm(log.wounded ~ log.killed))
abline(mod3,col="red",lwd=2)
legend("topleft",legend=expression(paste("Wounded=0.82+0.83 killed"," ",R^2==0.84)))

Figure 10.19: 受傷與死亡人數散佈圖2

以下將介紹三種轉換自變數或者依變數的模型。有興趣的同學可參考Kenneth Benoit的說明。

10.4 linear-log模型

log有以下特性：

1. log(a)的指的是以\(e\)為底對a取自然對數。\(e\)是數學常數，約等於2.718282(\(e\approx 2.718282\))。
2. log(a)的正式定義是積分，公式為：

\[\rm{log}(a)=\int_{1}^{a}\frac{1}{x}dx\]

• 把上面的公式寫成R的自訂函數，然後利用\(\tt{integrate}\)求1到a之間的面積，等於是log(a)，例如a=3，代入積分公式如下：

h=function(x){1/x}
a=3
integrate(h, 1, 3)
## 1.099 with absolute error < 7.6e-14
log(3)
## [1] 1.099

• 請注意\(e\)=exp(1)=2.718282。
  
• log(e)=log(exp(1))=log(2.718282)=1

h=function(x){1/x}
a=exp(1)
integrate(h, 1, a)
## 1 with absolute error < 1.1e-14

對X取log，也就是log(X)，稱為linear-log 模型，寫成

\[\hat{Y}=\hat{\beta}_{0}+\hat{\beta}_{1}\rm{log}X\]

我們模擬自變數x以及誤差項e，然後模擬y為取自然對數的x以及e的線性組合。畫圖10.20對照x有無取自然對數的散佈圖：

x <- seq(0.1,5,length.out = 100)
set.seed(1)
e <- rnorm(100, mean = 0, sd = 0.2)
y <- 1.2 + 0.5 * log(x) + e
llm1 <- lm(y ~ x)
llm2 <- lm(y ~ log(x))
par(mfrow=c(1,2))
car::crPlot(llm1, variable="x")
car::crPlot(llm2, variable="log(x)")

Figure 10.20: linear-log模型

因為\(\rm{log}(e)=1\)，所以\(\rm{log}X\)增加一單位可寫成logX+1，過程如下：

\[\begin{align*} \rm{log}X+1=\rm{log}X+\rm{log}e=\rm{log}(eX) \end{align*}\]

因此，解釋linear-log模型可說當X乘以2.718，Y變動\(\hat{\beta}_{1}\)單位。

X乘以2.72等於X增加1.72倍，或者增加172%。那麼如果X增加10%，Y變動了多少\(\hat{\beta}_{1}\)單位？

如果X本身（不是log(x)）從1增加\(m\)百分比，可以寫成\(\frac{100}{100}+\frac{m}{100}\)，加上\(\hat{\beta}\)，Y的變動為\(\hat{\beta}_{1}\times \rm{log}(1+m\%)\)。如果\(m=10\)，Y的變動為\(\rm{log}(1.1)\times\hat{\beta}_{1}=0.095\hat{\beta}_{1}\)。如果當\(m=1\)，log(1+1%)=log(1.01)\(\approx 0.01\)。所以當X增加1%，Y變動\(\hat{\beta}_{1}/100\)。

數學上可以證明對\(y=\beta\hspace{.2em} \rm{log}x\)的x微分，\(\frac{d}{dx}y=\beta\frac{1}{x}dx\)

因此在linear-log模型，當X增加1%，Y增加\(\hat{\beta}_{1}/100\)，等於\(100\times\hat{\beta}_{1}\%\)，也就是\(\frac{1}{100}\)的倒數。

• 表10.1顯示\(\beta_{0}\)以及\(\beta_{1}\)，我們分別乘以(1, 1.01)，(1, 1.02)，(1, 1.03)，得到的\(\hat{Y}\)分別比當x=1時，多了0.005，也就是\(\beta_{1}/100=0.5/100\)，然後是多了0.01，也就是\(2\times \beta_{1}/100=1/100\)，然後是\(3\times \beta_{1}/100=1.5/100\)。

stargazer::stargazer(llm2,  title='linear-log迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 10.1: linear-log迴歸模型

	Dependent variable:
	y
log(x)	0.500***
	(0.021)
Constant	1.222***
	(0.023)
Observations	100
R2	0.850
Adjusted R2	0.849
Residual Std. Error	0.181 (df = 98)
F Statistic	556.100*** (df = 1; 98)
Note:	p<0.1; p<0.05; p<0.01

cat('X=1, Y=',sum(as.vector(llm2$coefficients)*c(1, 1)),'\n')
## X=1, Y= 1.722
x_1=sum(as.vector(llm2$coefficients)*c(1, 1))
cat('X+0.01,Y+',sum(as.vector(llm2$coefficients)*c(1, 1.01))-x_1,'\n')
## X+0.01,Y+ 0.004998
cat('X+0.02,Y+',sum(as.vector(llm2$coefficients)*c(1, 1.02))-x_1,'\n')
## X+0.02,Y+ 0.009996
cat('X+0.03,Y+',sum(as.vector(llm2$coefficients)*c(1, 1.03))-x_1,'\n')
## X+0.03,Y+ 0.01499

10.5 log-linear模型

如果對Y取log， 稱為log-linear模型。\(\rm{log}Y=\hat{\beta}_{0}+\hat{\beta}_{1}X\)

在log-linear模型，當自變數X增加1個單位，\(\rm{log}Y\)增加\(\hat{\beta}_{1}\)單位。

因為：

\[\mathrm{log}(e^{a})=a\]

• 所以：

\[\hat{\beta}_{1}=\mathrm{log}(e^{\hat{\beta}_{1}})\]

• log(Y)增加\(\hat{\beta}_{1}\)可以寫成\(\rm{log}(e^{\hat{\beta}_{1}}Y)\)。這是因為：

\[\begin{align*} \rm{log}Y+\hat{\beta}_{1}&=\rm{log}(Y)+\rm{log}(e^{\hat{\beta}_{1}})\\ &=\mathrm{log}(e^{\hat{\beta}_{1}}Y) \end{align*}\]

• 用一組方程式來計算X增加1個單位時Y的預期變動：

\[\begin{align*} \mathrm{log}Y_{1} & =\hat{\beta}_{0}+\hat{\beta}_{1}X_{1} \tag{10.2} \end{align*}\] \[\begin{align*} \mathrm{log}Y_{2} & =\hat{\beta}_{0}+\hat{\beta}_{1}(X_{1}+1) \tag{10.3} \end{align*}\]

• 方程式(10.3)-方程式(10.2)

\[\begin{align*} \hat{\beta}_{1} & =\mathrm{log}Y_{2}-\mathrm{log}Y_{1}\\ & =\mathrm{log}\frac{Y_{2}}{Y_{1}}\\ \frac{Y_{2}}{Y_{1}}=e^{\hat{\beta}_{1}}\\ Y_{2} = Y_{1}\times e^{\hat{\beta}_{1}} \end{align*}\]

因此在log-linear模型，當X增加1單位，Y乘上\(e^{\hat{\beta}_{1}}\)倍。當X增加\(m\)單位，Y變為\(m\times e^{\hat{\beta}_{1}}\)。

如果\(\hat{\beta}_{1}\)很小(<0.69)，當X增加1單位，Y乘上\(e^{\hat{\beta}_{1}}\)倍，也就是增加\(e^{\hat{\beta}_{1}}-1\)單位，或者是\(100 \times (e^{\hat{\beta}_{1}}-1)\%\)

\(\hat{\beta}_{1}\)解釋成當X變動1單位，Y乘上\(e^{\hat{\beta}_{1}}\)，等於Y的預期變動量。

結論：\(\hat{\beta}_{1}\)可解釋成X變動1個百分比時，Y增加\(e^{1.01\times \hat{\beta}_1}\)倍。

想像Y是各個國家貨幣兌換美金，X是工業化程度。Y可以取log。

10.6 log-log模型

如果X,Y各取對數，稱為log-log模型。

\[ \mathrm{log}Y=\hat{\beta}_{0}+\hat{\beta}_{1}\mathrm{log}X \]

對上面的方程式的左右兩邊進行微分

\[\begin{align*} \partial\mathrm{log}Y_{1} & =\partial\hat{\beta}_{1}\mathrm{log}X_{1} \end{align*}\]

• 因為\(\partial \mathrm{log}X=\frac{1}{X}\partial X\)

• 所以\(\hat{\beta}_{1}=\mathcal{\frac{\partial Y/Y}{\partial X/X}}\)

\(\partial Y/Y\)代表Y的變動百分比(percent change)，\(\frac{\partial Y/Y}{\partial X/X}\)代表兩者的比值(odds ratio)，又稱為elasticity。

結論是：X變動1個百分點，則Y變動\(\hat{\beta}_{1}\)百分點。

10.7 範例

10.7.1 單迴歸模型

\(\tt{carData::Leinhardt}\)有income跟Infant mortality兩個變數，我們先檢視這兩個變數的分佈是否呈現線性分佈如圖10.21：

library(carData)
ggpubr::ggscatter(data=Leinhardt,x='income', y='infant',  color = "#22AF11", fill = "white", add = "loess",shape=4,
  cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 2000, label.sep = "\n"))

Figure 10.21: 所得與嬰兒死亡率散佈圖1

• 圖10.22顯示，大部分殘差值集中0的水平線的右邊，也就是自變數與依變數並不是線性的關係（自變數變動1單位、依變數固定變動若干單位）。

M1<-lm(infant~income, data=Leinhardt)
Mres<-data.frame(residual=M1$residuals,fitted=M1$fitted.values)
G<-ggpubr::ggscatter(data=Mres,x='fitted', y='residual',  color = "#0F1EDC")
G + geom_hline(yintercept = 0, linetype='dashed') 

Figure 10.22: 所得與嬰兒死亡率迴歸模型的殘差值與預測值散佈圖

我們用\(\tt{log(x,a)}\)轉換兩個變數，其中a是對數的底，轉換後成線性分佈如圖10.22：

Leinhardt$log.income<-log(Leinhardt$income,10)
Leinhardt$log.infant<-log(Leinhardt$infant,10)
ggpubr::ggscatter(data=Leinhardt,x='log.income', y='log.infant',  color = "#22AF11", add = "reg.line",shape=4,
  cor.coef = TRUE, 
   cor.coeff.args = list(method = "pearson", label.x = 2.8, label.sep = "\n"))

Figure 10.23: 所得與嬰兒死亡率散佈圖2

• 接下來我們估計這個新變數構成的迴歸模型：

LM1<-lm(log.infant~log.income,data=Leinhardt)
stargazer::stargazer(LM1,  title='所得與嬰兒死亡率',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 10.2: 所得與嬰兒死亡率

	Dependent variable:
	log.infant
log.income	-0.512***
	(0.051)
Constant	3.103***
	(0.137)
Observations	101
R2	0.502
Adjusted R2	0.497
Residual Std. Error	0.298 (df = 99)
F Statistic	99.840*** (df = 1; 99)
Note:	p<0.1; p<0.05; p<0.01

• 圖10.24顯示，大部分殘差值在0的水平線隨機地變動，也就是殘差值並沒有包含我們沒估計到的資訊。

LM<-data.frame(residual=LM1$residuals,fitted=LM1$fitted.values)
G<-ggpubr::ggscatter(data=LM,x='fitted', y='residual',  color = "#AA33C1")
G + geom_hline(yintercept = 0, linetype='dashed') 

Figure 10.24: 所得與嬰兒死亡率迴歸模型的殘差值與預測值散佈圖

根據表10.2，我們可以得到結論為收入每增加\(1\%\)，平均而言，嬰兒死亡率下降\(.51\%\)。兩者呈負相關的關係。相關資料可以參考William Jacoby的說明。

10.7.2 複迴歸模型

如果是複迴歸模型，也就是有超過一個自變數，我們也可以轉換變數。

Penn State University的Eberly College of Science提供線上的統計課程資料，其中一筆是房地產的資料，我們想解釋房屋的銷售價格，自變數是房子的面積以及是否有冷氣，加上兩者的交互作用。模型寫成：

\[Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}++\beta_{3}X_{1}X_{2}\]

• 其中\(X_{1}\)是面積、\(X_{2}\)是有無冷氣。
• 我們估計線性迴歸模型之後，得到的預測值以及殘差值的散佈圖如圖10.25，顯示殘差值並非隨機分佈。

est<-here::here('data','realestate.txt')
dat<-read.table(est)
dat<-dat[-1,]
dat<-dat%>%dplyr::select(V1,V2,V5)%>%
  mutate_all(as.numeric)
EM<-lm(V1~V2+V5+V2:V5, dat)
EMR<-data.frame(residual=EM$residuals,fitted=EM$fitted.values)
G<-ggpubr::ggscatter(data=EMR,x='fitted', y='residual',  color = "#22CCDC")
G + geom_hline(yintercept = 0, linetype='dashed') 

Figure 10.25: 房價迴歸模型的殘差值與預測值散佈圖

• 注意我們用\(\tt{mutate\_all(as.numeric)}\)把三個變數從字串轉為數字。

接著，我們把V1，V2用\(\tt{mutate\_at(x,log)}\)轉換為自然對數，估計線性模型，圖10.26，顯示殘差值已經是隨機分佈。

est<-here::here('data','realestate.txt')
dat<-read.table(est)
dat<-dat[-1,]
dat<-dat%>%dplyr::select(V1,V2,V5)%>%
  mutate_all(as.numeric) %>%
  mutate_at(c('V1','V2'), log)
EM<-lm(V1~V2+V5+V2:V5, dat)
EMR<-data.frame(residual=EM$residuals,fitted=EM$fitted.values)
G<-ggpubr::ggscatter(data=EMR,x='fitted', y='residual',  color = "#BEAADC")
G + geom_hline(yintercept = 0, linetype='dashed') 

Figure 10.26: 轉換後房價迴歸模型的殘差值與預測值散佈圖

最後呈現估計結果在表10.3：

est<-here::here('data','realestate.txt')
dat<-read.table(est)
dat<-dat[-1,]
dat<-dat%>%dplyr::select(V1,V2,V5)%>%
  mutate_all(as.numeric) %>%
  mutate_at(c('V1','V2'), log)
EM<-lm(V1~V2+V5+V2:V5, dat)
stargazer::stargazer(EM,  title='房價、面積、冷氣',covariate.labels = c('Square footabe','Air condition','Square footage X Air condition'), column.labels = c('Price'),
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 10.3: 房價、面積、冷氣

	Dependent variable:
	V1
	Price
Square footabe	0.198*
	(0.101)
Air condition	-0.849***
	(0.280)
Square footage X Air condition	0.270***
	(0.060)
Constant	2.871***
	(0.457)
Observations	521
R2	0.584
Adjusted R2	0.582
Residual Std. Error	0.575 (df = 517)
F Statistic	242.200*** (df = 3; 517)
Note:	p<0.1; p<0.05; p<0.01

• 結果顯示面積與面積與空調的交互作用都會影響房價。如果沒有空調，面積增加1個百分點，價格提高0.97個百分點。如果有空調，面積增加1個百分點，價格提高1.2個百分點。

11 結論

複迴歸模型涉及更多自變數，因此需要聯合檢定否所有迴歸係數等於0。

複迴歸模型可能有部分自變數之間有交互作用，需要納入考慮交互作用變數以及變數本身，也就是一次項以及交乘項。

依變數與自變數之間如果不是線性關係，可以考慮自然對數轉換，讓殘差值成隨機分佈。

自變數或者依變數轉換後有log-linear, linear-log, log-log等模型，需注意如何詮釋當X變動1個單位，Y的預期變動程度。對於三種模型的解釋有興趣的同學可參考這個網頁。

---

12 作業

1. 請寫語法計算F值，公式如下。並請用以下語法產生的資料測試。

\[F=\frac{MSR}{MSE}\] \[\sum (\hat{Y}-\bar{Y})^2=SSR\] \[MSR=\frac{SSR}{k}\] \[\sum (Y-\hat{Y})^2=SSE\] \[MSE=\frac{SSE}{n-k-1}\]

set.seed(02139)
X1 <- rnorm(200, 0, 1); X2 <- rnorm(200, 0, 0.5); X3 <- rnorm(200, 0, 0.75)
Y<-c(sample(X1, 130), sample(X2, 20), sample(X3, 30), rnorm(20, 0, 0.75))

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

3 回答現代統計學第二版的15.6題。

4. \(\tt{pscl::admit}\)這筆資料中有GRE的字彙分數(gre.verbal)、GRE的數學分數(gre.quant)以及性別等三個變數。請建立數學分數對字彙分數以及性別的迴歸方程式，並且解釋迴歸係數之含意。

5. 請使用現代統計學第二版的dataet 5.csv檔案（已經上傳到政大moodle），回答以下兩題：

1. 分別估計聯電與台積電的每個月報酬率變化對大盤股價變化率的迴歸方程式，報酬率變化率等於：

\[\Delta R=\frac{R_{t}-R_{t-1}}{R_{t-1}}\]

• 大盤股價變化率則為：

\[\Delta S=\frac{S_{t}-S_{t-1}}{S_{t-1}}\]

• 方程式為：

\[\rm{\Delta R}=\beta_{0}+\beta_{1}\Delta S+\epsilon\]

2. 合併聯電與台積電的資料，再估計每個月報酬率變化對大盤股價變化率的迴歸方程式如下，並且指出新的迴歸係數等於前面迴歸方程式的組合。（提示：假設台積電編碼為1，代入D=1可以得到台積電模型的迴歸係數）

\[\rm{\Delta R}=\alpha+\gamma\Delta S+\delta_{1}\times D+\delta_{2}(\Delta S\times D)+\epsilon\]

13 更新日期

## 最後更新日期 05/19/2022