在R軟件中實(shí)現(xiàn)回歸估計(jì)

1、3.1 最小二乘法有三種方式可以實(shí)現(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)#針對(duì)于擬合后的模型計(jì)算方差協(xié)方差矩陣回歸中允許使用諸如log()和sqrt()等相對(duì)復(fù)雜的項(xiàng)目作為自變量，但如果設(shè)計(jì)幕，就必須先計(jì)算，然后才能回歸；或者使用I(),它可以在對(duì)公式估值前強(qiáng)制完成計(jì)算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)如果我們要進(jìn)行無常數(shù)項(xiàng)，一般在回歸中引

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

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

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

6、驗(yàn)。根據(jù)我們對(duì)模型的設(shè)定，它既可以用一般的方法或調(diào)整后的協(xié)方差矩陣進(jìn)行F或Wald檢驗(yàn)。為了進(jìn)行假設(shè)檢驗(yàn)，我們必須構(gòu)造一個(gè)假設(shè)矩陣和一個(gè)右手邊的向量。例如，如果我們有一個(gè)包括常數(shù)項(xiàng)的五個(gè)參數(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得到的，默認(rèn)狀態(tài)下將會(huì)進(jìn)行F檢驗(yàn)。如果是由glm得到的，取而代之的將是Kai方檢驗(yàn)。檢驗(yàn)的類型可以通過type進(jìn)行修改。同樣，如果我們想利用異方差或自相關(guān)穩(wěn)健標(biāo)準(zhǔn)誤進(jìn)行檢驗(yàn)，我們既可以通過設(shè)定white.adjust=TRUE來使用white標(biāo)準(zhǔn)誤，也可以利用vcov計(jì)算我們自己的協(xié)方差矩陣。例如，如果我們想使用上述的Newey-West修正協(xié)方差矩陣，我們可以進(jìn)行如下的設(shè)定：linear.hypothesis(unrestricted,hm,rhs,vcov=NeweyWest(unrestricted)注意：設(shè)定white.adjust=TRUE將會(huì)通過提高white估計(jì)量的精度來修正

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

9、型。如果回歸中的變量屬于這種情況，必要的虛擬變量將會(huì)被自動(dòng)生成。例如，如果我們要進(jìn)行個(gè)人電腦的使用(pc)對(duì)公司雇員數(shù)(emple)和每一種狀態(tài)的虛擬變量(其中state是兩個(gè)縮寫字母的向量)，我們可以簡單的進(jìn)行如下的回歸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)變成一個(gè)因子，可以簡單的使用factor()命令。甚至可以在回歸中間使用這個(gè)命令。例如，如果我們想做同樣的回歸，但是要區(qū)分以數(shù)字編碼代表的不同區(qū)域，我們

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

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

14、以及u(i)較大的系數(shù)都是非常難于處理的。一個(gè)更一般的方法是通過timedemeaning(翻譯不好，請(qǐng)大家?guī)椭?的方法將每個(gè)變量的固定效應(yīng)剔除(所謂內(nèi)部within估計(jì)量)。則上述方程變?yōu)椋?參照原文)多數(shù)的統(tǒng)計(jì)軟件(例如stata的xtreg命令)都使用這種方法處理固定效應(yīng)模型。使用R,你可以手工對(duì)自變量和因變量進(jìn)行timedemean如果d是一個(gè)包含year,firm,profit和predictor的數(shù)據(jù)框，同時(shí)我們希望timedemean每一個(gè)公司的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)要注意的是，回歸中所報(bào)告的標(biāo)準(zhǔn)誤偏低。lm()報(bào)告的方差使用的公式為而真正的固定效應(yīng)回歸中的標(biāo)準(zhǔn)誤使用的公式為(參照原文)對(duì)于T較小的樣本，二者之間存在較為顯著的差異。另一種并不常用的方法是用一階差分，公式為

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

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

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

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

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

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

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

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