Modelling patient waiting time from an aggregated dataset with brms

Angelo_D_Ambrosio · July 6, 2019, 10:44am

Hello!
I have a dateset organized this way:
$ Service: "chrc: 13 cases.“Neurology visit”: 63, “Ginecology visit”: 64, "Gastr…
$ Provider: "chrc: 9 cases.“Provider 1”: 101, “Provider 2”: 104, “Provider 3”: 85, "Provid…
$ Class: “chrc: 4 cases.“A”: 203, “B”: 207, “C”: 214, “D”: 167”
$ N.patients: “nmrc: mean: 342 ± 1 041, median: 47 (IQR: 11 - 214)”
$ W.time: “nmrc: mean: 27.1 ± 29.7, median: 18.4 (IQR: 8.9 - 33)”
$ Year: “fctr: 2 cases.“1”: 413, “2”: 378”
$ Months: “nmrc: 2 cases.“4”: 378, “12”: 413”
Plus variable Context, defined as paste(Service, Provider, Class)

The data show for each context (combinations of Service, Provider and Class) and for two consecutive years (2018, 2019), the number of served patients and the average waiting time. To note, the number of months is different between the years, being 12 months for the first and 4 for the second.

Our aim is to estimate, for determined groups of contexts, which combinations of number of patients and average waiting time for the remaining 8 months of 2019 will predict that the average waiting time for the whole 2019 for the same contexts is below certain thresholds. As an alternative, which is the maximum waiting time per patient allowed to keep the year average below the threshold.
In order to do this I need to simulate the patients for the remaining months of 2019 and their waiting time.

Here’s how I would do it, but I’m totally open to suggestions:

Jointly model the number of patients and the average waiting time for each context via the formula:
bf(N.patients ~ Service + Provider + Class + Year + Months + (1 | Context), family = poisson()) + bf(W.time ~ Service + Provider + Class + Year + Months + weight(N.patients) + (1 | Context), sigma ~ Service + Provider + Class + Year + Months + weight(N.patients) + (1 | Context), family = lognormal()); is it well specified? I think I need to put some strong priors, because when I tried with just the Poisson part of the model with Cauchy(0, 25) on the coefficients, the convergence was really bad; any suggestions?
Predict the n. of patients and average waiting time for each context in the remaining 8 months of 2019. The outcome of this prediction is the number of patients, isn’t it? Or a distribution of the lambda parameter of the Poisson? I expect an output in which, for every context there is a distribution of patients counts.
For each simulation I assign a waiting time to the single patients per context using a lognormal distribution with the predicted W.Time and sigma as parameters. What about exponential distribution instead? I guess that step 2 and 3 cannot be done directly together with step 1 in one brms call, isn’t it?
I prepare a grid in which for each combination of average W.time and N.patients (for the remaining parts of the year) there is the probability that the global waiting time (for the whole year) is below the thresholds. Instead, I don’t really know how to setup the alternative version, in which the whole year compliance with the thresholds is linked to keeping all patient waiting times below a context-specific limit; any suggestions?

The aim is to give a tool to health operators that tells them that if they stay below the values shown in the grids, there is an high chance that at the end of the year the waiting time goals will be achieved. Between the two metods for computi the grid presented in (4), the first is easier to implement, but the second is easier for the health operators to follow.

paul.buerkner · July 6, 2019, 10:44am

Please use the Interfaces - brms tag for brms related questions. Otherwise they will likely end up being overlooked. I changed manually this time.

Angelo_D_Ambrosio · July 6, 2019, 7:32pm

Ok thank you! Didn’t know about it!

Angelo_D_Ambrosio · July 10, 2019, 5:18pm

Updates on my work:
I tried some models, but they don’t converge adequately, running very slow (~5-10 h), with low effective samples and many divergent transitions.

Before going to my results let’s repeat the data-generating hypothesis: for each context, as defined by crossing 3 variables (Service, Provider, Class), I got a number of patients, each waiting a certain amount of days. In the data, I have the number of patients and their average waiting time, with two observation per cluster (two years, one per year). Each cluster was observed for a different number of months (this of course influence the number of patients).
Since the biggest problems are with the Waiting Times, I’ll describe the data and the results using the Class variable, which is the one with the highest impact on this outcome.

This is the distribution in the original data:

