################################################# ## 線性迴歸(Linear Regression) ## Ref: John Verzani, Chapter 3.4 & Chapter 10 ## Date: 2011.12.06 ## Instructor: 劉正山 ################################################# ## 從這一講開始,我們將從統計學進入計量經濟學(econometrics)的領域,但 ## 所學過的統計學中母體、假設檢定、顯著性等概念還是會繼續派上用場。 ## 迴歸的基本觀念: ## 一、理論、假設、因果關係、相關分析與迴歸分析 ## 二、線性模型(或作模式)的結構 ## 三、用來估計參數的最小平方法(ordinal least squares, OLS) ## 四、殘差項(the error term): 預設為相互獨立且呈常態分佈 ## 五、總變異(Total sum of squared deviations, TSS), ## 可解釋變異(Regression sum of squared deviations, RSS), ## 隨機變異(Error sum of squared deviations, ESS) ## 六、線性迴歸中的母體參數:alpha, beta, residuals(error term) ## 七、迴歸結果的判讀:迴歸係數beta的t檢定 & R^2. ## X1->Y ## X2->Y ## Y = a + b*X1 模式/模型 ## H0: X1對Y無影響 (b=0) ## H1: X1對Y有影響 (b!=0) ## 迴歸的進階觀念: ## 一、變數的挑選與指定錯誤(specification error) ## 二、控制變數(control variables) ## 三、虛擬變數(dummy variables) ## 四、模型的適配(the fitting of the model) ## 五、交叉變數(interaction variables) ## 六、共線性問題及其他進行線性迴歸分析的條件 ##============== ## 1. 簡單線性迴歸 ##============== ## 第一步:要建立模型或假設你一定要有理論和文獻。 ## 第二步:依據你的假設來選擇變數。取得變數的資料後,先檢視變數的常態性, ## 再來探討變數之間的關聯性. ## 資料檔:Lewis-Beck (1995), p.53-67 # load("/home/ashan/teaching/statistics/data/wgcoll.rda") load("d:/teaching/statistics/data/wgcoll.rda") attach(wgc) fix(wgc) ## Y = a + bX ## "~" 以左為依變數,以右為自變數 plot(aa ~ pe) ## academic average & parents' education cor(aa,pe) ## correlation coefficient ##注意:只有當兩個變數都是連續 ## 型變數的時候才能使用相關係數,否則就只能使用CrossTable() ## 第三步:進行迴歸分析、試畫迴歸線、檢查是否有極端個案 model1 <- lm(aa ~ pe) model1 plot(aa~pe, cex=10*sqrt(cooks.distance(model1))) # 用 Cook's distance的 # 概念來找出影響力最大 # 的個案 (Verzani, # p.287) abline(h=mean(aa), lty=2) # 畫水平線 abline(v=mean(pe), lty=3) # 畫垂直線 abline(model1, data=wgc, lty=1, col="red") with(wgc, identify(aa ~ pe, labels=id)) #手動點出1個極端值(左下角) abline(lm(aa ~ pe, data=wgc, subset=c(-8)), lty=2, col="blue") # 在圖上加圖示 labels <- c("all", "-8") legend(locator(1), legend=labels, lty=1:2, col=c("red","blue")) ## 第四步:迴歸結果的判讀 summary(model1) ## 這個模型有解釋力嗎?(F檢定的虛無假設是(所有)自變數的迴歸參數為零) ## 這個模型解釋了多少依變數的變異量? ## R^2 = RSS/TSS (即:可被迴歸解釋的變異/總變異)。R^2就是兩個變數相關 ## 係數的平方。 cor(aa, pe)^2 ##在複迴歸中我們常用adjusted-R^2,它考量了自變數的增加對模式解釋力產生 ##的負面影響。 ##==================== ## 2. 複迴歸/多變數迴歸 ##==================== ## 為方便學習起見,先說明多變數回歸的程式碼; ## 在實際研究時,應先進行第三節所提到的檢查。 ## 2.1 基本語法和基本檢查 ## 注意:你是不是清楚了加進控制變數的邏輯了呢?控制變項不宜隨便加入模型 ## 中(當然,依理論而建構的模型裡有的變數也不宜輕易刪除),必須是當在理 ## 論和常理上這個變數會影響到依變數和自變數的情形下,才將這個變數加入作 ## 為控制變數。在迴歸模型裡,控制變數的作用與自變數並無根本上的差異,增 ## 加或減少控制變數非隨性而為、隨意增加,而是應該以理論、常理為基礎,把 ## 控制變數的數目控制在合理的範圍內。 ## 如果居住地(community type: 0=urban, 1=rural)對學生成績有影響呢? model2 <- lm(aa ~ pe + c, data=wgc) summary(model2) plot(fitted(model2), resid(model2)) ## 虛擬變數的使用:兩個虛擬變數可以表示三種名目變數的情境 ## 學生成就動機(student motivation)對學生成績有影響嗎? sm1 <- ifelse(sm==1,1,0) sm2 <- ifelse(sm==2,1,0) model3 <- update(model2, .~.+ sm, data=wgc) # model3 <- lm(aa ~ pe + c + sm1 + sm2, data=wgc) summary(model3) plot(fitted(model3), resid(model3)) ## 交叉變數 ## 如果家長的教育程度其實是受到其居住地的影響呢? ## (the effect of X1 dependents on X2) model4 <- lm(aa ~ pe*c, data=wgc) ## model4 <- lm(aa ~ pe + c + pe:c, data=wgc) summary(model4) ## 2.2 模型選擇(model selection) ## Akaike's information criterion (AIC): 愈小愈好。 AIC(model1) AIC(model2) AIC(model3) AIC(model4) ## 寶寶出生重量的迴歸模型(p.303, Ex 10.5) library(UsingR) data(babies) subbabies = subset(babies, gestation < 350 & age < 99 & ht <99 & wt1 <999 & dage < 99 & dht <99 & dwt <999) res.lm= lm(wt ~ gestation + age + ht + wt1 + dage + dht + dwt, data = subbabies) summary(res.lm) plot(fitted(res.lm), resid(res.lm)) plot(res.lm) # 注意:R需要你按Enter鍵來切換四張圖;也注意一下261這個 # outlier babies[260:263,] # the 261st (id=4604) has a very short gestation # period, so we think about dropping it. subbabies2 = subset(subbabies, id!=4604) res.lm2 = lm(wt ~ gestation + age + ht + wt1 + dage + dht + dwt, data = subbabies2) plot(res.lm2) ## 課堂練習: 請問性別對成績有影響嗎?先用cor檢查一下。如果這個變數加進了 ## model3會有什麼影響嗎?請用lm和update指令做個model5來印證你的說法。 model5 <- update(model3, .~. + g) summary(model5) ## 補充閱讀: ##------------------------------------------------------------------ ## 迴歸係數的標準化(standardized regression coefficients):若研究者要 ## 進行迴歸參數之間的比較,則必須先將迴歸係數標準化。因此,未標準化之前 ## 的迴歸係數不宜拿來進行比較。 ## Introduction: ## http://en.wikipedia.org/wiki/Standardized_coefficient Use R: ## http://tolstoy.newcastle.edu.au/R/help/06/08/33638.html John Fox's ## heuristic answer: ## http://finzi.psych.upenn.edu/R/Rhelp02a/archive/6038.html John ## Fox's perspective about using src: ## http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92774.html ## ------------------------------------------------------------------- ##=============================== ## 3. 變異量分析 ##============================== ## 嚴格來說,建立多變數線性模型之前必須先滿足以下幾個條件: ## 一、自變數分布呈線性(multivariate normality): 即自變數的分佈得(幾 ## 近)呈常態分配. ## 二、變異數齊一性(homoscedasticity; equality of variance): 每個樣本 ## 的殘差(即樣本觀察值到迴歸線的距離)變異量(幾近)相等,亦即,每個樣 ## 本點的殘差值不會因所有自變數的值放大而跟著放大或縮小; ## 三、任兩組殘差之間的獨立性: 每個殘差值之間相互獨立,不應該發生 ## autocorrelation of errors的情形; ## 四、自變數彼此之間的獨立性:要盡量避免共線性(multicollinearity)的 ## 發生,因為變數之間的高度共線性會導致迴歸參數在估計時的不穩定。 ## 五、殘差項與自變數之間的獨立性:即殘差項與自變數無關。 ## 六、依變數必須是連續變數:線性迴歸不適用於依變數為順序變數或名目變數的情形 ## 當變數個數大於兩個的時候,顯然我們無法用2D的圖形來作相關性分析;這時 ## 我們得透過對模型變異量的分析來判斷這些自變數是不是符合以上的條件。 ## 3.1 檢視所有的自變數整體來看是不是呈常態分佈(multivariate ## normality) hist(resid(model3)) plot(density(resid(model3))) ## a better way to examine the normality of the residuals (p.286) qqnorm(resid(model3)) qqline(resid(model3)) ## 當分佈不似常態時,必須分別檢視每個變數的常態性以及可能的極端觀察值 ## 如果是偏態出現問題,則必須使用變數的轉換(見本檔案最後一節的補充) ## 3.2 檢查殘差變異量是否相等 plot(resid(model3) ~ fitted(model3)) ## 我們希望看到水平線,水平線表示模型所解釋的變異量與未被解釋到的殘差變 ## 異量之間互不相關。如果這個預設不成立,則應使用weighted least ## squares (請自行延申學習)。 ## 你可以用R一次畫出所有可以用於分析資料變異的圖形: par(mfrow=c(2,2)) plot(model3) ## 3.3 檢查每個樣本之間相互獨立(或說是殘差值之間相互獨立) library(car) durbin.watson(model3) ## 注意:這個"Durbin-Watson"檢定的虛無假設是:殘差值之間相互獨立 ## 3.4 檢視共線性是否發生 library(car) vif(model3) # variance inflation factor; 若值大於4則表示變數間有高度共 # 線性的情形 ## 較傳統的作法是將所有自變數分別當作依變數去找出每個線性模型的R^2,再拿 ## 這每個R^2值去和原本模型(即依變數在等式左邊,自變數都在等式右邊)產 ## 生的R^2相比。如果你發現前者有個R^2 大於 後者的R^2,則表示該(罝於等 ## 式左邊的)自變數與某個其它自變數之間有高度的共線性。 ## 解決共線性問題的方法:1、增加樣本數;2、整合變數(將多個相關變數合併 ## 為一個);3、以圖表呈現多個模型(含原始模型和不含造成高度共線性的變 ## 數的模型)供讀者自行比對。 detach(wgc) ##===================== ## 4. 變數的轉換(進階補充) ##===================== ## data transformation (see also MLB, 1995, p.67-71) library(UsingR) attach(kid.weights) plot(weight ~ height, data=kid.weights, main="Before Transforming Height") m.before <- lm(weight ~ height, data=kid.weights) m.before abline(m.before) with(kid.weights, identify(height, weight, labels=age)) plot(weight ~ I(height^2), data=kid.weights, main="After Transforming Hieight") m.after <- lm(weight ~ I(height^2), data=kid.weights) m.after abline(m.after) with(kid.weights, identify(I(height^2), weight, labels=age)) detach(kid.weights) ################# ## new stuff ################ ## Analyzing Subsets # example1: # summary(lm(sr[ddpi >= 3] ~ dpi[ddpi >= 3], data = LS)) # example2: # summary(lm(sr ~ dpi + pop15 + pop75, subset = (ddpi >= 3), data = LS)) ## GLM poisson distribution # summary(glm(discoveries ~ year, family = poisson(log)))