【原】R语言异方差回归模型建模：用误差方差解释异方差

原文链接：http:///?p=10207

---

在社会科学中将OLS估计应用于回归模型时，其中的一个假设是同方差，我更喜欢常误差方差。这意味着误差方差没有系统的模式，这意味着该模型在所有预测级别上都同样差。

异方差性是同方差性的补充，不会使OLS产生偏差。如果您不像社会科学中的大多数人那样关心p值，那么异方差性可能不是问题。

计量经济学家已经开发出各种各样的异方差一致性标准误差，因此他们可以继续应用OLS，同时调整非恒定误差方差。这些更正的Wikipedia页面列出了这些替代标准错误所使用的许多名称。

我们提供了似然函数，并且两个函数都将找到使似然最大化的参数估计。

让我们来看一个简单的例子：

首先，我从均值3和标准差1.5的正态分布中提取500个观测值，并将其保存到数据集中：

dat <- data.frame(y = rnorm(n = 500, mean = 3, sd = 1.5))


样本的平均值和标准偏差为：

mean(dat$y)[1] 2.999048
sd(dat$y)[1] 1.462059




我也可以这样问这个问题，正态分布，均值和标准差的哪些参数可以最大程度地提高观察到的变量的可能性？

m.sd <- mle2(y ~ dnorm(mean = a, sd = exp(b)), data = dat, start = list(a = rnorm(1), b = rnorm(1)))




在上面的语法中，R变量y的平均值是一个常数a，而y的标准偏差是一个常数b。标准差取幂，确保它永远不会为负数。我们提供初始值，因此它可以在收敛到使可能性最大化的值之前开始估算。随机数足以满足初始值。

m.sd
Call:mle2(minuslogl = y ~ dnorm(mean = a, sd = exp(b)), start = list(a = rnorm(1), b = rnorm(1)), data = dat)
Coefficients: a b2.9990478 0.3788449
Log-likelihood: -898.89




系数a非常类似于数据的平均值。必须对系数b取幂，以获得标准偏差：

exp(coef(m.sd)[2]) b1.460596




这类似于我们上面获得的标准偏差。上面的语法演示的另一个有趣的事实是lm()类似的函数coef()，summary()并且可以在mle2()对象上使用。

我们上面执行的最大似然估计类似于使用OLS估计的仅截距回归模型：

coef(lm(y ~ 1, dat))(Intercept) 2.999048
sigma(lm(y ~ 1, dat))[1] 1.462059




截距是数据的平均值，残留标准偏差是标准偏差。

异方差回归模型

考虑以下研究。我们分配了两组，一个是治疗组，一个是30个人，另一个是对照组，每个是100个人，与治疗组相匹配的是决定结果的协变量。因此，我们对治疗效果感兴趣，并让我们假设一个简单的均值差就足够了。碰巧，这种治疗方法除了有效之外，还具有均质作用，例如，受试者被洗脑后对结果的改善更好。以下数据集应符合上述方案：

有100名参与者的治疗状态为0（对照组），平均值为0，标准差为1。有30名参与者的治疗状态为1（治疗组），平均值为0.3，标准值为1，偏差0.25。

这种情况显然违反了同方差假设，但是，我们继续对治疗效果进行OLS估计：




Call:
Residuals: Min 1Q Median 3Q Max-2.8734 -0.5055 -0.0287 0.4231 3.4097
Coefficients: Estimate Std. Error t value Pr(>|t|)(Intercept) 0.03386 0.09298 0.364 0.716treat 0.21733 0.19355 1.123 0.264
Residual standard error: 0.9298 on 128 degrees of freedomMultiple R-squared: 0.009754, Adjusted R-squared: 0.002018F-statistic: 1.261 on 1 and 128 DF, p-value: 0.2636




治疗效果为0.22，无统计学意义，p = 0.26p=.26在一个αα的.05级。但是我们知道方差不是同方差的，因为我们创建了数据，并且残差对拟合值的简单诊断图证实了这一点：

首先，我记录一下重新创建OLS模型：

在此函数中，我为结果的平均值创建一个模型，该模型是截距的函数b_int，以及治疗预测因子的系数b_treat。标准偏差还是一个指数常数。该模型将等效于线性模型。

但是，我们知道方差不是恒定的，而是两组不同。我们可以将标准偏差指定为组的函数：

