График гетероскедастических остатков, не совпадающий с графиком Пиньеро и Бейтса - PullRequest
Новая
       16

График гетероскедастических остатков, не совпадающий с графиком Пиньеро и Бейтса

0 голосов
/ 30 декабря 2018

У меня возникли проблемы при построении графика стандартизированных остатков в сравнении с ковариатным сопоставлением графика, показанного в Пиньеро и Бейтсе. Модели смешанных эффектов в 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)

График не выглядит значительным улучшением на первом графике, поскольку вертикальный разброс остаточных значений все еще выглядит значительно выше при более высоких давлениях.

Сюжет для второй улучшенной модели в Пиньеро и Бейтсе, однако.показывает аналогичный вертикальный разброс остатков на всех уровнях давления, что имеет смысл, учитывая, что гетероскедастичность была явно учтена в этой улучшенной модели.

Что я делаю неправильно?

1 Ответ

0 голосов
/ 30 декабря 2018

Там, где вы ошиблись, говорится, что 

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 --- необязательная символьная строка, указывающая тип используемых остатков.Если «ответ», используются «необработанные» остатки (наблюдаемые - установлены);иначе, если «Пирсон», используются стандартизированные остатки (необработанные остатки, деленные на соответствующие стандартные ошибки);иначе, если "нормализовано", используются нормализованные остатки (стандартизованные остатки, предварительно умноженные на обратный квадратный корень из оцененной матрицы корреляции ошибок).Используется частичное сопоставление аргументов, поэтому необходимо указать только первый символ.По умолчанию используется «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"
        }
    }
}
...