r - 异方差残差图与 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"
        }
    }
}