-
Notifications
You must be signed in to change notification settings - Fork 6
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Cohort based model #221
Comments
This comment was marked as outdated.
This comment was marked as outdated.
Work on this is here: https://github.com/epinowcast/primarycensoreddist |
Steps here IMO:
|
Using > summary(fit_direct)
Family: lognormal
Links: mu = identity; sigma = identity
Formula: delay_daily ~ 1
Data: sim_obs (Number of observations: 500)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.77 0.02 1.73 1.82 1.00 3457 2386
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.49 0.02 0.46 0.52 1.00 3250 2531
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
> summary(fit_direct_weighted)
Family: lognormal
Links: mu = identity; sigma = identity
Formula: delay | weights(n) ~ 1
Data: cohort_obs (Number of observations: 21)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.77 0.02 1.73 1.82 1.00 3542 2778
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.49 0.02 0.46 0.52 1.00 3094 2518
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1). |
int check_for_analytical {
if (dist_id == 2 && primary_id == 1) return 1; // Gamma delay with Uniform primary
if (dist_id == 1 && primary_id == 1) return 1; // Lognormal delay with Uniform primary
if (dist_id == 3 && primary_id == 1) return 1; // Weibull delay with Uniform primary
return 0; // No analytical solution for other combinations
} So this allows |
For the
Notes here:
|
for (n in 1:N) {
target += weights[n] * (primarycensored_lognormal_uniform_lcdf_lpdf(Y[n] | mu[n], sigma, pwindow));
} Have this issue. Even if define Also not currently putting |
This comment was marked as outdated.
This comment was marked as outdated.
|
This comment was marked as outdated.
This comment was marked as outdated.
@seabbs @SamuelBrand1 tough ask but I don't suppose either of you have any guesses about what might be going wrong here to cause compilation issues? The Stan code is:
And the compilation error is: Compiling Stan program... /Users/adamhowes/.cmdstan/cmdstan-2.35.0/stan/src/stan/model/model_base_crtp.hpp:93:50: note: in instantiation of function template specialization 'model_6fa90462cf3e16ffa43f0361963d4387_model_namespace::model_6fa90462cf3e16ffa43f0361963d4387_model::log_prob<false, false, double>' requested here /var/folders/4r/hkp4v9fn3wx044_hk8qnsjxw0000gn/T/RtmptN5cBH/model-111232eafcdd.hpp:1799:47: error: non-constant-expression cannot be narrowed from type 'size_t' (aka 'unsigned long') to 'int' in initializer list [-Wc++11-narrowing] 4 errors generated. make: *** [/var/folders/4r/hkp4v9fn3wx044_hk8qnsjxw0000gn/T/RtmptN5cBH/model-111232eafcdd.o] Error 1 Error: An error occured during compilation! See the message above for more information. |
Looks like |
Not sure what that means. So I guess what we want is to input empty which I mentioned in an issue (epinowcast/primarycensored#170) doesn't seem to work. So I guess to debug I should try different ways to specify Perhaps having a zero length vector into Stan is just not going to work? Wonder if you have vignettes over at Edit: I can see you have:
|
That might be the problem; worth checking |
I think you can't create params inside your wrapper function as it looks like it wants it to be defined in the model block? Seems like a weird edge case but a place to start |
I think to use the custom family functionality of |
Good news! My
Suggests to me it's some kind of install / compiler issue I have (rather than an issue with the model). I have tried:
The result of: > cmdstanr::cmdstan_make_local()
[1] "CXXFLAGS += -Wno-deprecated-declarations" suggests my Maybe something similar to acfr/snark#117 |
I'm a bit confused about how to solve this. I'd be interested @seabbs or @SamuelBrand1 if you try running the code in https://github.com/epinowcast/epidist/pull/426/files#diff-7080e507dc8e86af57fa62d8ef779fa1a108c00ff88a09b526df205512881bca and it might work for you if you are able to run the vignette in |
* Add cohort model template * Fix to previous commit (name as cohort_model) * Generate simulated cohort data * Add unweighted and weighted direct models * Thinking about custom family for pcd function * functions component of stanvars * Add transformed parameters for cohort model * Progress on implementing PCD model * Set q as vector * Get rid of params input * This would work, apart from it's the CDF. For the PMF need to import many primarycensored functions.. * Almost working with "import all functions" strategy * Wrap up on attempt * Rename to marginal model * Move towards single wrap function with others imported * Just running into some C++ errors.. * This doesn't change anything * Move to marginal model name and lint * Create aggregate data inside model conversion function for now * First draft on moving marginal_model into functions * Run document * Tests working up to valid Stan code * Regex version of marginal model * Use prep_marginal_obs * Improve assert for marginal model * Add pkgdown and document * Clean up scratch implementation * remove scratch file * update data format, formula, and family * update stan code * basic working version * add transformm data methods * start exploring the ebola example * add a helper to find meaningful relative_obs_times * get the full ebola vignette working with new variable requirements * improve return messages * update approx vignette * get ebole vignette passing by checking pp and related inputs * add marginal model integration tests * expand post processing tests * add marginal model * use the right transform data keyword * add ... pass through to make constructors correct * fix .summarise_n_by_formula test so error message is as expected * drop not required .row_id * check using ... properly * make the progress messages prettier for reducing data complexit: * check post process tests again * add a test for the specific transform data method * add some tests for the generic transform data method * put transform data tests in the correct folder * add a news update * change vignette language to talk about marginal model * update the FAQ to use the marginal variables * call it transformed_data not trans_data * change the error message to make it clear its a epidist limitation * update stan docs * Update NEWS.md * Update NEWS.md Co-authored-by: Adam Howes <[email protected]> * Update NEWS.md Co-authored-by: Adam Howes <[email protected]> * Update inst/stan/latent_model/functions.stan Co-authored-by: Adam Howes <[email protected]> * Update setup.R --------- Co-authored-by: Sam <[email protected]>
Done in #426 |
It would be nice to offer a model that scales with unique delays and strata versus by individual delay data points.
One approach for this is estimating a correct double censored PMF (by integrating the primary event distribution and the delay distribution to get a CDF and then taking the differences to get the secondary event PMF) for each observation time and then weighting the log pmf for each delay by the number of observed for that delay + obs time combination.
A julia prototype for this method is here: CDCgov/Rt-without-renewal#386 (comment)
For the simple single strata case this is quite simple but as strata become more numerous (and arbitrary) it becomes more of an engineering challenge. I plan to do this in the simple case for stan in
EpiNow2
(epiforecasts/EpiNow2#350) which should hopefully serve as a guide.We then need to think what this looks like within a brms implementation. I think it might involve injecting transformed parameters. This could be made simpler if brms has a cohort based model which it may have.
For anything beyond daily censoring this approach becomes quite complicated and so would need some more thought.
The text was updated successfully, but these errors were encountered: