在R软件中实现回归估计

1、3.1 最小二乘法有三种方式可以实现最小二乘法的简单线性回归，假设数据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()只能得出回归系数，要想得到更为详尽的回归信息，应该将结果作为数据保存或者使用“拟合模型"(fittedmodel)result<-lm(salaryage+exper+age*exper,data=byu)summary(result)myresid<-result$resid#

2、获得残差vcov(result)#针对于拟合后的模型计算方差协方差矩阵回归中允许使用诸如log()和sqrt()等相对复杂的项目作为自变量，但如果设计幕，就必须先计算，然后才能回归；或者使用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)如果我们要进行无常数项，一般在回归中引

3、入0result<-lm(smokes0+male+female,data=smokerdata)3.2 从回归中提取统计量一些统计量和参数都被存储在lm或者summary中output<-summary(result)SSR<-deviance(result)#残差平方和；(另一种方法：RSquared<-output$r.squared)LL<-logLik(result)#对数似然统计量DegreesOfFreedomv-result$df#自由度Yhat<-result$fitted.values#拟合值向量Residv-result$residua

4、lss<-output$sigma#误差标准差的估计值(假设同方差)CovMatrix<-sA2*output$cov#系数的方差一协方差矩阵(与vcov()同)3.3 异方差及相关问题3.3.1 异方差的BreuschPagan检验为了检验异方差是否存在，我们可以用lmtest包中的BreuschPagan!佥验。或者利用car包中的ncv.test()函数。二者工作的原理都是相同的。在回归之后，我们可以对拟合的模型采用bptest()函数unrestrictedv-lm(z-x)bptest(unrestricted)这将得到检验的“学生化的"(studentized

