-
Notifications
You must be signed in to change notification settings - Fork 2
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
Validate ABC-SMC estimation for compartmental (SIR) models #42
Comments
This is more of a minor annoyance, but I'm unsure about the format of variable
Looking into
This code was derived from In I figured when we run the |
This is a bug in rcolgem that I fixed a while back but didn't make it into this release for some reason. We need to replace: if (!is.na(sampleStates)) with if (all(!is.na(sampleStates)) This should result in a condition of length 1. We could make this change to our fork of rcolgem and direct people to that version. |
I've simulated a target tree, and was able to calculate kernel distances for a varying parameter.
The traceback error:
I noticed that when I calculated kernel distances for varying
One possible reason for this may be how I've set my prior distribution for |
and then checks
births = (ODE expr) therefore m = 1 This is okay for the sir.nondynamic model but not for the rest, since I'm including compartment R in the ODE expressions as well as compartment E for the seir model. This changes length(x0) to either 3 or 4. Before:
to:
I would also have to change births to |
We don't want to incorporate the recovered compartment in the rcolgem expressions, so change this: deaths <- rbind(c('parms$gamma * I - parms$mu * R')) to this: deaths <- rbind(c('parms$gamma * I')) and exclude |
Update on the following error:
Parameter (
fgy[[1]], # time axis of ODE solution
fgy[[2]], # births
fgy[[3]], # migrations
fgy[[4]], # deme sizes
fgy[[5]] # ODE solution
)
Plot for `x=fgy[[2]], y=fgy[[1]]`
I changed |
There's something funny going on with the kernel scores here.
Can you please check whether they're getting normalized correctly? I've run into a similar problem before, I'll dig around and make an issue.
… On Jun 28, 2017, at 12:46 PM, Tammy Ng ***@***.***> wrote:
Individual validation of each parameter:
Target tree parameters (true values):
theta <- c(t.end=50, N=5000, beta=0.1, gamma=1/520, mu=1/3640, alpha=0)
Configuration settings of priors:
> config$priors
$t.end
[
1] "rnorm(n=1,mean=50,sd=5)"
$N
[
1] "rnorm(n=1,mean=5000,sd=100)"
$beta
[
1] "rgamma(n=1,shape=50,rate=1)"
$gamma
[
1] "rgamma(n=1,shape=20,rate=1)"
$mu
[
1] "rgamma(n=1,shape=20,rate=1)"
$alpha
[
1] "rnorm(n=1,mean=1,sd=0.1)"
Kernel distances for varying N, using seq(1000, 10000, 1000) # (from, to, step)
Kernel distances for varying beta using seq(0.01, 0.36, 0.025)
Kernel distances for varying gamma using seq(0.001, 0.08, 0.0025)
Kernel distances for varying mu using seq(0.0001, 0.01, 0.0005)
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub, or mute the thread.
|
Okay, now that I understand these are kernel distances I think we need to redo the analyses for I also expect some of these parameters to be confounded, like |
I've set gamma to 1/520. and mu to 1/3640. Should I switch the values? |
Those target values are fine, but I would refine the range for simulations to evaluate against the target.
… On Jun 28, 2017, at 5:00 PM, Tammy Ng ***@***.***> wrote:
I've set gamma to 1/520. and mu to 1/3640. Should I switch the values?
norm.mode is set to 'NONE'
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub <#42 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABDtUDzSCvsGtuCPaaVu9Yk03J4rOvTNks5sIr73gaJpZM4MDQzY>.
|
Rosemary's paper focuses on a very different kind of model (network based epidemic).
However there are certainly similarities - it is a susceptible-infected process on a network. Actually her model found some signal for N, only I (number of infected) and N are confounded.
At this point I think we should switch out emphasis to bringing in other measures of tree shape, to see if we can improve the performance of ABC-SMC. Let's plan it out at our meeting Monday.
… On Jun 29, 2017, at 10:48 AM, Tammy Ng ***@***.***> wrote:
Kernel distances for varying N, for nsim=100 and seq(1000, 10000, 1000)
In Rosemary's paper, it says that parameter N and m had little to no identifiability with ABC. Is this for a different kind of scenario?
Kernel distances for beta for nsim=100 and seq(0.01, 0.28, 0.025)
Kernel distances for gamma for nsim=100 and seq(0.0005, 0.03, 0.0015)
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub, or mute the thread.
|
By the way, plots are much improved, thanks for refining the range on beta and gamma. |
Yep, I'll look for a theta that reproduces the negative probability error. |
Try using the following approach to examine the distribution of kernel score over two model parameters: plot(x, y, cex=sqrt(z)) # where x and y are two parameters and z is kernel score/distance You can also use the |
Well, the trick is that you need to vary both beta and N - I suggest doing a "grid search" where you examine all pairwise combinations of values {beta} x {N}.
… On Jul 5, 2017, at 9:27 AM, Tammy Ng ***@***.***> wrote:
Am I doing this correctly?
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub, or mute the thread.
|
I have started a grid search script in |
There are a few things to try here:
|
There's a bug in births <- rbind(c('parms$beta * S * I / (S+I) - parms$gamma * I')) should be births <- rbind(c('parms$beta * S * I / (S+I)')) The only place that |
Correcting the
|
The parameter values for some actually didn't make sense. In a previous comment, I've set the target tree with theta <- c(t.end=50, N=5000, beta=0.1, gamma=1/520, mu=1/3640, alpha=0) Configuration settings of priors > config$priors
$beta
[1] "rgamma(n=1,shape=50,rate=1)"
$gamma
[1] "rgamma(n=1,shape=20,rate=1)"
$mu
[1] "rgamma(n=1,shape=20,rate=1)" And now that I've changed them to: > config$priors
$beta
[1] "rgamma(n=1,shape=1,rate=5)"
$gamma
[1] "rgamma(n=1,shape=1,rate=10)"
$mu
[1] "rgamma(n=1,shape=0.2,rate=50)" I can now |
Ran with parameters in commit 4d37ca4 executed with 500 particles finished in approximately 12.39 hours Trajectories of the mean estimates of parameters in In the plots it looks like there's still room to converge for parameters beta, gamma, and mu. But if I lower my accept rate, the runtime would take a while longer. |
Check whether collecting samples with invalid values of |
@gtng92 observed that fixing to actual |
From Issue #120, I decided to implement the idea of incorporating a dummy tree upon an error encountered in True Params: theta <- c(t.end=200, N=10000, beta=0.1, gamma=0.2, mu=0.01, alpha=5) Grid-search Params: t.end <- seq(50, 400, 50)
N <- seq(7000, 14000, 1000)
beta <- seq(0.01, 0.36, 0.05)
gamma <- seq(0.01, 0.36, 0.05)
mu <- seq(0.001, 0.071, 0.01)
alpha <- 5 Average of top 10 avg.theta <- c(t.end=215, N=10900, beta=0.12, gamma=0.235, mu=0.03451, alpha=5) Average of top 20 avg.theta <- c(t.end=217.5, N=10277.78, beta=0.115, gamma=0.23, mu=0.0365, alpha=5) |
Ok, careful about this though -- if there was no real difference among trees with respect to kernel score, one could arbitrarily pick 10 or 20 trees and the mean parameter estimates would still be around these values based on the way you've designed the grid. It will be more convincing when we see the actual distribution of distances over the grid. Thanks! |
I think we're in danger of fishing for results here. Let's see what distances are like using other metrics. Thanks! |
Deprecated by #123 |
Just like kernel-ABC validation
The text was updated successfully, but these errors were encountered: