R 编程/线性模型

标准线性模型

在本节中，我们介绍了通过普通最小二乘法 (OLS) 估计的标准线性模型的估计函数。异方差和内生性将在下面讨论。主要的估计函数是lm().

伪数据模拟

我们首先生成一个伪数据集，以便没有异方差、没有内生性和误差项之间没有相关性。因此，普通最小二乘估计量是无偏且有效的。我们选择一个包含两个变量的模型，并将所有系数设置为 1。

$y_{i}=1+x_{1,i}+x_{2,i}+u_{i}$

> N <- 1000
> u <- rnorm(N)
> x1 <- rnorm(N)
> x2 <- 1 + x1 + rnorm(N)
> y <- 1 + x1 + x2 + u
> df <- data.frame(y,x1,x2)

最小二乘估计

• 估计简单线性模型的标准函数是 lm()。
• lsfit()执行最小二乘程序，但输出没有以时尚的方式格式化。
• ols()(Design) 是另一种选择。

我们使用以下方法估计模型lm(). 我们将结果存储在fit中，并使用以下方法打印结果summary()这是标准函数。

> fit <- lm(y ~ x1 + x2, data = df)
> summary(fit)

有一些替代方案可以显示结果。

• display()在 arm 包中是其中之一。
• coefplot()(arm) 绘制了带置信区间的估计系数。这是展示结果的好方法。
• mtable()在 memisc 包中可以将一组回归的结果显示在同一个表中。

> library("arm")
> display(fit)
> coefplot(fit)

fit是一个对象列表。您可以通过键入以下内容查看这些对象的列表names(fit). 我们也可以将函数应用于fit.

我们可以使用以下方法获取估计系数fit$coeff或coef(fit).

> fit$coeff
(Intercept)          x1          x2 
  1.2026522   0.8427403   1.5146775
> coef(fit)
(Intercept)          x1          x2 
     0.7541      1.7844      0.7222 
> output <- summary(fit)
> coef(output) 
             Estimate Std. Error  t value    Pr(>|t|)
(Intercept) 1.1945847  0.2298888 5.196359 0.001258035
x1          0.6458170  0.3423214 1.886581 0.101182585
x2          0.6175165  0.2083628 2.963660 0.020995713

se.coef()(arm) 返回估计系数的标准误差。

拟合值向量可以通过以下方式返回fit$fitted, fitted(fit)或predict()函数。该predict()函数还返回预测的标准误差和置信区间。

 
> fit$fitted
> fitted(fit)

残差向量

> fit$resid
> residuals(fit)

自由度

> fit$df

置信区间

我们可以使用以下方法获取置信区间confint()或conf.intervals()在 alr3 包中。

> confint(fit, level = .9)
                   5 %     95 %
(Intercept) -0.7263261 1.200079
x1          -0.5724022 1.909924
x2           0.6185011 2.475079
> confint(fit, level = .95)
                 2.5 %   97.5 %
(Intercept) -0.9652970 1.439050
x1          -0.8803353 2.217858
x2           0.3881923 2.705388
> confint(fit, level = .99)
                 0.5 %   99.5 %
(Intercept) -1.5422587 2.016012
x1          -1.6237963 2.961319
x2          -0.1678559 3.261436
> library(alr3)
> conf.intervals(fit)
                 2.5 %   97.5 %
(Intercept) -0.9652970 1.439050
x1          -0.8803353 2.217858
x2           0.3881923 2.705388

检验

coeftest()(lmtest) 对系数执行学生 t 检验和 z 检验。

> library("lmtest")
> coeftest(fit) # t-test
> coeftest(fit,df=Inf) # z-test (for large samples)

linear.hypothesis()(car) 对线性假设执行有限样本 F 检验，或使用 $\mathrm {X} ^{2}$ 统计量执行渐近 Wald 检验。

> library("car")
> linear.hypothesis(fit,"x1 = x2") # tests Beta1 = Beta2
> linear.hypothesis(fit,c("(Intercept)", "x1","x2"),rep(1,3)) # Tests  Beta0 = Beta1 = Beta2 = 1
> linear.hypothesis(fit,c("(Intercept)", "x1","x2"),rep(0,3)) # Tests  Beta0 = Beta1 = Beta2 = 0
> linear.hypothesis(fit,c("x1","x2"),rep(0,2)) # Tests Beta1 = Beta2 = 0

