Heteroscedastic residuals plot not matching plot shown in Pinhiero and Bates

user_4727699 0
Viewed 288

I am having trouble making a plot of standarised residuals versus a covariate match a plot shown in Pinhiero and Bates Mixed-Effects Models in S and S-Plus. The model being plotted is a general formulation of a nonlinear mixed-effects model, contained in the nlme package

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

When we plot the standardised residuals versus the transmembrane pressure in this model

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

The resulting plot shows evidence of heteroscedasticity across different pressures. Thus we fit a new model with a power variance function to account for this.

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

Which is clearly superior to the first model

anova(fm1Dial.gnls, fm2Dial.gnls)

However when we plot the standarised residuals versus transmembrane pressure for the new improved model

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

The plot doesn't look like much of an improvement on the first plot, with the vertical spread of the residuals still appearing to be much higher at higher pressures.

The plot for the second improved model in Pinhiero and Bates, however. shows a similar vertical spread of residuals at all levels of pressure, which makes sense given that the heteroscedasticity has explicitly been accounted for in this improved model.

What am I doing wrong?

nlme r
1 Answers

Where you were wrong is saying that

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

are standardized residuals while they indeed are not. You correctly found that 

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

or, more completely,

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

gives the same plot as in the book and those are standardized residuals.

?residuals.gnls explains quite a bit:

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".

From this description we also see why choosing type as "normalized" and "pearson" gives the same result: the former option would take into account the dependence structure of the error, but since we only relaxed the homoskedasticity assumption, we still have no dependence. That is also evident in nlme:::residuals.gnls in

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"
        }
    }
}
user_13205350
Related