outbreaker2
outbreaker2 copied to clipboard
interhosp to do
- [x] potential_colonised movement (Gibbs sampling?) @mirjamlaager
- [x] poisson_scale movement @PascalCrepey
- [x] poisson_scale prior (gamma?) @PascalCrepey
- [x] setting up create_param for new parameters @finlaycampbell
- [x] setting up create_config for new parameters @finlaycampbell
- [ ] first run of outbreaker on dataset @ClementMassonnaud
Github apparently doesn't let you assign issues to users that don't have push access, so I've just tagged you in the relevant issues for now. Tomorrow I'll move everything to my fork and then give you push access to that, then you don't have to do pull requests every time!
So I've gotten a first version working - it appears to be converging, which is a good start. We still need to do a bunch of validation to make sure things are behaving as expected, but I tested the cpp likelihoods against R implementations and they gave the same results. Below is a run on a sample dataset:
dates <- seq.int(0, 100, by = 25)
w <- dnorm(1:50, 25, 10)
n_cases <- as.integer(runif(6, 0, 25))
hosp_matrix <- matrix(runif(36), 6, 6)
diag(hosp_matrix) <- 0
colnames(hosp_matrix) <- rownames(hosp_matrix) <- 1:6
data <- outbreaker_data(dates = dates,
w_dens = w,
n_cases = n_cases,
hosp_matrix = hosp_matrix)
config <- create_config(init_potential_colonised = rep(2, 5),
sd_potential_colonised = 5,
prior_poisson_scale = c(5, 2),
pb = TRUE,
find_import = FALSE,
data = data)
res <- outbreaker(data, config)
plot(res)
@PascalCrepey @ClementMassonnaud @mirjamlaager
Great ! I quickly looked at the code, you're not estimating the "scale of the Poisson" for now, right ? next step is trying it on the simulated data ! ;-)
Ah you're right! Now estimating poisson_scale as of 83ddf91ef53db63d13fb79bd4c2e03f1080f6b21
When testing on datasets created by @ClementMassonnaud, the value of poisson_scale keep on rising with the number of iterations. As we don't have priors on its real value, we are doing some tests with @PascalCrepey to see the influence on results of fixing it at 1.
Yes that makes sense, I think this was always going to be a parameter that we either fix or put very strong priors on - it's really just making our assumptions about the relationship between observed and unobserved cases explicit. However, a continuously rising parameter value suggests the MCMC is either not converging or we have a bug - I would imagine it's the former, and that the proposal distribution for poisson_scale is too narrow. I would suggest increasing the value of sd_poisson_scale from 0.1 by doing create_config(sd_poisson_scale = 2)
or whatever value works. If the parameter is simply poorly defined then it should then be jumping up and down across a wide range of values, not steadily increasing or decreasing.