Как добавить доверительные интервалы для нулевой завышенной отрицательной биномиальной функции в ggplot2 R? - PullRequest
Новая
       134

Как добавить доверительные интервалы для нулевой завышенной отрицательной биномиальной функции в ggplot2 R?

3 голосов
/ 19 сентября 2019

У меня есть набор данных, который смоделирован как отрицательный бином с нулевым раздуванием со смешанными эффектами.Я хочу получить доверительные интервалы из прогнозов модели и построить график среднего значения и доверительных интервалов.Я попытался построить среднее значение модели.Может кто-нибудь сообщить мне, если это правильный способ сделать это?Я не знаю, как построить доверительные интервалы от моделей на ggplot2.Я хотел бы построить прогнозируемое среднее значение данных вместе с доверительным интервалом.Моя основная попытка на графике была сделана ниже в коде. 

library(pscl)
library(lmtest)

df <- data.frame(
      fertilizer = c("N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P","N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P"), 
      level = c("low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","low"), 
      repro = c(0,90,2,4,0,80,1,90,2,33,56,0,99,100,66,80,1,0,2,33,0,0,1,2,90,5,2,2,5,8,0,1,90,2,4,66,0,0,0,0,1,2,90,5,2,5,8,55)
    )    

    model <- formula(repro ∼ fertilizer + level | fertilizer * level)
    modelzinb <- zeroinfl(model, dist = "negbin", link = "logit",data =df)
    summary(modelzinb)

    df$predict <- predict(modelzinb)

    ggplot(df, aes(x=fertilizer, y=predict, color = fertilizer)) + theme_bw() +  stat_summary(aes(color = fertilizer),fun.y = mean, geom = "point", size = 4, position = position_dodge(0.1)) +
      scale_x_discrete("Fertlizer") +
      facet_wrap(.~level)

1 Ответ

1 голос
/ 20 сентября 2019

Я не уверен что вы хотели бы построить.Но я могу показать вам, как рассчитать доверительные интервалы прогнозирования для вашей отрицательно-биномиальной модели с нулевым раздувом.

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

  1. Давайте начнем с подгонки модели к вашим данным. 

    library(pscl)
model <- formula(repro ~ fertilizer + level | fertilizer * level)
modelzinb <- zeroinfl(model, dist = "negbin", link = "logit", data = df)

  2. Мы можем использовать predict, чтобы получить модельные оценки для ответа.

    resp <- predict(modelzinb)

    Обратите внимание, что в вашей модели с нулевым раздувом NB predict.zeroinfl (по умолчанию) возвращает оценочный средний ответ по шкале наблюдаемых подсчетов repro.

    Относительно достоверностиВ интервалах (CI) проблема заключается в том, что precict.zeroinfl не имеет аргумента interval, который позволяет напрямую рассчитывать CI.Примечание: кажется, что pscl::zeroinfl использовал для включения этой функции, см. Документацию для версии 0.54 .Возможно, стоило бы связаться с разработчиком пакетов по этому поводу.

    Суть в том, что мы должны найти другой способ расчета прогнозных КИ.Для этого мы используем начальную загрузку.Библиотека R boot предоставляет все необходимые функции и инструменты для этого.

  3. Для начала мы можем использовать функцию boot::boot для генерации загрузочных копий предсказанных ответов. boot::boot нужна функция, которая генерирует интересующее количество (здесь прогнозируемый ответ) на основе данных.Итак, сначала мы определим такую ​​функцию stat.В этом конкретном случае stat должен принимать два аргумента: первый аргумент - это исходные данные, второй аргумент - вектор индексов строк (которые определяют случайную выборку начальной загрузки данных).Затем функция будет использовать загруженные данные, чтобы соответствовать модели, а затем использовать модель для прогнозирования среднего отклика на основе полных данных. 

    stat <- function(df, inds) {
    model <- formula(repro ~ fertilizer + level | fertilizer * level)
    predict(
        zeroinfl(model, dist = "negbin", link = "logit", data = df[inds, ]),
        newdata = df)
}

    Чтобы обеспечить воспроизводимость результатов, мы устанавливаем фиксированное случайное значениеЗаполнить и сгенерировать 100 выборок начальной загрузки для оценки среднего ответа 

    set.seed(2018)
