Update a zero-inflated Poisson model to adjust model predictions

Question

I am trying to model out how a clinical metric declines over time with various therapies.

I'm a bit new to R and statistics, so appreciate the patience and help.

I have two data sets - the first a longer-term view of the clinical metric (x-axis is in days):

The second a shorter term view of a different trial (here the x-axis is weeks, both y-axis are the same clinical metric):

I've fitted a zero-inflated Poisson, which by testing with vuong() shows that the zero-inflated Poisson is a better fit.

Below is my code and additional information for the model followed by 2 issues I am having:

1) Distribution of BVAS scores showing zero-inflation

2) code for zero-inflated Poisson

summary(RAVE_CYC_PR3_reg2 <- zeroinfl(BVAS ~ Days, data = RAVE_CYC_PR3))  

vuong(RAVE_CYC_PR3_reg, RAVE_CYC_PR3_reg2)

dput(coef(RAVE_CYC_PR3_reg2, "count"))
dput(coef(RAVE_CYC_PR3_reg2, "zero"))

f <- function(data, i) {
  require(pscl)
  m <- zeroinfl(BVAS ~ Days, data = data[i,],
                start = list(count = c(1.839, -0.008), zero = c(-0.078, 0.006)))
                as.vector(t(do.call(rbind, coef(summary(m)))[,1:2]))
}

set.seed(10)

res <- boot(RAVE_CYC_PR3, f, R = 1200, parallel = "snow", ncpus = 4)

print(res)

# basic parameter estimates with percentile and bias adjusted CIs
parms <- t(sapply(c(1, 3, 5, 7), function(i) {
  out <- boot.ci(res, index = c(i, i + 1), type = c("perc", "bca"))
  with(out, c(Est = t0, pLL = percent[4], pUL = percent[5],
              bcaLL = bca[4], bcaLL = bca[5]))
}))


# add row names
row.names(parms) <- names(coef(RAVE_CYC_PR3_reg2))

parms

confint(RAVE_CYC_PR3_reg2)

## exponentiated parameter estimates with percentile and bias adjusted CIs
expparms <- t(sapply(c(1, 3, 5, 7), function(i) {
  out <- boot.ci(res, index = c(i, i + 1), type = c("perc", "bca"), h = exp)
  with(out, c(Est = t0, pLL = percent[4], pUL = percent[5],
              bcaLL = bca[4], bcaLL = bca[5]))
}))

The model seems to fit the data relatively well:

However, I have 2 issues:

1) When I apply my model to predict scores over time I get an output that is well above the clinical range [0-36]

predict_at <- seq(0, 365, 30)
means <- exp(predict(RAVE_CYC_PR3_reg2, newdata = list(Days = predict_at)))
simulated_data <- map(means, ~rpois(500, .x))
simd <- data_frame(times = predict_at, 
                   y = simulated_data) %>% unnest()

lined = data_frame(times = predict_at, y  = means)

simd %>%
  ggplot(aes(times, y)) +
  geom_point(color = "blue", alpha = 0.15) +
  geom_line(data=lined, color="red")

I'm surprised that the output is so far from any of the initial data, and I'd be curious how to better bound the predictions

2) I'd like to use the second image posted (with the dark red lines) to use this data to update the model so that the final model can accommodate a view on how an alternative therapy may more rapidly reduce scores to 0, while still having some real data beyond the 12 week period. From a statistical point, I'm not sure what is the most appropriate way to do this - do I just include this additional data when fitting the model? Can I weight this new data to a greater degree that the initial data? Lastly, I'm still quite new to R so the implementation of adding data with a particular weight is very unclear to me.

Thank you all for the input, truly appreciated.

Answer 1

A couple things:

1) If your data are bounded for good reason (e.g. your data is a rating, and 36 is the highest the rating can go), then Poisson regression is not appropriate. You should use something like Beta Regression or logistic regression, and transform your outcome to a number between 0 and 1 (1 being the largest the outcome can be). This can be done by dividing the observation by the max (so y/36).

2) Now this is sounding like you want to do inference, which is a lot different than just simulation. Are you interested in examining the effect of an alternative therapy as compared to standard of care ?