文章出处：http://cos.name/cn/topic/3272/page/2
3、Cross Sectional Regression
3.1 最小二乘法
有三种方式可以实现最小二乘法的简单线性回归，假设数据byu
（1）lm(byu$salary ~ byu$age + byu$exper)
（2）lm (salary ~ age + exper, data= byu)
（3）attach(byu)
lm(salary~age+exper)
lm()只能得出回归系数，要想得到更为详尽的回归信息，应该将结果作为数据保存或者使用“拟合模型”（fitted model）
result<-lm(salary~age+ exper + age*exper, data=byu)
summary(result)
myresid<-result$resid #获得残差
vcov(result) ＃针对于拟合后的模型计算方差－协方差矩阵
回归中允许使用诸如log（）和sqrt（）等相对复杂的项目作为自变量，但如果设计幂，就必须先计算，然后才能回归；或者使用I（），它可以在对公式估值前强制完成计算
salary$agesq <- (salary$age)^2
result <- lm(salary ~ age + agesq + log(exper) + age*log(exper),data=byu)
或者
result <- lm(salary ~ age + I(age^2) + log(exper) + age*log(exper),data=byu)
如果我们要进行无常数项，一般在回归中引入0
result <- lm(smokes ~ 0 + male + female ,data=smokerdata)
3.2 从回归中提取统计量
一些统计量和参数都被存储在lm或者summary中
output <- summary (result)
SSR<- deviance（result）＃残差平方和；（另一种方法：RSquared <- output$r.squared)
LL<- logLik(result) #对数似然统计量
DegreesOfFreedom<- result$df #自由度
Yhat<- result$fitted.values #拟合值向量
Resid<- result$residuals
s<- output$sigma #误差标准差的估计值（假设同方差）
CovMatrix <- s^2*output$cov #系数的方差－协方差矩阵（与vcov()同）
3.3 异方差及相关问题
3.3.1 异方差的Breusch－Pagan检验
为了检验异方差是否存在，我们可以用lmtest包中的Breusch－Pagan检验。或者利用car包中的ncv.test()函数。二者工作的原理都是相同的。在回归之后，我们可以对拟合的模型采用bptest（）函数
unrestricted<- lm(z~x)
bptest(unrestricted)
这将得到检验的“学生化的”（studentized）结果。如果为了保持与其他软件结论的一致性（包括ncv.test())，我们可以设置studentize=FALSE
3.3.2 异方差（自相关）稳健性协方差矩阵
在存在异方差的情况下，ols的估计是无偏的，但是所得到的关于β系数方差的估计则是不正确的。为了计算异方差一致性的协方差矩阵，我们可以利用car包中的hccm（）函数，而不是vcov()。
sandwich包中的vcovHC()命令可以实现同样的功能。同时利用vcovHAC()或者NeweyWest（）函数可以进行异方差和自相关稳健性Newey—West估计。
####################
还有一章节，special regression，这个可能是计量中用的比较多的回归，比如logit，probit，还有tobit回归模型，以及M估计，稳健方法等等。
另外关于一般的lm方法由于用到的是ols或者称为lse方法，所以不存在log likelihood的说法，所以求出来的Loglik值可能没有意义的说。
另外用splus算过probit和logit回归，然后使用loglik函数得到的值竟然是错误的，手动演算一下跟书上一样，参考武德里奇的《计量》一书中离散选择模型第一五章的第一个例子。
####################
3.4 线性假设检验（Wald和F）
car包中提供了一个函数可以自动的进行线性假设检验。根据我们对模型的设定，它既可以用一般的方法或调整后的协方差矩阵进行F或Wald检验。为了进行假设检验，我们必须构造一个假设矩阵和一个右手边的向量。例如，如果我们有一个包括常数项的五个参数的模型，并且我们的零假设如下（见原文）
则假设的矩阵和右手边的向量将为如下的形式（见原文）
我们可以用如下的命令加以实现
unrestricted<-lm(y~x1+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是由lm得到的，默认状态下将会进行F检验。如果是由glm得到的，取而代之的将是Kai方检验。检验的类型可以通过type进行修改。
同样，如果我们想利用异方差或自相关稳健标准误进行检验，我们既可以通过设定white.adjust=TRUE来使用white标准误，也可以利用vcov计算我们自己的协方差矩阵。例如，如果我们想使用上述的Newey-West修正协方差矩阵，我们可以进行如下的设定：
linear.hypothesis(unrestricted, hm, rhs, vcov=NeweyWest(unrestricted))
注意：设定white.adjust=TRUE将会通过提高white估计量的精度来修正异方差；如果要使用经典的white估计量，我们可以设定white.adjust="hc0"
3.5 加权和广义最小二乘法
你可以通过使用带有权重的lm()来进行加权最小二乘
result<-lm(smokes~0+male+female, data=smokerdata, weights=myweights)
广义最小二乘法可以通过MASS包中的lm.gls()命令实现。它将包含一个函数、权重矩阵和一个数据框（可选）。
glm（）命令也为使用其他高级线性回归方法提供了渠道，详见帮助文件。
3.6 带有因子（Factors）或组别（Groups）的模型
在R中对于定性的因子有特定的数据类型。如果回归中的变量属于这种情况，必要的虚拟变量将会被自动生成。例如，如果我们要进行个人电脑的使用（pc）对公司雇员数（emple）和每一种状态的虚拟变量（其中state是两个缩写字母的向量），我们可以简单的进行如下的回归
summary(lm(pc~emple+state)
Call:
lm(formula = pc ~ emple + state)
Residuals:
Min 1Q Median 3Q Max
-1.7543 -0.5505 0.3512 0.4272 0.5904
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.572e-01 6.769e-02 8.232 <2e-16 ***
emple 1.459e-04 1.083e-05 13.475 <2e-16 ***
stateAL -4.774e-03 7.382e-02 -0.065 0.948
stateAR 2.249e-02 8.004e-02 0.281 0.779
stateAZ -7.023e-02 7.580e-02 -0.926 0.354
stateDE 1.521e-01 1.107e-01 1.375 0.169
...
stateFL -4.573e-02 7.136e-02 -0.641 0.522
stateWY 1.200e-01 1.041e-01 1.153 0.249
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4877 on 9948 degrees of freedom
Multiple R-Squared: 0.02451, Adjusted R-squared: 0.01951
F-statistic: 4.902 on 51 and 9948 DF, p-value: < 2.2e-16
为了将数据（序列string型或者数值型）转变成一个因子，可以简单的使用factor()命令。甚至可以在回归中间使用这个命令。例如，如果我们想做同样的回归，但是要区分以数字编码代表的不同区域，我们可以用如下的命令
myout<-lm(pc~emple+factor(naics6))
这一命令将naic6转换成了因子，生成了相应的虚拟变量，进而进行标准的回归。
####################
4. 特殊回归
4.1 固定和随机效应模型
注意：这里所用的“固定”和“随机”的概念与计量经济学家通常使用的概念相同。
经济学中，固定和随机效应是用来解释面板数据（panel data）模型的截距项中的截面（cross sectional）差异的。令i表示截面指标（或表示数据是有组别的），t为时间指标（或为组别差异指标）。一个标准的固定效应模型可以写作（参照原文）
本质上说，每一个个体都有不同的非随时间变化的截距项。通常我们感兴趣的是β，而不是ui。随机效应模型有同样的方程，但是相对固定效应模型而言，它有附加的限制：个体特殊的效应与解释变量x（it）不相关，即E[uiXit] = 0。从计量经济学的角度讲，这只是一个在固定效应模型的基础上，附有更加严格的限制条件的模型（它允许“效应”与外生变量相关）。
4.1.1 固定效应
在截面数据的维数不大的情况下，进行固定效应估计的简单方法是在每一个个体中加入一个虚拟变量，即将截面指标视为一个因子。如果指标能够在样本中将个体识别出来，则有
lm(y~factor(index)+x)
这个回归可以进行固定效应估计并能够正确的报告β的标准误。但不幸的是，在样本中存在很多个体的情况下，我们不再关心固定效应的值。因此在这种情况下，lm（）的结果以及u（i）较大的系数都是非常难于处理的。
一个更一般的方法是通过time demeaning（翻译不好，请大家帮助）的方法将每个变量的固定效应剔除（所谓内部within估计量）。则上述方程变为：（参照原文）
多数的统计软件（例如stata的xtreg命令）都使用这种方法处理固定效应模型。使用R，你可以手工对自变量和因变量进行time demean。如果d是一个包含 year, firm, profit和predictor的数据框，同时我们希望time demean每一个公司的profit和predictor，我们可以使用如下的命令：
> g<-d
> for (i in unique(d$firm)){
+ timemean<-mean(d[d$firm==i,])
+ g$profit[d$firm==i]<-d$profit[d$firm==i]-timemean["profit"]
+ g$predictor[d$firm==i]<-d$predictor[d$firm==i]-timemean["predictor"]
+ }
> output<-lm(profit~predictor,data=g)
要注意的是，回归中所报告的标准误偏低。lm()报告的方差使用的公式为
而真正的固定效应回归中的标准误使用的公式为（参照原文）
对于T较小的样本，二者之间存在较为显著的差异。
另一种并不常用的方法是用一阶差分，公式为（参照原文）
其手工计算方法可以参照上述的within estimator
4.1.2 随机效应模型
nlme包包括了在线性或非线性数据框中进行随机效应回归的方法（而不是固定效应，本文只涉及对“固定效应”项的统计解释）。
假设存在如下的以sic3编码为随机效应的线性模型
我们可以用下列命令进行拟合
lme(ldsal~lemp+ldnpt,random=~1|sic3)
一般而言，在随机参数模型中，随机效应被置于竖线之后。波浪线和竖线之间的1表示随机效应是一个截距项。如果随机效应是一个截距项和一个外生变量，我们应当将该变量与1放在一起。例如：
lme(ldsal~lemp+ldnpt,random=~1+lemp|sic3)
对应于
对于非线性随机效应模型而言，只是将lme()替换为nlme()就可以了。
#####################
4.2 定性相应模型（Qualitative Response)
4.2.1 Logit/Probit
有很多种方法可以在R中实现logit和probit回归。最简单的方法是使用glm()命令及相应的选项。
h<-glm(c~y,family=binomial(link="logit"))
对于probit回归，将logit替换为probit即可。glm()函数的输出结果与lm()的类似，因此可以使用summary()命令加以分析。为了从中提取出对数似然统计量，可以使用logLik()命令
> logLik(h)
‘log Lik.’ -337.2659 (df=1)
4.2.2 多项式Logit （Multinormial Logit）
在nnet包中有一个multinom()函数可以用来进行多项式Logit估计。为了使用这一函数，首先应该将因变量转化成一个因子向量（包括所有情况），并且使用诸如正态回归这样的语法（syntax）。如果我们的因子以虚拟变量的形式被储存，则我们可以利用十进制数字的特征为所有的组合赋予唯一的因子。假如我们的因子变量是pc,inetacc和iapp，那么
> g <- pc*1 + inetacc*10 + iapp*100
> multinom(factor(g)~pc.subsidy+inet.subsidy+iapp.subsidy+emple+msamissing)
这样，我们利用所有因子的组合得到了一个多项式Logit。
多项式Probit模型的一个特征就是常被ill conditioned（如何译？）。一个解决的方法是使用MNP包中的mnp（）命令进行马尔可夫链蒙特卡罗模拟。
4.2.3 Ordered Logit/Probit
MASS包中有一个polr()函数可以进行ordered logit或probit回归。如果Sat表示一个顺序（ordered）的因子向量，则
house.plr <- polr(Sat ~ Infl + Type + Cont, method="probit")
4.3 Tobit和阶段（Censored）回归
我们用survival包来估计存在截断数据变量的模型。使用的函数是survreg(),其中将因变量作为一个Surv对象。假设我们要进行y对x和z的回归，但是大量的y数据是左侧截断的，并用0将其替换。
result <- survreg(Surv(y,y>0,type='left') ~ x + z, dist='gaussian')
第二点需要注意的是无论观测值是否是截断的，都可以使用Surv()函数（1表示是可以被观测的；0表示数据是截断的）。第三点需要说明的是，数据在哪一侧被截断。既然本例中数据在分布的低尾（lower tail）处被截断，我们使用left。dist选项对于survreg是必需的，这样才能得到一个经典的Tobit模型。
4.4 分位数回归
最小二乘回归方法可以估测因变量对自变量的条件期望。拟合值即是条件均值的估计。如果我们不想得到条件均值，而想估计预期条件中位数或其他分位数的话，我们可以使用quantreg包中的rq()命令。语法与lm()基本相同，除了我们要使用一个介于0和1之间的分位数tau。默认的情况下，tau＝.5，即为中位数，另一个名字是最小绝对偏差回归（least absolute deviation regression)
4.5 稳健性回归——M估计量
对于一些数据集，奇异值对最小二乘回归线的影响远远超出了我们的预想。一个解决的办法是使用包括残差平方和（对应于最小化L2）在内的一些方法求极小值并以此作为目标方程。通常的选择是使用绝对离差和（L1）和Huber法——一种将L1和L2混合的方法。R使用MASS包中的rlm()进行稳健性回归。语法与lm()相同除了它允许选择最小化作为目标方程。进行这种选择可以使用参数psi。可能的选项包括psi.huber, psi.hampel, psi.bisquare。
为了进行psi函数的定制，我们写了一个函数，如果deriv=0，函数返回ψ(x)/x；如果deriv=1，返回ψ′(x)/x。This function than then be passed to rlm() using the psi parameter.＃不清楚函数内容及语意。
4.6 非线性最小二乘
有些时候，经济中的模型并不是线性的。R可以进行如下形式的广义最小二乘
注意，残差项必须是附加在函数形式上的。如果不是，则必须进行相应的变换以达到这种形式。R中进行非线性最小二乘的函数是nls()，其语法与lm()相同。考虑如下的非线性例子：
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包。
4.7单一结构方程的两阶段最小二乘
为了实现单一方程的两阶段最小二乘，最简单的方法是使用sem包中的tsls()命令。如果我们想考察在控制了婚姻状况的情况下，教育对工资的影响，但是考虑到educ可能是内生的，则我们可能会使用motheduc和fatheduc作为工具变量进行回归
library(sem)
outputof2sls<-tsls(lwage~educ+married,~married+motheduc+fatheduc)
第一点需要说明的是我们进行估计的是一个结构性方程；第二，波浪线后面需要附加结构性方程中的工具变量和外生变量，这些条件的满足需要下面的2SLS估计量满足
输出结果可以用summary()函数和其他ols函数进行分析。注意，既然这个命令输出的是两阶段最小二乘的对象，包括标准误在内的描述性统计量都是正确的。如果我们用实际的两阶段方法进行估计，则输出的标准误是存伪的。
4.8 方程组
用于处理方程组（包括工具变量、两阶段最小二乘、似不相关回归[seemingly unrelated regression]
和变量)的命令包含在systemfit包中。注意，在R中，一个方程（包括标题）就恰恰是一种数据类型。因此，我们可以通过通常的赋值运算来构造方程模型列表以及相应的标签列表。
> demand <- q ~ p + d
> supply <- q ~ p + f + a
> system <- list(demand,supply)
> labels <- list("demand","supply")
4.8.1 似无关回归
一旦我们建立了系统和标签，我们就可以使用systemfit()函数结合SUR选项来识别似无关回归。
> resultsur <- systemfit("SUR",system,labels)
4.8.2 方程组的两阶段最小二乘
可以使用工具变量来对上述的方程组进行两阶段最小二乘。我们创建一个模型对象（不包含左侧）来表示我们要使用的工具变量和2SLS选项
> inst <- ~ d + f + a
> result2sls <- systemfit("2SLS",system,labels,inst)
与之相类似，我们还会遇到三阶段最小二乘、加权两阶段最小二乘和相应的其他模型。
####################
about fixed- and random-effect model, can refer to plm package.
about multinomial logit model, can refer to mlogit package.
about tobit model, refer to sampleSelection package.
####################
颇值得一提的是，回归中最紧要看residuals plot，
example(plot.lm)
####################
3、Cross Sectional Regression
3.1 最小二乘法
有三种方式可以实现最小二乘法的简单线性回归，假设数据byu
（1）lm(byu$salary ~ byu$age + byu$exper)
（2）lm (salary ~ age + exper, data= byu)
（3）attach(byu)
lm(salary~age+exper)
lm()只能得出回归系数，要想得到更为详尽的回归信息，应该将结果作为数据保存或者使用“拟合模型”（fitted model）
result<-lm(salary~age+ exper + age*exper, data=byu)
summary(result)
myresid<-result$resid #获得残差
vcov(result) ＃针对于拟合后的模型计算方差－协方差矩阵
回归中允许使用诸如log（）和sqrt（）等相对复杂的项目作为自变量，但如果设计幂，就必须先计算，然后才能回归；或者使用I（），它可以在对公式估值前强制完成计算
salary$agesq <- (salary$age)^2
result <- lm(salary ~ age + agesq + log(exper) + age*log(exper),data=byu)
或者
result <- lm(salary ~ age + I(age^2) + log(exper) + age*log(exper),data=byu)
如果我们要进行无常数项，一般在回归中引入0
result <- lm(smokes ~ 0 + male + female ,data=smokerdata)
3.2 从回归中提取统计量
一些统计量和参数都被存储在lm或者summary中
output <- summary (result)
SSR<- deviance（result）＃残差平方和；（另一种方法：RSquared <- output$r.squared)
LL<- logLik(result) #对数似然统计量
DegreesOfFreedom<- result$df #自由度
Yhat<- result$fitted.values #拟合值向量
Resid<- result$residuals
s<- output$sigma #误差标准差的估计值（假设同方差）
CovMatrix <- s^2*output$cov #系数的方差－协方差矩阵（与vcov()同）
3.3 异方差及相关问题
3.3.1 异方差的Breusch－Pagan检验
为了检验异方差是否存在，我们可以用lmtest包中的Breusch－Pagan检验。或者利用car包中的ncv.test()函数。二者工作的原理都是相同的。在回归之后，我们可以对拟合的模型采用bptest（）函数
unrestricted<- lm(z~x)
bptest(unrestricted)
这将得到检验的“学生化的”（studentized）结果。如果为了保持与其他软件结论的一致性（包括ncv.test())，我们可以设置studentize=FALSE
3.3.2 异方差（自相关）稳健性协方差矩阵
在存在异方差的情况下，ols的估计是无偏的，但是所得到的关于β系数方差的估计则是不正确的。为了计算异方差一致性的协方差矩阵，我们可以利用car包中的hccm（）函数，而不是vcov()。
sandwich包中的vcovHC()命令可以实现同样的功能。同时利用vcovHAC()或者NeweyWest（）函数可以进行异方差和自相关稳健性Newey—West估计。
####################
还有一章节，special regression，这个可能是计量中用的比较多的回归，比如logit，probit，还有tobit回归模型，以及M估计，稳健方法等等。
另外关于一般的lm方法由于用到的是ols或者称为lse方法，所以不存在log likelihood的说法，所以求出来的Loglik值可能没有意义的说。
另外用splus算过probit和logit回归，然后使用loglik函数得到的值竟然是错误的，手动演算一下跟书上一样，参考武德里奇的《计量》一书中离散选择模型第一五章的第一个例子。
####################
3.4 线性假设检验（Wald和F）
car包中提供了一个函数可以自动的进行线性假设检验。根据我们对模型的设定，它既可以用一般的方法或调整后的协方差矩阵进行F或Wald检验。为了进行假设检验，我们必须构造一个假设矩阵和一个右手边的向量。例如，如果我们有一个包括常数项的五个参数的模型，并且我们的零假设如下（见原文）
则假设的矩阵和右手边的向量将为如下的形式（见原文）
我们可以用如下的命令加以实现
unrestricted<-lm(y~x1+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是由lm得到的，默认状态下将会进行F检验。如果是由glm得到的，取而代之的将是Kai方检验。检验的类型可以通过type进行修改。
同样，如果我们想利用异方差或自相关稳健标准误进行检验，我们既可以通过设定white.adjust=TRUE来使用white标准误，也可以利用vcov计算我们自己的协方差矩阵。例如，如果我们想使用上述的Newey-West修正协方差矩阵，我们可以进行如下的设定：
linear.hypothesis(unrestricted, hm, rhs, vcov=NeweyWest(unrestricted))
注意：设定white.adjust=TRUE将会通过提高white估计量的精度来修正异方差；如果要使用经典的white估计量，我们可以设定white.adjust="hc0"
3.5 加权和广义最小二乘法
你可以通过使用带有权重的lm()来进行加权最小二乘
result<-lm(smokes~0+male+female, data=smokerdata, weights=myweights)
广义最小二乘法可以通过MASS包中的lm.gls()命令实现。它将包含一个函数、权重矩阵和一个数据框（可选）。
glm（）命令也为使用其他高级线性回归方法提供了渠道，详见帮助文件。
3.6 带有因子（Factors）或组别（Groups）的模型
在R中对于定性的因子有特定的数据类型。如果回归中的变量属于这种情况，必要的虚拟变量将会被自动生成。例如，如果我们要进行个人电脑的使用（pc）对公司雇员数（emple）和每一种状态的虚拟变量（其中state是两个缩写字母的向量），我们可以简单的进行如下的回归
summary(lm(pc~emple+state)
Call:
lm(formula = pc ~ emple + state)
Residuals:
Min 1Q Median 3Q Max
-1.7543 -0.5505 0.3512 0.4272 0.5904
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.572e-01 6.769e-02 8.232 <2e-16 ***
emple 1.459e-04 1.083e-05 13.475 <2e-16 ***
stateAL -4.774e-03 7.382e-02 -0.065 0.948
stateAR 2.249e-02 8.004e-02 0.281 0.779
stateAZ -7.023e-02 7.580e-02 -0.926 0.354
stateDE 1.521e-01 1.107e-01 1.375 0.169
...
stateFL -4.573e-02 7.136e-02 -0.641 0.522
stateWY 1.200e-01 1.041e-01 1.153 0.249
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4877 on 9948 degrees of freedom
Multiple R-Squared: 0.02451, Adjusted R-squared: 0.01951
F-statistic: 4.902 on 51 and 9948 DF, p-value: < 2.2e-16
为了将数据（序列string型或者数值型）转变成一个因子，可以简单的使用factor()命令。甚至可以在回归中间使用这个命令。例如，如果我们想做同样的回归，但是要区分以数字编码代表的不同区域，我们可以用如下的命令
myout<-lm(pc~emple+factor(naics6))
这一命令将naic6转换成了因子，生成了相应的虚拟变量，进而进行标准的回归。
####################
4. 特殊回归
4.1 固定和随机效应模型
注意：这里所用的“固定”和“随机”的概念与计量经济学家通常使用的概念相同。
经济学中，固定和随机效应是用来解释面板数据（panel data）模型的截距项中的截面（cross sectional）差异的。令i表示截面指标（或表示数据是有组别的），t为时间指标（或为组别差异指标）。一个标准的固定效应模型可以写作（参照原文）
本质上说，每一个个体都有不同的非随时间变化的截距项。通常我们感兴趣的是β，而不是ui。随机效应模型有同样的方程，但是相对固定效应模型而言，它有附加的限制：个体特殊的效应与解释变量x（it）不相关，即E[uiXit] = 0。从计量经济学的角度讲，这只是一个在固定效应模型的基础上，附有更加严格的限制条件的模型（它允许“效应”与外生变量相关）。
4.1.1 固定效应
在截面数据的维数不大的情况下，进行固定效应估计的简单方法是在每一个个体中加入一个虚拟变量，即将截面指标视为一个因子。如果指标能够在样本中将个体识别出来，则有
lm(y~factor(index)+x)
这个回归可以进行固定效应估计并能够正确的报告β的标准误。但不幸的是，在样本中存在很多个体的情况下，我们不再关心固定效应的值。因此在这种情况下，lm（）的结果以及u（i）较大的系数都是非常难于处理的。
一个更一般的方法是通过time demeaning（翻译不好，请大家帮助）的方法将每个变量的固定效应剔除（所谓内部within估计量）。则上述方程变为：（参照原文）
多数的统计软件（例如stata的xtreg命令）都使用这种方法处理固定效应模型。使用R，你可以手工对自变量和因变量进行time demean。如果d是一个包含 year, firm, profit和predictor的数据框，同时我们希望time demean每一个公司的profit和predictor，我们可以使用如下的命令：
> g<-d
> for (i in unique(d$firm)){
+ timemean<-mean(d[d$firm==i,])
+ g$profit[d$firm==i]<-d$profit[d$firm==i]-timemean["profit"]
+ g$predictor[d$firm==i]<-d$predictor[d$firm==i]-timemean["predictor"]
+ }
> output<-lm(profit~predictor,data=g)
要注意的是，回归中所报告的标准误偏低。lm()报告的方差使用的公式为
而真正的固定效应回归中的标准误使用的公式为（参照原文）
对于T较小的样本，二者之间存在较为显著的差异。
另一种并不常用的方法是用一阶差分，公式为（参照原文）
其手工计算方法可以参照上述的within estimator
4.1.2 随机效应模型
nlme包包括了在线性或非线性数据框中进行随机效应回归的方法（而不是固定效应，本文只涉及对“固定效应”项的统计解释）。
假设存在如下的以sic3编码为随机效应的线性模型
我们可以用下列命令进行拟合
lme(ldsal~lemp+ldnpt,random=~1|sic3)
一般而言，在随机参数模型中，随机效应被置于竖线之后。波浪线和竖线之间的1表示随机效应是一个截距项。如果随机效应是一个截距项和一个外生变量，我们应当将该变量与1放在一起。例如：
lme(ldsal~lemp+ldnpt,random=~1+lemp|sic3)
对应于
对于非线性随机效应模型而言，只是将lme()替换为nlme()就可以了。
#####################
4.2 定性相应模型（Qualitative Response)
4.2.1 Logit/Probit
有很多种方法可以在R中实现logit和probit回归。最简单的方法是使用glm()命令及相应的选项。
h<-glm(c~y,family=binomial(link="logit"))
对于probit回归，将logit替换为probit即可。glm()函数的输出结果与lm()的类似，因此可以使用summary()命令加以分析。为了从中提取出对数似然统计量，可以使用logLik()命令
> logLik(h)
‘log Lik.’ -337.2659 (df=1)
4.2.2 多项式Logit （Multinormial Logit）
在nnet包中有一个multinom()函数可以用来进行多项式Logit估计。为了使用这一函数，首先应该将因变量转化成一个因子向量（包括所有情况），并且使用诸如正态回归这样的语法（syntax）。如果我们的因子以虚拟变量的形式被储存，则我们可以利用十进制数字的特征为所有的组合赋予唯一的因子。假如我们的因子变量是pc,inetacc和iapp，那么
> g <- pc*1 + inetacc*10 + iapp*100
> multinom(factor(g)~pc.subsidy+inet.subsidy+iapp.subsidy+emple+msamissing)
这样，我们利用所有因子的组合得到了一个多项式Logit。
多项式Probit模型的一个特征就是常被ill conditioned（如何译？）。一个解决的方法是使用MNP包中的mnp（）命令进行马尔可夫链蒙特卡罗模拟。
4.2.3 Ordered Logit/Probit
MASS包中有一个polr()函数可以进行ordered logit或probit回归。如果Sat表示一个顺序（ordered）的因子向量，则
house.plr <- polr(Sat ~ Infl + Type + Cont, method="probit")
4.3 Tobit和阶段（Censored）回归
我们用survival包来估计存在截断数据变量的模型。使用的函数是survreg(),其中将因变量作为一个Surv对象。假设我们要进行y对x和z的回归，但是大量的y数据是左侧截断的，并用0将其替换。
result <- survreg(Surv(y,y>0,type='left') ~ x + z, dist='gaussian')
第二点需要注意的是无论观测值是否是截断的，都可以使用Surv()函数（1表示是可以被观测的；0表示数据是截断的）。第三点需要说明的是，数据在哪一侧被截断。既然本例中数据在分布的低尾（lower tail）处被截断，我们使用left。dist选项对于survreg是必需的，这样才能得到一个经典的Tobit模型。
4.4 分位数回归
最小二乘回归方法可以估测因变量对自变量的条件期望。拟合值即是条件均值的估计。如果我们不想得到条件均值，而想估计预期条件中位数或其他分位数的话，我们可以使用quantreg包中的rq()命令。语法与lm()基本相同，除了我们要使用一个介于0和1之间的分位数tau。默认的情况下，tau＝.5，即为中位数，另一个名字是最小绝对偏差回归（least absolute deviation regression)
4.5 稳健性回归——M估计量
对于一些数据集，奇异值对最小二乘回归线的影响远远超出了我们的预想。一个解决的办法是使用包括残差平方和（对应于最小化L2）在内的一些方法求极小值并以此作为目标方程。通常的选择是使用绝对离差和（L1）和Huber法——一种将L1和L2混合的方法。R使用MASS包中的rlm()进行稳健性回归。语法与lm()相同除了它允许选择最小化作为目标方程。进行这种选择可以使用参数psi。可能的选项包括psi.huber, psi.hampel, psi.bisquare。
为了进行psi函数的定制，我们写了一个函数，如果deriv=0，函数返回ψ(x)/x；如果deriv=1，返回ψ′(x)/x。This function than then be passed to rlm() using the psi parameter.＃不清楚函数内容及语意。
4.6 非线性最小二乘
有些时候，经济中的模型并不是线性的。R可以进行如下形式的广义最小二乘
注意，残差项必须是附加在函数形式上的。如果不是，则必须进行相应的变换以达到这种形式。R中进行非线性最小二乘的函数是nls()，其语法与lm()相同。考虑如下的非线性例子：
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包。
4.7单一结构方程的两阶段最小二乘
为了实现单一方程的两阶段最小二乘，最简单的方法是使用sem包中的tsls()命令。如果我们想考察在控制了婚姻状况的情况下，教育对工资的影响，但是考虑到educ可能是内生的，则我们可能会使用motheduc和fatheduc作为工具变量进行回归
library(sem)
outputof2sls<-tsls(lwage~educ+married,~married+motheduc+fatheduc)
第一点需要说明的是我们进行估计的是一个结构性方程；第二，波浪线后面需要附加结构性方程中的工具变量和外生变量，这些条件的满足需要下面的2SLS估计量满足
输出结果可以用summary()函数和其他ols函数进行分析。注意，既然这个命令输出的是两阶段最小二乘的对象，包括标准误在内的描述性统计量都是正确的。如果我们用实际的两阶段方法进行估计，则输出的标准误是存伪的。
4.8 方程组
用于处理方程组（包括工具变量、两阶段最小二乘、似不相关回归[seemingly unrelated regression]
和变量)的命令包含在systemfit包中。注意，在R中，一个方程（包括标题）就恰恰是一种数据类型。因此，我们可以通过通常的赋值运算来构造方程模型列表以及相应的标签列表。
> demand <- q ~ p + d
> supply <- q ~ p + f + a
> system <- list(demand,supply)
> labels <- list("demand","supply")
4.8.1 似无关回归
一旦我们建立了系统和标签，我们就可以使用systemfit()函数结合SUR选项来识别似无关回归。
> resultsur <- systemfit("SUR",system,labels)
4.8.2 方程组的两阶段最小二乘
可以使用工具变量来对上述的方程组进行两阶段最小二乘。我们创建一个模型对象（不包含左侧）来表示我们要使用的工具变量和2SLS选项
> inst <- ~ d + f + a
> result2sls <- systemfit("2SLS",system,labels,inst)
与之相类似，我们还会遇到三阶段最小二乘、加权两阶段最小二乘和相应的其他模型。
####################
about fixed- and random-effect model, can refer to plm package.
about multinomial logit model, can refer to mlogit package.
about tobit model, refer to sampleSelection package.
####################
颇值得一提的是，回归中最紧要看residuals plot，
example(plot.lm)
####################
扫一扫
188
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
