replace latent model with the new version

[parksw3]

I wrote latent_truncation_censoring_adjusted_delay_zero, which allows zero delays for the latent model. It seems to be working (test/zero.R) and just as fast as the original model when we don't have zero delays (test/zero_comparespeed.R).

Think we should replace latent_truncation_censoring_adjusted_delay with the new one if we think both points above are true. But didn't want to do that before getting reviewed.

[parksw3]

I guess this is what the pull request is for...

[seabbs]

This is really nice!

I was a little concerned about settings where we might have overlapping primary and secondary observation windows and generally distressed by for loops so hacked around on this a bit (see main). I don't know if there is much between the implementations so happy to go with whichever. I think speed-wise your for loop-based approach is fractionally faster and maybe easier to follow?

[parksw3]

I like avoiding loops but I'm confused about this (I'm at school and don't have my stan-running laptop machine with me right now so can't test numerically but):

pwindow = pwindow_raw .* (to_vector(vreal2) + swindow .* to_vector(vreal4));

vreal4 is woverlap which is either 0 or 1. So if there's overlap, woverlap = 1, then pwindow gets larger?

Consider a case where ptime_upr = 1, ptime_lwr = 0, stime_upr=0, stime_lwr=0 such that there is an exact overlap. Then,

as.numeric((stime_lwr - ptime_lwr + 1) <= (ptime_upr - ptime_lwr)) = 1

And so pwindow can be greater than 1? Did you mean to do something like:

pwindow = pwindow_raw .* (to_vector(vreal2) .* (1-to_vector(vreal4))) + swindow .* to_vector(vreal4));

[parksw3]

I think there may be some examples where both of our approaches might fail. For example, what happens if 0 < p < 1 and 0 < s < 2. If a sample is at stime = 1.5 then ptime is going to be between 0 and 1.5. So maybe better to start the inference with ptime and then work out stime?

I do think that

pwindow = pwindow_raw .* (to_vector(vreal2) .* (1-to_vector(vreal4))) + swindow .* to_vector(vreal4));

this might be a bit clunky and hard to read though. And don't see much disadvantage in using loops.

[seabbs]

haha if you thought that was clunky you should look at what I just pushed. Some really pointless optimisation happening here.

[seabbs]

Consider a case where ptime_upr = 1, ptime_lwr = 0, stime_upr=0, stime_lwr=0 such that there is an exact overlap. Then,

So for this case I guess we need to check what stime_upr is.

[parksw3]

Consider a case where ptime_upr = 1, ptime_lwr = 0, stime_upr=0, stime_lwr=0 such that there is an exact overlap. Then,

Sorry, typo. I meant stime_upr = 1. I was just trying to make a point that

pwindow = pwindow_raw .* (to_vector(vreal2) + swindow .* to_vector(vreal4));

this doesn't seem right.

haha if you thought that was clunky you should look at what I just pushed.

Just looking at it now and trying to parse what's going on...

[seabbs]

If a sample is at stime = 1.5 then ptime is going to be between 0 and 1.5. So maybe better to start the inference with ptime and then work out time?

So this is what redefining the upper bound of ptime was meant to deal with.

[seabbs]

Hmm done some more redefining but concluding it may be time for me to go and have a cup of tea...

[parksw3]

A few typos in the code. Making suggestions here.

So this is what redefining the upper bound of ptime was meant to deal with.

I'm not sure that this is a good idea. Right now, we have

pwindow_upr := ifelse(
          (stime_lwr - ptime_upr) <= -1, 
          stime_lwr - ptime_lwr,
          ptime_upr - ptime_lwr
        )

Consider a case in which ptime_lwr=0, ptime_upr=1, stime_lwr=0, and stime_upr=2. Then, pwindow_upr = stime_lwr - ptime_lwr = 0. So here, I think you mean to write stime_upr - ptime_lwr instead of stime_lwr - ptime_lwr.