5、)结果。如果为了保持与其他软件结论的一致性(包括ncv.test(),我们可以设置studentize=FALSE3.3.2 异方差(自相关)稳健性协方差矩阵在存在异方差的情况下，ols的估计是无偏的，但是所得到的关于B系数方差的估计则是不正确的。为了计算异方差一致性的协方差矩阵，我们可以利用car包中的hccm()函数，而不是vcov()。sandwich包中的vcovHC()命令可以实现同样的功能。同时利用vcovHAC()或者NeweyWest()函数可以进行异方差和自相关稳健性Newey-West估计。3.4 线性假设检验(Wald和F)car包中提供了一个函数可以自动的进行线性假设检

6、验。根据我们对模型的设定，它既可以用一般的方法或调整后的协方差矩阵进行F或Wald检验。为了进行假设检验，我们必须构造一个假设矩阵和一个右手边的向量。例如，如果我们有一个包括常数项的五个参数的模型，并且我们的零假设如下(见原文)则假设的矩阵和右手边的向量将为如下的形式(见原文)我们可以用如下的命令加以实现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得到的，默认状态下将会进行F检验。如果是由glm得到的，取而代之的将是Kai方检验。检验的类型可以通过type进行修改。同样，如果我们想利用异方差或自相关稳健标准误进行检验，我们既可以通过设定white.adjust=TRUE来使用white标准误，也可以利用vcov计算我们自己的协方差矩阵。例如，如果我们想使用上述的Newey-West修正协方差矩阵，我们可以进行如下的设定：linear.hypothesis(unrestricted,hm,rhs,vcov=NeweyWest(unrestricted)注意：设定white.adjust=TRUE将会通过提高white估计量的精度来修正

8、异方差；如果要使用经典的white估计量，我们可以设定white.adjust="hc0"3.5加权和广义最小二乘法你可以通过使用带有权重的lm()来进行加权最小二乘result<-lm(smokes0+male+female,data=smokerdata,weights=myweights)广义最小二乘法可以通过MASSfe中的lm.gls()命令实现。它将包含一个函数、权重矩阵和一个数据框(可选)。glm()命令也为使用其他高级线性回归方法提供了渠道，详见帮助文件。> 带有因子(Factors)或组别(Groups)的模型在R中对于定性的因子有特定的数据类

9、型。如果回归中的变量属于这种情况，必要的虚拟变量将会被自动生成。例如，如果我们要进行个人电脑的使用(pc)对公司雇员数(emple)和每一种状态的虚拟变量(其中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为了将数据(序列string型或者数值型)转变成一个因子，可以简单的使用factor()命令。甚至可以在回归中间使用这个命令。例如，如果我们想做同样的回归，但是要区分以数字编码代表的不同区域，我们

12、可以用如下的命令myout<-lm(pcemple+factor(naics6)这一命令将naic6转换成了因子，生成了相应的虚拟变量，进而进行标准的回归。4.特殊回归固定和随机效应模型注意：这里所用的“固定”和“随机”的概念与计量经济学家通常使用的概念相同。经济学中，固定和随机效应是用来解释面板数据(paneldata)模型的截距项中的截面(crosssectional)差异的。令i表示截面指标(或表示数据是有组别的)，t为时间指标(或为组别差异指标)。一个标准的固定效应模型可以写作(参照原文)本质上说，每一个个体都有不同的非随时间变化的截距项。通常我们感兴趣的是B,而不是ui。随机效

13、应模型有同样的方程，但是相对固定效应模型而言，它有附加的限制：个体特殊的效应与解释变量x(it)不相关，即EuiXit=0从计量经济学的角度讲，这只是一个在固定效应模型的基础上，附有更加严格的限制条件的模型(它允许“效应”与外生变量相关)。固定效应在截面数据的维数不大的情况下，进行固定效应估计的简单方法是在每一个个体中加入一个虚拟变量，即将截面指标视为一个因子。如果指标能够在样本中将个体识别出来，则有lm(yfactor(index)+x)这个回归可以进行固定效应估计并能够正确的报告B的标准误。但不幸的是,在样本中存在很多个体的情况下，我们不再关心固定效应的值。因此在这种情况下，lm()的结果

14、以及u(i)较大的系数都是非常难于处理的。一个更一般的方法是通过timedemeaning(翻译不好，请大家帮助)的方法将每个变量的固定效应剔除(所谓内部within估计量)。则上述方程变为：(参照原文)多数的统计软件(例如stata的xtreg命令)都使用这种方法处理固定效应模型。使用R,你可以手工对自变量和因变量进行timedemean如果d是一个包含year,firm,profit和predictor的数据框，同时我们希望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()报告的方差使用的公式为而真正的固定效应回归中的标准误使用的公式为(参照原文)对于T较小的样本，二者之间存在较为显著的差异。另一种并不常用的方法是用一阶差分，公式为

16、(参照原文)其手工计算方法可以参照上述的withinestimator随机效应模型nlme包包括了在线性或非线性数据框中进行随机效应回归的方法(而不是固定效应，本文只涉及对“固定效应”项的统计解释)。假设存在如下的以sic3编码为随机效应的线性模型我们可以用下列命令进行拟合lme(ldsallemp+ldnpt,random=1|sic3)一般而言，在随机参数模型中，随机效应被置于竖线之后。波浪线和竖线之间的1表示随机效应是一个截距项。如果随机效应是一个截距项和一个外生变量，我们应当将该变量与1放在一起。例如：lme(ldsallemp+ldnpt,random=1+lemp|sic3)对应于

17、对于非线性随机效应模型而言，只是将lme()替换为nlme()就可以了。定性相应模型(QualitativeResponse)Logit/Probit有很多种方法可以在R中实现logit和probit回归。最简单的方法是使用glm()命令及相应的选项。h<-glm(cy,family=binomial(link="logit")对于probit回归，将logit替换为probit即可。glm()函数的输出结果与lm()的类似，因此可以使用summary。命令加以分析。为了从中提取出对数似然统计量，可以使用logLik()命令>logLik(h)'logL

18、ik.'-337.2659(df=1)多项式Logit(MultinormialLogit)在nnet包中有一个multinom()函数可以用来进行多项式Logit估计。为了使用这一函数，首先应该将因变量转化成一个因子向量(包括所有情况)，并且使用诸如正态回归这样的语法(syntax)。如果我们的因子以虚拟变量的形式被储存，则我们可以利用十进制数字的特征为所有的组合赋予唯一的因子。假如我们的因子变量是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()函数可以进行orderedlogit或probit回归。如果Sat表示一个顺序(ordered)的因子向量，则house.plr<-polr(SatInfl+Type+Cont,method="probit")Tobit和阶

20、段(Censored)回归我们用survival包来估计存在截断数据变量的模型。使用的函数是survreg(),其中将因变量作为一个Surv对象。假设我们要进行y对x和z的回归，但是大量的y数据是左侧截断的，并用0将其替换。result<-survreg(Surv(y,y>0,type='left')x+z,dist='gaussian')第二点需要注意的是无论观测值是否是截断的，都可以使用Surv()函数(1表示是可以被观测的；0表示数据是截断的)。第三点需要说明的是，数据在哪一侧被截断。既然本例中数据在分布的低尾(lowertail)处被截断，我

21、们使用left。dist选项又t于survreg是必需的，这样才能得到一个经典的Tobit模型。分位数回归最小二乘回归方法可以估测因变量对自变量的条件期望。拟合值即是条件均值的估计。如果我们不想得到条件均值，而想估计预期条件中位数或其他分位数的话，我们可以使用quantreg包中的rq()命令。语法与lm()基本相同，除了我们要使用一个介于0和1之间的分位数tau。默认的情况下，tau=.5,即为中位数，另一个名字是最小绝对偏差回归(leastabsolutedeviationregression)4.5稳健性回归M估计量对于一些数据集，奇异值对最小二乘回归线的影响远远超出了我们的预想。一个解

22、决的办法是使用包括残差平方和(对应于最小化L2)在内的一些方法求极小值并以此作为目标方程。通常的选择是使用绝对离差和(L1)和Huber法一一一种将L1和L2混合的方法。R使用MASSfe中的rlm()进行稳健性回归。语法与lm()相同除了它允许选择最小化作为目标方程。进行这种选择可以使用参数psio可能的选项包括psi.huber,psi.hampel,psi.bisquare。为了进行psi函数的定制，我们写了一个函数，如果deriv=0,函数返回巾(x)/x；如果deriv=1,返回也(x)/x°Thisfunctionthanthenbepassedtorlm()usingt

23、hepsiparameter.#不清楚函数内容及语意。非线性最小二乘有些时候，经济中的模型并不是线性的。R可以进行如下形式的广义最小二注意，残差项必须是附加在函数形式上的。如果不是，则必须进行相应的变换以达到这种形式。R中进行非线性最小二乘的函数是nls(),其语法与lm()相同c考虑如下的非线性例子：nls()用来估计第二个方程的第一个部分。方程的全部内容都需要被指定，包括参数。R要求指定待估参数的初始值。result<-nls(log(Y)-log(1+exp(a*X1+b*X2),start=list(a=1,b=1),data=mydata)a和b的估计值被存放于nls的对象中，称作result。估计结果可以用summary()命令进行浏览。在高版本的R中，nls()命令是基本包中的一部分，而在低版本中，则必须加载nls包。单一结构方程的两阶段最小二乘为了实现单一方程的两阶段最小二乘，最简单的方法是使用sem包中的tsls()命令。如果我们想考察在控制了蜡姻状况的情况下，教育对工资的影响，但是考虑到educ可能是内生的，则我们可能会使用motheduc和fatheduc作为工具变量进行回归library(sem)outputof2sl