异方差残差图与 Pinhiero 和 Bates 中显示的图不匹配

Question

我在制作标准化残差图与协变量图匹配 Pinhiero 和 Bates S 和 S-Plus 中的混合效应模型 中显示的图时遇到了问题。正在绘制的模型是非线性混合效应模型的一般公式，包含在 nlme 包

中

library(nlme)
options(contrasts = c("contr.helmert", "contr.poly"))
fm1Dial.gnls <- gnls(rate ~ SSasympOff(pressure, Asym, lrc, c0),
                     data = Dialyzer,
                     params = list(Asym + lrc ~ QB, c0 ~ 1),
                     start = c(53.6, 8.6, 0.51, -0.26, 0.225))

当我们在这个模型中绘制标准化残差与跨膜压力的关系图时

plot(fm1Dial.gnls, resid(.) ~ pressure, abline = 0)

结果图显示了不同压力下的异方差性证据。因此，我们拟合了一个具有幂方差函数的新模型来解决这个问题。

fm2Dial.gnls <- update(fm1Dial.gnls, weights = varPower(form = ~ pressure))

明显优于第一个模型

anova(fm1Dial.gnls, fm2Dial.gnls)

然而，当我们绘制新改进模型的标准化残差与跨膜压力的关系图时

plot(fm2Dial.gnls, resid(.) ~ pressure, abline = 0)

该图看起来与第一个图相比没有太大改进，残差的垂直分布在较高压力下仍然显得高得多。

然而，Pinhiero 和 Bates 中第二个改进模型的情节。显示了在所有压力水平下类似的残差垂直分布，鉴于此改进模型中已明确考虑了异方差性，这是有道理的。

我做错了什么？

Answer 1

你说错了

plot(fm2Dial.gnls, resid(.) ~ pressure, abline = 0)

是标准化残差，但实际上不是。您正确地发现

plot(fm2Dial.gnls, resid(., type = "p") ~ pressure, abline = 0)

或者，更完整地说，

plot(fm2Dial.gnls, resid(., type = "pearson") ~ pressure, abline = 0)

给出了与书中相同的图，这些是标准化残差。

?residuals.gnls 解释了很多：

type --- an optional character string specifying the type of residuals to be used. If "response", the "raw" residuals (observed - fitted) are used; else, if "pearson", the standardized residuals (raw residuals divided by the corresponding standard errors) are used; else, if "normalized", the normalized residuals (standardized residuals pre-multiplied by the inverse square-root factor of the estimated error correlation matrix) are used. Partial matching of arguments is used, so only the first character needs to be provided. Defaults to "response".

从这个描述中我们还可以看出为什么选择 type 作为 "normalized" 和 "pearson" 会得到相同的结果：前一个选项会考虑 依赖性 结构错误，但由于我们只是放宽了同方差假设，我们仍然没有依赖性。这在 nlme:::residuals.gnls 和

中也很明显

if (type != "response") {
    val <- val/attr(val, "std")
    lab <- "Standardized residuals"
    if (type == "normalized") {
        if (!is.null(cSt <- object$modelStruct$corStruct)) {
            val <- recalc(cSt, list(Xy = as.matrix(val)))$Xy[, 
              1]
            lab <- "Normalized residuals"
        }
    }
}

异方差残差图与 Pinhiero 和 Bates 中显示的图不匹配

Heteroscedastic residuals plot not matching plot shown in Pinhiero and Bates

r

nlme