在R軟件中實現(xiàn)回歸估計_第1頁
在R軟件中實現(xiàn)回歸估計_第2頁
在R軟件中實現(xiàn)回歸估計_第3頁
在R軟件中實現(xiàn)回歸估計_第4頁
在R軟件中實現(xiàn)回歸估計_第5頁
免費預(yù)覽已結(jié)束,剩余1頁可下載查看

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進行舉報或認領(lǐng)

文檔簡介

1、3.1 最小二乘法有三種方式可以實現(xiàn)最小二乘法的簡單線性回歸,假設(shè)數(shù)據(jù)byu1.1.1 lm(byu$salarybyu$age+byu$exper)2.2.2 lm(salaryage+exper,data=byu)3.3.3 attach(byu)lm(salaryage+exper)lm()只能得出回歸系數(shù),要想得到更為詳盡的回歸信息,應(yīng)該將結(jié)果作為數(shù)據(jù)保存或者使用“擬合模型"(fittedmodel)result<-lm(salaryage+exper+age*exper,data=byu)summary(result)myresid<-result$resid#

2、獲得殘差vcov(result)#針對于擬合后的模型計算方差協(xié)方差矩陣回歸中允許使用諸如log()和sqrt()等相對復(fù)雜的項目作為自變量,但如果設(shè)計幕,就必須先計算,然后才能回歸;或者使用I(),它可以在對公式估值前強制完成計算salary$agesq<-(salary$age)A2result<-lm(salaryage+agesq+log(exper)+age*log(exper),data=byu)或者result<-lm(salaryage+I(ageA2)+log(exper)+age*log(exper),data=byu)如果我們要進行無常數(shù)項,一般在回歸中引

3、入0result<-lm(smokes0+male+female,data=smokerdata)3.2 從回歸中提取統(tǒng)計量一些統(tǒng)計量和參數(shù)都被存儲在lm或者summary中output<-summary(result)SSR<-deviance(result)#殘差平方和;(另一種方法:RSquared<-output$r.squared)LL<-logLik(result)#對數(shù)似然統(tǒng)計量DegreesOfFreedomv-result$df#自由度Yhat<-result$fitted.values#擬合值向量Residv-result$residua

4、lss<-output$sigma#誤差標準差的估計值(假設(shè)同方差)CovMatrix<-sA2*output$cov#系數(shù)的方差一協(xié)方差矩陣(與vcov()同)3.3 異方差及相關(guān)問題3.3.1 異方差的BreuschPagan檢驗為了檢驗異方差是否存在,我們可以用lmtest包中的BreuschPagan!僉驗?;蛘呃胏ar包中的ncv.test()函數(shù)。二者工作的原理都是相同的。在回歸之后,我們可以對擬合的模型采用bptest()函數(shù)unrestrictedv-lm(z-x)bptest(unrestricted)這將得到檢驗的“學(xué)生化的"(studentized