Class  .category  descr                                               range        
B      N.patients mean: 127 ± 263, median: 37 (IQR: 12 - 112)         1 - 2020
B      W.times     mean: 16.3 ± 11.4, median: 13.1 (IQR: 9.57 - 20)    1.33 - 89.5
D      N.patients mean: 129 ± 296, median: 38 (IQR: 13.5 - 113)       1 - 3220   
D      W.times     mean: 29.6 ± 20, median: 25.4 (IQR: 16.4 - 36.6)    1.14 - 127 
P      N.patients mean: 1 001 ± 1 807, median: 346 (IQR: 127 - 1 066) 2 - 13900  
P      W.times     mean: 49 ± 39.8, median: 34.1 (IQR: 20.3 - 66.9)    0.25 - 202 
U      N.patients mean: 25.6 ± 51.5, median: 8 (IQR: 3 - 20.5)        1 - 370  
U      W.times     mean: 9.78 ± 21.1, median: 4.06 (IQR: 2.18 - 9.73)  0.067 - 235

Here’s the first model, with only a fixed effect for Months while the combination of the other factors is treated as a random intercept:

 Family: MV(poisson, lognormal) 
  Links: mu = log
         mu = identity; sigma = log 
Formula: N.patients ~ Months + (Months | ID | Context) 
         W.time ~ Months + (Months | ID | Context) 
         sigma ~ Months + (Months | ID | Context)
   Data: TdA (Number of observations: 778) 
Samples: 16 chains, each with iter = 10000; warmup = 5000; thin = 2;
         total post-warmup samples = 40000
 Priors: c(set_prior('cauchy(0, 10)', class = 'Intercept'), set_prior("cauchy(0, 2.5)", class = 'b') )
Warning message: There were 1386 divergent transitions after warmup.

Here are the eff. sample ratios by parameter type

Par         Par.class   descr                                                         range
W.times     b           mean: 0.127 ± 0.16, median: 0.0665 (IQR: 0.0215 - 0.172)      0.017; 0.36 
W.times     cor         mean: 0.114 ± 0.058, median: 0.137 (IQR: 0.0684 - 0.155)      0.036; 0.17 
W.times     r           mean: 0.205 ± 0.164, median: 0.149 (IQR: 0.0728 - 0.303)      0.015; 0.71 
W.times     sd          mean: 0.0317 ± 0.0123, median: 0.0294 (IQR: 0.0261 - 0.0349)  0.019; 0.049
lp          lp          mean: 0.0138 ± NA, median: 0.0138 (IQR: 0.0138 - 0.0138)      0.014; 0.014
N.patients  b           mean: 0.159 ± 0.175, median: 0.159 (IQR: 0.097 - 0.221)       0.035; 0.28 
N.patients  cor         mean: 0.132 ± 0.0751, median: 0.111 (IQR: 0.0716 - 0.159)     0.046; 0.29 
N.patients  r           mean: 0.279 ± 0.192, median: 0.228 (IQR: 0.112 - 0.429)       0.037; 0.75 
N.patients  sd          mean: 0.182 ± 0.162, median: 0.182 (IQR: 0.125 - 0.239)       0.067; 0.3

And Rhat

Par         Par.class   descr                                                   range
W.times     b           mean: 1.01 ± 0.0145, median: 1.01 (IQR: 1 - 1.02)       1; 1.03 
W.times     cor         mean: 1.01 ± 0.00455, median: 1 (IQR: 1 - 1.01)         1; 1.01 
W.times     r           mean: 1 ± 0.00436, median: 1 (IQR: 1 - 1.01)            1; 1.03 
W.times     sd          mean: 1.02 ± 0.00687, median: 1.02 (IQR: 1.01 - 1.02)   1.01; 1.03
lp          lp          mean: 1.04 ± NA, median: 1.04 (IQR: 1.04 - 1.04)        1.04; 1.04
N.patients  b           mean: 1 ± 0.00604, median: 1 (IQR: 1 - 1.01)            1; 1.01 
N.patients  cor         mean: 1 ± 0.00315, median: 1 (IQR: 1 - 1.01)            1; 1.01 
N.patients  r           mean: 1 ± 0.00196, median: 1 (IQR: 1 - 1)               1; 1.01 
N.patients  sd          mean: 1 ± 0.00306, median: 1 (IQR: 1 - 1)               1; 1

Finally, here’s the predictive distribution by class and variable created starting using add_predicted_draws. I would also like get access to the predicted sigma estimates but I don’t know how. Here I considered as independent all the draws for all original data single observations (so these are not aggregated).

Class   .category       descr                                                                                 range
B       N.patients      mean: 127 ± 262, median: 38 (IQR: 12 - 109)                                           0; 2300
B       W.times         mean: 54 726 ± 74 961 320, median: 13.7 (IQR: 8.93 - 20.5)                            9.5 / 10^30; 1.55 x 10^11
D       N.patients      mean: 129 ± 296, median: 37 (IQR: 13 - 113)                                           0; 3570
D       W.times         mean: 12 214 220 076 910 ± 35 146 404 053 632 668, median: 24.4 (IQR: 14.7 - 37.4)    1.4 / 10^15; 1.01 x 10^20
P       N.patients      mean: 1 000 ± 1 803, median: 352 (IQR: 131 - 1 052)                                   0; 14700
P       W.times         mean: 5 785 573 ± 16 909 344 272, median: 34.9 (IQR: 20.2 - 66.1)                     2.5 / 10^12; 4.95 x 10^13
U       N.patients      mean: 47.2 ± 294, median: 9 (IQR: 4 - 26)                                             0; 31400
U       W.times         mean: *some huge number, practically Inf*, median: 5.69 (IQR: 2.68 - 12.3)            4.8 / 10^65; 1.43 x 10^99

It’s possible to notice some problems with W.times variance estimates, both for the low n.eff ratio that from the totally out-of-range predicted values.

I tried then two more models, one adding prior(cauchy(0,2), class = "sd", group = 'Context', resp = 'N.patients'), prior(cauchy(0,2), class = "sd", group = 'Context', resp = 'W.times') in addition to the prior of the previous model, and one with also prior(lkj(2), class = "cor"). Both models perform worse, with a median of 40% lower n.eff ratios and the one without prior(lkj(2), class = "cor") being slightly worse. The model with the lkj prior is better than the first model regarding the lp parameter, as expected (n.eff mean: 0.0206 ± NA, median: 0.0206, IQR: 0.0206 - 0.0206). Regarding the predictive efficiency, the means and medians of the posterior estimates on N.patients were more or less similar to those of the first model, and so were the medians of the estimates on W.times. Instead the means and sd of W.times became even more extreme with this two models. Also the divergent transitions were almost doubled.

My opinion is that having just two observation per cluster, Stan cannot explore well W.times. Especially for the U class, which has high variability due to fewer patients. I tried to add something like weights(N.patients) in the W.times and sigma formula, but it gave me error ($ operator is invalid for atomic vectors). Then I thought that being a multivariate model, the inverse relationship between N.patients and sigma should have been already modeled. I thought about reparametrization with non-centered random effects, but I read that brms already do that.

Finally I tried a model in which I add Provider, Class and Service as fixed effect in addition to Context, hoping that those variables would provide some cross-cluster information (especially at the class level) and shrink such extreme estimates; but the results were even worse, with almost all n.eff ratio below 1%, longer computing time, ten times more divergent transitions (13607) and similar predictive efficiency (extreme W.times draws).

I’m running out of options here. Should I put a prior or a bound on W.times predictions? but I don’t know how.
Should I remove extreme draws manually? It’s the model totally wrong? In that case, can you help me to write it better?

Thanks

Emmanuel_Charpentier · July 10, 2019, 6:36pm

I note :

This is a red flag…

MegBallard · July 11, 2019, 5:36am

My first instinct is to ask do you expect 2018 to be much different from 2019.

You don’t need to predict anything for 2018 because you have all that data. So you want to know what predictors are necessary to adequately represent the wait times observed.

I would start by just trying to make THAT model.

And note that provider number you have 9 values/levels/cases there. But you have 13 services. So unless some of those providers are on more than one service (which is suspect is not the case) you actually have up to 117 providers (but probably not that many). they should be specified by a provider_id (for example: service_provider)

So I personally would start at one of the following:

wait ~ n.patients + Provider_id
wait ~ n.patients + Service
wait ~ n.patients + Class

And you may or may not want to include one or all of the categorical variables as group level effects instead and maybe with varying slopes.

And you could add month (if you have this information, ie the observations of patients and wait times in Jan, Feb, Mar… ) as pop level predictors. This would allow you to look at 2018 and then use the specific remaining months (summer, fall, winter) which probably incur different numbers of patients and therefore waiting time to help predict how many patients you can expect in the end of the year and how that would impact waiting time.

But I could, as I always say, be totally wrong.

Though I can nearly guarantee some of the model problems are due to mis-specification of levels of provider (and Context already encodes that information into the formula so you do not need to include them individually necessarily). Class of visit should also mean the same thing across all groups. if it does you can include it separately as it’s own group level predictor. If it doesn’t you should specify each “class” within its Service, or provider, specifically

Angelo_D_Ambrosio · July 11, 2019, 7:27am

I know, that’s what I’m trying to solve!

Angelo_D_Ambrosio · July 11, 2019, 8:37am

In this case I’m assuming that there is no year trend.
As you say N.patients and W.times do depend from the interrelation of Provider, Service and Class. Each Provider can offer only certain services, so the actual combinations are 110, and not each service is offered in all the classes. In total we have 416 combinations of Provider, Service and Classes, which I joined in a variable called Context (repeated in almost all cases twice, once per year, for a total of 791 observations). I put them together because the influence of these 3 variable is not simply additive (both in theory and as verified in my data -> larger prediction error).

So the models were always Outcome ~ Months + (Months | ID | Context) with Outcome being N.patients, W.times and sigma, in a multivariate fashion (therefore the ID to join them).

My fear is that using just Context as grouping variable would remove the cross-context regularizing effect that could have provided using them separately (eg if a Context in class U has high W.Times variance due to few patients, it could take information from other U classes). But I tried also Outcome ~ Provider + Class + Service + Months + (Months | ID | Context) and the situation got worse (~13000 divergent transitions, n.eff rate < 1%). Tried also no random effect but Provider * Class * Service and it takes forever.
Does defining the grouping part as (Months | ID | Provider / Service / Class) change anything or is it equivalent to just using Context?

N.patients do influences W.times positively, but not so much (Services with bigger numbers have more personal and are scaled to face larger amount of patients), while instead influences wildly sigma (fewer patients means larger variability in the measured waiting times, due to outliers).
The number of months instead influences a lot N.patients, obviously, but not so much W.times mean and variability, if not through N.patients itself.
Here’s the W.times x N.patients distribution, with also log transformation (W.times = Attesa, N.Patients = N.prest)

I plan to have data on the single months soon, but first I need to specify the model in a way to make it efficient and fast to compute.

If I don’t solve it fast I will have to resolve to a MLE plus random generation solution, which is not ideal (lack of control on priors, regularization, no multivariate and sigma modeling).

SMALL UPDATE:
trying Outcome ~ Months + (Months | Context) with glmer() gives me error because the number of random effects is higher than the number of observations.

MegBallard · July 11, 2019, 10:34pm

so you are trying to model number of patients AND waiting times, rather than just waiting times? (the first post looked like you were trying ultimately to predict waiting times.

I’m only partially familiar with multivariate regression. The responses may need to be on similar scales to make model work more efficiently (I know this is true for numerical predictors).

I was assuming “provider” refers to people (doctors, nurses) and not offices full of people… that may be my misunderstanding. My husband is a physician and here a provider would be a person. And the departments they work for are “services” (he is on Infectious Disease service). So he would be 1 of 8 providers in his service. If that is how your data is set up (and each service has different doctors), than labeling level explicitly is highly recommended. That way when you have nested categorical data you can’t make mistakes when writing the formula. The nesting structure in your formula above seems reasonable but if you have each level labeled so the nesting is “explicit” you can be confident in your choices and can usually write the formula as though the levels are crossed. Assuming you have looked at GLMM FAQ.

So if your data represents 15 providers (doctors/nurses) you should have 15 levels.

The (Month | Context) specification in the glmer could be reduced to (1|Context) to see if that helps. Adding the slope requires more observations per level of context than I think your model contains.

I usually suggest starting as basically as you can with a model and build it up piecewise when you are having trouble with it.

If I have misunderstood all of this I apologize for wasting your time! If I haven’t I’ll check back and can hopefully help!

Angelo_D_Ambrosio · July 12, 2019, 9:18am

Hello, probably is a misunderstanding due to translation from Italian.
Provider is a healthcare facility, like a clinic. Service is the specific healthcare procedure, eg in the case of infectious diseases, could be a visit or a microbiology test. Class is the urgency class (eg. this procedure is urgent, ie < 2 days wait, or can be delayed n days). I called the crossing of Provider, Service and Class as Contexts.
The aim is to model to how many patients (N.patients) will be served a given service in a given healthcare facility (Provider), in a certain amount of time (months), given their urgency class. For example, how many urgent electrocardiogram exams will be performed at the clinic of Springfield in 8 months? Given this number I want to predict the distribution of individual waiting times, starting from the per-context average waiting time I have in my data. Of course, the less patients you have, the more variability in the mean waiting times is to be expected.

Eventually, for each Context, I want to have a distribution of number of patients seen with their respective waiting time, for a number of months of observation.