library(boot)
res <- boot(df, stat, R = 100)

    Выходной объект res представляет собой list, который содержит оценку среднего ответа для полных данных в элементе t0 (проверьте с помощью all.equal(res$t0, predict(modelzinb))) и расчетный средний отклик для выборок начальной загрузки в элементе t (который представляет собой матрицу измерения R x nrow(df), где R - количество выборок начальной загрузки).

  4. Все, что осталось сделать, - это рассчитать доверительный интервал из оцененных средних ответов на основе моделей, установленных на выборках начальной загрузки.Для этого мы можем использовать удобную функцию boot::boot.ci.Функция позволяет рассчитывать различные типы КИ (базовый, нормальный, BCa и т. Д.).Здесь мы используем "basic" в демонстрационных целях.Я не утверждаю, что это лучший выбор.

    boot::boot.ci принимает аргумент index, который соответствует вводу оцененного среднего вектора ответа.Фактические интервалы хранятся в последних 2 столбцах (столбцы 4 и 5) матрицы, которая хранится как элемент list возвращаемого объекта boot.ci. 

    CI <- setNames(as.data.frame(t(sapply(1:nrow(df), function(row)
    boot.ci(res, conf = 0.95, type = "basic", index = row)$basic[, 4:5]))),
    c("lower", "upper"))

  5. Сейчасмы закончили, и мы можем связать столбцы CI с исходными данными, включая предсказанные средние значения ответа. 

    df_all <- cbind.data.frame(df, response = predict(modelzinb), CI)
#   fertilizer level repro response      lower    upper
#1           N   low     0 29.72057   5.876731 46.48165
#2           N   low    90 29.72057   5.876731 46.48165
#3           N  high     2 21.99345 -15.228956 38.86421
#4           N  high     4 21.99345 -15.228956 38.86421
#5           N   low     0 29.72057   5.876731 46.48165
#6           N   low    80 29.72057   5.876731 46.48165
#7           N  high     1 21.99345 -15.228956 38.86421
#8           N  high    90 21.99345 -15.228956 38.86421
#9           N   low     2 29.72057   5.876731 46.48165
#10          N   low    33 29.72057   5.876731 46.48165
#11          N  high    56 21.99345 -15.228956 38.86421
#12          N  high     0 21.99345 -15.228956 38.86421
#13          P   low    99 24.22626  -9.225827 46.17656
#14          P   low   100 24.22626  -9.225827 46.17656
#15          P  high    66 22.71826   2.595246 39.88333
#16          P  high    80 22.71826   2.595246 39.88333
#17          P   low     1 24.22626  -9.225827 46.17656
#18          P   low     0 24.22626  -9.225827 46.17656
#19          P  high     2 22.71826   2.595246 39.88333
#20          P  high    33 22.71826   2.595246 39.88333
#21          P   low     0 24.22626  -9.225827 46.17656
#22          P   low     0 24.22626  -9.225827 46.17656
#23          P  high     1 22.71826   2.595246 39.88333
#24          P  high     2 22.71826   2.595246 39.88333
#25          N   low    90 29.72057   5.876731 46.48165
#26          N   low     5 29.72057   5.876731 46.48165
#27          N  high     2 21.99345 -15.228956 38.86421
#28          N  high     2 21.99345 -15.228956 38.86421
#29          N   low     5 29.72057   5.876731 46.48165
#30          N   low     8 29.72057   5.876731 46.48165
#31          N  high     0 21.99345 -15.228956 38.86421
#32          N  high     1 21.99345 -15.228956 38.86421
#33          N   low    90 29.72057   5.876731 46.48165
#34          N   low     2 29.72057   5.876731 46.48165
#35          N  high     4 21.99345 -15.228956 38.86421
#36          N  high    66 21.99345 -15.228956 38.86421
#37          P   low     0 24.22626  -9.225827 46.17656
#38          P   low     0 24.22626  -9.225827 46.17656
#39          P  high     0 22.71826   2.595246 39.88333
#40          P  high     0 22.71826   2.595246 39.88333
#41          P   low     1 24.22626  -9.225827 46.17656
#42          P   low     2 24.22626  -9.225827 46.17656
#43          P  high    90 22.71826   2.595246 39.88333
#44          P  high     5 22.71826   2.595246 39.88333
#45          P   low     2 24.22626  -9.225827 46.17656
#46          P   low     5 24.22626  -9.225827 46.17656
#47          P  high     8 22.71826   2.595246 39.88333
#48          P   low    55 24.22626  -9.225827 46.17656
#...

Образец данных 

df <- data.frame(
    fertilizer = c("N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P","N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P"),
    level = c("low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","low"),
    repro = c(0,90,2,4,0,80,1,90,2,33,56,0,99,100,66,80,1,0,2,33,0,0,1,2,90,5,2,2,5,8,0,1,90,2,4,66,0,0,0,0,1,2,90,5,2,5,8,55))
...