And here's the reason why redefining the upper bound is bad. In the same case, we have pwindow_upr = stime_lwr - ptime_lwr = 2. But there's no reason for us to sample between 1 and 2 for ptime since we know that it's between 0 and 1.

But actually, sampling ptime first doesn't help either because we can have ptime_lwr=0, ptime_upr=2, stime_lwr=1, and stime_upr=2. Same problem...

Also:

pwindow[woverlap] = (to_vector(vreal2[woverlap]) +  swindow[woverlap]) .*
      pwindow_raw[woverlap];

Think this can be simplified:

pwindow[woverlap] = (swindow[woverlap]) .*
      pwindow_raw[woverlap];

right? (sorry to keep suggesting without making changes myself. Can't test any of this on the current laptop so don't want to break the code without knowing.)

[seabbs]

I think there may be some examples where both of our approaches might fail. For example, what happens if 0 < p < 1 and 0 < s < 2. If a sample is at stime = 1.5 then ptime is going to be between 0 and 1.5. So maybe better to start the inference with ptime and then work out time?

Yes agree. I have reformulated along these lines but as I currently have it currently there are some rather large sampling problems that may be on me or may indicate more general issues with the shape of the posterior.

In terms of posterior shape I am not clear there is any information when we are assuming a uniform prior so perhaps we should be making a censoring correction for these values vs trying to model explicitly.

I think as written though it clears up your other points. For some of the simplifications you suggested the reason to keep them (they are gone now but similar ones exist in this implementation) ifsfor the cases when we have very oddly sized censoring windows. That being said if we are going to keep this complexity we perhaps may need to explore it at least a little (though I am not sure that is a good idea).

[seabbs]

In conclusion - I have broken everything 😆

In the simple daily censoring case where there is overlap as I have written it everything cancels down so that we are just looking at pwindow ~ Unif(0,1) and pwindow + (1 - pwindow) Unif(0,1) so I guess it isn't surprising that this is an issue.

      swindow[woverlap] = pwindow_overlap + (
          to_vector(vreal3[woverlap]) - pwindow_overlap
        ) .* swindow_raw[overlap];

[parksw3]

I went through your code and spent some time thinking about this general problem. I feel like the best way to proceed is to make a simplifying assumption so that we don't have to try to cover all possible cases. As I have explained above, it seems like it would be impossible to cover all possible cases, even for very simple examples:

And here's the reason why redefining the upper bound is bad. In the same case, we have pwindow_upr = stime_lwr - ptime_lwr = 2. But there's no reason for us to sample between 1 and 2 for ptime since we know that it's between 0 and 1.

But actually, sampling ptime first doesn't help either because we can have ptime_lwr=0, ptime_upr=2, stime_lwr=1, and stime_upr=2. Same problem... [which requires sampling stime first]

So maybe we shouldn't try to generalize. Instead, in case there's any overlap, what if we assume that ptime_upr should be equal to stime_upr by assumption? For example, we can imagine these synthetic data:

• case 1: ptime_lwr = 0, ptime_upr = 2, stime_lwr = 1, stime_upr = 5. If we know that stime is between 1 and 5, why should we necessarily limit that ptime_upr = 2? It would make more sense to allow for the possibility that ptime can be greater than 2, and we should be setting ptime_upr = 5 by assumption.

• case 2:

  pwindow_upr = stime_lwr - ptime_lwr = 2. But there's no reason for us to sample between 1 and 2 for ptime since we know that it's between 0 and 1.

I'm going to take back what I said before. If we think stime is between 0 and 2, I don't think it doesn't make sense to consider the possibility that ptime can be greater than 1. In the idealized daily censoring, it obviously doesn't make sense. But in real settings, it doesn't seem like a terrible assumption.

I'm not sure if I'm making sense but this assumption simplifies a lot of the things... I think this gets back us to the very beginning, but I think it may be better to keep things simple, fast, and robust rather than to be too general...

[seabbs]

Just went to run this and I am seeing the following when running targets.

