# 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.

• Please edit the question to say more about the nature of the data and its censoring. It sounds like this might be handled by some of the tools of survival analysis that allow for frailty/random-effect terms, but there aren't enough details in the question to say for sure. Please provide extra information by editing the question, as comments are easy to overlook and can be deleted.
– EdM
Oct 25, 2022 at 15:49
• Thanks for the update. Is the censored outcome the pendulum height with impact mass treated as a covariate, or some impact "intensity" that takes impact mass into account? Are there other covariates in your model, or is your model just intended to describe your data efficiently among a group of participants who aren't otherwise distinguished by covariate values? Again, please provide that information by editing the question.
– EdM
Oct 25, 2022 at 19:58
• Also, in terms of taking care of the repeated measures, do you just want to adjust the coefficient (co)variance estimates to take the clustering into account (GEE-like sandwich, as with a cluster term in the R survival package), or do you want to model the random effect itself? How many subjects are involved?
– EdM
Oct 25, 2022 at 20:23
• My problem is that I do not know how I can model the distribution of the censored data by also taking the repeated measures adequately into account. I assume that the answer is simple, but obviously not for me. I've just gone blank. Oct 25, 2022 at 20:40
• It is probably better for your paper to ask a statistician to collaborate instead of asking on Stack Exchange. And that way they can be an author on the paper they helped to write. Oct 25, 2022 at 21:29

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.

• Thanks, @EdM. I tried the third approach and came to the conclusion that it worked best for my problem. A convincing advantage for me was that it was more or less an upgrade to what I have already done. I found this blog post that helped me lot dive a deeper into the idea and its implementation. Thanks for your help. Much appreciated. Nov 1, 2022 at 20:35