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):

enter image description here

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

enter image description here

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

enter image description here

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) {
  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]))


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


# 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))



## 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:

enter image description here

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.

enter image description here


1 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 ?

  • $\begingroup$ Yes, exactly - I'm looking to compare an alternative therapy to the standard of care, however, I only have ~14 weeks of data for the alternative therapy while I have a broader data set for the standard of care. Can you provide some brief guidance as to how the prior approach should change given this is a rating and I am trying to infer how this metric would change over time? $\endgroup$ Jan 30, 2019 at 21:32
  • $\begingroup$ By comparing the two therapies, you have now entered the realm of inference, which is a little tricker than simulation. From your data it looks like you follow multiple patients over time, which probably calls for a longitudinal analysis. I can't say more without knowing about your study design, goals, data collection, etc. My advice is to make friends with a statistician at your institution. Sorry I can't be more help. $\endgroup$ Jan 30, 2019 at 22:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.