Maximum likelihood estimation with censored data that include repeated measures I have received the review for my latest journal article in which I present an experimental study with human subjects. In this study, we have loaded our subjects with impact loads created by a pendulum. The study's objective was to calculate biomechanical limits from the experimental impact data that can be used for safe physical interactions with machines like collaborative robots.
The impacts' intensity was gradually increased until the subjects' felt a pain that has exceeded their personal pain-tolerance thresholds (PTT). During the tests, we have increased the impact intensity (i.e., maximum impact force) by increasing the pendulum's starting height. The height could only be adjusted in discrete steps. This condition is the reason why our impact forces measured in the tests are interval-censored. The true but unknown PTT lies somewhere between the force that was measured in the final and penultimate impact test. Some test results are, however, left-censored (when the lowest adjustable impact intensity causes pain beyond the PTT) or right-censored (when the highest adjustable impact intensity causes no pain beyond the PTT).
Whenever a subject's PPT was exceeded, we have repeated the tests with an higher impact mass. Altogether, we have increased the impact mass two times so that we received three results from each subject. Our data stem therefore from repeated measures.
The input variables of each test were the impact velocity (depends on the height; potential energy to kinetic energy) and impact mass. The output variable was the maximum impact force as measured by a force sensor that was mounted on the pendulum's front.
I have already created a linear-mixed model to find out that the impact mass is not a significant covariate. I did not analyse the effect of the impact velocity since the relation between this variable and the maximum impact force must be linear according to the following physical model
$F = \sqrt{ m c } v$
where $F$ is the max. impact force, $m$ the impact mass, $c$ the stiffness of the soft tissue under load, and $v$ the impact velocity. This model was confirmed by various articles about impact experiments.
After this analysis, I have estimated parameters for a log-logit CDF using MLE (with R using fitdistcen package). Based on the inverted CDF, I was then able to calculate the desired limits for a specific quantile.
One of the reviewers insisted that my data include repeated measures (which is true) and that I need to take care of it when using MLE. Unfortunately, I could not find a procedure that factors in censored data and data  from repeated measures at the same time. After days of searching for a solution, I feel completely lost. Has anyone an idea how I can treat the repeated measures?
EDIT 1 Added more background information and context to my question as advised by EdM.
EDIT 2 Added further information about the variables in my model as advised by EdM.

The main problem here is that the observations within each subject are not independent. Thus the standard errors of your estimates of the parameters of your model will tend to be too narrow.
A well respected way to "take care of" repeated measures in models fit by maximum likelihood is to keep the same point estimates of coefficients but take the clustering of observations within subjects into account. One way this is done is by calculating a robust, "sandwich" version of the coefficient covariance matrix.
The survreg() function of the R survival package provides a cluster argument through which you can specify which observations come from which individuals to get a robust coefficient covariance matrix. That function can accept Surv() values as outcomes that cover all the types of censoring that you have, and provides your "loglogistic" distribution family as an option.
An alternative cluster implementation that might work immediately would be to use the tools of the sandwich package on the object returned by the fitdistcen() function and a specification of subject identifiers in a cluster argument. I don't know if the sandwich tools will work on that type of model object, however.
Another approach is cluster bootstrapping. Instead of repeated sampling with replacement from the individual observations as in standard bootstrapping, you resample with replacement from the subjects, keeping all of each selected subject's individual observations in as many copies as the individual is represented in the bootstrap sample. That mimics the process of repeated sampling of subjects from the underlying population. You model each cluster-bootstrap sample as you did the full data set, keeping the corresponding coefficient covariance matrix returned by the vcov() function from the fitdistcen object. Do a few hundred such samples, and report the average of the coefficient covariance matrices. That represents the estimated error due to differences among the individuals in the population rather than the covariances weighted toward within-individual distributions returned in your original model.
Fitting a fully parametric random-effects model to this type of censored data is beyond my intellectual PTT. If that's what you want, I suspect that you will be best served by Bayesian modeling. The brms package allows for censored outcomes via its addition-terms. Section 25.4 of this document illustrates how to handle censored data in brms.