另请参阅waldtest()(lmtest) 用于嵌套模型。

方差分析

我们还可以使用以下方法进行方差分析anova().

> anova(fit)

模型搜索和信息准则

> # Akaike Information Criteria
> AIC(fit)
[1] 26.72857
> # Bayesian Information Criteria
> AIC(fit,k=log(N))
[1] 27.93891

stats4 包含AIC()和BIC()函数

> library(stats4)
> ?BIC
> lm1 <- lm(Fertility ~ . , data = swiss)
> AIC(lm1)
[1] 326.0716
> BIC(lm1)
[1] 339.0226

该step()函数使用 Akaike 信息准则执行模型搜索。

> N <- 10^3
> u <- rnorm(N)
> x1 <- rnorm(N)
> x2 <- rnorm(N) + x1
> x3 <- rnorm(N)
> y <- 1+ x1 + x2 + u
> fit <- lm(y~x1+x2 + x3)
> step.fit <- step(fit)

Zelig

• 该方法也在 Zelig 中得到支持

> N <- 1000
> u <- rnorm(N)
> x <- rnorm(N)
> y <- 1 + x + u
> mydat <- data.frame(y,x)
> z.out <- zelig(y ~  x, model = "ls", data = mydat)
> x.out <- setx(z.out, x = 10)
> s.out <- sim(z.out, x.out)
> summary(s.out)

贝叶斯估计

• MCMCregress()(MCMCpack)
• BLR()(BLR)

> N <- 1000
> u <- rnorm(N)
> x <- rnorm(N)
> y <- 1 + x + u
> mydat <- data.frame(y,x)
> 
> posterior <- MCMCregress(y ~ x, data = mydat)
> summary(posterior)
> plot(posterior)

异方差

• 请参阅 lmtest 和 sandwich 包。
• gls()(nlme) 计算广义最小二乘估计量。
• 请参阅 Mahmood Arai 的“使用 R 进行聚类稳健标准误差”(pdf)。他建议使用两个函数来计算聚类稳健标准误差。clx()允许进行单向聚类，而mclx()允许进行双向聚类。可以使用以下命令加载它们source("http://people.su.se/~ma/clmclx.R").

> N <- 10 # 10 people
> T <- 5 # 5 times
> id <- rep(1:N,T)
> f <- rep(rnorm(N),T) # is individual specific
> u <- rnorm(N*T)
> x1 <- rnorm(N*T) 
> x2 <- rnorm(N*T) + x1
> y <- 1 + x1 + x2 + f + u
> fit <- lm(y ~ x1 + x2 )
> source("http://people.su.se/~ma/clmclx.R")
> clx(fit, 1, id)

稳健性

库克距离

>library(car)
> cookd(fit)
           1            2            3            4            5 
0.0006205008 0.0643213760 0.2574810866 1.2128206779 0.2295047699 
           6            7            8            9           10 
0.3130578329 0.0003365221 0.0671830241 0.0048474954 0.0714255871

影响图

> influence.plot(fit)

杠杆图

> leverage.plot(fit,term.name=x1)
> leverage.plot(fit,term.name=x2)

邦费罗尼异常值检验

> outlier.test(fit)

max|rstudent| = 2.907674, degrees of freedom = 6,
unadjusted p = 0.02706231, Bonferroni p = 0.2706231

Observation: 3

另请参阅outlier.t.test()在 alr3 包中。

• inf.index()在 alr3 包中计算所有稳健性统计量（库克距离、学生化残差、异常值检验等）
• rlm()执行稳健估计

• 请参阅 UCLA 示例
• 另请参阅 robustbase 包

工具变量

• ivreg()在 AER 包中
• tsls()在 sem 包中。
• 也可以使用gmm()命令在 gmm 包中。请参阅矩量法以获取示例。

伪数据模拟

我们首先模拟一个假数据集，其中x与u相关，z与u独立，x与z相关。因此x是y的内生解释变量，z是x的有效工具变量。

