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.