专栏首页拓端tecdatR语言异方差回归模型建模:用误差方差解释异方差
原创

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

原文链接:http://tecdat.cn/?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.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)))

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

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

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

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

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

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

coef(lm(y ~ 1, dat))(Intercept)   2.999048sigma(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.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

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

首先,我记录一下重新创建OLS模型:

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

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

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

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

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

治疗效果大致相同,但现在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      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

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

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

我来说两句

0 条评论
登录 后参与评论

相关文章

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

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

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

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

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

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

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

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

    用户7886150
  • java中的lastIndexOf( )函数是什么意思

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

    黑泽君
  • 二进制的骚操作,省字段,省带宽,提效率...

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

    Java识堂
  • VUE开发一个组件——移动端弹出层(IOS版)

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

    Javanx
  • 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...

    逍遥之
  • 数据库事务总结

    SuperHeroes

扫码关注云+社区

领取腾讯云代金券