Mfcttrf

I have a set of healthcare providers serving patients. In a given amount of time, a specific provider can see only a certain amount of patients, depending on the medical procedure and other variables.

A patient, since the booking of the procedure, is seen by the doctor after a certain amount of days, again specific to provider, procedure etc...

I have a data set in which for each provider, for a given amount of time, I know how many patients have been seen and what was the average waiting time.
I want to simulate more patients for each provider and their individual waiting times.

I assumed that the number of patients seen in the given time can be modeled as a $sim Poisson(lambda)$ with lambda depending on provider and procedure characteristics and amount of time.

I modeled the average waiting time as $sim lognormal(mu_global, sigma)$ with parameters as function of the same variables of before plus the log of n.patients.

Finally, I'm modeling the simulated new patients waiting time as $sim Gamma(mu/theta,theta)$ with $mu$ predicted from the model above and $theta$ chosen using domain knowledge since I don't have past information on individual waiting times.

I would like to know if I choose the right distributions given the problem at hand.

It is best to look at the distribution of waiting times for a particular provider. My first thought would be that if the process is anything like a queueing If process that the distribution should be nearly exponential. So I would check to see if the sample mean and standard deviation are approximately equal. If so, I would look to see if an empirical CDF (ECDF) of the data
roughly fits the CDF of $mathsfExp(textrate = lambda),$ where $mu = 1/lambda$ is estimated as $hat mu = bar X,$ the sample mean.

Only if that doesn't seem to work well, would I pursue fitting the data to a gamma distribution. This is also a plausible possibility, partly because the sum of $k$ exponential waiting times (of the same rate) is gamma-distributed with shape parameter $k.$ [If you use a 'distribution ID' program, a 'gamma' distribution will almost always win out over 'exponential' because the family $mathsfExp$ is a sub-family of $mathsfGamma.$]

Exponential data. As an initial example, let's pursue data randomly sampled from $mathsfExp(textrate = 1/10),$ so that the average waiting time is (a perhaps optimistic) 10 days. Suppose we
have waiting times for $n = 500$ patients.

set.seed(714) # for reproducibility
x = rexp(500, 0.1)
mean(x); sd(x)
[1] 9.909112
[1] 10.36662

So the sample mean and SD are about the same. In practice, I would not know
the true rate $lambda$ so I will estimate it as $hat lambda = 1/9.9 = 0.101.$

In the plot below, the boxplot shows many high 'outliers', as is typical of
an exponential sample. The Density function of $mathsfExp(0.101),$ is a
reasonable fit to the histogram of the data. Also (usually more revealing),
the ECDF plot of the sample is well-approximated by the CDF of this distribution. [The ECDF is a 'stairstep' plot that jumps up by $1/500$ at each of the $500$ observed
values.]

par(mfrow=c(1,3))
boxplot(x, col="skyblue2", pch=19, main="Boxplot of DATA")
hist(x, prob=T, br=20, col="skyblue2", main="Histogram with EXP(.101) Density")
 rug(x); curve(dexp(x, .101), add=T, col="red")
plot(ecdf(x), main="ECDF with EXP(.101) CDF")
 curve(pexp(x, .101), add=T, col="red")
par(mfrow=c(1,1))

These favorable results are hardly surprising because data were sampled from an exponential population. If real data performs as well, then you could simulate additional data from a similar population with R code rexp(n, 0.101), where parameter n is the desired number of simulated values.

However, you must realize that you are not gaining additional information about actual patient waiting times by doing that. All the 'information' you have is given by the sample and the assumption that the population is exponential.

Gamma data. If the exponential model does not seem to fit, perhaps the next step is to assume that data are gamma-distributed, to estimate the parameters, and make similar plots to see if you get a better fit to the data.
[Several pages on this site and online discuss estimation of gamma
parameters; one recent page discusses both MMEs and MLEs.]

Just to see what happens if we have gamma data and try to fit an exponential model, I repeat the simulation above, but using x = rgamma(500, 2, .2).

set.seed(714) # for reproducibility
x = rgamma(500, 2, 0.2)
mean(x); var(x)
[1] 10.62662
[1] 59.49749

Pretending that these data are exponential and estimating the rate as
$hat lambda = 0.1062,$ R code similar to that above gives the following
graphs--with noticeably unsatisfactory fits.

Using a gamma model with method-of-moments estimators (MMEs) from the link above, I estimate
the shape parameter as $hat alpha = 1.90$ and the rate parameter as
$hat lambda = 0.178.$ [Maximum likelihood estimators (MLEs) would likely be
more accurate, but MMEs are good enough to give an idea how the graphing procedure works.]