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

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

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

`mean(dat\$y)[1] 2.999048sd(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)))`

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

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

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

## 异方差回归模型

`Call:Residuals:    Min      1Q  Median      3Q     Max-2.8734 -0.5055 -0.0287  0.4231  3.4097Coefficients:            Estimate Std. Error t value Pr(>|t|)(Intercept)  0.03386    0.09298   0.364    0.716treat        0.21733    0.19355   1.123    0.264Residual 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`

`Maximum likelihood estimationCall:(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`

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

`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`

`par(mfrow = c(1, 1))`

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

`(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      126Conditional 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 ‘ ’ 1Dispersion 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`

0 条评论

• ### R语言中使用线性模型、回归决策树自动组合特征因子水平

没有定义一个（连续的）因变量，没有定义一个连续的协变量，也没有定义一个分类变量，此处有十个级别。我们可以使用

• ### Python贝叶斯回归分析住房负担能力数据集

我想研究如何使用pymc3在贝叶斯框架内进行线性回归。根据从数据中学到的知识进行推断。

• ### 用Python的Numpy求解线性方程组

解决线性方程组的最终目标是找到未知变量的值。这是带有两个未知变量的线性方程组的示例：

• ### java的StringBuffer可变字符串的追加及修改、查找

-----------java的StringBuffer可变字符串的追加及修改--------------

• ### java中的lastIndexOf( )函数是什么意思

String中的lastIndexOf方法，是获取要搜索的字符、字符串最后次出现的位置。

• ### 二进制的骚操作，省字段，省带宽，提效率...

上一个礼拜和一个同事对接口，前端同事问我是不是接口文档写错了，一个订单的异常标签有多个，不应该返回一个数组吗？为啥只返回了一个数字。

• ### VUE开发一个组件——移动端弹出层（IOS版）

再造一轮，vue移动端弹出层，包括confrim询问框，tips提示框，popup选择器等。

• ### WhylSon：在Why3中证明你的 Michelson 智能合约的重要性（CS PL）

本文介绍了 WhylSon，这是一个用 Michelson 编写的智能合约的推理验证工具，Michelson 是 Tezos 区块链的底层语言。WhylSon ...

• ### Power Query中数据分割函数详解(4)

Table.SplitColumn(table as table, sourceColumn as text,splitter as function,opti...