> N <- 1000
> z <- rnorm(N)
> u <- rnorm(N) 
> x <- 1 + z + u + rnorm(N) # x is correlated with the error term u (endogeneity) and the instrument z
> y <- 1 + x + u

两阶段最小二乘法

然后我们使用OLS估计模型（lm()）并使用z作为x的工具变量进行IV估计。

> ols <- lm(y ~ x)
> summary(ols) # ols are biased
> library("AER")
> iv <- ivreg(y ~ x | z)
> summary(iv) # IV estimates are unbiased
> library("sem")
> iv2 <- tsls(y  ~ x, instruments = ~ z)
> summary(iv2)
> library("gmm")
> iv3 <- gmm(y ~ x, z)
> summary(iv3)

我们绘制结果 

> plot(y ~ x, col = "gray")
> abline(a  = 1,b = 1, lty = 1, col = 1, lwd = 2)
> abline(ols,  lty = 2, col = 2 , lwd = 2)
> abline(iv, lty = 3, col = 3, lwd = 2)
> legend("topleft", legend = c("True values","OLS","IV"), col = 1:3, lwd = rep(2,3), lty = 1:3)

面板数据

plm()(plm) 实现标准的随机效应、固定效应、差分方法。它类似于Stata的xtreg命令。

请注意，plm输出与xtable()和mtable()不兼容，无法用于出版物质量输出。

• lme4和gee实现随机效应和多级模型。
• 另请参见BayesPanel

随机效应模型

为了实现随机效应模型，我们生成一个包含1000个观测值和5个时间段的假数据集。

> N <- 1000
> T <- 5
> library(mvtnorm)
> sig <- diag(rep(1,T))
> r <- rmvnorm(N, sigma = sig)
> wide <- data.frame(id = 1:N,f = rnorm(N), u = r)
> long <- reshape(wide, varying = list(3:7), v.names = "u", direction = "long", timevar = "year")
> long$x1 <- 1 + rnorm(N*T) 
> long$x2 <- 1 + rnorm(N*T) + long$x1
> long$y <- 1 + long$x1 + long$x2 + long$f + long$u
> head(long[order(long$id),])

我们使用plm()函数和model = "random"选项估计随机效应模型。

> library("plm")
> panel <- plm.data(long, index = c("id","year"))
> # panel <- pdata.frame(long,c("id","year"))
> eq <- y ~ x1 + x2
> re <- plm(eq, model = "random", data=panel)
> summary(re)

固定效应模型

对于固定效应模型，我们生成一个假数据集，并将固定效应f与协变量相关联 

> N <- 1000
> T <- 5
> library(mvtnorm)
> sig <- diag(rep(1,T))
> r <- rmvnorm(N, sigma = sig)
> wide <- data.frame(id = 1:N,f = rnorm(N), u = r)
> long <- reshape(wide, varying = list(3:7), v.names = "u", direction = "long", timevar = "year")
> long$x1 <- 1 + rnorm(N*T) + long$f
> long$x2 <- 1 + rnorm(N*T) + long$x1
> long$y <- 1 + long$x1 + long$x2 + long$f + long$u
> head(long[order(long$id),])

我们首先使用plm.data()将我们的数据转换为plm数据框。我们使用plm()估计固定模型，model = "within"作为选项。然后，我们比较固定效应模型和随机效应模型的估计结果，并进行豪斯曼检验。最后，我们绘制固定效应的密度图。

> library("plm")
> panel <- plm.data(long, index = c("id","year"))
> #panel <- pdata.frame(long,c("id","year"))
> eq <- y ~ x1 + x2
> fe <- plm(eq, model = "within", data=panel)
> summary(fe)
> re <- plm(eq, model = "random", data=panel)
> summary(re)
> phtest(fe, re)
> plot(density(fixef(fe)))
> rug(fixef(fe))

动态面板数据

• pgmm()(plm) 实现Arellano Bond估计过程。它类似于Stata的xtabond2命令。

联立方程模型

对于[:w:Simultaneous_equations_model|联立方程模型]，需要以下软件包 

• sem软件包
• systemfit软件包

外部链接

• 回归系数的可视化