原文鏈接：http://tecdat.cn/?p=10207

?

在社會(huì)科學(xué)中將OLS估計(jì)應(yīng)用于回歸模型時(shí)，其中的一個(gè)假設(shè)是同方差，我更喜歡常誤差方差。這意味著誤差方差沒有系統(tǒng)的模式，這意味著該模型在所有預(yù)測(cè)級(jí)別上都同樣差。

異方差性是同方差性的補(bǔ)充，不會(huì)使OLS產(chǎn)生偏差。如果您不像社會(huì)科學(xué)中的大多數(shù)人那樣關(guān)心p值，那么異方差性可能不是問題。

計(jì)量經(jīng)濟(jì)學(xué)家已經(jīng)開發(fā)出各種各樣的異方差一致性標(biāo)準(zhǔn)誤差，因此他們可以繼續(xù)應(yīng)用OLS，同時(shí)調(diào)整非恒定誤差方差。這些更正的Wikipedia頁面列出了這些替代標(biāo)準(zhǔn)錯(cuò)誤所使用的許多名稱。

我們提供了似然函數(shù)，并且兩個(gè)函數(shù)都將找到使似然最大化的參數(shù)估計(jì)。

讓我們來看一個(gè)簡單的例子：

首先，我從均值3和標(biāo)準(zhǔn)差1.5的正態(tài)分布中提取500個(gè)觀測(cè)值，并將其保存到數(shù)據(jù)集中：

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

樣本的平均值和標(biāo)準(zhǔn)偏差為：

1. mean(dat$y)

2. [1] 2.999048

3. 


4. sd(dat$y)

5. [1] 1.462059

我也可以這樣問這個(gè)問題，正態(tài)分布，均值和標(biāo)準(zhǔn)差的哪些參數(shù)可以最大程度地提高觀察到的變量的可能性？

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

2. start = list(a = rnorm(1), b = rnorm(1)))

在上面的語法中，R變量y的平均值是一個(gè)常數(shù)a，而y的標(biāo)準(zhǔn)偏差是一個(gè)常數(shù)b。標(biāo)準(zhǔn)差取冪，確保它永遠(yuǎn)不會(huì)為負(fù)數(shù)。我們提供初始值，因此它可以在收斂到使可能性最大化的值之前開始估算。隨機(jī)數(shù)足以滿足初始值。

1. m.sd

2. 


3. Call:

4. mle2(minuslogl = y ~ dnorm(mean = a, sd = exp(b)), start = list(a = rnorm(1),

5. b = rnorm(1)), data = dat)

6. 


7. Coefficients:

8. a ? ? ? ? b

9. 2.9990478 0.3788449

10. 


11. Log-likelihood: -898.89

系數(shù)a非常類似于數(shù)據(jù)的平均值。必須對(duì)系數(shù)b取冪，以獲得標(biāo)準(zhǔn)偏差：

1. exp(coef(m.sd)[2])

2. b

3. 1.460596

這類似于我們上面獲得的標(biāo)準(zhǔn)偏差。上面的語法演示的另一個(gè)有趣的事實(shí)是lm()類似的函數(shù)coef()，summary()并且可以在mle2()對(duì)象上使用。

我們上面執(zhí)行的最大似然估計(jì)類似于使用OLS估計(jì)的僅截距回歸模型：

1. coef(lm(y ~ 1, dat))

2. (Intercept)

3. 2.999048

4. 


5. sigma(lm(y ~ 1, dat))

6. [1] 1.462059

截距是數(shù)據(jù)的平均值，殘留標(biāo)準(zhǔn)偏差是標(biāo)準(zhǔn)偏差。

異方差回歸模型

考慮以下研究。我們分配了兩組，一個(gè)是治療組，一個(gè)是30個(gè)人，另一個(gè)是對(duì)照組，每個(gè)是100個(gè)人，與治療組相匹配的是決定結(jié)果的協(xié)變量。因此，我們對(duì)治療效果感興趣，并讓我們假設(shè)一個(gè)簡單的均值差就足夠了。碰巧，這種治療方法除了有效之外，還具有均質(zhì)作用，例如，受試者被洗腦后對(duì)結(jié)果的改善更好。以下數(shù)據(jù)集應(yīng)符合上述方案：

有100名參與者的治療狀態(tài)為0（對(duì)照組），平均值為0，標(biāo)準(zhǔn)差為1。有30名參與者的治療狀態(tài)為1（治療組），平均值為0.3，標(biāo)準(zhǔn)值為1，偏差0.25。

這種情況顯然違反了同方差假設(shè)，但是，我們繼續(xù)對(duì)治療效果進(jìn)行OLS估計(jì)：

1. 


2. Call:

3. 


4. Residuals:

5. Min ? ? ?1Q ?Median ? ? ?3Q ? ? Max

6. -2.8734 -0.5055 -0.0287 ?0.4231 ?3.4097

7. 


8. Coefficients:

9. Estimate Std. Error t value Pr(>|t|)

10. (Intercept) ?0.03386 ? ?0.09298 ? 0.364 ? ?0.716

11. treat ? ? ? ?0.21733 ? ?0.19355 ? 1.123 ? ?0.264

12. 


13. Residual standard error: 0.9298 on 128 degrees of freedom

14. Multiple R-squared: ?0.009754, Adjusted R-squared: ?0.002018

15. F-statistic: 1.261 on 1 and 128 DF, ?p-value: 0.2636

治療效果為0.22，無統(tǒng)計(jì)學(xué)意義，p = 0.26p=.26在一個(gè)αα的.05級(jí)。但是我們知道方差不是同方差的，因?yàn)槲覀儎?chuàng)建了數(shù)據(jù)，并且殘差對(duì)擬合值的簡單診斷圖證實(shí)了這一點(diǎn)：

首先，我記錄一下重新創(chuàng)建OLS模型：

在此函數(shù)中，我為結(jié)果的平均值創(chuàng)建一個(gè)模型，該模型是截距的函數(shù)b_int，以及治療預(yù)測(cè)因子的系數(shù)b_treat。標(biāo)準(zhǔn)偏差還是一個(gè)指數(shù)常數(shù)。該模型將等效于線性模型。

但是，我們知道方差不是恒定的，而是兩組不同。我們可以將標(biāo)準(zhǔn)偏差指定為組的函數(shù)：

在此，我們?yōu)闃?biāo)準(zhǔn)差指定了一個(gè)模型，該模型作為截距的函數(shù)s_int，代表控制組，并且與該截距的偏差為s_treat。

我們可以做得更好。我們可以利用系數(shù)從OLS模型作為初始值b_int和b_treat。運(yùn)行模型：

1. 


2. 


3. Maximum likelihood estimation

4. 


5. Call:

6. (minuslogl = y ~ dnorm(mean = b_int + b_treat * treat, sd = exp(s_int +

7. s_treat * treat)), start = list(b_int = coef(m.ols)[1], b_treat = coef(m.ols)[2],

8. s_int = rnorm(1), s_treat = rnorm(1)))

9. 


10. Coefficients:

11. Estimate Std. Error ?z value ? Pr(z)

12. b_int ? ?0.033862 ? 0.104470 ? 0.3241 0.74584

13. b_treat ?0.217334 ? 0.112249 ? 1.9362 0.05285 .

14. s_int ? ?0.043731 ? 0.070711 ? 0.6184 0.53628

15. s_treat -1.535894 ? 0.147196 -10.4344 < 2e-16 ***

16. ---

17. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

18. 


19. -2 log L: 288.1408

治療效果大致相同，但現(xiàn)在p值為.053。遠(yuǎn)小于假設(shè)為純正方差分析的0.26。b_treat變量的精度要高得多，因?yàn)榇颂幍臉?biāo)準(zhǔn)誤差.11小于.19。

標(biāo)準(zhǔn)差模型建議標(biāo)準(zhǔn)差為：

1. exp(coef(m.het)[3])

2. 


3. s_int

4. 1.044701

對(duì)照組和1.045：

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

2. 


3. s_int

4. 0.2248858

.22為治療組。這些值接近我們所知道的模擬值。我們可以確認(rèn)樣本統(tǒng)計(jì)數(shù)據(jù)為：

1. 


2. treat ? ? ? ? y

3. 1 ? ? 0 1.0499657