✖ error branch standata_truncation_adjusted_a7cd01ec
Error: Truncation bounds are invalid: lb >= ub

✔ skip branch standata_truncation_adjusted_b1647d3c
• start branch standata_truncation_adjusted_ef56afa9
• record workspace standata_truncation_adjusted_8f0e551d
✖ error branch standata_truncation_adjusted_8f0e551d
Error: Some responses are outside of the truncation bounds.

Got to drop off and do other things now but will have a look at this later if time and if not tomorrow.

[parksw3]

Hmmm, I didn't get this error when I ran it... I can also check when I get home today. I'm now at the office so don't have access to stan...

[parksw3]

Hmm, actually I just tried running my files and started seeing the above error. So maybe I missed it before.

But I don't see why we would get this error... It looks like the error must be coming from here:

trunc(lb = 1e-3, ub = censored_obs_time)

when lb >= ub. But we would only get this if censored_obs_time = 0. Currently, we have

censored_obs_time := obs_at - (ptime_lwr + (ptime_upr - ptime_lwr) / 2)

And so we can get censored_obs_time <= 0 if the midpoint between ptime_lwr and ptime_upr is greater than obs_at. This should never happen if we're dealing with a one day width but can happen if we're dealing with more than one day width. But this would also require ptime_upr > obs_at in which case the data point should've been truncated anyway? more digging to do...

[parksw3]

Running

min(sapply(list_observations, function(x) min(x$censored_obs_time)))

gives 0.5 so we should never get ub < lb... what's going on!!!!

[parksw3]

@seabbs ignore all my messages above and just read this one lol. Sorry for spamming.

Ahhhh I found it! Something wrong with Ebola data. Looking at one of the workspaces, we see that delay_daily > censored_obs_time in some cases.

      id onset_date  test_date ptime stime ptime_daily ptime_lwr ptime_upr stime_daily stime_lwr stime_upr delay_daily delay_lwr delay_upr
   <int>     <Date>     <Date> <num> <num>       <num>     <num>     <num>       <num>     <num>     <num>       <num>     <num>     <num>
1:  1280 2014-09-09 2014-09-15   114   120         114       114       115         120       120       121           6         5         7
   obs_at obs_time censored_obs_time censored scenario  obs_type tar_group sample_size        data_type
    <num>    <num>             <num>   <char>   <fctr>    <char>     <int>       <num>           <char>
1:    120        6               5.5 interval 120 days real-time         1         200 ebola_case_study

The problem is here:

ptime = as.numeric(onset_date - min(onset_date)),
stime = as.numeric(test_date - min(onset_date))

where we're taking day 120 as stime for this data even though the true stime is probably between 120 and 121. So when we truncate the data at t=120, this should've been truncated if we knew that true stime was between 120 and 121.

Two solutions:

Solution 1: Add 0.5 to ptime and stime. These ptime and stime are placeholders anyway so it doesn't hurt to add 0.5, I think?

ptime = as.numeric(onset_date - min(onset_date))+0.5,
stime = as.numeric(test_date - min(onset_date))+0.5

Solution 2: Do a better job with filter_obs_by_obs_time. Right now, we have;

DT(stime <= obs_at)

which is clearly not working for discrete cases. One thing we could consider is:

DT(floor(stime) + 1 <= obs_at)

This would clear out all the continuous time cases as it did before. And also help us get around the discrete cases as well.

It might be worth implementing both? I think the first problem definitely needs to be fixed. Less clear on the second one. I think keeping stime <= obs_at technically is correct and makes sense, and I'm not too sure if we want to introduce a work around for special cases... I've commented where we can introduce these changes but wanted to discuss before implementing them.

Also, it seems like it takes around ~2 mins for me to compile each model so running to target file might not be bad at all for me. But let's resolve this first.

[seabbs]

I think (🤞🏻) these issues are now as resolved as they are getting. We can reopen if/when that is not actually the case and then (hopefully) finally resolve.