Interval censoring: missing time2

[Lukas Van Riel]

Hi,

I'm exploring whether the INLAjoint package would be suitable to model my dataset as a multistate model without absorbing states. For this I generated a small mock dataset. I succeeded running the model when coding the transitions as observed (event=1), but when I try to run it it as interval censored (event=3), It get the following error:

Error in inla.surv.1(time, event, time2, truncation, cure, .special = .special) :
  'time2' has to be present for interval censored data or in-interval events

I initiate time2 when I create the survival objects, so I'm not sure how to resolve this. I have included a script that should reproduce the error. Do you have a suggestion how to resolve this?

Thank you

Lukas

[Denis Rustand]

Hi Lukas,

Your code seems to run fine for me, are you using the latest GitHub version of the package? You can get it here (see section "Install"): github.com/DenisRustand/INLAjoint

I'll update the CRAN version soon.

Best,

Denis Rustand

he...@r-inla.org

[Lukas Van Riel]

Hi Denis,

Thank you for looking into my code. Updating did resolve the issue, sorry for not reporting back earlier.

I'm reaching out about a new issue related to interval censoring in reversible multi-state models. Please let me know if I should open a new thread instead of continuing here.

For context: I'm trying to model forest plot transitions between classes over ~60 years with approximately decadal observations. To validate my approach before applying it to my real 9-class dataset, I've created a minimal 3-state example with known parameters.

The problem is this: With exact transition times, INLAjoint recovers all parameters (shape, intercepts, covariates). However, when I introduce interval censoring, simulating artificial observation times to mimic my real dataset, Weibull shape parameters are systematically higher than their true values while intercepts are biased towards lower values (more negative). This pattern is consistent across all 6 transitions in my 3-state model. Covariate coefficients show a mixed response, most are retrieved within the credible interval but some lie outside significantly, sometimes even with a sign switch. 

I've tried several things to resolve this: Including/excluding individual ID as frailty, collapsing observation intervals to minimize redundancy (while retaining all information), reducing transition rates to minimize indirect paths, running gap time interpretations instead of clock-forward (although this is tricky in the case of interval censoring). None of these approaches eliminated the shape parameter bias, though they did fix some covariate sign issues related to misclassified indirect transitions.

I'm not sure what else I could try. I keep coming back to the same questions, I was hoping you can help me with some of them:

• Does INLA support using the weibullsurv survival specification for multi-state models with interval censoring in competing risks format, or should I use something different? Ideally the output of my model would be transformed to transition probabilities afterwards so I would prefer to use a parametric approach.

• Could it be related to how INLA(joint) interprets time intervals in multi-row-per-individual data? Does it treat intervals as gap-time durations rather than calendar time windows? 

• Are there known limitations for interval-censored competing risks models that might explain systematic shape parameter bias while covariate effects remain unbiased (maybe something about trading off scale and shape parameters)?

• Do you know of any working examples of interval-censored multi-state models I could use as reference? I’m looking for the correct way to implement this but I do not find a case that exactly matches mine?

In case it could be helpful, I attached a minimal working example that shows the encountered bias. It runs in a few minutes on my personal machine.

Thank you for your time,

Lukas Van Riel

Op donderdag 20 februari 2025 om 15:05:10 UTC-5 schreef Denis Rustand:

[Denis Rustand]

Hi Lukas,

I think the main issue is that you do not account for delayed entry in the at-risk for each transition, while you did it for the "exact times" model. For example if we look at the first individual and the first transition. This individual is observed from time 25 to 34 and right censored, this should be a left-truncated observation. Then this individual has an interval censored event between time 34 and 49, this is fine no need for truncation there as the boundaries are defined by the interval censoring. In the dataset you shared, when you build the survival outcome for each transition (lines 408-414), you should add this truncation for the lines that are not interval censored and do not start at time 0. You also need to properly set the right censoring time for these observations as it uses the "time" argument and not "time2".

In the attached code, I added this truncation and the results look better, although not perfectly aligned. I guess this may be related to the loss of information from the interval censoring and/or the data transformation. For example the first individual is entering the at-risk at time 20 in the "exact times" data while he enters at 25 in the "interval-censored" data, I am not sure why but this could contribute to the differences in the results.

Best,
Denis Rustand

[Lukas Van Riel]

Hi Denis,

Thank you for your swift response and for helping to solve the problem. 

Running your suggestion indeed improves the estimates for the intercepts and hyperparameters. However, it seems to recover the covariate coefficients less reliably. Is there anything else I could try to add to the dataset or model specification to improve this? 

I found indications that the covariate issue might be because of the existence of indirect paths (due to the long observation intervals). However, decreasing transition rates to decrease indirect paths did make the results for intercepts and hyperparameters worse.

Regarding the loss of information due to interval censoring: wouldn't INLA widen its credible intervals to reflect its uncertainty in the case of interval censoring? Or is that not how it works?

Do you think INLA(joint) is suitable to model this kind of multi-state model? Do you have any knowledge of a multi-state model with interval censoring that was modelled with INLA?

Happy to share any code if that would be helpful.

Thank you,

Lukas

Op dinsdag 28 oktober 2025 om 05:45:54 UTC-4 schreef Denis Rustand:

[Denis Rustand]

Hi Lukas,

Overall the main issue with interval censoring in multi-state models is the entry time in the at-risk period. If you consider an illness-death model where an individual is interval censored for the transition from state 1 (Healthy) to state 2 (Sick), the entry in the "at-risk" period for the transition from state 2 (Sick) to state 3 (Death) is also interval-censored, which is not something INLAjoint can handle at the moment and I am not aware of the alternative software that can handle this.

In the data you shared, as I mentioned in my previous message, for the exact data you have different entry times compared to the interval-censored (i.e. 20 vs. 25 for first individual). This should not be different as we need to know the exact entry time for both approaches and they should be the same, here the difference is likely introducing the bias you observe. You can try to use the true values for entry times to see if it fixes the bias.

When you don't have the true value from data simulation, I guess one option is to impute the entry time but probably not optimal, likelihood approaches may exist but are not implemented in INLA at the moment. Maybe sampling in the intervals and Bayesian averaging of the resulting models could work too... I am not very familiar with this.

Best,
Denis Rustand

[Lukas Van Riel]

Hi Denis,

Thank you very much for your suggestions. 

I implemented some of them, and it seems it is possible (in my case) to estimate the exact transition time well enough so that it approaches the results of using the exact transition times. However, I still notice a bias (mostly in the shape parameters and intercepts) that seems to worsen with increasing sample size. Creating the data set with perfectly symmetrical transition hazards between all states did not solve the issue, which I believe means the problem does not lie with indirect paths. It could point towards me running the model incorrectly. 

I am creating the survival objects as follows now (same as before, except the imputed times are used as truncation for the interval censored rows): 

Surv_obs.list <- list() 
for(i in seq_along(data_obs_list)) {     
# Initialize truncation   
data_obs_list[[i]]$truncation <- data_obs_list[[i]]$imputed_entry     

# Override for right-censored: use Tstart instead   
idx_censored <- data_obs_list[[i]]$status == 0   
data_obs_list[[i]]$truncation[idx_censored] <- data_obs_list[[i]]$Tstart[idx_censored]     
# Adjust time vector   
time_vec <- data_obs_list[[i]]$Tstart   
time_vec[idx_censored] <- data_obs_list[[i]]$Tstop[idx_censored]     

Surv_obs.list[[i]] <- inla.surv(     
time = time_vec,    
 time2 = data_obs_list[[i]]$Tstop,     
event = data_obs_list[[i]]$status,     
truncation = data_obs_list[[i]]$truncation 
) } 

With the imputed_entry the calculated or exact transition time. 

Is this how truncation for interval censored transitions should be implemented? I also tried to include this truncation for the right censored rows, but this seemed to worsen the bias.  
Another issue could be that I do not use a frailty term on the ID currently, so perhaps the bias arises from some double counting. Could this be the problem? Is there something besides the id="id" argument in the joint function that could solve this? 

Thank you very much for taking the time to help me out,

Lukas

Op woensdag 29 oktober 2025 om 04:46:23 UTC-4 schreef Denis Rustand:

[Denis Rustand]

Hi Lukas,

After carefully checking all the code and Weibull implementation, I do not see any issue on the estimation part. Truncation is used to define the entry in the at-risk period, if one individual is in state 1, he's not at risk for the transition 2->3 so the time he spends in state 1 should not be included in the model for 2->3. If he transitions to state 2 after 3 units of time, left-truncation should be 3 for the 2->3 model, regardless if he has an interval-censored event or if he is right censored later on. So in case of right censored observation, the data is: truncation = entry time, time = censoring time, status = 1. For interval censored event, the data is: truncation = entry time, time = lower bound for interval censored event, time2 = upper bound, event = 3.

There are two implementations of Weibull in INLA, the default for INLAjoint is the standard parametrization used elsewhere but it is possible to switch to the alternative one which is more stable numerically. Check inla.doc("weibullsurv") in R for details on these two parametrizations. I ran a couple simulations for a simple illness-death model where death is interval censored and noticed more stable results with variant 1 instead of the default variant 0:

Starting with variant 0 and increasing the (fixed and constant) size of the intervals for interval-censored events:
Variant 0, follow-up range: 50, sample size: 3000
Interval censoring of size: 1
Shape:
     True       Est         SD  Lower_95 Upper_95
[2,]  0.8 0.8835395 0.05715800 0.7751403 1.000080

Not looking so good when the interval range is increased:
Variant 0, follow-up range: 50, sample size: 3000
Interval censoring of size: 3
Shape:
     True       Est         SD  Lower_95 Upper_95
[2,]  0.8 0.9362208 0.06584126 0.8067366 1.065382

Increasing sample size leads to better results here:
Variant 0, follow-up range: 50, sample size 10000
Interval censoring of size: 3
Shape:
     True       Est         SD  Lower_95 Upper_95
[2,]  0.8 0.8386719 0.03531101 0.7706673 0.909669

Switch to variant 1:
Variant 1, follow-up range: 50, sample size 3000
Interval censoring of size: 1
Shape:
     True       Est         SD  Lower_95 Upper_95
[2,]  0.8 0.8314387 0.02541669 0.7823051 0.882369

Looks better than variant 0:
Variant 1, follow-up range: 50, sample size 3000
Interval censoring of size: 3
Shape:
     True       Est         SD  Lower_95  Upper_95
[2,]  0.8 0.8426661 0.02912759 0.7864788 0.9011487

Still good with much larger interval compared to variant 0 (verified with multiple seeds):

Variant 1, follow-up range: 50, sample size 3000
Interval censoring of size: 5
Shape:
     True       Est         SD  Lower_95  Upper_95
[2,]  0.8 0.8522674 0.03337600 0.7881691 0.9195598

Here the results start looking bad:
Variant 1, follow-up range: 50, sample size 3000
Interval censoring of size: 10
Shape:
     True      Est         SD Lower_95 Upper_95
[2,]  0.8 1.111211 0.04959942 1.015864 1.211107

Increasing sample size does not help, we reached the limit of my simulation design (need to adapt the design, e.g. longer follow-up):
Variant 1, follow-up range: 50, sample size 10000
Interval censoring of size: 10
Shape:
     True      Est         SD Lower_95 Upper_95
[2,]  0.8 1.098923 0.02675153 1.046942 1.152265

You are likely in this situation where the simulated data does not inform properly all the parameters, in my case it's because the length of the interval I'm trying to simulate is too large compared to the fixed length of the follow-up but it could come from something else in the design.

To switch from variant 0 to variant 1, add "variant=1" in control options. You can also change the priors, you need to run the model with argument "run=F" to build the model without calling INLA, then we can switch priors before running INLA with joint.run():
fit <- joint(
  formSurv = list(s12 ~ 1, s13 ~ 1, s23 ~ 1),
  basRisk = rep("weibullsurv", 3),
  dataSurv = list(d12, d13, d23),
  id = "id",
  control = list(config=TRUE, variant=1), run=F # added variant=1 and run=F

)

# Define priors manually for each parameter of each transition:

fit$.args$control.family[[1]]$hyper <- list(alpha=list(param=5)) # shape prior for 1->2
fit$.args$control.family[[2]]$hyper <- list(alpha=list(param=5)) # shape prior for 1->3
fit$.args$control.family[[3]]$hyper <- list(alpha=list(param=5)) # shape prior for 2->3
fit$.args$control.fixed$mean <- list(Intercept_S1=0, # mean of log scale for 1->2
                                     Intercept_S2=0, # mean of log scale for 1->3
                                     Intercept_S3=0) # mean of log scale for 2->3
fit$.args$control.fixed$prec <- list(Intercept_S1=0.01, # prec of log scale for 1->2
                                     Intercept_S2=0.01, # prec of log scale for 1->3
                                     Intercept_S3=0.01) # prec of log scale for 2->3
runM <- joint.run(fit) # run INLA with modified priors
summary(runM)

You can try with variant 1 and to adapt the priors to see how they impact your results. You can learn more about the Weibull prior distributions in INLA here:
https://www.sciencedirect.com/science/article/abs/pii/S0167715221000602

I don't know about the frailty, you should not include it if you don't simulate one. The bias you observe is likely the result of the simulation design.

Best,
Denis Rustand

[Lukas Van Riel]

Hi Denis,

Thank you for your response and for taking the time to look into it with your simulation analysis. I've implemented your recommendations and ran systematic tests, but unfortunately the bias problem persists with variant 1 and various priors on shapes and intercepts. 

I tested the following on my artificial multi-state model: 

• Sample size analysis: I fitted the model with 1000, 5000 and 10 000 individuals. The bias increases with sample size. 
• Interval width: I fitted models with 60 years of follow-up, for both 12 year and 4 year average observation intervals. The 4 year intervals performed worse for all sample sizes. 
• Prior sensitivity analysis: I varied the priors for the intercepts and shape parameters. This did not significantly impact the results.

The pattern of bias increasing with N is highly unusual and suggests something more fundamental. However, I can not find any model misspecification after implementing all your suggestions. Your simulation showed stable bias at different N, but mine shows systematic worsening. Worse bias in the case of 4 year intervals could also point to the same sample size effect since the follow-up was the same in both cases.

Do you think the interval length could have anything to do with an apparent bias-sample size correlation? Or could there be an interaction between multi-state structure, interval censoring, and truncation that becomes more problematic with larger N?

I have added the code I ran in case that would be helpful. I recognize this is taking considerable time and really appreciate your willingness to help me work through it.

Thank you,

Lukas

Op zaterdag 15 november 2025 om 03:25:57 UTC-5 schreef Denis Rustand:

[Denis Rustand]

Hi Lukas,

As shown in my previous message: yes, the interval size will create bias if assumed too large compared to the rest of the simulation design. Another possible source for the bias is your reversible transitions, when combined with wide intervals, we might miss some unobserved back-and-forth events (I think the msm package addresses some issues related to this).

Anyway, I still think the bias comes from the simulation design rather than INLA/INLAjoint. I made a simulated example for an illness-death model where both the illness and death events are interval censored (but exact entry time in the at-risk for transition to death from illness is known). The results are looking fine (code attached to replicate):

 Parameter True   Estimate         SD
    alpha1  1.5  1.4950706 0.02678937
    alpha2  0.8  0.8326571 0.02821338
    alpha3  1.2  1.1918505 0.03738116
     beta1 -2.5 -2.4946629 0.01504283
     beta2 -3.5 -3.4068396 0.05285284
     beta3 -1.8 -1.7915562 0.05118753

You should try to simplify your code until you identify the origin of the bias (or start from a simple code that works and add layers until the bias appears).

Best,
Denis Rustand