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

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