Оценить SE для всех уровней факторов с нулевой инфляцией модели

Question

У меня довольно сложная модель ZINB. Я попытался повторить основную структуру того, что я пытаюсь сделать:

MyDat<-cbind.data.frame(fac1 = rep(c("A","B","C","D"),10), 
     fac2=c(rep("X",20),rep("Y",20)), 
     offset=c(runif(20, 50,60),runif(20,150,165)), 
     fac3=rep(c(rep("a1",4),rep("a2",4),rep("a3",4),rep("a4",4),rep("a5",4)),2), 
     Y=c(0,0,0,1,0,0,11,10,0,0,0,5,0,0,0,35,60,0,0,0,0,2,0,0,16,0,0,0,0,0,3,88,0,0,0,0,0,0,27,0)) 

f<-formula(Y~fac1+ offset(log(offset))|fac3+ fac2) 
ZINB <-zeroinfl(f, dist = "negbin",link = "logit", data = MyDat) 
summary(ZINB) 

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

Вот результат:

Call: 
zeroinfl(formula = f, data = MyDat, dist = "negbin", link = "logit") 

Pearson residuals: 
     Min  1Q Median  3Q  Max 
-0.418748 -0.338875 -0.265109 -0.001566 2.682920 

Count model coefficients (negbin with log link): 
      Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.7192  0.9220 -1.865 0.062239 . 
fac1B  -4.4161  1.4700 -3.004 0.002663 ** 
fac1C  -1.2008  1.2896 -0.931 0.351778  
fac1D   0.1928  1.3003 0.148 0.882157  
Log(theta) -1.7349  0.4558 -3.806 0.000141 *** 

Zero-inflation model coefficients (binomial with logit link): 
      Estimate Std. Error z value Pr(>|z|) 
(Intercept) -11.5899 210.8434 -0.055 0.956 
fac3a2  -0.4775  2.4608 -0.194 0.846 
fac3a3  -11.2284 427.5200 -0.026 0.979 
fac3a4  10.7771 210.8056 0.051 0.959 
fac3a5  -0.3135  2.3358 -0.134 0.893 
fac2Y  11.8292 210.8298 0.056 0.955 
--- 
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 0.1764 
Number of iterations in BFGS optimization: 76 
Log-likelihood: -63.82 on 11 Df 

Я консультировался документы и статистика книг и форумов, но я до сих пор не знаю, как представить эту информацию. То, что я действительно хочу, это штриховой график, показывающий эффекты на оси Y и на 4 уровнях на X.

Если я правильно понял, уровень A fac1 в настоящее время установлен на 0 и является моим эталонным уровнем (пожалуйста исправьте меня, если я ошибаюсь здесь). Итак, я могу сделать график из 4 уровней (включая уровень А как ноль). Это не кажется идеальным. Мне бы очень хотелось иметь 95% CI для всех уровней.

Я также могу использовать функцию прогнозирования, однако pred.zeroinfl не дает оценок ошибок, и я не уверен, как интерпретировать эффект смещения.

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

Ниже приведен код, и сюжет, чтобы создать предсказанные значения:

MyDat$phat<-predict(ZINB, type="response") 
MyDat$phat_os<-MyDat$phat/MyDat$offset 

plot(phat~fac1, MyDat) 

Predictions plot

ли самонастройки пути пойти? Я пробовал это и сталкивался со всеми неприятностями, потому что я не уверен, что это необходимо.

Заранее благодарю вас, и, пожалуйста, пройдите на меня, если я сделаю глупый недосмотр/допущение. Я все еще учусь, но эта статистика немного от меня доходит.

Answer 1

Для начала вы можете рассчитать коэффициенты модели с их доверительными интервалами. Пакет arm имеет функцию coefplot, но у него нет метода для моделей zeroinfl, поэтому я создал простой график коэффициентов ниже, используя ggplot2. Метод predict для моделей zeroinfl не дает доверительных интервалов для прогнозов, но this answer to a question on CrossValidated показывает, как построить загрузочные доверительные интервалы для моделей zeroinfl.

Что касается уровней fac1: A является опорным уровнем, поэтому коэффициенты для других уровней относятся к fac1 = "A".

library(pscl) 
library(ggplot2) 

MyDat<-cbind.data.frame(fac1 = rep(c("A","B","C","D"),10), 
         fac2=c(rep("X",20),rep("Y",20)), 
         offset=c(runif(20, 50,60),runif(20,150,165)), 
         fac3=rep(c(rep("a1",4),rep("a2",4),rep("a3",4),rep("a4",4),rep("a5",4)),2), 
         Y=c(0,0,0,1,0,0,11,10,0,0,0,5,0,0,0,35,60,0,0,0,0,2,0,0,16,0,0,0,0,0,3,88,0,0,0,0,0,0,27,0)) 

f<-formula(Y ~ fac1 + offset(log(offset))|fac3 + fac2) 
ZINB <-zeroinfl(f, dist = "negbin",link = "logit", data = MyDat) 

# Extract coefficients and standard errors from model summary 
coefs = as.data.frame(summary(ZINB)$coefficients$count[,1:2]) 
names(coefs)[2] = "se" 
coefs$vars = rownames(coefs) 

# Coefficient plot 
ggplot(coefs, aes(vars, Estimate)) + 
    geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") + 
    geom_errorbar(aes(ymin=Estimate - 1.96*se, ymax=Estimate + 1.96*se), 
       lwd=1, colour="red", width=0) + 
    geom_errorbar(aes(ymin=Estimate - se, ymax=Estimate + se), 
       lwd=2.5, colour="blue", width=0) + 
    geom_point(size=4, pch=21, fill="yellow") + 
    theme_bw() 

И вот как выглядит сюжет.