Как получить стандартные ошибки из регрессии данных подсчета с нулевым раздувом? [закрыто]

9

Следующий код

PredictNew <- predict (glm.fit, newdata = Predict, X1 =X1, Y1= Y1, 
                       type = "response", se.fit = TRUE)

выдает 3 столбца data.frame--PredictNew, подогнанные значения, стандартные ошибки и остаточный член масштаба.

Идеально ... Однако, используя модель, оснащенную zeroinfl {pscl}:

PredictNew <- predict (zeroinfl.fit, newdata = Predict, X1 =X1, Y1= Y1, 
                       type = "response", se.fit = TRUE)

или

PredictNew <- predict (zeroinfl.fit, newdata = Predict, X1 =X1, Y1= Y1, 
                       type = "response", se.fit = TRUE, MC = 2500, conf = .95))

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

(Код несколько упрощен, у меня на самом деле четыре переменные и смещение - никаких проб с predict.glmи se.fit = TRUEпроизводящих SE).

KalahariKev
источник
5
Посмотрите эту ветку на R-Help: stat.ethz.ch/pipermail/r-help/2008-De December/thread.html#182806 (в частности, сообщение от Ахима Цейлиса, который предоставляет код для выполнения того, что, я думаю, вы пытаюсь сделать). Не похоже, что в эту predict()функцию были встроены стандартные ошибки zeroinfl().
Смиллиг
Спасибо, этот код, казалось, дал довольно разумные результаты. Другие должны заметить, что параметр Forex () в новой функции zeroinfl.predict для se.fit = TRUE был изменен на se = TRUE, чтобы извлечь прогнозируемые интервалы и se
KalahariKev

Ответы:

4

Насколько мне известно, predictметод для получения результатов zeroinflне включает стандартные ошибки. Если ваша цель - построить доверительные интервалы, одной из привлекательных альтернатив является использование начальной загрузки. Я считаю привлекательным, потому что у начальной загрузки есть потенциал, чтобы быть более устойчивым (с потерей эффективности, если все предположения для SE выполнены).

Вот примерный код для того, чтобы делать то, что вы хотите. Это не будет работать точно, но, надеюсь, вы сможете внести необходимые исправления.

## load boot package
require(boot)
## output coefficients from your original model
## these can be used as starting values for your bootstrap model
## to help speed up convergence and the bootstrap
dput(round(coef(zeroinfl.fit, "count"), 3))
dput(round(coef(zeroinfl.fit, "zero"), 3))

## function to pass to the boot function to fit your model
## needs to take data, an index (as the second argument!) and your new data
f <- function(data, i, newdata) {
  require(pscl)
  m <- zeroinfl(count ~ child + camper | persons, data = data[i, ], start = list(count = c(1.598, -1.0428, 0.834), zero = c(1.297, -0.564)))
  mparams <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
  yhat <- predict(m, newdata, type = "response")
  return(c(mparams, yhat))    
}

## set the seed and do the bootstrap, make sure to set your number of cpus
## note this requires a fairly recent version of R
set.seed(10)
res <- boot(dat, f, R = 1200, newdata = Predict, parallel = "snow", ncpus = 4)

## get the bootstrapped percentile CIs
## the 10 here is because in my initial example, there were 10 parameters before predicted values
yhat <- t(sapply(10 + (1:nrow(Predict)), function(i) {
  out <- boot.ci(res, index = i, type = c("perc"))
  with(out, c(Est = t0, pLL = percent[4], pUL = percent[5]))
}))

## merge CIs with predicted values
Predict<- cbind(Predict, yhat)

Я нарисовал этот код на двух написанных мной страницах: одна - параметры начальной загрузки из регрессии Пуассона с zeroinfl нулевым раздувом с нулевым раздувом Пуассона, а другая - как получить доверительные интервалы с начальной загрузкой для прогнозируемых значений из отрицательной биномиальной модели с нулевым усечением , Надеемся, что в совокупности это даст вам достаточное количество примеров, чтобы заставить его работать с прогнозируемыми значениями из нулевого раздувания Пуассона. Вы также можете получить некоторые графические идеи :)

Джошуа
источник
Я пытался адаптировать ваш код для усеченной до нуля отрицательной биномиальной модели в пакете VGAM, но получил ошибку. Должен ли я создать новый вопрос здесь в резюме и ссылку здесь? Буду очень признателен за вашу помощь с этим. В частности, это ошибка , я получаю: Error in X.vlm.save %*% coefstart : non-conformable arguments.
Рафаэль К