在此，我们为标准差指定了一个模型，该模型作为截距的函数s_int，代表控制组，并且与该截距的偏差为s_treat。

我们可以做得更好。我们可以利用系数从OLS模型作为初始值b_int和b_treat。运行模型：





Maximum likelihood estimation
Call:(minuslogl = y ~ dnorm(mean = b_int + b_treat * treat, sd = exp(s_int + s_treat * treat)), start = list(b_int = coef(m.ols)[1], b_treat = coef(m.ols)[2], s_int = rnorm(1), s_treat = rnorm(1)))
Coefficients: Estimate Std. Error z value Pr(z)b_int 0.033862 0.104470 0.3241 0.74584 b_treat 0.217334 0.112249 1.9362 0.05285 .s_int 0.043731 0.070711 0.6184 0.53628 s_treat -1.535894 0.147196 -10.4344 < 2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-2 log L: 288.1408




治疗效果大致相同，但现在p值为.053。远小于假设为纯正方差分析的0.26。b_treat变量的精度要高得多，因为此处的标准误差.11小于.19。

标准差模型建议标准差为：

exp(coef(m.het)[3])
s_int1.044701




对照组和1.045：

exp(coef(m.het)[3] + coef(m.het)[4])
s_int0.2248858




.22为治疗组。这些值接近我们所知道的模拟值。我们可以确认样本统计数据为：




treat y1 0 1.04996572 1 0.2287307




在没有异方差且允许异方差的情况下，也可以轻松地对模型进行模型比较：




Likelihood Ratio TestsModel 1: m.mle, y~dnorm(mean=b_int+b_treat*treat,sd=exp(s1))Model 2: m.het, y~dnorm(mean=b_int+b_treat*treat,sd=exp(s_int+s_treat*treat)) Tot Df Deviance Chisq Df Pr(>Chisq)1 3 347.98 2 4 288.14 59.841 1 1.028e-14 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1




似然比测试建议我们改进了模型，χ 2（1 ）= 59.81 ，p < 0.001χ2（1个）=59.81，p<.001。

因此，我们可以确认在此单个示例中对方差建模可以提高精度。当影响为零并且我们具有异方差性时，很容易编写一个将异方差MLE与OLS估计进行比较的仿真代码。

我从上面对代码进行了更改，方法是给治疗组的平均值为零，以使两组之间没有均值差。我重复了该过程500次，从OLS及其p值中节省了治疗效果，从异方差MLE及其p值中节省了治疗效果。

然后，我绘制结果：




par(mfrow = c(1, 1))




OLS和异方差性MLE的治疗效果相似。但是，当null为true时，异方差MLE模型的p值表现得更好。如果null为true，则可以期望p值均匀分布。OLS迭代的p值堆叠在高端。

这次，我重复此过程，使治疗组的平均值为0.15，因此零效果的null假设为假。 

治疗效果再次具有相同的分布。然而，与OLS相比，异方差MLE的p值要小得多，异方差MLE具有更大的统计功效来检测治疗效果。

---

首先，为负对数可能性指定一个函数，然后将此函数传递给MLE。




(minuslogl = ll, start = list(b_int = rnorm(1), b_treat = rnorm(1), s_int = rnorm(1), s_treat = rnorm(1)))
Coefficients: Estimate Std. Error z value Pr(z)b_int 0.033862 0.104470 0.3241 0.74584 b_treat 0.217334 0.112249 1.9362 0.05285 .s_int 0.043733 0.070711 0.6185 0.53626 s_treat -1.535893 0.147196 -10.4343 < 2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-2 log L: 288.1408


Family: gaussian ( identity )Formula: y ~ treatDispersion: ~treatData: dat
AIC BIC logLik deviance df.resid 296.1 307.6 -144.1 288.1 126

Conditional model: Estimate Std. Error z value Pr(>|z|)(Intercept) 0.03386 0.10447 0.324 0.7458 treat 0.21733 0.11225 1.936 0.0528 .---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Dispersion model: Estimate Std. Error z value Pr(>|z|)(Intercept) 0.08746 0.14142 0.618 0.536 treat -3.07179 0.29439 -10.434 <2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1




在这种情况下，离散度在对数方差的范围内，因此必须取平方的指数对数方差平方根才能检索上述的组标准差。