4. 2 ? ? 1 0.2287307

在沒有異方差且允許異方差的情況下，也可以輕松地對(duì)模型進(jìn)行模型比較：

1. 


2. Likelihood Ratio Tests

3. Model 1: m.mle, y~dnorm(mean=b_int+b_treat*treat,sd=exp(s1))

4. Model 2: m.het, y~dnorm(mean=b_int+b_treat*treat,sd=exp(s_int+s_treat*treat))

5. Tot Df Deviance ?Chisq Df Pr(>Chisq)

6. 1 ? ? ?3 ? 347.98

7. 2 ? ? ?4 ? 288.14 59.841 ?1 ?1.028e-14 ***

8. ---

9. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

似然比測(cè)試建議我們改進(jìn)了模型，χ 2（1 ）= 59.81 ，p < 0.001χ2（1個(gè)）=59.81，p<.001。

因此，我們可以確認(rèn)在此單個(gè)示例中對(duì)方差建?？梢蕴岣呔?。當(dāng)影響為零并且我們具有異方差性時(shí)，很容易編寫一個(gè)將異方差MLE與OLS估計(jì)進(jìn)行比較的仿真代碼。

我從上面對(duì)代碼進(jìn)行了更改，方法是給治療組的平均值為零，以使兩組之間沒有均值差。我重復(fù)了該過程500次，從OLS及其p值中節(jié)省了治療效果，從異方差MLE及其p值中節(jié)省了治療效果。

然后，我繪制結(jié)果：

1. 


2. par(mfrow = c(1, 1))

OLS和異方差性MLE的治療效果相似。但是，當(dāng)null為true時(shí)，異方差MLE模型的p值表現(xiàn)得更好。如果null為true，則可以期望p值均勻分布。OLS迭代的p值堆疊在高端。

這次，我重復(fù)此過程，使治療組的平均值為0.15，因此零效果的null假設(shè)為假。?

治療效果再次具有相同的分布。然而，與OLS相比，異方差MLE的p值要小得多，異方差MLE具有更大的統(tǒng)計(jì)功效來檢測(cè)治療效果。

首先，為負(fù)對(duì)數(shù)可能性指定一個(gè)函數(shù)，然后將此函數(shù)傳遞給MLE。

1. 


2. (minuslogl = ll, start = list(b_int = rnorm(1), b_treat = rnorm(1),

3. s_int = rnorm(1), s_treat = rnorm(1)))

4. 


5. Coefficients:

6. Estimate Std. Error ?z value ? Pr(z)

7. b_int ? ?0.033862 ? 0.104470 ? 0.3241 0.74584

8. b_treat ?0.217334 ? 0.112249 ? 1.9362 0.05285 .

9. s_int ? ?0.043733 ? 0.070711 ? 0.6185 0.53626

10. s_treat -1.535893 ? 0.147196 -10.4343 < 2e-16 ***

11. ---

12. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

13. 


14. -2 log L: 288.1408

?

1. 


2. Family: gaussian ?( identity )

3. Formula: ? ? ? ? ?y ~ treat

4. Dispersion: ? ? ? ? ~treat

5. Data: dat

6. 


7. AIC ? ? ?BIC ? logLik deviance df.resid

8. 296.1 ? ?307.6 ? -144.1 ? ?288.1 ? ? ?126

9. 


10. 


11. Conditional model:

12. Estimate Std. Error z value Pr(>|z|)

13. (Intercept) ?0.03386 ? ?0.10447 ? 0.324 ? 0.7458

14. treat ? ? ? ?0.21733 ? ?0.11225 ? 1.936 ? 0.0528 .

15. ---

16. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

17. 


18. Dispersion model:

19. Estimate Std. Error z value Pr(>|z|)

20. (Intercept) ?0.08746 ? ?0.14142 ? 0.618 ? ?0.536

21. treat ? ? ? -3.07179 ? ?0.29439 -10.434 ? <2e-16 ***

22. ---

23. Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

在這種情況下，離散度在對(duì)數(shù)方差的范圍內(nèi)，因此必須取平方的指數(shù)對(duì)數(shù)方差平方根才能檢索上述的組標(biāo)準(zhǔn)差。