5、)結(jié)果。如果為了保持與其他軟件結(jié)論的一致性(包括ncv.test(),我們可以設(shè)置studentize=FALSE3.3.2 異方差(自相關(guān))穩(wěn)健性協(xié)方差矩陣在存在異方差的情況下,ols的估計是無偏的,但是所得到的關(guān)于B系數(shù)方差的估計則是不正確的。為了計算異方差一致性的協(xié)方差矩陣,我們可以利用car包中的hccm()函數(shù),而不是vcov()。sandwich包中的vcovHC()命令可以實現(xiàn)同樣的功能。同時利用vcovHAC()或者NeweyWest()函數(shù)可以進行異方差和自相關(guān)穩(wěn)健性Newey-West估計。3.4 線性假設(shè)檢驗(Wald和F)car包中提供了一個函數(shù)可以自動的進行線性假設(shè)檢

6、驗。根據(jù)我們對模型的設(shè)定,它既可以用一般的方法或調(diào)整后的協(xié)方差矩陣進行F或Wald檢驗。為了進行假設(shè)檢驗,我們必須構(gòu)造一個假設(shè)矩陣和一個右手邊的向量。例如,如果我們有一個包括常數(shù)項的五個參數(shù)的模型,并且我們的零假設(shè)如下(見原文)則假設(shè)的矩陣和右手邊的向量將為如下的形式(見原文)我們可以用如下的命令加以實現(xiàn)unrestricted<-lm(yx1+x2+x3+x4)rhs<-c(0,1)hm<-rbind<-(c(1,0,0,0,0),c(0,0,1,1,0)linear.hypothesis(unrestricted,hm,rhs)注意:如果unrestricted是由

7、lm得到的,默認狀態(tài)下將會進行F檢驗。如果是由glm得到的,取而代之的將是Kai方檢驗。檢驗的類型可以通過type進行修改。同樣,如果我們想利用異方差或自相關(guān)穩(wěn)健標準誤進行檢驗,我們既可以通過設(shè)定white.adjust=TRUE來使用white標準誤,也可以利用vcov計算我們自己的協(xié)方差矩陣。例如,如果我們想使用上述的Newey-West修正協(xié)方差矩陣,我們可以進行如下的設(shè)定:linear.hypothesis(unrestricted,hm,rhs,vcov=NeweyWest(unrestricted)注意:設(shè)定white.adjust=TRUE將會通過提高white估計量的精度來修正

8、異方差;如果要使用經(jīng)典的white估計量,我們可以設(shè)定white.adjust="hc0"3.5加權(quán)和廣義最小二乘法你可以通過使用帶有權(quán)重的lm()來進行加權(quán)最小二乘result<-lm(smokes0+male+female,data=smokerdata,weights=myweights)廣義最小二乘法可以通過MASSfe中的lm.gls()命令實現(xiàn)。它將包含一個函數(shù)、權(quán)重矩陣和一個數(shù)據(jù)框(可選)。glm()命令也為使用其他高級線性回歸方法提供了渠道,詳見幫助文件。> 帶有因子(Factors)或組別(Groups)的模型在R中對于定性的因子有特定的數(shù)據(jù)類

9、型。如果回歸中的變量屬于這種情況,必要的虛擬變量將會被自動生成。例如,如果我們要進行個人電腦的使用(pc)對公司雇員數(shù)(emple)和每一種狀態(tài)的虛擬變量(其中state是兩個縮寫字母的向量),我們可以簡單的進行如下的回歸summary(lm(pcemple+state)Call:lm(formula=pcemple+state)Residuals:Min1QMedian3QMax-1.7543-0.55050.35120.42720.5904Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)5.572e-016.769e-02

10、8.232<2e-16*emple1.459e-041.083e-0513.475<2e-16*stateAL-4.774e-037.382e-02-0.0650.948stateAR2.249e-028.004e-020.2810.779stateAZ-7.023e-027.580e-02-0.9260.354stateDE1.521e-011.107e-011.3750.169.stateFL-4.573e-027.136e-02-0.6410.522stateWY1.200e-011.041e-011.1530.249Signif.codes:0'*'0.00

11、1'*'0.01'*'0.05'.0.1''1Residualstandarderror:0.4877on9948degreesoffreedomMultipleR-Squared:0.02451,AdjustedR-squared:0.01951F-statistic:4.902on51and9948DF,p-value:<2.2e-16為了將數(shù)據(jù)(序列string型或者數(shù)值型)轉(zhuǎn)變成一個因子,可以簡單的使用factor()命令。甚至可以在回歸中間使用這個命令。例如,如果我們想做同樣的回歸,但是要區(qū)分以數(shù)字編碼代表的不同區(qū)域,我們

12、可以用如下的命令myout<-lm(pcemple+factor(naics6)這一命令將naic6轉(zhuǎn)換成了因子,生成了相應(yīng)的虛擬變量,進而進行標準的回歸。4.特殊回歸固定和隨機效應(yīng)模型注意:這里所用的“固定”和“隨機”的概念與計量經(jīng)濟學(xué)家通常使用的概念相同。經(jīng)濟學(xué)中,固定和隨機效應(yīng)是用來解釋面板數(shù)據(jù)(paneldata)模型的截距項中的截面(crosssectional)差異的。令i表示截面指標(或表示數(shù)據(jù)是有組別的),t為時間指標(或為組別差異指標)。一個標準的固定效應(yīng)模型可以寫作(參照原文)本質(zhì)上說,每一個個體都有不同的非隨時間變化的截距項。通常我們感興趣的是B,而不是ui。隨機效

13、應(yīng)模型有同樣的方程,但是相對固定效應(yīng)模型而言,它有附加的限制:個體特殊的效應(yīng)與解釋變量x(it)不相關(guān),即EuiXit=0從計量經(jīng)濟學(xué)的角度講,這只是一個在固定效應(yīng)模型的基礎(chǔ)上,附有更加嚴格的限制條件的模型(它允許“效應(yīng)”與外生變量相關(guān))。固定效應(yīng)在截面數(shù)據(jù)的維數(shù)不大的情況下,進行固定效應(yīng)估計的簡單方法是在每一個個體中加入一個虛擬變量,即將截面指標視為一個因子。如果指標能夠在樣本中將個體識別出來,則有l(wèi)m(yfactor(index)+x)這個回歸可以進行固定效應(yīng)估計并能夠正確的報告B的標準誤。但不幸的是,在樣本中存在很多個體的情況下,我們不再關(guān)心固定效應(yīng)的值。因此在這種情況下,lm()的結(jié)果

14、以及u(i)較大的系數(shù)都是非常難于處理的。一個更一般的方法是通過timedemeaning(翻譯不好,請大家?guī)椭?的方法將每個變量的固定效應(yīng)剔除(所謂內(nèi)部within估計量)。則上述方程變?yōu)椋?參照原文)多數(shù)的統(tǒng)計軟件(例如stata的xtreg命令)都使用這種方法處理固定效應(yīng)模型。使用R,你可以手工對自變量和因變量進行timedemean如果d是一個包含year,firm,profit和predictor的數(shù)據(jù)框,同時我們希望timedemean每一個公司的profit和predictor,我們可以使用如下的命令:g<-dfor(iinunique(d$firm)+timemean&l

15、t;-mean(dd$firm=i,)+g$profitd$firm=i<-d$profitd$firm=i-timemean"profit"+g$predictord$firm=i<-d$predictord$firm=i-timemean"predictor"+output<-lm(profitpredictor,data=g)要注意的是,回歸中所報告的標準誤偏低。lm()報告的方差使用的公式為而真正的固定效應(yīng)回歸中的標準誤使用的公式為(參照原文)對于T較小的樣本,二者之間存在較為顯著的差異。另一種并不常用的方法是用一階差分,公式為

16、(參照原文)其手工計算方法可以參照上述的withinestimator隨機效應(yīng)模型nlme包包括了在線性或非線性數(shù)據(jù)框中進行隨機效應(yīng)回歸的方法(而不是固定效應(yīng),本文只涉及對“固定效應(yīng)”項的統(tǒng)計解釋)。假設(shè)存在如下的以sic3編碼為隨機效應(yīng)的線性模型我們可以用下列命令進行擬合lme(ldsallemp+ldnpt,random=1|sic3)一般而言,在隨機參數(shù)模型中,隨機效應(yīng)被置于豎線之后。波浪線和豎線之間的1表示隨機效應(yīng)是一個截距項。如果隨機效應(yīng)是一個截距項和一個外生變量,我們應(yīng)當將該變量與1放在一起。例如:lme(ldsallemp+ldnpt,random=1+lemp|sic3)對應(yīng)于

17、對于非線性隨機效應(yīng)模型而言,只是將lme()替換為nlme()就可以了。定性相應(yīng)模型(QualitativeResponse)Logit/Probit有很多種方法可以在R中實現(xiàn)logit和probit回歸。最簡單的方法是使用glm()命令及相應(yīng)的選項。h<-glm(cy,family=binomial(link="logit")對于probit回歸,將logit替換為probit即可。glm()函數(shù)的輸出結(jié)果與lm()的類似,因此可以使用summary。命令加以分析。為了從中提取出對數(shù)似然統(tǒng)計量,可以使用logLik()命令>logLik(h)'logL

18、ik.'-337.2659(df=1)多項式Logit(MultinormialLogit)在nnet包中有一個multinom()函數(shù)可以用來進行多項式Logit估計。為了使用這一函數(shù),首先應(yīng)該將因變量轉(zhuǎn)化成一個因子向量(包括所有情況),并且使用諸如正態(tài)回歸這樣的語法(syntax)。如果我們的因子以虛擬變量的形式被儲存,則我們可以利用十進制數(shù)字的特征為所有的組合賦予唯一的因子。假如我們的因子變量是pc,inetacc和iapp,那么>g<-pc*1+inetacc*10+iapp*100>multinom(factor(g)pc.subsidy+inet.subs

19、idy+iapp.subsidy+emple+msamissing)這樣,我們利用所有因子的組合得到了一個多項式Logit。多項式Probit模型的一個特征就是常被illconditioned(如何譯?)。一個解決的方法是使用MNFfe中的mnp()命令進行馬爾可夫鏈蒙特卡羅模擬。OrderedLogit/ProbitMASS包中有一個polr()函數(shù)可以進行orderedlogit或probit回歸。如果Sat表示一個順序(ordered)的因子向量,則house.plr<-polr(SatInfl+Type+Cont,method="probit")Tobit和階

20、段(Censored)回歸我們用survival包來估計存在截斷數(shù)據(jù)變量的模型。使用的函數(shù)是survreg(),其中將因變量作為一個Surv對象。假設(shè)我們要進行y對x和z的回歸,但是大量的y數(shù)據(jù)是左側(cè)截斷的,并用0將其替換。result<-survreg(Surv(y,y>0,type='left')x+z,dist='gaussian')第二點需要注意的是無論觀測值是否是截斷的,都可以使用Surv()函數(shù)(1表示是可以被觀測的;0表示數(shù)據(jù)是截斷的)。第三點需要說明的是,數(shù)據(jù)在哪一側(cè)被截斷。既然本例中數(shù)據(jù)在分布的低尾(lowertail)處被截斷,我

21、們使用left。dist選項又t于survreg是必需的,這樣才能得到一個經(jīng)典的Tobit模型。分位數(shù)回歸最小二乘回歸方法可以估測因變量對自變量的條件期望。擬合值即是條件均值的估計。如果我們不想得到條件均值,而想估計預(yù)期條件中位數(shù)或其他分位數(shù)的話,我們可以使用quantreg包中的rq()命令。語法與lm()基本相同,除了我們要使用一個介于0和1之間的分位數(shù)tau。默認的情況下,tau=.5,即為中位數(shù),另一個名字是最小絕對偏差回歸(leastabsolutedeviationregression)4.5穩(wěn)健性回歸M估計量對于一些數(shù)據(jù)集,奇異值對最小二乘回歸線的影響遠遠超出了我們的預(yù)想。一個解

22、決的辦法是使用包括殘差平方和(對應(yīng)于最小化L2)在內(nèi)的一些方法求極小值并以此作為目標方程。通常的選擇是使用絕對離差和(L1)和Huber法一一一種將L1和L2混合的方法。R使用MASSfe中的rlm()進行穩(wěn)健性回歸。語法與lm()相同除了它允許選擇最小化作為目標方程。進行這種選擇可以使用參數(shù)psio可能的選項包括psi.huber,psi.hampel,psi.bisquare。為了進行psi函數(shù)的定制,我們寫了一個函數(shù),如果deriv=0,函數(shù)返回巾(x)/x;如果deriv=1,返回也(x)/x°Thisfunctionthanthenbepassedtorlm()usingt

23、hepsiparameter.#不清楚函數(shù)內(nèi)容及語意。非線性最小二乘有些時候,經(jīng)濟中的模型并不是線性的。R可以進行如下形式的廣義最小二注意,殘差項必須是附加在函數(shù)形式上的。如果不是,則必須進行相應(yīng)的變換以達到這種形式。R中進行非線性最小二乘的函數(shù)是nls(),其語法與lm()相同c考慮如下的非線性例子:nls()用來估計第二個方程的第一個部分。方程的全部內(nèi)容都需要被指定,包括參數(shù)。R要求指定待估參數(shù)的初始值。result<-nls(log(Y)-log(1+exp(a*X1+b*X2),start=list(a=1,b=1),data=mydata)a和b的估計值被存放于nls的對象中,稱作result。估計結(jié)果可以用summary()命令進行瀏覽。在高版本的R中,nls()命令是基本包中的一部分,而在低版本中,則必須加載nls包。單一結(jié)構(gòu)方程的兩階段最小二乘為了實現(xiàn)單一方程的兩階段最小二乘,最簡單的方法是使用sem包中的tsls()命令。如果我們想考察在控制了蠟姻狀況的情況下,教育對工資的影響,但是考慮到educ可能是內(nèi)生的,則我們可能會使用motheduc和fatheduc作為工具變量進行回歸library(sem)outputof2sl

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負責。
  • 6. 下載文件中如有侵權(quán)或不適當內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論