Почему ggplot2 95% ДИ и прогноз 95% ДИ, рассчитанные вручную, отличаются? - PullRequest
Новая
       130

Почему ggplot2 95% ДИ и прогноз 95% ДИ, рассчитанные вручную, отличаются?

0 голосов
/ 01 февраля 2019

Я хотел бы знать, почему при расчете 95% -ных доверительных полос из линейной модели смешанных эффектов ggplot2 создает более узкие полосы, чем при ручном вычислении, например, следуя здесь методу Бена Болкера доверительные интервалы при предсказаниях .То есть, дает ли ggplot2 неточное представление модели?

Вот воспроизводимый пример, использующий набор данных sleepstudy (модифицированный, чтобы быть структурно похожим на df, над которым я работаю): 

data("sleepstudy") # load dataset 
height <- seq(165, 185, length.out = 18) # create vector called height
Treatment <- rep(c("Control", "Drug"), 9) # create vector called treatment
Subject <- levels(sleepstudy$Subject) # get vector of Subject
ht.subject <- data.frame(height, Subject, Treatment) 
sleepstudy <- dplyr::left_join(sleepstudy, ht.subject, by="Subject") # Append df so that each subject has its own height and treatment
sleepstudy$Treatment <- as.factor(sleepstudy$Treatment)

Создание модели, добавление прогнозов к исходному df и построение графика 

m.sleep <- lmer(Reaction ~ Treatment*height + (1 + Days|Subject), data=sleepstudy)
sleepstudy$pred <- predict(m.sleep)
ggplot(sleepstudy, aes(height, pred, col=Treatment)) + geom_smooth(method="lm")[2]

Расчет доверительных интервалов по методу Болкера 

newdf <- expand.grid(height=seq(165, 185, 1),
                   Treatment=c("Control","Drug"))
newdf$Reaction <- predict(m.sleep, newdf, re.form=NA) 
modmat <- model.matrix(terms(m.sleep), newdf)
pvar1 <- diag(modmat %*% tcrossprod(vcov(m.sleep), modmat))
tvar1 <- pvar1+VarCorr(m.sleep)$Subject[1]
cmult <- 1.96

newdf <- data.frame(newdf
,plo = newdf$Reaction-cmult*sqrt(pvar1)
,phi = newdf$Reaction+cmult*sqrt(pvar1)
,tlo = newdf$Reaction-cmult*sqrt(tvar1)
,thi = newdf$Reaction+cmult*sqrt(tvar1))

# plot confidence intervals
ggplot(newdf, aes(x=height, y=Reaction, colour=Treatment)) + 
geom_point() +
geom_ribbon(aes(ymin=plo, ymax=phi, fill=Treatment), alpha=0.4)[2]

1 Ответ

0 голосов
/ 08 февраля 2019

С некоторыми изменениями это кажется последовательным.Доверительные интервалы действительно больше, но не намного больше.Имейте в виду, что ggplot подходит для очень другой модели;он подбирает отдельные линейные (не линейные смешанные) модели с помощью обработки, которая игнорирует (1) повторные измерения и (2) влияние дня.

Кажется странным подгонять модель со случайными наклонами, но без уровня населениянаклон (например, здесь ), поэтому я добавил фиксированный эффект Days: 

m.sleep <- lmer(Reaction ~ Treatment*height + Days +
                (1 + Days|Subject),
                data=sleepstudy)

Я немного реорганизовал код построения: 

theme_set(theme_bw())
gg0 <- ggplot(sleepstudy, aes(height, colour=Treatment)) +
    geom_point(aes(y=Reaction))+
    geom_smooth(aes(y=pred), method="lm")
  • Если вы хотите вычислить доверительные интервалы (которые были бы сопоставимы с тем, что делает lm() / ggplot2), то вам, вероятно, следует не добавить VarCorr(m.sleep)$Subject[1] к дисперсии (tvar1 переменная из примера часто задаваемых вопросов предназначена для создания интервалов прогнозирования , а не доверительных интервалов ...)
  • , так как у меня было Days в модели выше, я добавилmean(sleepstudy$Days) к кадру данных прогноза. 
newdf <- expand.grid(height=seq(165, 185, 1),
                     Treatment=c("Control","Drug"),
                     Days=mean(sleepstudy$Days))
newdf$Reaction <- newdf$pred <- predict(m.sleep, newdf, re.form=NA) 
modmat <- model.matrix(terms(m.sleep), newdf)
pvar1 <- diag(modmat %*% tcrossprod(vcov(m.sleep), modmat))
tvar1 <- pvar1
cmult <- 1.96

newdf <- data.frame(newdf
,plo = newdf$Reaction-cmult*sqrt(pvar1)
,phi = newdf$Reaction+cmult*sqrt(pvar1)
,tlo = newdf$Reaction-cmult*sqrt(tvar1)
,thi = newdf$Reaction+cmult*sqrt(tvar1))

gg0 + 
    geom_point(data=newdf,aes(y=Reaction)) +
    geom_ribbon(data=newdf,
                aes(ymin=plo, ymax=phi, fill=Treatment), alpha=0.4,
                colour=NA)

enter image description here

Сравнение с оценочными уклонами и стандартными ошибками: 

m0 <- lm(Reaction~height*Treatment,sleepstudy)
ff <- function(m) {
    print(coef(summary(m))[-1,c("Estimate","Std. Error")],digits=2)
}

> ff(m0)
##                      Estimate Std. Error
## height                   -0.3       0.94
## TreatmentDrug          -602.2     234.01
## height:TreatmentDrug      3.5       1.34

ff(m.sleep)
##                      Estimate Std. Error
## TreatmentDrug          -55.03      425.3
## height                   0.41        1.7
## Days                    10.47        1.5
## TreatmentDrug:height     0.33        2.4

Это выглядит последовательно / примерно правильно: смешанная модель дает большую стандартную ошибкуrs для склона относительно высоты и высоты: взаимодействие лечения.(Основные эффекты TreatmentDrug выглядят сумасшедшими, потому что они ожидаемые эффекты лечения при height==0 ...)

В качестве перекрестной проверки я могу получить аналогичные ответы от sjPlot::plot_model() ... 

library(sjPlot)
plot_model(m.sleep, type="pred", terms=c("height","Treatment"))

enter image description here

...