# Modélisation de données longitudinales : protocoles robustes dans les essais cliniques, interactions bactériophages/bactéries

Jérémy Seurat

► To cite this version:

Jérémy Seurat. Modélisation de données longitudinales : protocoles robustes dans les essais cliniques, interactions bactériophages/bactéries. Méthodologie [stat.ME]. Université de Paris, 2021. Français.
NNT: . tel-04452525

HAL Id: tel-04452525

https://theses.hal.science/tel-04452525

Submitted on 12 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Université de Paris

École doctorale 393 Pierre Louis de santé publique

Épidémiologie et sciences de l'information biomédicale

UMR 1137 - IAME - Infection, Antimicrobiens, Modélisation, Évolution,
INSERM, Université de Paris

Équipe BIPID – Modélisation biostatistique, pharmacométrie, et investigation
clinique en maladies infectieuses

Modélisation de données longitudinales : protocoles
robustes dans les essais cliniques, interactions
bactériophages/bactéries

Par Jérémy Seurat

Thèse de doctorat de biostatistiques et biomathématiques

Dirigée par le Pr. France Mentré, co-encadrée par le Dr. Thu Thuy Nguyen

Présentée et soutenue publiquement le 8 mars 2021

Devant un jury composé de :

Florence Gattacceca,
Nicolas Grégoire,
Vincent Jullien,
Estelle Kuhn,
France Mentré,
Thu Thuy Nguyen,
Luc Pronzato,

MCU, Université d'Aix-Marseille,
MCU, Université de Poitiers,
PU-PH, Université Sorbonne Paris Nord,
DR, Université Paris-Saclay,
PU-PH, Université de Paris,
Dr., Novartis,
DR, Université Côte d'Azur,

examenateur
rapporteur
examenateur
rapporteur
directeur de thèse
co-encadrant
examenateur
Remerciements

Cette thèse fut remplie de moments scientifiquement et humainement enrichissants. Je remercie
les personnes rencontrées sur le chemin menant au diplôme.

A Estelle Kuhn et Nicolas Grégoire qui m'ont fait l'honneur d'être rapporteurs de ma thèse,
pour leur lecture attentive du manuscrit et leurs suggestions.

A Florence Gattacceca, Vincent Jullien et Luc Pronzato qui m'ont fait l'honneur d'examiner ce
travail et de faire partie de mon jury de thèse.

A France, qui a accepté de diriger cette thèse. Ta rigueur, ta volonté de transmettre, de faire
progresser dans les différents aspects scientifiques et du monde de la recherche sont une source
d'inspiration. Je te remercie aussi d'avoir pu être si disponible, y compris ces derniers mois,
d'avoir fait preuve de pragmatisme et d'humanité tout au long de ces années de doctorat.

A Thu Thuy, qui a accepté de m'encadrer en stage de Master 2 puis de co-encadrer cette thèse.
Je te remercie pour ton extrême gentillesse, ton adaptabilité, tes connaissances et je te remercie
surtout pour ta faculté à les expliquer. Cela a été essentiel pour bâtir les fondations de ma thèse.

A Jérémie, qui m'a encadré sur le travail de modélisation des données de phages, merci pour ta disponibilité, le partage de tes connaissances et ton esprit critique.

Aux personnes ayant participé et aidé aux différents travaux menés durant ma thèse, en
particulier Yuxin, Hervé, Ivana, Romain, Pascal Girard, Kosalaram Goteti, Raphaëlle et
Laurent.

Aux gestionnaires de IAME, merci beaucoup Houda de me faciliter la vie depuis mon arrivée
dans l'unité ainsi que pour ton dynamisme au labo et en dehors, merci aussi à Tiphanie, Myriam
et Vincent.

A l'équipe informatique/Catibiomed, sans qui aucun run Monolix ne serait possible, merci
Hervé (encore !), Lionel, Karim, Saadia et François.

Aux étudiants ou chercheurs de IAME, passés ou présents, en particulier à François, Charles,
Julie, Jimmy, Emmanuelle, Cédric, Niccolo, Tram, Florence, Emilie, Morgane (les deux),
Pauline, Marc, Antonio, Vincent, Simon, Nathan, Claire, Marilou, Coralie, Selma, Mathieu,
Manel, Jinju, Guillaume, Aurélien, Marion, Mélanie, Alexandra, Drifa, Arthur, Anne,
Emmanuelle, Dara, Alice, Lucas, Nadia, Sylvie, Anthony, Nadège, Gaspard et Donia. Je vous
remercie chaleureusement pour l'émulation et pour nos moments de franche camaraderie.
A mes amis (les vrais, ceux qui se reconnaitrons) d'enfance, du lycée, de pharma ou d'ailleurs,
merci pour votre soutien et pour tous les moments partagés et à venir.

A ma famille, en particulier à mes chers parents. Merci pour votre éducation, votre amour, votre
générosité et votre soutien sans faille.

A Laure-Hélène, merci d'ensoleiller ma vie.
# **Production scientifique dans le cadre de la thèse**

## **Articles**

### Publications acceptées

**Seurat J**, Nguyen TT, Mentré F. Robust designs accounting for model uncertainty in longitudinal studies with binary outcomes. Stat Methods Med Res. 2020 Mar;29(3):934–52.

**Seurat J**, Girard P, Goteti K, Mentré F. Comparison of various Phase I combination therapy designs in oncology for evaluation of early tumor shrinkage using simulations. CPT Pharmacometrics Syst Pharmacol. 2020 Dec;9(12):686-694.

**Seurat J***, Tang Y*, Mentré F, Nguyen TT. Finding optimal design in nonlinear mixed effect models using multiplicative algorithms. Comput Methods Programs Biomed. In Press.

### Publication soumise

Delattre R*, **Seurat J***, Haddad F, Nguyen TT, Gaborieau B, Kane R, Ricard JD, Debarbieux L, Guedj J. Deployment of phage therapy: from *in vivo* to *in silico* models for the treatment of bacterial pneumonia.

## **Communications orales**

Loingeville F*, **Seurat J***, Nguyen TT, Rivière MK, Mentré F. Using Hamiltonian Monte-Carlo to design longitudinal studies with discrete outcomes accounting for parameter and model uncertainties. *Population Optimum Design of Experiments* (PODE). Bâle, Suisse. 2017.

**Seurat J**, Nguyen TT, Mentré F. Robust designs accounting for model uncertainty in longitudinal studies with binary outcomes. *GDR « Statistiques et Santé »*. Paris, France. 2019.

<footnote>* : co-premier auteur</footnote>

<footnote>* : co-présentateur</footnote>
Présentations affichées

**Seurat J**, Nguyen TT, Mentré F. Robust designs accounting for model uncertainty in longitudinal studies with binary outcomes. *Design of experiments* (DoE). Marseille, France. 2018.

**Seurat J**, Nguyen TT, Mentré F. Robust designs accounting for model uncertainty in longitudinal studies with binary outcomes. *Population Approach Group Europe* (PAGE). Montreux, Suisse. 2018.

**Seurat J**, Girard P, Sharma VD, Goteti K, Mentré F. Comparison by clinical trial simulation of various Phase I combination therapy designs in oncology for evaluation of early tumor shrinkage. *Population Approach Group Europe* (PAGE). Stockholm, Suède. 2019.

**Seurat J**, Tang Y, Nguyen TT, Mentré F, Le Nagard H, on behalf of the PFIM group. New version of PFIM for optimal design in nonlinear mixed effects models using R S4. *Population Approach Group Europe* (PAGE). Stockholm, Suède. 2019.
## Liste des principales abréviations

AGQ : *Adaptive Gaussian quadrature* (Quadrature de Gauss adaptative)

AIC : *Akaike information criteria* (Critère d'information d'Akaike)

AUC : *Area under the curve* (Aire sous la courbe)

CFU : *Colony-forming unit* (Unité formant colonie)

CTS : *Clinical trial simulation* (Simulation d'essais cliniques)

EBE : *Empirical Bayes estimates* (Estimation empirique bayésienne)

ETS : *Early tumor shrinkage* (Décroissance tumorale précoce)

FIM : *Fisher information Matrix* (Matrice d'information de Fisher)

FO : *First order* (Premier ordre)

FW : *Fedorov-Wynn*

HMC : *Hamiltonian Monte-Carlo* (Monte-Carlo Hamiltonien)

MC : *Monte-Carlo*

MNLEM : *Modèle non-linéaire à effets mixtes*

NPDE : *Normalized predictive distribution error* (Résidus simulés normalisés décorrélés)

PD : *Pharmacodynamie*

PFU : *Plaque-forming unit* (Unité formant plaque)

PK : *Pharmacokinetics* (Pharmacocinétique)

QRMC : *Quasi random Monte-Carlo* (Monte-Carlo quasiment aléatoire)

RSE : *Relative standard error* (Erreur standard relative)

SAEM : *Stochastic Approximation Expectation Maximization*

VPC : *Visual predictive check* (Vérification visuelle des prédictions)
# Sommaire

## PREMIÈRE PARTIE : INTRODUCTION .................................................................................... 10

### 1.1. MODÈLES NON-LINÉAIRES À EFFETS MIXTES .................................................................... 10

#### 1.1.1. MODÈLES POUR DONNÉES CONTINUES........................................................................ 13

#### 1.1.2. MODÈLES POUR DONNÉES NON CONTINUES.................................................................... 14

#### 1.1.3. MÉTHODES ET LOGICIELS D'ESTIMATION .................................................................... 16

#### 1.1.4. SÉLECTION ET ÉVALUATION DE MODÈLES .................................................................... 17

#### 1.1.5. MODÉLISATION ET ONCOLOGIE.................................................................................. 19

### 1.2. EVALUATION ET OPTIMISATION DE PROTOCOLES DANS LES MODÈLES MIXTES ........................................ 20

#### 1.2.1. MATRICE D'INFORMATION DE FISHER ........................................................................ 21

##### 1.2.1.1. Méthodes d'évaluations de la FIM........................................................................ 23

##### 1.2.1.2. Critères d'optimalité ................................................................................................. 24

##### 1.2.1.3. Incertitude du modèle et des paramètres ................................................................. 25

##### 1.2.1.4. Logiciels et algorithmes d'optimisation ................................................................. 26

#### 1.2.2. SIMULATIONS D'ESSAIS CLINIQUES ........................................................................ 27

### 1.3. PHAGOTHÉRAPIE ................................................................................................................ 28

#### 1.3.1. PHAGES ET PHAGOTHÉRAPIE ................................................................................ 29

##### 1.3.1.1. Généralités à propos des bactériophages................................................................. 29

##### 1.3.1.2. Phagothérapie........................................................................................................ 31

#### 1.3.2. MODÉLISATION DE L'INTERACTION PHAGE/BACTÉRIE.................................................... 33

##### 1.3.2.1. Modélisation de la cinétique bactérienne ................................................................. 34

##### 1.3.2.2. Modélisation de la dynamique virale .................................................................... 37

##### 1.3.2.3. Modèles mathématiques existants pour la description de l'interaction phage/bactérie.. 39

### 1.4. OBJECTIFS DE LA THÈSE ........................................................................................................ 40

## DEUXIÈME PARTIE : EVALUATION ET OPTIMISATION DE PROTOCOLES DANS LES MODÈLES MIXTES .................................................................................................................... 42

### 2.1. PROTOCOLES ROBUSTES TENANT COMPTE DE L'INCERTITUDE DU MODÈLE POUR RÉPONSE BINAIRE RÉPÉTÉE .... 42

#### 2.1.1. RÉSUMÉ ................................................................................................................................. 42

#### 2.1.2. ARTICLE 1 (PUBLIÉ)........................................................................................................ 44
## 2.2. ALGORITHMES POUR L'OPTIMISATION DE PROTOCOLES.................................................. 68

### 2.2.1. ALGORITHMES MULTIPLICATIFS.................................................................................. 68

2.2.1.1. Résumé .................................................................................................................. 68

2.2.1.2. Article 2 (soumis) .................................................................................................... 69

### 2.2.2. CONTRIBUTION À LA NOUVELLE VERSION DE PFIM.................................................. 121

## 2.3. COMPARAISON DE PROTOCOLES POUR ÉVALUER LA DÉCROISSANCE TUMORALE PRÉCOCE.................. 122

### 2.3.1. RÉSUMÉ .................................................................................................................. 122

### 2.3.2. ARTICLE 3 (ACCEPTÉ) .................................................................................................. 124

# TROISIÈME PARTIE : MODÉLISATION DE L'INTERACTION BACTÉRIOPHAGES/BACTÉRIES POUR LE

## TRAITEMENT DE LA PNEUMONIE .................................................................................. 138

### 3.1. RÉSUMÉ.................................................................................................................. 138

### 3.2. ARTICLE 4 (SOUMIS) .................................................................................................. 140

# QUATRIÈME PARTIE : DISCUSSION .................................................................................. 177

### 4.1. PROTOCOLES EN MODÉLISATION .................................................................................. 177

#### 4.1.1. PROTOCOLES ROBUSTES PRENANT EN COMPTE L'INCERTITUDE DU MODÈLE .................................................. 177

#### 4.1.2. ALGORITHMES MULTIPLICATIFS.................................................................................. 179

#### 4.1.3. COMPARAISON DE PROTOCOLES POUR ÉVALUER LA DÉCROISSANCE TUMORALE PRÉCOCE .................................... 183

### 4.2. MODÉLISATION DE L'INTERACTION BACTÉRIOPHAGES/BACTÉRIES POUR LE TRAITEMENT DE LA PNEUMONIE ... 185

### 4.3. CONCLUSION GÉNÉRALE.................................................................................................. 188

# BIBLIOGRAPHIE.................................................................................................................. 190

# ANNEXE........................................................................................................................................ 203
# Première partie : Introduction

Les travaux de recherche composant ce doctorat ont reposé sur l'utilisation des modèles non-linéaires à effets mixtes (MNLEM), outils adaptés à l'analyse de données longitudinales. Ces travaux sont décomposés en deux parties. La première thématique concerne l'évaluation et l'optimisation de protocoles d'études longitudinales en amont de l'analyse par modélisation. Les travaux méthodologiques en planification d'essais constituent un champ de recherche important. Le choix du protocole est essentiel pour rendre les résultats d'analyse par MNLEM précis et donc permettre des inférences fiables. La seconde thématique est plus appliquée et utilise les MNLEM pour analyser des données précliniques de phagothérapie. La phagothérapie est une alternative à l'antibiothérapie pour le traitement des infections bactériennes, elle consiste en l'administration de bactériophages. Les bactériophages sont des virus infectant uniquement les bactéries. L'émergence de la résistance bactérienne aux antibiotiques est de plus en plus préoccupante, mener des études autour de la compréhension des mécanismes et des solutions pour parer à l'antibiorésistance est donc d'importance majeure. Il demeure nécessaire de bien comprendre et de quantifier les interactions entre bactéries et bactériophages.

Dans la première partie de cette introduction, nous décrivons les MNLEM. Les méthodes actuellement utilisées pour l'évaluation et l'optimisation de protocoles pour des essais analysés par MNLEM sont décrites dans la deuxième partie. Ensuite, la phagothérapie est présentée dans la troisième partie. Dans cette troisième partie, il est aussi présenté l'état de l'art des principaux modèles existants en bactériologie et en virologie qu'il est nécessaire de connaître pour la modélisation de données d'études précliniques sur la phagothérapie. Les objectifs de la thèse sont présentés en fin d'introduction.

## 1.1. Modèles non-linéaires à effets mixtes

Les modèles mathématiques visent à représenter les données issues du monde réel (Witelski and Bowen 2015). Ils peuvent aussi expliquer les mécanismes aboutissant aux observations ou être utilisés dans un objectif de prédiction. Ils peuvent être appliqués dans divers domaines scientifiques (physique, santé, sciences sociales, etc.) et se traduisent sous forme d'équations. Ces équations peuvent être sous forme de solutions analytiques. Elles peuvent aussi être sous la forme d'équations différentielles, notamment utilisées dans les systèmes dynamiques qui permettent de décrire l'évolution temporelle de différentes variables observées. Ce dernier type
de modèle est souvent utilisé dans le domaine de la santé, par exemple pour décrire l'évolution d'une épidémie dans une population ou encore pour décrire l'effet d'un médicament sur la progression d'une maladie à l'échelle individuelle. Souvent, ils visent à représenter les phénomènes le plus mécanistiquement possible, c'est à dire décrivant au mieux les processus biologiques sous-jacents. Le défi dans ce type d'approche est d'estimer les paramètres du modèle, compte tenu des données dont on dispose. En effet, un modèle plus complexe et donc avec davantage de paramètres peut généralement mieux caractériser les phénomènes réels. Néanmoins, si les informations apportées par les données ne sont pas suffisantes, l'estimation des paramètres est imprécise et il est fréquent de ne pas pouvoir identifier ces paramètres. Lors du choix du modèle et de sa paramétrisation, il est donc nécessaire de trouver un compromis entre la représentation fidèle des phénomènes biologiques qui sont complexes et la simplification de ces phénomènes, en fonction des connaissances et données à disposition.

Les modèles mixtes sont particulièrement adaptés à l'analyse de données longitudinales, c'est-
à-dire les données recueillies à différents temps lors du suivi de chaque individu d'une
population (Pinheiro and Bates 2000). Ce type d'approche permet de compenser le manque
d'information pour un individu, en utilisant l'ensemble des données de tous les individus d'une
étude. Les observations de tous les individus sont classiquement décrites par un même modèle
mathématique structural. Les modèles mixtes permettent d'estimer, d'une part, les
caractéristiques de population c'est-à-dire les paramètres moyens (effets fixes), et d'autre part
la variabilité inter- voire intra-individuelle (effets aléatoires), ce qui permet d'obtenir des
paramètres pouvant varier d'un individu à un autre et d'une occasion à une autre,
respectivement. La variabilité peut aussi être « expliquée » par des covariables catégorielles
(e.g. le génotype) ou continues (e.g. l'âge). Enfin la différence entre une observation et la
prédiction du modèle est l'erreur résiduelle, c'est-à-dire la part non expliquée par le modèle.
Cette erreur résiduelle peut comporter, entre autres, l'erreur de mesure et un mauvais choix de
modèle.

Les MNLEM (Davidian and Giltinan 2017) sont fréquemment appliqués en pharmacométrie, l'étude de la pharmacologie quantitative qui comprend la pharmacocinétique (PK) et la pharmacodynamie (PD) (Ette and Williams 2007). La PK qui étudie le devenir d'une substance active contenue dans un médicament après son administration dans l'organisme, c'est-à-dire la concentration (le plus souvent plasmatique) d'un médicament au cours du temps. Les données de PK ont été analysées par MNLEM à partir des années 1970 (Sheiner et al. 1972). Les modèles permettent aussi de décrire la PD, c'est-à-dire l'effet du médicament sur l'organisme. Il est ainsi
possible de mettre en relation la concentration d'un médicament, son effet et le temps par l'utilisation d'un modèle PKPD. La PD peut alors être dépendante de la concentration courante de médicament (réponse directe), un délai peut aussi être pris en compte entre la concentration et l'effet (réponse indirecte). En outre, la PD peut être dépendante de l'exposition au médicament (mesurée par l'aire sous la courbe ou AUC) ou encore du pic de concentration. La pharmacométrie est aujourd'hui employée au cours des différentes étapes du développement d'un médicament, et dans toutes les aires thérapeutiques (Marshall et al. 2019).

Les MNLEM ont de nombreuses applications en pharmacométrie et sont de plus en plus employés pour aider à la prise de décision suite aux études cliniques (Kimko and Pinheiro 2015). De manière générale, le fait de prendre en compte l'ensemble des données longitudinales conduit à une plus grande puissance de tests. Premièrement, ils permettent de mettre en évidence des effets de covariables, par test de Wald ou de rapport de vraisemblance, et ces effets de covariables peuvent avoir un impact clinique. Lorsque des effets de covariables sont mis en évidence, cela peut conduire à des adaptations de posologie en fonction des caractéristiques d'un patient avant même le début de la thérapie, de manière à éviter des doses de médicament inefficaces ou, au contraire, toxiques. Par exemple, les effets de 4 covariables (âge, poids, sexe et créatininémie, marqueur de la fonction rénale) ont été mis en évidence sur la clairance du carboplatine, un anticancéreux (Chatelut et al. 1995). Une fois la thérapie débutée, l'adaptation de la posologie individuelle peut être réalisée. A partir d'un modèle et de quelques mesures chez un patient, il est possible d'estimer ses paramètres individuels, généralement obtenus à partir des EBE (empirical Bayes estimates), et ainsi définir la posologie adaptée pour atteindre la concentration, l'exposition voire l'effet cible. L'intérêt de cette méthode d'ajustement de posologie a par exemple été montré pour la vancomycine, un antibiotique (Berthaud et al. 2019). Une autre utilisation possible des modèles mixtes est la réalisation de tests statistiques basés sur les prédicitions individuelles, c'est-à-dire à partir des EBE. A partir d'un MNLEM, il est aussi possible d'effectuer des simulations étudiant des conditions différentes de l'étude à partir de laquelle le modèle a été construit. Par exemple, dans le cadre de la maladie à virus Ebola, des posologies (doses de charge et d'entretien) pour de futures études ont été recommandées par simulations (Madelain et al. 2017). Il a aussi été montré l'impact du moment de l'initiation du traitement sur le risque de décès (Madelain et al. 2018). De manière plus générale, les MNLEM sont souvent utilisés pour la sélection de dose durant le développement pharmaceutique (Bonate 2000; Miller 2005). Enfin, les MNLEM sont utiles pour évaluer et optimiser les protocoles d'études (section 2.).
### 1.1.1. Modèles pour données continues

Les MNLEM, notamment en pharmacométrie, sont le plus souvent employés pour décrire des données continues, par exemple la concentration d'un médicament, la glycémie, la taille tumorale ou encore la charge bactérienne.

Pour décrire un vecteur de données continues $y_i$ (de longueur $n_i$) d'un individu $i$ ($i = 1, ..., N$), un MNLEM peut s'écrire de la manière suivante (équation 1) :

$$y_i = f(\theta_i, \xi_i) + \varepsilon_i, \quad (1)$$

où la fonction $f$ est le modèle structurel, $\theta_i$ est le vecteur des paramètres individuels, $\xi_i$ est le protocole élémentaire (temps de mesure, schéma d'administration, etc., voir section 1.2.) et $\varepsilon_i$ le vecteur d'erreurs résiduelles individuel ($\varepsilon_i = (\varepsilon_{i1}, ..., \varepsilon_{in_i})^T$).

Dans ce type de modèle, $f$ est donc une fonction mathématique commune à l'ensemble des individus. Les différences entre les individus sont décrites par les paramètres individuels $\theta_i$, qui peuvent être donnés par la fonction $g(\mu, b_i, z_i, \beta)$ avec $\mu$ le vecteur de paramètres à effets fixes, $b_i$ le vecteur d'effets aléatoires, $z_i$ le vecteur des éventuelles covariables et $\beta$ le vecteur d'effets de covariables. Supposons d'abord un modèle simple sans effet de covariable. La fonction $g$ est le plus fréquemment additive ($\theta_i = \mu + b_i$) ou exponentielle ($\theta_i = \mu \times e^{b_i}$). Le vecteur d'effets aléatoires $b_i$ suit une loi normale d'espérance 0 et de matrice de variance-covariance $\Omega$. Lorsque $g$ est additive, la distribution de $\theta_i$ est donc normale, lorsque $g$ est exponentielle, la distribution de $\theta_i$ est log-normale. D'autres distributions sont possibles pour $\theta_i$ (logit-normale, probit-normale, multimodale, etc.) mais ne seront pas détaillées ici. Les effets aléatoires peuvent caractériser la variabilité inter-individuelle (BSV) ou la variabilité inter-occasion ou intra-individuelle (WSV), c'est-à-dire la variation des paramètres entre différentes périodes (Karlsson and Sheiner 1993).

Les effets de covariables peuvent expliquer une partie de la variabilité observée. Il existe deux grands types de covariables : catégorielle (par exemple le sexe, un variant génétique ou l'échelle de performance ECOG (Eastern Cooperative Oncology Group) à l'inclusion) ou continues (par exemple l'âge ou la masse corporelle). Pour une covariable catégorielle, il est commun de considérer une classe de référence et d'estimer l'effet des autres classes sur un paramètre. Les paramètres individuels peuvent ainsi s'écrire $\theta_i = \mu + b_i + \beta z_i$, avec $z_i = 0$ pour la classe de référence et $z_i = 1$ pour l'autre classe, pour une covariable binaire. Il est aussi possible de rassembler des individus de différentes classes afin de gagner en puissance statistique. Un
modèle de covariable peut s'écrire de différentes manières, notamment pour les covariables continues. En plus des effets de covariables à l'inclusion, l'effet de covariables temps-dépendantes peut être étudié, par exemple, la variation de la créatininémie d'un individu depuis le début de l'étude (Wählby et al. 2004).

Le vecteur de résidus est noté $\varepsilon_i$ et est supposé suivre une loi normale de moyenne 0 et de variance $\Sigma_i$, avec $\Sigma_i$ qui est souvent une matrice diagonale de taille $n_i$, définie telle que dans l'équation 2 :

$$ \Sigma_i = \mathrm{diag} \left( \sigma_{\mathrm{inter}} + \sigma_{\mathrm{slope}} \times f(\theta_i, \xi_i) \right)^2, \quad (2) $$

où $\sigma_{inter}$ et $\sigma_{slope}$ sont respectivement les paramètres d'erreur additive et proportionnelle. Des modèles d'erreur plus complexes peuvent aussi être utilisés.

Pour les modèles de données continues, le vecteur de paramètres estimés $\psi$, de longueur $P$, est composé de chaque effet fixe $\mu$, de chaque effet de covariable $\beta$, des éléments de la matrice de variance-covariance $\Omega$ et des paramètres d'erreur $\sigma_{inter}$ et $\sigma_{slope}$.

### 1.1.2. Modèles pour données non continues

De plus en plus, les MNLEM sont utilisés pour analyser des données longitudinales non continues. C'est-à-dire les données binaires, scores cliniques, données de comptage, temps d'événement, etc.

Les données binaires répétées analysées par MNLEM sont au centre de l'un des travaux de ce doctorat. Dans différents domaines, la réponse peut être définie de manière binaire. Les critères peuvent être de type oui/non pour définir, entre autres, la guérison, la présence d'un symptôme ou d'un effet indésirable. Par exemple, lors d'un test diagnostique du virus de l'immunodéficience humaine, la charge virale est détectable si supérieure à 200 copies par mL ou indétectable sinon (Rodger et al. 2019). Le portage bactérien peut être dichotomisé de manière similaire et mesuré à différents temps. Par exemple, dans une étude sur *Streptococcus pneumoniae*, des écouvillons ont été réalisés chez des enfants à différents temps, et la relation entre la probabilité de portage et différentes variables telles que l'âge ont été étudiées (Syrjänen et al. 2001). Les études sur la maladie à SARS-CoV2 2019 rapportent aussi des critères binaires pour le diagnostic, la rémission ou la séroconversion, bien que la dichotomisation de l'information en positif/négatif peut conduire à une perte de puissance statistique (Blain et al. 2020; Fragkou et al. 2020). Lorsque les données binaires sont recueillies de manière
longitudinale, leur analyse nécessite l'utilisation de modèles mixtes linéaires généralisés tels que les modèles logistiques qui peuvent s'écrire, pour une donnée binaire $y_{ij}$ pour l'individu $i$ et l'observation $j$, au temps de mesure $t_{ij}$ (équation 3) :

$$ \operatorname{logit} \left( p(y_{ij} = 1 | b_i) \right) = f(\theta_i, t_{ij}). \qquad (3) $$

Par exemple, pour un lien logit-linéaire avec le temps, la probabilité d'observer la modalité $y_{ij} = 1$ (équation 4) s'écrit :

$$ p(y_{ij} = 1 | b_i) = \frac{\exp(\theta_{1i} + \theta_{2i} \times t_{ij})}{1 + \exp(\theta_{1i} + \theta_{2i} \times t_{ij})}, \qquad (4) $$

où $\theta_{1i}$ et $\theta_{2i}$ sont respectivement les paramètres d'origine et de pente de l'individu $i$.

D'autres liens logit peuvent être utilisés (logarithmique, quadratique, exponentiel, etc.). Comme pour les modèles de données continues, les paramètres individuels $\theta_i$ sont donnés par la fonction $g(\mu, b_i, z_i, \beta)$ avec $\mu$ le vecteur de paramètres à effets fixes, $b_i$ le vecteur d'effets aléatoires, $z_i$ le vecteur des éventuelles covariables et $\beta$ le vecteur d'effets de covariables.

Pour les modèles de données binaires répétées, le vecteur de paramètres estimés $\psi$, de longueur $P$, est composé de chaque effet fixe $\mu$, de chaque effet de covariable $\beta$ et des éléments de la matrice de variance-covariance $\Omega$.

Au cours de cette thèse, des données de survenue d'événement ont aussi été analysées. Différents types d'analyse de ce type de données existent, les plus connues étant l'estimation non paramétrique selon Kaplan et Meier (Kaplan and Meier 1958) et semi-paramétrique selon Cox (Cox 1972). Dans les modèles de Cox (modèles à risques proportionnels), le risque instantané $h_i(t)$ pour l'individu $i$ au temps $t$, qui correspond à la probabilité de survenue d'un événement sur un intervalle de temps sachant que l'événement ne s'est pas encore produit, s'écrit (équation 5) :

$$ h_i(t | z_i) = h_0(t) \exp(\beta \times z_i), \qquad (5) $$

avec $h_0(t)$ le risque de base et $\beta$ l'effet de la covariable $z_i$.

La probabilité de survie de l'individu $i$ $S_i(t)$ jusqu'à $t$ est alors décrite par l'équation 6 :

$$ S_i(t | z_i) = e^{-\int_0^t h_i(s | z_i) ds}. \qquad (6) $$

Il est aussi possible d'étudier la survenue d'un événement par approche paramétrique, avec $h_0$ une fonction de risque de base dépendant de paramètres. Enfin, le lien entre un biomarqueur longitudinal et le risque instantané peut être établi. Lorsque cet événement est le décès, ce type
d'analyse s'appelle la modélisation conjointe (Tsiatis et al. 1995; Wu et al. 2012). Dans ce type d'analyse, les données longitudinales du biomarqueur sont décrites par un modèle, par exemple, un MNLEM, et la fonction de risque instantané peut s'écrire (équation 7) :

$$h_i(t|F_i(\theta_i, t), \theta_s) = h_0(t) \exp(\beta \times f(\theta_i, t)), \quad (7)$$

avec $F_i(\theta_i, t)$ l'historique de la prédiction individuelle jusqu'au temps $t$, $\theta_s$ les paramètres du modèle de survie (pour le risque de base et le lien avec le marqueur longitudinal), $h_0(t)$ la fonction décrivant le risque de base et $f(\theta_i, t)$ la prédiction individuelle de la variable longitudinale au temps $t$. D'autres fonctions de lien entre le risque instantané de décès et la variable longitudinale sont possibles, par exemple un lien non linéaire (Madelain et al. 2018). Pour les modèles conjoints, le vecteur de paramètres estimés $\psi$, de longueur $P$, est composé des paramètres du modèle longitudinal (i.e. pour données continues, voir section 1.1.1.) et des paramètres de survie.

Les modèles conjoints présentent certains avantages. Ils permettent de prédire la probabilité de survie future d'un individu de manière dynamique, à partir du modèle et de la trajectoire du marqueur longitudinal d'intérêt chez cet individu (Desmée et al. 2015; Riglet et al. 2020). Il a par exemple été montré que l'utilisation des données longitudinales dans un modèle en oncologie permet de meilleures prédictions que celle des données à l'inclusion uniquement (Tardivon et al. 2019). En outre, l'estimation des paramètres longitudinaux peut être biaisée lorsqu'un nombre conséquent d'individus sortent rapidement de l'étude de manière non aléatoire (décès ou à échec de la thérapie par exemple). Dans ce cas, pour un modèle qui ne tient pas compte de la survie (donc de la sortie d'étude), les valeurs de paramètres peuvent être influencées par les données des individus survivants, qui ne sont pas nécessairement représentatifs de l'ensemble de la population étudiée. Les modèles conjoints permettent de corriger ce biais (Rizopoulos 2012; Desmée et al. 2015).

### 1.1.3. Méthodes et logiciels d'estimation

La majorité des méthodes d'estimation de paramètres dans les MNLEM sont fondées sur le maximum de vraisemblance. Cette vraisemblance, dans le cas de modèles non linéaires, n'a pas de solution analytique. La vraisemblance $L$ des observations $y_i$ d'un individu $i$ s'écrit (équation 8) :

$$L(y_i|\psi) = \int_{b_i} p(y_i|b_i, \psi)p(b_i|\psi)db_i, \quad (8)$$
avec $p(y_i|b_i, \psi)$ la distribution de probabilité conditionnelle (PDF) des $y_i$ sachant $b_i$ et $p(b_i|\psi)$ la PDF de $b_i$. Il existe différents algorithmes permettant d'estimer les paramètres des MNLEM.
Le premier algorithme a été développé au début des années 1970 par Sheiner. Cet algorithme, nommé *First Order* (FO) repose sur une linéarisation du modèle autour de l'espérance des effets aléatoires (Sheiner et al. 1972). Des extensions de cette méthode ont été proposées ensuite : par exemple, l'algorithme *First Order Conditional Estimation* (FOCE) se base sur la linéarisation du modèle autour des estimations individuelles des effets aléatoires (Lindstrom and Bates 1990). Ces deux méthodes sont intégrées au logiciel NONMEM (https://www iconplc.com/innovation/nonmem/, Bauer 2019a, 2019b). Une autre méthode populaire pour l'estimation des paramètres des MNLEM est l'algorithme *Stochastic Approximation Expectation Maximisation* (SAEM) couplé aux chaînes de Markov de Monte carlo (MCMC) (Kuhn and Lavielle 2005). A chaque itération, une première phase (« Expectation ») consiste à simuler des paramètres individuels par MCMC puis à estimer la vraisemblance. La seconde phase (« Maximisation ») consiste à déterminer les paramètres maximisant la vraisemblance. Les paramètres sont alors mis à jour pour l'itération suivante jusqu'à convergence. Cet algorithme est intégré aux logiciels NONMEM, MONOLIX (http://lixoft.com/products/monolix/) et au package R *saemix* (Comets et al. 2017). L'algorithme SAEM est généralement plus efficace que FOCE dans le cas de données éparses et de modèles complexes (Bauer 2019b).

Dans le cadre de ces travaux de thèse, les paramètres des NLMEM ont été estimés avec l'algorithme SAEM du logiciel MONOLIX. Cet algorithme permet de fournir des estimations pour des modèles décrivant différents types de données : continues, discrètes ou de survie.

### 1.1.4. Sélection et évaluation de modèles

Les méthodes de sélection de modèle sont une thématique majeure en biostatistiques et pharmacométrie. Usuellement en pharmacométrie, la comparaison de modèles se fait selon le critère d'Akaike (AIC) ou de Schwartz (BIC) qui sont fondés sur la vraisemblance (Delattre et al. 2014). Un modèle avec un AIC ou BIC plus faible est plus parcimonieux. La précision d'estimation des paramètres doit aussi être vérifiée. L'évaluation de la précision d'estimation des paramètres se fait en analysant l'erreur standard relative (RSE). Les RSE des paramètres sont obtenues en divisant l'erreur standard (SE) par l'estimation du paramètre, avec les SE des paramètres obtenues à partir de la racine carrée des éléments diagonaux de l'inverse de la matrice d'information de Fisher. Bien que ce critère ne soit pas clairement établi dans la
littérature, on considère qu’un paramètre est précisément estimé si la RSE ne dépasse pas 30 % pour les paramètres à effets fixes et 50 % pour les variances des effets aléatoires.

Au cours de la construction d'un modèle et avant d'utiliser un modèle, il est nécessaire de l'évaluer selon différents critères tels que l'adéquation entre données et prédicitions du modèle et, si possible, la vérification des hypothèses émises. Récemment, des standards ont été proposés concernant l'évaluation de modèles, en particulier au moyen d'outils graphiques. Un tutoriel a été proposé par la société internationale de la pharmacométrie (ISoP) pour l'évaluation de modèles en pharmacométrie pour les données continues (Nguyen et al. 2017). Des recommandations ont aussi été émises par la FDA (Food and Drug Administration), l'agence réglementaire américaine, pour l'évaluation de modèles pharmacocinétiques de population (Food and Drug Administration 2019). Les principaux graphiques d'évaluation sont les observations en fonction des prédicitions (individuelles ou moyennes), qui doivent être également répartis de part et d'autre de la droite d'identité, et les résidus (individuels, de population, pondérés, etc.) en fonction du temps ou des prédicitions qui doivent être répartis de part et d'autre de l'axe y = 0. Les graphiques individuels, c'est-à-dire la représentation des prédicitions et des observations en fonction du temps pour chaque individu, peuvent aussi être utilisés à des fins diagnostiques. Pour une évaluation avancée du modèle, des méthodes utilisant des simulations sont nécessaires. Les vérifications visuelles des prédicitions (VPC, Visual Predictive Check) permettent de comparer les simulations basées sur le modèle aux observations pour différents percentiles (Yano et al. 2001). Les résidus simulés normalisés décorrélés (NPDE, Normalized Prediction Distribution Error) en fonction du temps ou des prédicitions, doivent être répartis également autour de l'axe y = 0 et leur distribution doit être proche d'une loi Normale d'espérance nulle et de variance 1 (Brendel et al. 2006).

Bien qu'il n'y ait pas de standard concernant l'évaluation de modèles pour données non continues répétées à ce jour, quelques outils sont disponibles. Par exemple, pour les modèles de données binaires, des résidus fondés sur une variable latente ont été proposés (Liu and Zhang 2018). Pour les modèles de données de temps d'événement tels que les modèles conjoints, il est possible de procéder à l'évaluation du modèle par VPC construits autour de l'estimateur de la survie de Kaplan-Meier (Friberg et al. 2009; Björnsson and Simonsson 2011) ou encore par NPDE (Cerou et al. 2018).
### 1.1.5. Modélisation et oncologie

Le choix de protocoles des essais cliniques en oncologie est un point critique, en particulier pour les phases précoces, où la réponse, la toxicité et la PK des anticancéreux sont évaluées. Dans ces études, il est de plus en plus recommandé d'utiliser des biomarqueurs plutôt que des critères tels que la survie des patients qui nécessitent des études longues (Seymour et al. 2013). Parmi ces biomarqueurs, la taille tumorale est un marqueur d'efficacité intéressant. Particulièrement, la décroissance précoce de la taille tumorale (ETS, Early Tumor Shrinkage), mesurable après quelques semaines de traitement, est significativement liée aux critères plus cliniques, couramment utilisés en oncologie tels que la survie ou la survie sans progression (Claret et al. 2009; Wang et al. 2009).

La littérature est particulièrement riche concernant les modèles d'évolution de la taille tumorale, la plupart sont décrits dans une revue (Ribba et al. 2014) et disponibles dans les librairies de modèles telles que DDMoRe (http://repository.ddmore.eu/). La taille tumorale est souvent mesurée comme la somme des plus longs diamètres des différentes lésions tumorales. Ces mesures sont néanmoins sujettes à une erreur de mesure potentiellement importante due à la variabilité entre les différents appareils d'imagerie et les opérateurs (Sullivan et al. 2013). Cela renforce l'utilité des MNLEM, qui prennent en compte ces erreurs de mesure au travers de l'estimation de l'erreur résiduelle.

Parmi les modèles les plus anciens, le modèle de Gompertz (équation 9) est utilisé pour décrire l'évolution naturelle de la taille tumorale en l'absence de traitement, c'est-à-dire une croissance exponentielle au début puis un ralentissement de la vitesse de croissance jusqu'à un plateau (Gompertz 1825). Il peut s'écrire sous cette forme :

$$ \frac{dTS}{dt} = KS \times TS \times \ln \left( \frac{TS}{TS_{max}} \right), \qquad (9) $$

avec TS la taille tumorale, t le temps, KS le paramètre de croissance tumorale et TSmax la taille tumorale maximale.

Simeoni a développé un modèle à partir de données chez la souris qui prend en compte les différents états cellulaires, et utilise donc des compartiments de transduction où les cellules tumorales prolifèrent de manière exponentielle et s'engagent vers les phases de mort cellulaires, notamment lorsqu'exposées à des traitements anticancéreux (Simeoni et al. 2004). Les modèles d'inhibition de la croissance tumorale (TGI, Tumor Growth Inhibition) peuvent prendre différentes formes : la croissance peut être linéaire, exponentielle, de type Gompertz ou autre
et la décroissance médiée par la concentration ou l'exposition médicamenteuse, avec ou sans résistance au traitement. Ces modèles peuvent s'appliquer à différents types de cancer et de thérapies. Par exemple, le modèle de Claret (équation 10) s'écrit :

$$ \frac{dTS}{dt} = KS \times TS - KD \times AUC_{X} \times e^{-KR \times t} \times TS, \quad (10) $$

avec TS la taille tumorale, t le temps, KS le paramètre de croissance tumorale, KD la décroissance tumorale liée à l'exposition ($AUC_X$) à l'anticancéreux X et KR la résistance à X. Ce modèle a été développé dans le cadre d'une étude sur le cancer colorectal (Claret et al. 2009) mais a aussi été utilisé pour des données d'étude sur le cancer thyroïdien (Lu et al. 2010).

La modélisation est aussi utile en phase 2 ou en phase 3, lorsque l'efficacité de la thérapie est souvent évaluée sur la survie ou la survie sans progression. En particulier, la méthodologie des modèles conjoints constitue un outil idéal pour l'analyse de ces données d'événement, en liant la pharmacocinétique du ou des médicaments, voire l'évolution d'un biomarqueur (par exemple, la taille tumorale (Tardivon et al. 2019) ou le PSA (Prostate Specific Antigen) dans le cancer de la prostate (Desmée et al. 2015)) au risque de décès.

En oncologie, des combinaisons de médicaments sont fréquemment l'objet d'essais cliniques. Lors de l'utilisation de MNLEM dans ce cadre, il est préférable de quantifier l'effet propre de chaque molécule impliquée ainsi que l'interaction entre les différentes molécules : antagonisme, additivité ou synergie, par exemple par le biais d'un modèle d'interaction pharmacodynamique globale (Wicha et al. 2017). Néanmoins, afin d'identifier tous ces éléments, il est nécessaire de disposer de données informatives, ce qui est possible si le protocole de l'étude a été élaboré judicieusement en amont.

## 1.2. Evaluation et optimisation de protocoles dans les modèles mixtes

Les modèles mixtes permettent d'analyser des données éparses, c'est-à-dire avec peu d'information par individu (nombre réduit de temps de mesures par exemple) par rapport au nombre de paramètres du modèle que l'on souhaite estimer. En amont de la réalisation d'une étude, son protocole doit être défini. En optimisation de protocoles, on s'intéresse aux méthodes permettant de récolter le maximum « d'information » dans le but de répondre aux objectifs d'une étude, en prenant en compte les contraintes imposées par celle-ci. Afin d'évaluer et d'optimiser un protocole avant analyse par MNLEM, les deux méthodes possibles sont
l'utilisation de la matrice d'information de Fisher (FIM) ou la simulation d'essais cliniques CTS (ou non-cliniques). Une première méthode d'évaluation de la FIM pour les modèles mixtes a été proposée à la fin des années 1990 (Mentré et al. 1997). Cette thématique de recherche a ensuite fait l'objet de nombreux développements (davantage détaillés en section 1.2.1.) : des travaux méthodologiques et appliqués; un congrès : le PODE (Population Optimum Design of Experiments, http://www.maths.qmul.ac.uk/~bb/PODE/PODE.html); des logiciels tels que PFIM, développé dans notre équipe de recherche.

Pour un individu *i* (*i* = 1, ..., *N*), le design élémentaire, noté ξᵢ, est défini par le nombre de mesures *n*ᵢ, par les temps de mesures (*t*ᵢ₁, ..., *t*ᵢ*n*ᵢ) et possiblement par d'autres éléments telle que la dose *d*ᵢ d'un médicament. Le design de population, noté Ξ, est défini par le nombre total d'individus *N* et par l'ensemble des designs élémentaires de chaque individu : Ξ = {ξ₁, ..., ξ_N}.
Les designs élémentaires sont souvent regroupés en *Q* groupes d'individus dans un design de population. Au sein d'un groupe *q* (*q* = 1, ..., *Q*), tous les individus *N*_q_ partagent alors le même protocole élémentaire ξ_q. Le protocole de population peut alors être défini tel que Ξ = {[ξ₁, *N*₁], ..., [ξ_Q, *N*_Q]}. Déterminer un protocole d'étude consiste à choisir le nombre d'individus, le nombre de groupes et d'individus par groupe, les temps de mesures, le nombre de temps de mesures ou la posologie (dose de charge, dose d'entretien, nombre de doses, intervalle de temps entre deux doses). Le choix du design influence directement la précision d'estimation des paramètres. En amont de la modélisation, le choix du design est donc crucial car mal choisi, il peut conduire à des études ininterprétables ou non concluantes (Jonsson et al. 1996). De plus, un bon choix de protocole peut aussi permettre de rendre une étude plus éthique et de réduire ses coûts, en ne mesurant la ou les variables d'intérêts qu'aux temps les plus informatifs, voire en réduisant le nombre d'individus au strict nécessaire pour répondre aux objectifs.

1.2.1. Matrice d'information de Fisher

Au début du 20ème siècle, Fisher a montré le lien existant entre l'information obtenue à partir
de données expérimentales et la matrice d'information de Fisher (Fisher 1922). L'approche
basée sur l'évaluation de la FIM (Atkinson et al. 2007; Fedorov and Leonov 2013; Pronzato
and Pázman 2013) est une alternative efficace à la méthode par CTS. Selon l'inégalité de
Cramér-Rao, l'inverse de la FIM est la borne inférieure de la matrice de variance-covariance
de tout estimateur non biaisé des paramètres. En d'autres termes, la racine carrée de l'inverse
des éléments diagonaux de la FIM constitue une prédiction des erreurs standard (SE) d'estimation des paramètres du modèle.

La FIM $\mathcal{M}$ d'un protocole élémentaire $\xi_i$ peut être exprimée de la manière suivante (équation 11) :

$$ \mathcal{M}(\psi, \xi_i) = E_{y_i} \left( \frac{\partial \log(L(y_i|\psi))}{\partial \psi} \frac{\partial \log(L(y_i|\psi))^T}{\partial \psi} \right), \quad (11) $$

avec $\psi$ le vecteur de paramètres à estimer, $y_i$ le vecteur d'observation de l'individu $i$ (dépendant notamment de $\xi_i$) et $L$ la vraisemblance.

La FIM d'un protocole de population $\Xi$ est en fait la somme des $N$ FIM élémentaires (équation 12) :

$$ \mathcal{M}(\psi, \Xi) = \sum_{i=1}^{N} \mathcal{M}(\psi, \xi_i). \quad (12) $$

Dans le cas particulier où tous les individus d'une étude partagent le même protocole élémentaire, la FIM de population est donc égale à $\mathcal{M}(\psi, \xi) \times N$.

La FIM est donc obtenue à partir du modèle et son vecteur de paramètres associés, pour un design de population. Il est ensuite possible de prédire la SE des paramètres d'un modèle à partir de la FIM. Cet outil permet donc d'éviter des protocoles conduisant à des estimations de paramètres trop peu précises pour être interprétées. Il est aussi possible de prédire la puissance du test de Wald pour détecter un effet de covariable catégorielle à partir de la FIM, la statistique du test de Wald étant $S_W(\beta) = \beta/SE(\beta)$ avec $\beta$ l'effet de covariable estimé (Retout et al. 2007). Par extension, il est possible de prédire le nombre de sujets nécessaires pour atteindre une certaine puissance de ce test de Wald (classiquement 80 ou 90 %). De même, la puissance du test d'équivalence entre deux groupes peut être prédite à partir de la FIM (Nguyen et al. 2012).

En plus de l'estimation de la matrice individuelle et de population, la matrice d'information bayésienne pour les effets aléatoires peut être évaluée (Mentré et al. 1995). A partir de cette matrice bayésienne, il est possible de prédire la régression vers la moyenne des paramètres individuels ou $\eta$-shrinkage (Combes et al. 2013).
1.2.1.1. Méthodes d'évaluations de la FIM

Comme énoncé dans la section 1.1.1.3., la vraisemblance, et donc la FIM n'ont pas de forme
analytique dans les MNLEM. Plusieurs méthodes ont été proposées pour évaluer la FIM. La
méthode FO (First Order), première méthode qui a été proposée pour évaluer la FIM dans les
MNLEM, utilise une linéarisation du premier ordre du modèle autour des effets aléatoires
(Mentré et al. 1997). Initialement proposée dans le cadre d'un modèle d'erreur additif, cette
méthode a été étendue à l'évaluation de la FIM pour MNLEM avec erreur proportionnelle (ou
combinée) (Retout et al. 2002). Cette méthode est très efficace en général et intégrée à la plupart
des logiciels d'optimisation de protocoles basés sur la FIM (Nyberg et al. 2015). Cependant, la
FIM prédite, après calcul par FO, peut présenter certaines limites en cas de variabilités élevées,
ou modèles très complexes, et particulièrement dans les modèles pour données discrètes. De
nombreuses méthodes alternatives ont donc été étudiées pour évaluer la FIM pour ces modèles
(Waite and Woods 2015). Des méthodes récentes ont été développées dans le but d'évaluer la
FIM pour des modèles de données continues et discrètes. Parmi celles-ci, des méthodes sont
fondées sur l'utilisation de simulations de Monte Carlo (MC), combinées à la quadrature de
Gauss adaptative (AGQ) ou au Monte Carlo Hamiltonien (HMC).

L'approche de Monte Carlo-quadrature de Gauss adaptative MC-AGQ a d'abord montré de
bonnes performances pour l'évaluation de la FIM dans le cadre de modèles présentant une non-
linéarité importante (Nguyen and Mentré 2014), puis pour des modèles de données discrètes :
binaires, binaires répétées, comptage, temps d'événements répétés (Ueckert and Mentré 2017).
En termes de temps de calcul, cette méthode est efficace pour des modèles avec peu d'effets
aléatoires mais présente un coût croissant avec le nombre d'effets aléatoires. Cette méthode a
donc été étendue très récemment : l'approche QRMC-AGQ (Quasi Random Monte Carlo-AGQ
permet de réduire le temps de calcul tout en conservant une évaluation correcte de la FIM
(Ueckert and Mentré 2015).

L'approche de Monte Carlo-Monte Carlo Hamiltonien MC-HMC a aussi montré de bons
résultats pour les modèles complexes tels que les modèles de données discrètes, en conservant
un temps de calcul acceptable avec un nombre important d'effets aléatoires (Riviere et al. 2016).
Cette approche est codée en R dans la librairie MIXFIM (Riviere and Mentré 2015), et nécessite
l'écriture de modèles en Stan, langage initialement développé pour les approches bayésiennes
(Stan Development Team 2016).
### 1.2.1.2. Critères d'optimalité

Quelle que soit la méthode utilisée pour évaluer la FIM, il est nécessaire de comparer les différents protocoles évalués selon un critère. Lorsque l'on cherche à optimiser un protocole, il est possible de comparer directement les SE prédites ou la puissance prédite du test de Wald, mais il est plus commun d'utiliser des critères d'optimalité calculés à partir de la FIM. Le calcul de la FIM nécessite de spécifier la structure du modèle et de supposer ses valeurs de paramètres. Les connaissances sur le modèle et les paramètres sont généralement obtenues à partir de précédentes études. L'ensemble des critères définis évoqués dans cette section sont utilisés dans le cadre de planification de protocole locale, c'est-à-dire avec un modèle et ses valeurs de paramètres définis.

De manière générale, maximiser l'information de Fisher revient à minimiser la variance autour de l'estimation de paramètres, et donc à maximiser la précision d'estimation. Le critère de D-optimalité $\varphi_D$, largement utilisé, est directement calculé à partir du déterminant de la FIM (Kiefer 1959). Le protocole pour lequel le critère D est maximisé est celui pour lequel la précision d'estimation attendue est la meilleure en considérant l'ensemble des paramètres $\psi$, dont le vecteur est de dimension P. Le critère D est défini par l'équation 13 :

$$ \varphi_D(\Xi) = \det(\mathcal{M}(\psi, \Xi))^{1/P}. \quad (13) $$

D'autres critères peuvent être utilisés (Atkinson et al. 2007). Par exemple, le critère A est calculé à partir de la trace de la FIM et le critère E à partir des valeurs propres. Le choix du critère d'optimalité peut dépendre de l'objectif de l'étude. Par exemple, lorsque l'on cherche à optimiser spécifiquement la précision d'estimation d'un paramètre ou d'un sous-ensemble de paramètres, le critère $D_s$ ($\varphi_{D_s}$) peut être choisi (Karlin and Studden 1966). Le calcul de ce critère (équation 14), proche de celui du critère D, vise à maximiser l'information sur un sous-ensemble de paramètres d'intérêt ($\psi_s$, de longueur $S$) :

$$ \varphi_{D_s}(\Xi) = \left( \frac{\det(\mathcal{M}(\psi, \Xi))}{\det(\mathcal{M}(\psi_T, \Xi))} \right)^{1/S}, \quad (14) $$

avec $\psi$ composé du vecteur de paramètres d'intérêt $\psi_s$ et du vecteur des autres paramètres $\psi_T$. Ce choix de critère peut être intéressant dans diverses situations, par exemple si l'objectif d'une étude est de détecter un effet de covariable par test de Wald, la puissance de ce test étant dépendante de la précision d'estimation de cet effet (voir section 1.2.1.). Cependant, négliger la précision d'estimation des autres paramètres d'un modèle est déconseillé. Un design $D_s$-
optimal peut conduire à des situations où les autres paramètres du modèle ne sont pas identifiables (Sölkner 1993). Il est donc préférable d'utiliser des critères mélangeant le critère $D_S$ à un autre critère. Par exemple, le critère $DD_S$ (Atkinson and Bogacka 1997) combine les critères $D$ et $D_S$, il s'écrit (équation 15) :

$$ \varphi_{DD_S}(\Xi) = \det(\mathcal{M}(\psi_T, \Xi))^{\frac{1-\alpha}{P-S}} \times \left( \frac{\det(\mathcal{M}(\psi, \Xi))}{\det(\mathcal{M}(\psi_T, \Xi))} \right)^{\frac{\alpha}{S}}, \quad (15) $$

où $\alpha$, compris entre 0 et 1, est un terme d'équilibre entre l'information sur les paramètres $\psi_T$ et sur les paramètres d'intérêt. Plus récemment, le critère $D_S$, combiné à un autre critère (le critère V, utilisé dans le but de minimiser la variance des prédictions), a montré son utilité pour le choix de doses à évaluer dans le cadre d'une combinaison de deux molécules (Papathanasiou et al. 2019).

Afin de comparer facilement deux protocoles, il est commun de calculer l'efficacité pour un critère $X$, défini comme le ratio de ce critère pour le protocole $\Xi$ par rapport à un protocole de référence $\Xi^*$, qui peut être le protocole optimal (équation 16) :

$$ E_X(\Xi) = \frac{\varphi_X(\Xi)}{\varphi_X(\Xi^*)}. \quad (16) $$

### 1.2.1.3. Incertitude du modèle et des paramètres

Dans le cas où le modèle final (c'est-à-dire celui effectivement utilisé pour décrire les données) et/ou les valeurs de paramètres sont différents du modèle spécifié a priori au moment du choix de protocole, le protocole choisi peut ne pas être assez informatif et conduire à des estimations peu précises de paramètres. Il est donc nécessaire de développer des approches de planification de protocole robuste, c'est-à-dire des méthodes prenant en compte l'incertitude sur le choix du modèle et des valeurs de paramètres. En cas d'incertitude, une première approche consiste à effectuer une analyse de sensibilité aux valeurs de paramètres, c'est-à-dire d'évaluer, voire optimiser un protocole avec différentes valeurs de paramètres. Des critères d'optimalité robustes, supposant une incertitude sur la valeur des paramètres, ont été proposés, par exemple le critère ED ou le critère HCD (Pronzato and Walter 1985; Foo and Duffull 2010). D'autres approches, plus récentes, ont été utilisées pour prendre en compte l'incertitude du modèle, dans le contexte où un protocole commun était choisi dans une étude d'un « cocktail » de médicaments (Nguyen et al. 2016). L'approche proposée dans ce travail était de supposer plusieurs modèles candidats et d'allouer un poids à chacun, sur le principe du *model averaging*
(Buckland et al. 1997; Buatois et al. 2018) et en utilisant la théorie de design optimal composé (Atkinson et al. 2007). Pour cette approche, différents modèles candidats doivent être définis et le critère robuste (ou composé) CD $\varphi_{CD}$ (équation 17) peut être utilisé afin de trouver le protocole optimal commun aux modèles proposés. Le calcul de ce critère est en fait le produit pondéré des critères D de chaque modèle :

$$ \varphi_{CD}(\Xi) = \prod_{m=1}^{M} (\varphi_{D,m}(\Xi))^{w_m}, \qquad (17) $$

avec m, un modèle candidat parmi les M modèles considérés et $w_m$ le poids attribué à m ($w_m$ tel que $0 \le w_m \le 1$ et $\sum_{m=1}^{M} w_m = 1$). Un $w_m$ plus proche de 1 signifie que le modèle m pourrait mieux décrire les données que les autres modèles candidats.

Appliqué aux modèles de Poisson pour données de comptage, un travail d'optimisation de protocole robuste, basé sur la FIM calculée par MC-HMC, a été conduit de manière à prendre en compte l'incertitude sur les paramètres pour cinq modèles différents, puis sur les paramètres et considérant l'incertitude du modèle à choisir parmi les cinq (Loingeville et al. 2019). La prise en compte de l'incertitude des valeurs de paramètres n'a pas eu d'impact sur le protocole optimal dans la plupart des cas. Néanmoins, le fait de pré-spécifier un modèle qui n'est pas le modèle finalement utilisé peut conduire à de mauvais choix de protocoles, conduisant à des D-efficacités inférieure à 30% par rapport au protocole optimal (pour la même quantité de données recueillies). L'utilisation du critère CD a conduit à un choix de doses robuste, le protocole CD-optimal conduisant à une D-efficacité d'au moins 80% pour les 5 modèles candidats. Concernant les travaux de cette thèse, ce type d'approche robuste par rapport au choix du modèle a donc été retenue (section 2.1.).

### 1.2.1.4. Logiciels et algorithmes d'optimisation

Le premier logiciel qui a été développé pour optimiser des protocoles basés sur la FIM de population est PFIM (Retout et al. 2001; Dumont et al. 2018), dont la version 4.0. est disponible sur R et en version interfacée. D'autres logiciels de design ont ensuite été développés, par exemple POPT (Duffull et al. 2005), PopED (Foracchia et al. 2004), PopDes (Gueorguieva et al. 2007) et PkStaMp (Aliev et al. 2012). Les différents logiciels précédemment cités ont été comparés en 2014 et ont montré des résultats proches en termes de FIM prédite par méthode FO pour modèles de données continues (Nyberg et al. 2015). En outre, ce travail a aussi montré que l'utilisation de la matrice de Fisher diagonale par blocs est plus proche des résultats
empiriques (obtenus par simulations) que la matrice de Fisher complète. Parmi les spécificités
de chaque logiciel au moment de cette comparaison, PopED permettait de prendre en compte
l'incertitude sur la valeur des paramètres et PFIM permettait de prendre en compte la FIM
observée d'une précédente étude (aujourd'hui, PopED a aussi cette fonctionnalité).

Les algorithmes développés pour l'optimisation de protocoles peuvent être séparés en deux
catégories : ceux qui permettent de trouver un protocole dans un espace continu et ceux utilisant
un espace de protocoles discret. Un algorithme utilisé pour une recherche de protocole optimal
dans un espace continu est l'algorithme Simplex (ou du simplexe), développé à partir de 1947
pour des fonctions linéaires par Dantzig (Nash 2000) puis étendus à des problèmes
d'optimisation non-linéaires (Nelder and Mead 1965). Le Simplex est l'algorithme intégré à
PFIM 4.0. pour optimiser un protocole dans un espace continu. D'autres algorithmes ont été
développés pour l'optimisation de protocoles pour MNLEM. Dans PopED, une combinaison
de deux algorithmes (un permettant une recherche aléatoire globale et un autre, fondé sur le
gradient stochastique, permettant d'affiner le résultat) est intégrée afin de chercher le protocole
optimal dans un espace continu (Foracchia et al. 2004). Dans ce même but, l'utilisation du recuit
simulé est une autre possibilité, il est notamment intégré à WinPOPT, la version interfacée de
POPT (Duffull et al. 2002).

L'optimisation dans un espace discret consiste à déterminer le ou les designs élémentaires les plus informatifs à partir des possibles éléments de designs pré-définis tels que temps de mesure ou doses possibles. Dans PFIM 4.0., l'algorithme de Fedorov-Wynn (FW) est utilisé pour ce type d'optimisation (Fedorov 1972; Wynn 1972). Dans PopED, une version modifiée de cet algorithme est intégrée. L'algorithme FW permet aussi de définir la proportion d'individus qu'il faut allouer aux protocoles élémentaires. Simplex et FW, tous deux présents dans PFIM 4.0., ont été comparés. Il a été montré que FW est moins sensible que Simplex aux valeurs initiales et que FW est plus performant pour les problèmes d'optimisation complexes (Retout et al. 2007). D'autres algorithmes existants pourraient être adaptés à l'optimisation de protocoles d'études analysées par modèles mixtes. C'est l'objectif de l'un des travaux de cette thèse que de développer un algorithme pour optimiser ces protocoles (section 2.2.1.).

1.2.2. Simulations d'essais cliniques

La méthode par simulations consiste à estimer les paramètres d'un modèle à partir de différents
jeux de données simulés à partir d'un modèle et du protocole à évaluer. Cette approche permet
de juger de l'informativité d'un protocole. Comme pour l'approche basée sur la FIM, elle
permet d'évaluer la précision d'estimation des paramètres, en calculant l'écart type empirique des valeurs estimées. Cette approche permet aussi de vérifier la justesse d'estimation d'un paramètre ; un protocole informatif permet de minimiser le biais d'estimation des paramètres, calculé comme l'écart moyen entre valeur estimée et valeur réelle. Il est aussi souvent reporté l'erreur quadratique moyenne qui mesure à la fois la précision et le biais d'estimation.

Cette approche possède différents types d'applications. Au début des années 1990, sa première utilité était de comparer le biais et la précision des paramètres entre différents protocoles en termes de nombre et d'allocation des temps de mesure en pharmacocinétique (al-Banna et al. 1990). Cette approche a aussi rapidement servi à évaluer différents régimes de doses et d'évaluer l'intérêt d'une augmentation de doses en termes de pourcentage d'individus qui atteignent un seuil d'efficacité désiré (Green and Duffull 2003). De manière similaire, les EBE obtenus après simulations peuvent être utilisés pour réaliser des tests et comparer les puissances ou les erreurs de type 1 obtenues. Un des travaux de recherche de ce doctorat utilise cette approche dans ce but, avec une application en oncologie (section 2.3.). L'approche par CTS est aussi utile pour vérifier les prédictions issues de l'approche basée par la FIM. Cela est notamment utile pour comparer les logiciels (Nyberg et al. 2015) ou pour développer de nouvelles méthodes d'évaluation de la FIM (Nguyen and Mentré 2014).

L'approche par CTS est néanmoins peu adaptée à l'optimisation de protocoles. En effet, pour un protocole à évaluer, après avoir simulé un grand nombre de jeu de données, il est nécessaire d'estimer les paramètres du modèle à partir de chaque jeu de données, ce qui est coûteux en temps de calcul. Par conséquent, elle ne permet d'évaluer qu'un nombre limité de protocoles et le risque avec cette méthode est de passer à côté d'autres protocoles bien plus informatifs.

En pratique, lorsqu'on optimise un protocole, la meilleure approche semble donc de trouver le protocole le plus informatif par une approche basée sur la FIM. Cependant, comme expliqué dans la section 1.2.1., la FIM est basée sur la vraisemblance, elle-même calculée par approximation dans le cadre des MNLEM. De plus, l'inverse de la FIM donne la borne inférieure de la matrice de variance-covariance de tout estimateur non biaisé des paramètres. Il est donc conseillé de vérifier par CTS que le protocole optimal permet effectivement d'atteindre l'objectif souhaité.

## 1.3. Phagothérapie

Les antibiotiques ont drastiquement diminué la mortalité associée aux maladies infectieuses.
Leur utilisation massive et parfois inadaptée, en thérapeutique humaine ou chez les animaux d'élevage a conduit à l'émergence de bactéries résistantes à ces molécules. Le rapport O'Neill, établi en 2014 alertait sur l'augmentation de l'antibiorésistance et ses conséquences (O'Neill 2014). Il est estimé que la résistance bactérienne est responsable d'environ 700 000 décès par an dans le monde. Dans ce rapport, il est aussi prédit qu'en 2050, la résistance bactérienne pourrait être la première cause de mortalité avec 10 millions de décès par an. Dans ce contexte, la résistance bactérienne est particulièrement surveillée pour certaines espèces bactériennes, c'est le cas d'Escherichia coli par exemple. Il est aussi recommandé de promouvoir une meilleure hygiène et un usage plus rationnel des antibiotiques, c'est-à-dire utiliser les antibiotiques les plus adaptés, et seulement en cas de nécessité, que ce soit en thérapeutique humaine ou dans l'industrie agro-alimentaire (McEwen and Collignon 2018). De plus, la recherche de nouvelles thérapies est nécessaire, que ce soit le développement de nouveaux antibiotiques ou l'exploration d'autres voies telles que les anticorps dirigés contre les toxines bactériennes ou la phagothérapie : l'administration de virus dirigés contre les bactéries.

### 1.3.1. Phages et phagothérapie

Les bactériophages sont des virus prédateurs des bactéries. Ils ont été découverts en plusieurs étapes. Hankin a été le premier à mettre en évidence l'existence naturelle d'un agent antibactérien dans l'eau du Gange et de la Jumna à la fin du XIXème siècle (Hankin 1896). En 1917, d'Hérelle a isolé et démontré la spécificité de ces virus pour les bactéries qu'il a nommé « bactériophages », puis a proposé de les utiliser en thérapeutique (« phagothérapie ») contre les dysenteries à Shigella (d'Herelle 1917). Supplantés par les antibiotiques tels que la pénicilline à partir des années 1940, ils sont alors peu utilisés à l'échelle mondiale, hormis dans certains instituts d'Europe de l'Est où la recherche et le développement de la phagothérapie a perduré (Sulakvelidze and Morris 2001). Aujourd'hui, à cause de l'émergence de la résistance bactérienne aux antibiotiques, un regain d'intérêt pour la phagothérapie est observé dans les pays occidentaux, notamment en utilisation compassionnelle, pour des patients infectés par des bactéries toto-résistantes, c'est-à-dire résistantes à tout antibiotique sur le marché (Patey et al. 2018).

#### 1.3.1.1. Généralités à propos des bactériophages

Les bactériophages infectent spécifiquement les bactéries. Il existe quatre types de cycles d'infection des bactéries par les phages (Figure 1). Le plus intéressant en thérapeutique est le
cycle lytique : après adsorption du phage à la surface de la bactérie, l'acide nucléique est transféré dans le cytoplasme. La transcription et la traduction de gènes se fait alors dans le but de synthétiser du matériel protéique et génomique viral plutôt que bactérien, afin de produire de nouveaux bactériophages. Vient ensuite la lyse bactérienne donc la libération des nouveaux virions qui infecteront d'autres bactéries. Parmi les autres cycles d'infection, la lysogénie et pseudolysogénie consistent en une possible lyse bactérienne après une ou plusieurs divisions bactériennes. Dans le cycle chronique, les bactéries deviennent productrices de phages sans lyse (Weinbauer 2004).

<figure><img src="image_2.png" /><figcaption>Figure 1. Les différents cycles d'infection des phages (d'après Weinbauer 2004).</figcaption></figure>

Les bactériophages sont présents dans tous les écosystèmes, notamment dans notre alimentation et donc notre microbiote intestinal (Letarov and Kulikov 2009). Il est estimé que le nombre de phages sur Terre est de $10^{30}$, soit en quantité plus importante que les bactéries. La diversité génétique des phages est aussi très importante, avec $10^8$ génomes différents (Rohwer 2003). Les bactériophages participent de manière importante à l'équilibre de la biomasse bactérienne, via une relation proie (les phages) – prédateur (les bactéries), on parle aussi de couples phage bactérie par exemple, le couple 536/536_P1 où 536 est la souche bactérienne *Escherichia coli* 536 et 536_P1 un phage particulièrement virulent vis-à-vis de cette souche. Comme dans tout système proie-prédateur, des mécanismes sont mis en jeu par les bactéries pour résister aux phages, et par les phages pour contrer ces résistances. Par ailleurs, il a été démontré que les
bactériophages ont une part de responsabilité dans l'extinction des épidémies bactériennes
(Faruque et al. 2005).

1.3.1.2. Phagothérapie

La préparation de phages avant administration chez l'animal ou l'homme nécessite le recueil du ou des phages d'intérêt, leur amplification et leur purification. Généralement, le spectre d'action d'un bactériophage est réduit. Il est souvent spécifique d'une espèce bactérienne et d'un nombre limité de souches au sein de cette espèce. Certains phages sont tout de même capable d'infecter différentes espèces (Jensen et al. 1998).

Cette particularité a un impact en thérapeutique. Lorsque l'infection est chronique avec un pronostic vital peu engagé, ou lorsque la bactérie responsable est identifiée, la phagothérapie peut être intéressante. Elle l'est d'autant plus si le lieu de l'infection est un tissu où les antibiotiques diffusent mal car les phages, contrairement aux médicaments chimiques, vont pouvoir se multiplier une fois sur le lieu de l'infection. De plus, leur spectre généralement étroit permet d'éviter la destruction des microbiotes (notamment intestinal) comme lors du traitement par antibiotiques. Dans le cas d'une infection aigüe avec pronostic vital engagé à court terme, lorsque la bactérie n'est pas identifiée ou lorsque l'infection est causée par plusieurs espèces bactériennes, il semble préférable, en première intention, d'utiliser des antibiotiques à large spectre (Wunderink and Waterer 2014).

De nombreuses preuves d'efficacité de la phagothérapie existent dans la littérature, notamment chez l'animal, tant en termes de réduction de la charge bactérienne que de mortalité (Caflisch et al. 2019). Cela a été montré pour différentes bactéries et pour divers types d'infection. Chez le lapin infecté par Vibrio cholerae, responsable du choléra, une réduction importante de la charge bactérienne et une absence de symptômes a été constatée pour les animaux traités par phages, contrairement aux non traités (Bhandare et al. 2019). Dans un projet précédent conduit dans le cadre d'une collaboration entre l'Institut Pasteur et notre unité de recherche, chez un modèle murin de pneumonie à E.coli 536, l'administration d'une dose de phage 536_P1 a fait aussi bien que 5 doses d'antibiotiques en termes de réduction de la charge bactérienne et de mortalité (Dufour et al. 2015). L'administration de plusieurs espèces de bactériophages permet d'étendre le spectre d'action de cette thérapie. Par exemple, l'étude d'un cocktail de trois phages chez la souris infectée par K.pneumoniae a montré un meilleur état de santé des souris traitées par le cocktail comparativement à celles traitées par l'un des trois phages (Gu et al. 2012). Chez l'homme, la plupart des études disponibles ont été faites il y a plusieurs décennies et n'ont pas
été conduites avec la rigueur méthodologique permettant de fournir la preuve d'efficacité clinique avec la rigueur actuelle (d'Herelle 1931; Tsulukidze 1941). Les publications d'essais cliniques récents sur la phagothérapie ne sont pas nombreuses pour plusieurs raisons. Les essais cliniques sont couteux et les laboratoires pharmaceutiques ont pour l'instant peu d'intérêt à développer la phagothérapie, les essais sont donc souvent menés par des équipes académiques. En outre, les autorités de santé ont tardé à établir un cadre pour faciliter la mise en place d'études cliniques (Furfaro et al. 2018; Górski et al. 2018). Certaines études cliniques récentes ont toutefois été menées et ont révélé une efficacité plus ou moins importante de la phagothérapie. Par exemple, chez des patients atteints d'otites chroniques à P.aeruginosa, il a été montré l'efficacité de l'administration d'une dose de phages sur les symptômes de la maladie et sur le compte de bactéries (Wright et al. 2009). En revanche, l'utilisation d'un cocktail de phages par voie orale n'a pas permis d'améliorer l'état clinique d'enfants présentant une diarrhée nécessitant une hospitalisation, par rapport au traitement de référence (Sarker et al. 2016). Par ailleurs, l'étude récente PhagoBurn n'a pas permis de conclure à l'efficacité de la phagothérapie chez les grands brûlés (Jault et al. 2019). Les raisons mises en avant pour expliquer ces échecs thérapeutiques sont un mauvais choix de schéma thérapeutique ou de critères d'inclusion. En effet, la phagothérapie semble plus pertinente pour les infections lorsque le pathogène présente une multi-résistance aux antibiotiques, et une sensibilité au bactériophage administré. Cependant, la phagothérapie est utilisée pour traiter les infections bactériennes chez l'homme depuis plusieurs décennies en Europe de l'Est, notamment dans l'Institut Eliava, situé en Géorgie (Kutateladze and Adamia 2008).

Les bactériophages présentent un profil de toxicité intéressant en thérapeutique. Ils n'entraînent pas de lésion au niveau des tissus dans lesquels ils se distribuent. Ils ne sont pas capables d'infanter les cellules eucaryotes et ne présentent donc pas de reprotoxicité ou de mutagénicité sur les cellules humaines. Bien que les phages puissent induire une réaction immunitaire aiguë en cas de mauvaise purification de la préparation, celle-ci est généralement réduite ou inexistante, les phages faisant partie de notre environnement en permanence. Chez un modèle murin d'infection pulmonaire, il a été montré que la réponse immunitaire innée n'est pas plus importante par phagothérapie comparée à l'antibiothérapie, et cela pour deux couples dont 536/536_P1 (Dufour et al. 2019). Certaines études suggèrent la possibilité d'une réponse adaptative avec notamment production d'anticorps neutralisants dirigés contre les phages (Ochs et al. 1992; Rubinstein et al. 2000). Les principaux effets toxiques observés lors d'administration de phages sont dus à la libération de toxines lors de la lyse bactérienne. Il a
néanmoins été montré que la libération d'endotoxines n'était pas supérieure à celle induite par
antibiothérapie (Matsuda et al. 2005). Cela a par exemple été montré pour le couple 536/536_P1
où l'exposition aux endotoxines était supérieure après traitement par céphalosporines ou par
imipénème par rapport à la phagothérapie (Dufour et al. 2017). Toutes ces preuves de non
toxicité des phages permettent un assouplissement progressif des conditions fixées par les
autorités de santé pour les études cliniques des bactériophages et une prescription facilitée en
usage compassionnel.

1.3.2. Modélisation de l'interaction phage/bactérie

Certains échecs des études de phagothérapie pourraient être causés par une mauvaise
compréhension de la relation entre le schéma d'administration des phages et leur efficacité
(Payne and Jansen 2001), conduisant à de mauvais choix de schémas d'administration des
phages. Par exemple, si la dose de phages est trop faible ou si la voie d'administration conduit
à une trop faible concentration de phages sur le lieu de l'infection, les bactéries peuvent
continuer à se multiplier et entraîner des lésions avant que les phages soient en quantité
suffisante pour contrôler la charge bactérienne. Le bactériophage étant un médicament, il est
nécessaire de bien comprendre les relations entre dose, concentration et efficacité. Ceci peut
être réalisé à l'aide de modèles mathématiques, plusieurs phénomènes sont alors à prendre en
compte.

De la même manière qu'une molécule classique, la PK propre des phages (absorption,
distribution et élimination) doit être étudiée. Cela se fait par administration d'une quantité de
phages connue chez l'individu non infecté et suivi de la concentration de phages dans le sang
et les tissus d'intérêt. Les bactériophages, contrairement aux molécules classiques, sont
capables de se multiplier sur le lieu d'une infection. Les processus spécifiques de la
multiplication virale sont donc à prendre en compte dans la manière de modéliser la dynamique
des phages. La capacité des phages à s'adsorber aux bactéries, le temps nécessaire pour lyser
une bactérie et le nombre de virus relargués lors de la lyse sont autant d'éléments à prendre en
compte dans le modèle. Ils peuvent être estimés à partir de données d'expériences *in vitro* et *in vivo*, où phages et bactéries sont mis en contact. La capacité des bactériophages à se multiplier
dépend de la probabilité de rencontre entre phages et bactéries et donc des concentrations
bactérienne et virale sur le lieu d'infection. Il est donc nécessaire de prendre en compte la
cinétique bactérienne. La multiplication « naturelle » des bactéries se quantifie idéalement en
l'absence de traitement, en suivant l'évolution de la charge bactérienne sur le lieu d'infection.
En outre, pour certains couples phages/bactéries, le phénomène de résistance doit aussi être intégré au modèle (Skurnik and Strauch 2006).

Afin d'émettre des recommandations de dose, de voie et de nombre d'administrations de bactériophages à partir d'un modèle, il est donc crucial d'intégrer à la fois la cinétique bactérienne, la PK des phages et les interactions entre phages et bactéries. De plus, ces éléments font partie de la liste des informations qu'il est recommandé de fournir dans les études cliniques sur la phagothérapie (Abedon et al. 2011).

Bien que les premières données d'expérimentation de phages datent d'il y a plus de 100 ans, la littérature concernant la modélisation des phages est peu dense, en particulier pour la modélisation de l'interaction phage/bactérie. Pour comprendre ce type de modèles, il est nécessaire de détailler certains standards de modélisation de cinétique bactérienne d'une part, et de dynamique virale d'autre part. Ici, nous n'évoquerons pas la modélisation de la propagation des infections au niveau populationnel mais seulement la modélisation intra-hôte.

### 1.3.2.1. Modélisation de la cinétique bactérienne

Les modèles mathématiques de cinétique bactérienne, fréquemment écrits sous la forme d'équations différentielles, permettent de représenter et de quantifier l'effet des médicaments, le plus fréquemment des antibiotiques, sur la dynamique bactérienne ainsi que sur l'apparition de résistances.

Dans la littérature, la description de cinétiques bactériennes par modélisation est retrouvée plus fréquemment *in vitro* (Mouton et al. 1997; Nielsen et al. 2007) qu'*in vivo* (Guo et al. 2011). Ceci est explicable par plusieurs facteurs. Premièrement, lorsque les infections bactériennes sont localisées dans un tissu spécifique, il est souvent nécessaire de sacrifier l'animal pour relever la charge bactérienne sur le site d'infection. Dans ce contexte, la récolte de données longitudinales et donc la modélisation de l'évolution de la charge bactérienne est complexe. Deuxièmement, les interactions entre le système immunitaire de l'hôte et les bactéries sont aussi à prendre en compte, mais les marqueurs de l'immunité ne sont pas toujours suivis. Enfin, la variabilité inter-expérience et inter-individuelle est généralement beaucoup plus importante *in vivo* que dans les études *in vitro*. L'utilisation de modèles mixtes est intéressante parce qu'ils permettent de prendre en compte la variabilité, même lorsque les données sont éparses, comme expliqué section 1.1.1.

Beaucoup de modèles de cinétique bactérienne différents ont été publiés, une revue datant de 2013 liste les modèles les plus retrouvés en pharmacométrie des antibactériens (Nielsen and
Friberg 2013). Nous allons présenter les principaux modèles de cinétique bactérienne permettant de comprendre les travaux de ce doctorat.

Dans un système *in vitro* où les bactéries ont des ressources et de l'espace pour se multiplier, leur cinétique peut être décrite par l'équation 18 :

$$ \frac{dB}{dt} = k_g B - k_d B, \qquad (18) $$

où B est la charge bactérienne, $k_g$ et $k_d$ sont respectivement les paramètres de croissance et de mort naturelle bactérienne. Ici, la condition initiale $B(0)$ peut être équivalente à l'inoculum bactérien choisi $I_B$ divisé par le volume de distribution estimé $V_B$ (cette condition initiale peut aussi, plus simplement, être estimée par le paramètre $B_0$, sans prendre en compte l'inoculum). Les paramètres $k_g$ et $k_d$ étant difficiles à identifier, l'écriture de l'équation de la cinétique bactérienne est alors fréquemment simplifiée (équation 19) :

$$ \frac{dB}{dt} = \alpha B, \qquad (19) $$

avec $\alpha$, la croissance « nette » bactérienne équivalente à $k_g - k_d$.

Dans la plupart des expériences l'espace de croissance bactérienne est délimité, il est alors observé un plateau $B_{max}$, la cinétique bactérienne s'écrit alors (équation 20) :

$$ \frac{dB}{dt} = \alpha B \left(1 - \frac{B}{B_{max}}\right). \qquad (20) $$

Il existe différentes manières d'écrire la perturbation de cette cinétique par un antibactérien ABx (Figure 2, équation 21) :

<figure><img src="image_3.png" /></figure>
Figure 2. Représentation schématique d'un modèle de cinétique bactérienne (d'après Friberg and Guedj 2019).

$$ \frac{dB}{dt} = \alpha B \left(1 - \frac{B}{B_{max}}\right) - ABx_{eff}B, \qquad (21) $$

où $ABx_{eff}$ est l'effet de cet antibactérien sur la cinétique bactérienne. Pour un modèle PKPD non-linéaire classique : $ABx_{eff} = \frac{E_{max}C}{EC_{50}+C}$ avec $E_{max}$ l'effet maximal et $EC_{50}$ la concentration $C$ de l'antibactérien $ABx$ qu'il est nécessaire d'atteindre pour obtenir la moitié de cet effet maximal. Dans certains cas, de manière moins mécanistique, cet effet peut être linéaire : $ABx_{eff} = \delta \times C$.

Lors de l'administration d'un médicament aussi, il est commun de pouvoir identifier l'apparition de souches résistantes, il est alors possible de décrire l'évolution de la charge bactérienne des deux sous-populations sensible (S, équation 22) et résistante (R, équation 23) :

$$ \frac{dS}{dt} = \alpha S \left(1 - \frac{S}{B_{max}}\right) - ABx_{eff}S \qquad (22) $$

$$ \frac{dR}{dt} = w\alpha R \left(1 - \frac{R}{B_{max}}\right), \qquad (23) $$

avec $B = R + S$, la charge bactérienne totale mesurée et $w$ un paramètre de valeur sélective de populations résistantes par rapport aux sensibles (c'est-à-dire le ratio entre la vitesse nette de croissance des résistants par rapport à celle des sensibles). Les conditions initiales de ce système $S(0)$ et $R(0)$ peuvent être estimées (respectivement par les paramètres $S_0$ et $R_0$), elles peuvent aussi être paramétrées telles que $S(0) = \frac{I_B}{V_B} \times (1 - f)$ et $R(0) = \frac{I_B}{V_B} \times f$, avec $f$ la fréquence de fluctuation, c'est-à-dire la proportion de bactéries résistantes en l'absence de traitement. Si les données à disposition le permettent, il est possible, à partir de modèles plus élaborés, de prendre en compte la pression de sélection due à l'antibactérien, c'est-à-dire le taux de mutation des souches sensibles afin qu'elles deviennent résistantes.

Lorsqu'on modélise la cinétique bactérienne *in vivo*, il peut être nécessaire de prendre en compte l'effet de l'immunité de l'hôte, en particulier l'immunité innée, médiée par les macrophages. Afin de prendre en compte cette composante dans le modèle, l'idéal est de suivre un ou des marqueurs de l'immunité de façon longitudinale, puis de considérer un système à plusieurs équations intégrant l'évolution de la charge bactérienne, celle de l'immunité et les interactions entre les deux (Roach et al. 2017). Lorsque aucun marqueur de l'immunité n'est
suivi durant l'étude, une façon simple de procéder est de considérer un temps nécessaire pour l'activation de l'immunité $t_I$ (équations 24 et 25) :

$$ \text{avant } t_I, \qquad \frac{dB}{dt} = \alpha B \left(1 - \frac{B}{B_{\text{max}}}\right) \qquad (24) $$

$$ \text{après } t_I, \qquad \frac{dB}{dt} = \alpha B \left(1 - \frac{B}{B_{\text{max}}}\right) - I_{\text{eff}}, \qquad (25) $$

avec $I_{\text{eff}}$ l'effet de l'immunité, qui peut s'écrire sous différentes formes.

Dans certaines études *in vivo*, si le décès est observé précocément pour un nombre conséquent d'individus, l'utilisation de modèles conjoints peut être un atout. En effet, la mortalité ici est liée à l'évolution de l'infection, les individus avec de hautes charges bactériennes ou un système immunitaire défaillant sont donc plus à risque de décès. Le risque est ici d'avoir une trop grande représentation des données d'individus survivants. Comme évoqué précédemment (voir section 1.1.1.2.), ceci peut biaiser l'estimation des paramètres d'un modèle incluant uniquement la description du marqueur longitudinal. L'utilisation d'un modèle conjoint (marqueur longitudinal-survie) permet de réduire ce biais.

### 1.3.2.2. Modélisation de la dynamique virale

Les bactériophages sont des virus dont les cellules cibles sont les bactéries. Afin de caractériser la dynamique des phages, l'utilisation des modèles classiques de dynamique virale, c'est-à-dire ceux développés pour les virus pathogènes, est nécessaire. Ces modèles ont initialement été développés pour des infections virales chroniques tels que le virus de l'immunodéficience humaine (Wei et al. 1995), puis ont été étendus à des infections aiguës telles que la grippe (Smith and Perelson 2011). Contrairement aux modèles de cinétique bactérienne, les modèles de dynamique virale prennent usuellement en compte la dynamique des cellules cibles des virus car les virus ont besoin de ces cellules pour se multiplier. La dynamique virale peut généralement être décrite par le système d'équations 26-28 (Perelson 2002) :

$$ \frac{dS}{dt} = \lambda - dS - \beta VS \qquad (26) $$

$$ \frac{dI}{dt} = \beta VS - \delta I \qquad (27) $$

$$ \frac{dV}{dt} = pI - cV. \qquad (28) $$
Ce modèle inclut à la fois les cellules cibles, aux stades non infectées (sensibles) S et infectées I, ainsi que la quantité de virus V. Les conditions initiales de ce système sont généralement $S(0) = \frac{\lambda}{d}$, $I(0) = 0$ et $V(0) = I_V$, avec $I_V$ l'inoculum viral (si le virus est inoculé à $t = 0$). Dans ce système, les cellules cibles sont naturellement générées à la vitesse $\lambda$ et subissent une élimination d'ordre 1 qui dépend du paramètre $d$. Le paramètre $\beta$ est ici le taux d'infectivité, qui caractérise la vitesse d'infection et la probabilité de rencontre entre les virus et les cellules cibles. Les cellules infectées sont ensuite éliminées par le virus à la vitesse $\delta$. De nouveaux virus sont produits à une vitesse $p$ par les cellules infectées, puis sont éliminés naturellement à la vitesse $c$.

Le compartiment de cellules infectées est fréquemment séparé en deux compartiments $I_1$ et $I_2$, avec les cellules $I_1$ qui sont latentes (non productrices) et $I_2$ productrices de virus (Figure 3, équations 29-32).

<figure><img src="image_4.png" /><figcaption>Figure 3. Représentation schématique d'un modèle de dynamique virale (d'après Friberg and Guedj 2019).</figcaption></figure>

$$ \frac{dS}{dt} = \lambda - dS - \beta VS \quad (29) $$

$$ \frac{dI_1}{dt} = \beta VS - eI_1 \quad (30) $$

$$ \frac{dI_2}{dt} = eI_1 - \delta I_2 \quad (31) $$

$$ \frac{dV}{dt} = pI_2 - cV. \quad (32) $$

Ce modèle permet de prendre en compte le délai nécessaire au virus pour « utiliser » la cellule cible avant que celle-ci puisse produire de nouveaux virions. Le paramètre $e$ caractérise ce
délai, mais l'estimation de ce paramètre est souvent difficile. Des expériences *in vitro* peuvent néanmoins aider à le déterminer.

D'autres modèles de dynamique virale plus élaborés ont été développés, incluant notamment des effets de traitement antiviraux, voire de cellules réfractaires. Enfin, comme en cinétique bactérienne, intégrer des marqueurs de la réponse immunitaire peut aider à expliquer les phénomènes observés *in vivo* (Friberg and Guedj 2019). Par exemple, dans le cadre d'une étude sur la maladie à virus Ebola chez le macaque, les lymphocytes T et trois cytokines ont été mesurés de manière longitudinale et intégrés au système d'équations. De plus, le lien entre charge virale, marqueurs immunologiques et survie a été montré par l'utilisation d'un modèle conjoint (Madelain et al. 2018).

### 1.3.2.3. Modèles mathématiques existants pour la description de l'interaction phage/bactérie

Les premiers modèles mathématiques des bactériophages reprenant les standards de cinétique bactérienne et de dynamique virale ont été développés à partir des années 2000 (Payne 2000). A partir d'un modèle composé d'un système de trois équations différentielles : bactéries non infectées S (équation 33), bactéries infectées I (équation 34) et phages V (équation 35) :

$$ \frac{dS}{dt} = \alpha S - \beta VS \quad (33) $$

$$ \frac{dI}{dt} = \alpha I + \beta VS - \delta I \quad (34) $$

$$ \frac{dV}{dt} = b\delta I - \beta VS - k_e V. \quad (35) $$

Dans ce modèle, les conditions initiales sont généralement $S(0) = S_0$, $I(0) = 0$ et $V(0) = 0$, lorsque les phages sont administrés après le début du suivi. Ce modèle a permis d'avancer l'importance du temps d'inoculation des phages pour le contrôle de l'infection bactérienne. Les paramètres clés caractérisant les interactions entre phages et bactéries sont le taux d'infectivité $\beta$ et la vitesse de lyse $\delta$. Le paramètre $b$, généralement absent en dynamique virale classique, caractérise le *burst size*, c'est-à-dire le nombre de virions libérés à la fin du cycle d'infection, lors de la lyse de la bactérie. Ce premier modèle mathématique constitue une base théorique pour la modélisation de l'interaction phage-bactérie. Il est bien adapté à la description d'expériences *in vitro* (Cairns et al. 2009). Il n'inclut cependant pas certaines complexités que l'on peut rencontrer lors de l'utilisation de phages *in vivo* (Weld et al. 2004). Quelques années
plus tard, un modèle étendait ces équations pour intégrer l'émergence de souches bactériennes résistantes au phage et une expansion bactérienne limitée par les ressources disponibles dans le système (Levin and Bull 2004). Certains modèles intègrent aussi la pression de sélection exercée par les bactériophages, par le biais d'un paramètre de taux de mutation des bactéries sensibles vers des bactéries résistantes (Cairns et al. 2009). Un modèle récent intègre aussi la dynamique de l'immunité de l'hôte et son effet sur la charge bactérienne (Roach et al. 2017). La prochaine étape pourrait être l'utilisation de modèles mixtes, de manière à intégrer et expliquer certaines variabilités observées *in vivo*. L'utilisation de modèles conjoints, c'est-à-dire liant la survie à la cinétique bactérienne ou aux marqueurs de l'immunité est aussi une piste à explorer. La modélisation permet une meilleure compréhension et la quantification de l'interaction phage-bactérie, et pourrait *in fine* aider au choix du schéma d'administration de la phagothérapie chez l'homme. Un objectif de ces travaux de thèse est d'utiliser les modèles mixtes afin de caractériser l'interaction entre bactériophages et bactéries (section 3.).

## 1.4. Objectifs de la thèse

La première partie de ce doctorat porte sur l'évaluation et l'optimisation de protocoles des études longitudinales analysées par modèles non linéaires à effets mixtes. Parmi les travaux de recherche réalisés, deux sont des développements de méthodes pour la planification de protocoles basée sur l'utilisation de la matrice d'information de Fisher et un autre travail porte sur la comparaison de protocoles par simulations, appliquée en oncologie.

Le premier travail s'inscrit dans la continuité de travaux de recherche menés au sein de notre équipe ; d'une part dans le développement de méthodes d'évaluation de l'information de Fisher adaptées aux modèles de données non continues (section 1.2.1.1.), et d'autre part dans le développement de méthodes de protocoles robustes, capables de prendre en compte l'incertitude du modèle et/ou des valeurs de paramètres (section 1.2.1.3.). Dans le cadre de l'analyse de données binaires répétées, notre objectif était donc de développer une méthode de planification prenant en compte l'incertitude du modèle, en évaluant la matrice d'information de Fisher par méthode de Monte Carlo-Monte Carlo Hamiltonien. Un objectif secondaire était d'évaluer la pertinence de l'approche développée, c'est pourquoi nous avons réalisé des simulations d'essais cliniques (section 2.1).

Nous avons vu précédemment la nécessité de développer de nouveaux algorithmes à l'optimisation de protocoles pour les modèles mixtes à partir de la matrice d'information de Fisher (section 1.2.1.4). C'est pourquoi l'objectif du deuxième travail de recherche était
d'adapter les algorithmes multiplicatifs (jusqu'ici appliqués aux régressions linéaires ou non-linéaires simples) à l'optimisation de protocole pour modèles mixtes. Un autre objectif était de l'inclure dans la prochaine version de PFIM (section 1.2.1.4). Ces travaux sont présentés dans la section 2.2.

Nous avons vu en section 1.1.5. que la modélisation de la taille tumorale mesurée de manière longitudinale était de plus en plus développée, en particulier dans le développement d'anticancéreux utilisés en combinaison. Nous avons d'autre part vu que les simulations d'essais cliniques étaient adaptées à la comparaison de différents protocoles d'études (section 1.2.2.). Le troisième travail de ce doctorat porte sur la comparaison de protocoles en phase précoce de développement clinique en oncologie. Pour cette application, l'objectif était de savoir quels étaient les avantages et inconvénients de différents protocoles dans le cadre d'une étude de combinaison d'anticancéreux, avec comme critère de décision la décroissance tumorale précoce. Par simulation d'essais cliniques, des protocoles à un bras (comparaison à une référence historique), deux bras (essai clinique randomisé classique) et quatre bras (avec plusieurs doses possibles pour un médicament) ont donc été comparés, principalement en termes de puissance statistique (et d'erreur de type I) pour détecter un effet de la combinaison par rapport au traitement de référence. Ce travail est présenté dans la section 2.3.

L'autre partie de ce doctorat porte sur la phagothérapie, pour laquelle nous avons vu que la modélisation pouvait permettre de mieux comprendre l'interaction entre bactériophages et bactéries (section 1.3.2.). Des études ont été conduites *in vitro* et *in vivo* sur la souche bactérienne *Escherichia coli* 536 et le phage 536_P1 afin de comprendre les interactions entre bactéries et phages, ainsi que le lien entre le schéma d'administration des phages et leur efficacité. Notre travail de recherche avait pour objectif de construire un modèle mathématique permettant de caractériser l'évolution de la concentration de phage et de la charge bactérienne et d'en estimer les paramètres. Ce travail est présenté dans la troisième partie du manuscrit.

Une discussion de l'ensemble des travaux de cette thèse est proposée en quatrième partie.
# Deuxième partie : Evaluation et optimisation de protocoles dans les modèles mixtes

Dans cette partie, nous présentons les travaux autour de l'évaluation et l'optimisation de protocoles d'études longitudinales analysées par modèles mixtes. En section 2.1. est abordé le développement d'une méthode de planification de protocole robuste prenant en compte l'incertitude dans le choix du modèle pour des données binaires répétées. La section 2.2. traite de l'extension d'un algorithme multiplicatif à l'optimisation de protocoles pour les modèles mixtes, ainsi que de son intégration à la version de PFIM actuellement en développement. Enfin, en section 2.3., nous proposons une comparaison de protocoles par simulations appliquée à l'étude d'une combinaison d'anticancéreux.

## 2.1. Protocoles robustes tenant compte de l'incertitude du modèle pour réponse binaire répétée

### 2.1.1. Résumé

Nous avons proposé une méthode pour la planification de protocole dans le cadre d'études avec réponse binaire répétée et analysée par modèles mixtes. Cette approche est basée sur la FIM évaluée par échantillonnage de Monte Carlo et de Monte Carlo Hamiltonien, méthode existante dans la librairie R MIXFIM (Riviere and Mentré 2015; Riviere et al. 2016). En planification de protocoles, il est nécessaire de spécifier le modèle *a priori*. Or, si celui-ci est très différent du modèle finalement utilisé pour décrire les données, le protocole choisi pourra être sous-optimal. Nous avons donc proposé une méthode robuste tenant compte de l'incertitude du modèle en s'appuyant sur plusieurs modèles candidats.

Nous avons appliqué cette approche dans le cadre d'une étude avec deux groupes d'individus (placebo & traité) et supposé des modèles logistiques mixtes. Chaque protocole était composé de quatre temps de mesures, dont un au début de l'étude et un à 12 mois, c'est-à-dire la fin de l'étude. Nous avons cherché à optimiser, de façon combinatoire, l'allocation des deux autres temps de mesures parmi les valeurs possibles suivantes : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 et 11 mois. Premièrement, nous avons réalisé une planification de protocole locale (non robuste) pour chacun des 4 modèles candidats et selon 3 critères d'optimalité différents. Le critère de D-optimalité, visant à assurer une précision globale de l'ensemble des paramètres, a été utilisé. Le
critère Ds a aussi été employé dans le but d'optimiser la précision d'estimation d'un sous-
ensemble de paramètres, dans notre cas, l'effet de la covariable « traitement ». Cet effet
traitement étant détecté par test de Wald, l'optimisation de sa précision d'estimation permet
d'augmenter la puissance de ce test. Enfin, nous avons proposé l'utilisation du critère DDs qui
permet d'assurer un compromis entre la précision d'estimation de l'ensemble des paramètres et
la puissance du test de Wald pour détecter l'effet traitement. Deuxièmement, afin de prendre en
compte l'incertitude du modèle, nous avons alloué un poids à chacun des modèles candidats
compris entre 0 et 1, poids qui représente l'aptitude que l'on donne à un modèle à être celui qui
sera utilisé pour décrire les données. Dans notre exemple, nous avons considéré le même poids
(0,25) pour chacun des quatre modèles candidats. Nous avons utilisé ces poids dans le calcul
des différents critères robustes ou composés (CD pour Compound D, CDs et CDDs) afin
d'optimiser les temps de mesures. Nous avons aussi utilisé la FIM afin de prédire la puissance
statistique moyenne avec ces modèles candidats, pour détecter l'effet traitement.

Les temps de mesures obtenus selon les modèles et selon les critères sont différents. Les D-,
Ds- et DDs-efficacités ont été calculées avec les différents protocoles. L'approche robuste
permet d'éviter des protocoles pouvant être peu efficaces sous certains modèles. Par exemple,
le protocole CD-optimal conduit à une D-efficacité supérieure à 0.8 avec n'importe lequel des
modèles candidats mais la planification locale, couplée à un mauvais choix de modèle a priori,
peut conduire à une D-efficacité inférieure à 0.65. Nous avons prédit le nombre de sujets
nécessaires (NSN) pour obtenir une puissance de 90 % permettant de détecter l'effet du
traitement avec le protocole robuste (CDDs-optimal) ξ<sub>CDDs</sub> = (0,4,11,12), un protocole
planifié avec une approche non robuste (DDs-optimal pour le modèle le plus simple) ξ<sub>DDs,1</sub> =
(0,2,11,12) et un protocole non optimisé (avec des temps de mesures équi-répartis) ξ<sub>ES</sub> =
(0,4,8,12) : ils étaient respectivement de 274, 320 et 358.

Cinq cents jeux de données ont été simulés sous chaque modèle candidat et avec les trois
protocoles différents : robuste, DDs-optimal, et non optimisé. Les paramètres ont ensuite été
estimés pour chaque jeu de données et le bon modèle (correspondant au modèle de simulation),
avec l'algorithme SAEM de Monolix 2016R1. Les résultats obtenus avec les trois protocoles
ont été comparés en termes de biais, de précision d'estimation des paramètres, d'erreur de type
I et de puissance. Les simulations ont confirmé que l'utilisation d'un protocole robuste permet
une précision d'estimation acceptable des paramètres ainsi qu'une puissance de tests supérieure
à 80 % avec chaque modèle.
Nous avons aussi appliqué ces méthodes à une étude réelle de données binaires répétées provenant de 4162 sujets de la cohorte EPESE (Established Populations for Epidemiologic Studies of the Elderly (Landerman et al. 2011)) et analysées par modèle logistique mixte. Nous avons montré que les temps de mesures et la taille de l'étude auraient pu être choisis de manière plus efficace afin de détecter des effets de covariables, par exemple en optimisant les temps de mesures selon le critère DDs.

La stratégie proposée dans ce travail est pertinente pour la planification de protocole d'études longitudinales en cas d'incertitude du modèle. Ce travail a fait l'objet d'un article publié dans *Statistical Methods in Medical Research*.

### 2.1.2. Article 1 (publié)

Robust designs accounting for model uncertainty in longitudinal studies for binary outcome,
publié dans *Statistical Methods in Medical Research*, 2020

Jérémy Seurat, Thu Thuy Nguyen, France Mentré
# Robust designs accounting for model uncertainty in longitudinal studies with binary outcomes

Jérémy Seurat, ◦ Thu Thuy Nguyen and France Mentré

Statistical Methods in Medical Research
2020, Vol. 29(3) 934-952

© The Author(s) 2019

Article reuse guidelines:

sagepub.com/journals-permissions

DOI: 10.1177/0962280219850588

journals.sagepub.com/home/smm

<figure><img src="image_5.png" /></figure>

SAGE

## Abstract

To optimize designs for longitudinal studies analyzed by mixed-effect models with binary outcomes, the Fisher information matrix can be used. Optimal design approaches, however, require a priori knowledge of the model. We aim to propose, for the first time, a robust design approach accounting for model uncertainty in longitudinal trials with two treatment groups, assuming mixed-effect logistic models. To optimize designs given one model, we compute several optimality criteria based on Fisher information matrix evaluated by the new approach based on Monte-Carlo/Hamiltonian Monte-Carlo. We propose to use the $DD_s$-optimality criterion, as it ensures a compromise between the precision of estimation of the parameters, and hence the Wald test power, and the overall precision of parameter estimation. To account for model uncertainty, we assume candidate models with their respective weights. We compute robust design across these models using compound $DD_s$-optimality. Using the Fisher information matrix, we propose to predict the average power over these models. Evaluating this approach by clinical trial simulations, we show that the robust design is efficient across all models, allowing one to achieve good power of test. The proposed design strategy is a new and relevant approach to design longitudinal studies with binary outcomes, accounting for model uncertainty.

## Keywords

Optimal design, compound optimality, Fisher information matrix, longitudinal binary data, nonlinear mixed effect models

## I Introduction

Nonlinear mixed effect models (NLMEMs) are frequently used in model-based drug development to analyze pharmacokinetic/pharmacodynamic data.¹ The use of NLMEMs (i.e. the population approach) allows the estimation of mean parameters, inter and intra-subject variabilities as well as covariate effects and is appropriate for exploiting the richness of longitudinal data.²⁻⁴ NLMEMs can compensate for the lack of individual information by borrowing the strength from the data in the whole population and therefore allow for precise parameter estimation even with sparse designs, where few samples are collected from each subject.

NLMEMs are increasingly used for the analysis of longitudinal clinical studies, for both continuous and discrete data (binary, count or time-to-event). Binary outcomes are frequently encountered to characterize the clinical response in different therapeutic areas such as infectiology (virological success, bacterial carrying), and can be used for responder/non-responder analysis as well as for toxicity analysis. For instance, the carriage of *Streptococcus pneumoniae* was studied in children between two months and two years, by collecting nasopharyngeal swabs and aspirates over 10 visits to detect the impact of determinants as age or health status on the proportion of positive bacterial samples.⁵ In another context, Mestre et al.⁶ recently proposed a binary outcome (corresponding to a minimal change in UDysRS Part III Impairment scale from baseline) to characterize the remission of dyskinesia in Parkinson's disease. The analysis of the longitudinal data in these kinds of outcomes required discrete response generalized linear mixed effect models (GLMMs) or discrete response NLMEMs and adequate estimation methods.

IAME, UMR 1137, INSERM, Université Paris Diderot, Paris, France

**Corresponding author:**

Jérémy Seurat, IAME, UMR 1137, INSERM, University Paris Diderot Sorbonne Paris Cité, 16, rue Henri Huchard, 75018 Paris, France.

Email: jeremy.seurat@inserm.fr
Before modeling, it is crucial to choose an appropriate design in order to obtain good precision of parameter estimates. Indeed, the informativeness of a dataset for parameter estimation depends on the design choice. A design in NLMEMs, called a population design, is composed of the number of elementary designs, the specification of each elementary design, and the associated number of subjects. In these settings, the term elementary design is used to describe a group of subjects with identical design characteristics such as the number of samples per subject and the allocation of informative times and doses. To evaluate and optimize designs, two approaches have been proposed. The first approach, based on clinical trial simulation (CTS), is very time-consuming and is therefore limited in term of designs that can be evaluated. Alternatively, the expected Fisher information matrix (FIM) can be used,⁷⁻⁹ as its inverse is defined as the lower bound of the variance-covariance matrix of any unbiased estimated parameters, according to the Cramér-Rao inequality. However, the FIM has no analytical form in NLMEMs and its computation, which requires multiple integrations, can be challenging. A first method to evaluate the FIM in NLMEMs, based on first order (FO) linearization of the model around the expectation of random effects, was proposed¹⁰ and implemented in several software programs.¹¹ Although efficient in general, FO presents limitations when used with complex models, large variability¹²,¹³ and when using longitudinal models for discrete data.¹³ Alternative to compute the FIM in continuous NLMEMs, without linearization, has been proposed using Monte Carlo (MC) integrations.¹⁴ There are also a large body of work on the FIM for GLMMs with random intercepts, using marginal quasi-likelihood, penalized quasi-likelihood, complete enumeration-based or MC methods.¹⁵,¹⁶ More recently, new methods have been developed to compute the FIM for both continuous and discrete NLMEMs using MC methods combined with adaptive Gaussian quadrature (AGQ)¹⁷ or Hamiltonian Monte Carlo (HMC).¹³ The latter approach (MC/HMC) was implemented in the R package *MIXFIM*,¹⁸ which uses functions written in the probabilistic language Stan,¹⁹ which was developed for Bayesian inference. This method, by efficiently drawing HMC samples and calculating partial derivatives of the log-likelihood, has been shown to be more suitable than MC/AGQ in complex models with a large number of random effects.¹³

From the expression of the FIM evaluated by these approaches, different optimality criteria can then be computed to optimize designs. The widely used D-optimality criterion consists in maximizing the determinant of the FIM, i.e. to optimize the precision of estimation for all model parameters. As an alternative to the D-optimality, the D<sub>S</sub>-optimality<sup>20</sup> can accommodate situations in which only a subset of the model parameters is of interest. This can be particularly useful to optimize the precision of a “treatment effect” parameter on the longitudinal evolution of biomarkers in clinical trials for example, which is directly linked to the power of the Wald test to detect this effect.²¹ However, D<sub>S</sub>-optimal designs can lead to problem in estimation of all the parameters of the model because of lack of experimental identifiability.²² To find a balance between optimizing the precision of the parameters to be tested and the precision all the model parameters, a mixture of D- and D<sub>S</sub>-optimality, the DD<sub>S</sub>-optimality, can be used.²⁰ However, these methods require a priori knowledge of the model and its parameter values, which can usually be obtained from previous studies but lead to designs that are locally optimal. If the final model is very different from the a priori assumed model, the design optimized with the wrongly assumed model might not be informative enough to precisely estimate parameters of the final model. Therefore, there is need to develop and evaluate new methods to optimize robust designs accounting for model uncertainty. In order to propose optimal designs across a set of candidate models, we make use of both the theory of compound optimal design<sup>7</sup> and the principle of model averaging.²³ The compound optimal design theory was previously used to propose informative designs for both estimation and model discrimination,<sup>24</sup> and recently to find a common design for cocktail of drugs.²⁵ Meanwhile, model averaging is a promising alternative approach to model selection when modeling data, which uses model weights (computed from AIC or BIC) to calculate a weighted average of the predictions.<sup>26-28</sup>

Combining these approaches, our main objective is to propose and evaluate a new methodological strategy for
robust designs in longitudinal studies with discrete outcomes taking into account model uncertainty, based on the
expected FIM in NLMEMs. We choose to use the MC/HMC method to compute the FIM without any
linearization. From expression of the FIM, we predict the power of the Wald test to detect a treatment effect
included in the model and to calculate the number of subjects needed for a required power. Our FIM-based
strategy aims to find designs that are both informative to achieve good average power over several models and to
precisely estimate the whole set of population parameters. We illustrate these methods through an example of
designing a trial with repeated binary outcomes, inspired from a previous work within the team,<sup>29</sup> considering
several candidate models which include a treatment effect on the longitudinal response. The relevance of
these approaches is evaluated by CTS in terms of bias and imprecision of parameter estimates and power of
the Wald test.
In Section 2, we detail the notations and present different optimality criteria for a given model or for taking into account model uncertainty. In Section 3, we show an application of these methods to design a longitudinal trial with binary outcomes. In Section 4, we present a simulation study to evaluate the relevance of the proposed design strategy. In Section 5, we show an application to a real study. Finally, we discuss the results and perspectives of this work in Section 6.

# 2 Methods

## 2.1 Basic concepts and notation

### 2.1.1 Population design

The elementary design $\xi_i$ for the individual $i$ ($i = 1, ..., N$) is defined by the number $n_i$ of samples and their allocation in time ($t_{i1}, ..., t_{in}$). The population design $\Xi = \{N, (\xi_1, ..., \xi_N)\}$ is defined by the number of individuals $N$ and the set of elementary designs to be performed in each individual $(\xi_1, ..., \xi_N)$, with a total number of observations $\sum_{i=1}^{N} n_i$.

### 2.1.2 NLMEMs for binary outcomes

This work considers a set of candidate NLMEMs $m = 1, ..., M$ for binary data. Considering a given logistic model $m$, the logit of the conditional probability for observation $y_{ij}$ of individual $i$ at sample $j$ ($j = 1, ..., n_i$) is written

$$ \operatorname{logit}(p(y_{ij} = 1|b_i)) = f_m(t_{ij}, g(\mu_m, b_i, z_i, \beta_m)) \quad (1) $$

In equation (1), $f_m$ is the function describing the logit-probability of $y_{ij}$ given individual random effects $b_i$, for a time $t_{ij}$ and a vector of subject-specific parameters modelled through $g$. The function $g$ depends on the vector of fixed effects $\mu_m$, $b_i$ as well as the vector of covariates $z_i$ and the vector of covariate effects $\beta_m$. The random effects are assumed to follow a normal distribution with mean zero and variance-covariance matrix $\Omega_m$ which accounts for the inter-individual variability, i.e. $b_i \sim N(0, \Omega_m)$. We denote by $\psi_m$ the vector of population parameters, i.e. $\psi_m = (\mu_m^T, \beta_m^T, \Omega_{m,u}^T)^T$, where $\Omega_{m,u}$ is a vector containing all unique elements of $\Omega_m$. Let $P_m$ denote the length of the vector $\psi_m$. Observations are usually assumed to be independent conditionally upon the random effects, i.e. the joint conditional probability for the vector of observations for individual $i$ ($y_i = (y_{i1}^T, ..., y_{in_i}^T)^T$) is

$$ p(y_i|b_i, \psi_m) = \prod_{j=1}^{n_i} \frac{\exp(f_m(t_{ij}, g(\mu_m, b_i, z_i, \beta_m)))}{1 + \exp(f_m(t_{ij}, g(\mu_m, b_i, z_i, \beta_m)))} $$

### 2.1.3 FIM

For a model $m$ and its parameter values $\psi_m$, if the $N$ subjects are independent, the population FIM, $\mathcal{M}(\psi_m, \Xi)$, is the sum of $N$ elementary FIMs, $\mathcal{M}(\psi_m, \xi_i)$, i.e. $\mathcal{M}(\psi_m, \Xi) = \sum_{i=1}^{N} \mathcal{M}(\psi_m, \xi_i)$. In this work, we assume the same elementary design for all individuals, i.e. $\xi_i = \xi$ for $i = 1, ..., N$, then $\Xi = \{N, \xi\}$ and $\mathcal{M}(\psi_m, \Xi) = N\mathcal{M}(\psi_m, \xi)$. The elementary FIM for $\psi_m$ can be expressed as

$$ \mathcal{M}(\psi_m, \xi) = E_y \left( \frac{\partial \log(L(y|\psi_m))}{\partial \psi_m} \frac{\partial \log(L(y|\psi_m))^T}{\partial \psi_m} \right) \quad (2) $$

where the likelihood $L$ of the observations vector $y$ of an individual (subscript $i$ is omitted for simplicity) is given by

$$ L(y|\psi_m) = \int_b p(y|b, \psi_m) p(b|\psi_m) db \quad (3) $$

with $p(y|b, \psi_m)$ the probability density function (p.d.f.) of $y$ given the random effects $b$, and $p(b|\psi_m)$ the p.d.f. of $b$. To evaluate the FIM in equation (2), one needs to compute two integrations: one over the observations $y$ and one in equation (3) over the distribution of random effects. We evaluate the former by MC integrations and the latter by HMC, as proposed in Riviere et al.,¹³ using the package *MIXFIM*¹⁸ in R v3.2.1. *MIXFIM* calls functions of the R package *rstan* written in the probabilistic language Stan¹⁹ to draw HMC samples and calculate partial derivatives of the log-likelihood.
## 2.2 Optimality criteria

The methods we propose to find robust designs are presented in this section. Different criteria to find an optimal design are presented, according to different purposes. For a given model, based on the determinant of the FIM, the D-optimality is used to obtain informative designs for a precise estimation of all parameters of a model, the D<sub>S</sub>-optimality accommodates situations in which only a subset of parameters is of interest, and the compound DD<sub>S</sub>-optimality offers a compromise between the D- and the D<sub>S</sub>-criteria. Averaging over several models and weighting each one, the compound CD-, CD<sub>S</sub>- and CDD<sub>S</sub>-optimalities accounts for model uncertainty.

### 2.2.1 Optimality criteria for a given model

To optimize the estimation precision for the whole set of population parameters $\psi_m$, we maximize the widely used<sup>7</sup> D-optimality criterion $\Phi_{D,m}(\Xi)$ (equation (4), Table 1).

Let $\psi_{S,m}$, a subset of $\psi_m$, be the vector of parameters of interest of length $S_m$. To optimize the estimation precision for the $S_m$ parameters, we use the D<sub>S</sub>-optimality $\Phi_{D_s,m}(\Xi)$ (equation (5), Table 1).<sup>7</sup> This criterion is particularly useful to find a design which maximizes the predicted power of the Wald test or minimizes the number of subjects needed to detect a significant treatment effect $\beta_m$, by focusing on the precision of estimation for $\beta_m$.

Let $\psi_{T,m}$ be the remaining parameters of $\psi_m$ beside $\psi_{S,m}$. To find a compromise between precision of $\psi_{S,m}$ and $\psi_{T,m}$, the compound DD<sub>S</sub>-optimality criterion,<sup>7,20</sup> $\Phi_{DD_s,m}(\Xi, \alpha_m)$ is used (equation (6), Table 1). $\alpha_m$ is a term between 0 and 1, which quantifies the interest in the estimation precision of $\psi_{S,m}$. When $\alpha_m = 1/P_m$, the DD<sub>S</sub>-criterion coincides with the D-criterion; thus, all the parameters of the model are of equal interest. When $\alpha_m = 1$, the compound DD<sub>S</sub>-criterion coincides with the D<sub>S</sub>-criterion. When $\alpha_m = 0$, the precision of the parameters $\psi_{T,m}$ is optimized. Maximizing the DD<sub>S</sub>-criterion with $1/P_m < \alpha_m < 1$ provides a compromise between the D- and the D<sub>S</sub>-criteria.

For any given optimality criterion $X$, the relative efficiency $E_{X,m}(\Xi)$ of a design $\Xi$, with respect to the $X$-optimal design $\Xi_{X,m}$ for model $m$, is computed as

$$E_{X,m}(\Xi) = \frac{\Phi_{X,m}(\Xi)}{\Phi_{X,m}(\Xi_{X,m})}$$

### 2.2.2 Robust optimality criteria accounting for model uncertainty

We propose to account for uncertainty in the model choice by computing optimality criteria based on the theory of compound optimality<sup>7</sup> and of model averaging.<sup>23</sup> We assume that a set of M candidate models is available, with each model associated to a weight $w_m$ ($\sum_{m=1}^M w_m = 1$) quantifying the balance between the M models. To find a common optimal design for the M models, according to the X-optimality, we maximize the weighted product of

<table><caption>**Table 1.** Various optimality criteria used depending on parameters of interest and accounting or not for model uncertainty.</caption><thead><tr><th>Parameters of interest</th><th>Given model m</th><th>Set of candidate models<br>m = 1, ..., M</th></tr></thead><tbody><tr><td>All the parameters ψ<sub>m</sub></td><td>D-optimality<br>Φ<sub>D,m</sub>(Ξ) = Det(M(ψ<sub>m</sub>, Ξ))<sup>1/P<sub>m</sub></sup> (4)</td><td>CD-optimality<sup>7,24,25</sup><br>Φ<sub>CD</sub>(Ξ) = Π<sup>M</sup><sub>m=1</sub> (Φ<sub>D,m</sub>(Ξ))<sup>w<sub>m</sub></sup> (7)</td></tr><tr><td>Subset of parameters ψ<sub>S,m</sub></td><td>D<sub>S</sub>-optimality<br>Φ<sub>D<sub>S</sub>,m</sub>(Ξ) = (Det(M(ψ<sub>T,m</sub>, Ξ)) / Det(M(ψ<sub>T,m</sub>, Ξ)))<sup>1/S<sub>m</sub></sup> (5)</td><td>CD<sub>S</sub>-optimality<br>Φ<sub>CD<sub>S</sub></sub>(Ξ) = Π<sup>M</sup><sub>m=1</sub> (Φ<sub>D<sub>S</sub>,m</sub>(Ξ))<sup>w<sub>m</sub></sup> (8)</td></tr><tr><td>Compromise between<br>ψ<sub>S,m</sub> and ψ<sub>T,m</sub></td><td>DD<sub>S</sub>-optimality<br>Φ<sub>DD<sub>S</sub>,m</sub>(Ξ, α<sub>m</sub>) = (Det(M(ψ<sub>T,m</sub>, Ξ)))<sup>1-α<sub>m</sub>/P<sub>m</sub>-α<sub>m</sub></sup><br>× (Det(M(ψ<sub>T,m</sub>, Ξ)))<sup>α<sub>m</sub>/S<sub>m</sub></sup> (6)</td><td>CDD<sub>S</sub>-optimality<br>Φ<sub>CDD<sub>S</sub></sub>(Ξ) = Π<sup>M</sup><sub>m=1</sub> (Φ<sub>DD<sub>S</sub>,m</sub>(Ξ))<sup>w<sub>m</sub></sup> (9)</td></tr></tbody></table>

Note: $\psi_m$ is the vector of population parameters (of size $P_m$) for a model $m$: $\psi_{S,m}$ the subset of parameters of interest (of size $S_m$) and $\psi_{T,m}$ which are not of interest. $M(\psi_m, \Xi)$ is the matrix containing the Fisher information for all the parameters $\psi_m$ and $M(\psi_{T,m}, \Xi)$ is obtained by truncation of $M(\psi_m, \Xi)$ and keeping only the rows and columns corresponding to $\psi_{T,m}$. $\alpha_m$ ($0 \le \alpha_m \le 1$) is a term which quantifies the balance between the D- and the D<sub>S</sub>-optimality criteria and expresses the interest in the precision of estimation for $\psi_{S,m}$. $w_m$ is the weight for the model $m$ ($\sum_{m=1}^M w_m = 1$).
efficiencies $\prod_{m=1}^{M} E_{X,m}^{w_m}$, with X being the D-, Ds- or DDs-optimality criterion. Maximizing the weighted product of efficiencies is equivalent to maximizing the weighted product of the criteria, hence the definition of compound optimality criteria,⁷ $\Phi_{CD}(\Xi)$, $\Phi_{CDs}(\Xi)$ or $\Phi_{CDDs}(\Xi, \alpha_m)$ as provided (equations (7) to (9), Table 1). Of note, it is also equivalent to minimizing the weighted sum of minus-log-criteria, that is a convex function to which the General Equivalence Theorem can be applied.³⁰

# 3 Application to design a longitudinal trial with binary response

In this section, we present how the methods introduced in Section 2 are applied to design a longitudinal study over one year, with two balanced treatment groups and binary response. Candidate NLMEMs including a treatment effect and describing the response probability over time are considered, inspired by Lestini et al.²⁹ and Ogungbenro and Aarons.³¹ We aim to determine the optimal location of a limited number of measurement times chosen among a discrete design composed of (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) in month unit, according to different optimality criteria. We investigate the predicted efficiency of the various optimal designs obtained as well as the power of the Wald test to detect a treatment effect over all the candidate models.

## 3.1 Models and parameters

The studied models describe binary responses ($y_{ij}=1$, responder or 0, non responder) recorded over time ($t_{ij}=0$ to 12 months) in two different treatment groups. The logit of response probability at time $t_{ij}$ is given by equation (1).

We denote $f_m(t_{ij}) = f_m(t_{ij}, g(\mu_m, b_i, z_i, \beta_m))$ and omit index $ij$ for simplicity in the followings. We consider four candidate models $f_m$ describing the evolution of the logit-probability of the response over time (Figure 1), which are the linear model (M1), loglinear model (M2), quadratic model (M3) and exponential model (M4).

$$M1: f_1(t) = \theta_1 + \theta_2(1 + \beta \times 1_T)t$$

$$M2: f_2(t) = \theta_1 + \theta_2(1 + \beta \times 1_T)\log(t + 1)$$

$$M3: f_3(t) = \theta_1 + \theta_2(1 + \beta \times 1_T)t^2$$

$$M4: f_4(t) = \theta_1 + \theta_2(1 + \beta \times 1_T)(\exp(\theta_3 t) - 1)$$

$1_T$ is a treatment group indicator variable (with $1_T = 0$ if control group, $1_T = 1$ if treated group). We denote $\beta$ the effect size of the treatment on parameter $\theta_2$. Model M1 and its parameters are inspired by Lestini et al.²⁹ Values of $\psi_1$ are given in Table 2, $\theta_1$ and $\theta_2$ follow a normal distribution:

$\theta_p = \mu_p + b_p$ where $b_p \sim N(0, \omega_p^2)$ for $p = \{1, 2\}$. Models M2 to M4 are alternative models with parameters normally distributed as in M1. Values of their fixed parameters are chosen to give the same mean value of the logit of the response probability as for M1 at the beginning ($t=0$ month) and the end ($t=12$ months) of the study in the two treatment groups. Similar treatment effect value and a standard deviation $\omega_2$ equal to 188% of $\mu_2$ for are assumed for all models. Assumed values of population parameters $\psi_m$ for each model $m$ ($m = 1, \dots, 4$) are given in Table 2.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Plots of the logit of the response probability ($p$) over 12 months in each treatment group for the four candidate models: M1 linear (in orange), M2 loglinear (in blue), M3 quadratic (in green) and M4 exponential (in purple), using the fixed effect values indicated in Table 2. Solid lines represent the control group and dashed lines the treated group.</figcaption></figure>
<table><caption>**Table 2.** Population parameter values $\psi_m$ for each candidate model $m$ ($m = 1, \dots, 4$): M1 is the linear model, M2 the loglinear model, M3 the quadratic model and M4 the exponential model.</caption><thead><tr><th></th><th>M1</th><th>M2</th><th>M3</th><th>M4</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>-2</td><td>-2</td><td>-2</td><td>-2</td></tr><tr><td>μ<sub>2</sub></td><td>0.09</td><td>0.42</td><td>7.50×10<sup>-3</sup></td><td>2.01×10<sup>-2</sup></td></tr><tr><td>μ<sub>3</sub></td><td>-</td><td>-</td><td>-</td><td>0.33</td></tr><tr><td>β</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>ω<sub>1</sub></td><td>0.70</td><td>0.70</td><td>0.70</td><td>0.70</td></tr><tr><td>ω<sub>2</sub></td><td>0.17</td><td>0.79</td><td>1.41×10<sup>-2</sup></td><td>3.79×10<sup>-2</sup></td></tr></tbody></table>

Note: $\mu$ are the fixed effects, $\omega$ the standard deviation of random effects and $\beta$ the effect of the treatment covariate.

## 3.2 Design optimization and power evaluation

We aim to propose an informative design with four measurement times identical in all subjects, denoted $\xi = \{t_1, t_2, t_3, t_4\}$, first, for each candidate model separately, then, by averaging over the four candidate models, given the parameters values $\psi_m$ given in Table 2.

### 3.2.1 Choice of MC-HMC samples for FIM evaluation

To perform design optimization for the presented models, we need to compute different optimality criteria, based on the FIM evaluated by MC/HMC. Since this is a stochastic approach, we can quantify uncertainty in the computation of the determinant of the FIM, according to the MC and HMC samples. For instance, the 95% uncertainty intervals of the D-criterion can be computed from the 2.5th and 97.5th percentiles of the distributions generated using non-parametric bootstrap. According to Riviere et al.,¹³ the FIM evaluation method implemented in *MIXFIM* performs well with 5000 Monte Carlo samples and 200 HMC samples, for a longitudinal logistic model for binary response (similar to our model M1) and a rich design. To ensure that these settings are appropriate for all candidate models even for sparser design (four measurement times), the convergence of the D-criterion is studied with respect to the MC samples increasing from 50 to 10,000 (see Supplementary Material, Figure S1), with a non-optimized equispaced design $\xi_{ES} = (0, 4, 8, 12)$. Two different configurations were studied: 1000 HMC samples and 1000 burn-in vs. 200 HMC and 500 burn-in. We found that the convergence plots of D-criterion shows a similar trend between the two configurations. Moreover, 5000 MC samples seem sufficient for a good convergence of the D-criterion. We therefore set the number of MC samples to 5000 and the number of HMC samples to 200 with 500 burn-in in *MIXFIM* to evaluate the FIM and to calculate the different optimality criteria with all the designs and the models of the study.

To ensure that these settings are satisfying, convergence of the D-criterion is also verified with an optimized design (Supplementary Material, Figure S2).

### 3.2.2 Design optimization

As in Lestini et al.,²⁹ we assumed that $t_1$ and $t_4$ are fixed to 0 and 12 months (i.e. the end of the study), only two times $t_2$ and $t_3$ are optimized among the following set: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 months. No repetition is considered, i.e. four measurement times must be different. Combinatorial optimization is performed for different optimality criteria. This corresponds to $\binom{11}{2} = 55$ possible elementary designs to be evaluated, and the optimal one is chosen according to each considered criterion. We considered $N = 100$ subjects equally distributed in the two treatment groups for design optimization, although this number has no influence on optimal sampling time or loss of efficiency in design comparisons, but only on the specific value of relative standard errors when reported.

#### 3.2.2.1 D-, Ds- and DDs-optimal designs for each model.
Using the D-, Ds- or DDs-optimality, the optimal location of measurement times are obtained and compared between the four models and between these optimality criteria. The relative D-, Ds- or DDs-efficiencies of each optimal design are also computed assuming each of the four models to evaluate the impact of model misspecification on design performance.

The Ds- and the DDs-criteria are computed with a particular interest on the treatment effect ($\psi_{S,m} = \beta$). The preliminary step to compute compound DDs-optimal designs is to find an appropriate $\alpha_m$ value for each model.
For each model, we therefore compute the DD<sub>S</sub>-optimal designs Ξ<sub>DD<sub>S</sub>,m</sub> for α<sub>m</sub> values from 0 to 1 by increment of 0.05. Then, we compare the DD<sub>S</sub>-optimal allocations of measurement times and the D- and D<sub>S</sub>-efficiencies of these DD<sub>S</sub>-optimal designs over α<sub>m</sub>. We consider that for a DD<sub>S</sub>-design to be satisfactory, α<sub>m</sub> should maximize E<sub>D,m</sub> × E<sub>D<sub>S</sub>,m</sub>. Moreover, both the D- and D<sub>S</sub>-efficiencies should be above 0.8 (other threshold values are not investigated in this work). The final DD<sub>S</sub>-optimal design is computed with the retained α<sub>m</sub> value.

3.2.2.2. Robust optimal designs. To propose a robust design accounting for model uncertainty, different compound criteria (CD-, CD<sub>S</sub>- and CDD<sub>S</sub>-optimality) are evaluated for a combination of the four candidate models. The same weight w<sub>m</sub> is assigned for each model (total uncertainty) i.e. w<sub>m</sub> = 1/4.

3.2.3 Expected average power using FIM

We aim to evaluate the average power of the Wald test for a design Ξ = {N, ξ} to detect a significant treatment effect over all the candidate models varying the total number of subject from 50 to 450. We also computed the number of subjects to achieve an average power of 0.9.

For that, first, we derive for each model *m* the standard error (SE) of *β* from the expected FIM and we compute the power π<sub>m</sub>(Ξ) to detect the treatment effect in that model, following what has been done in Retout et al.<sup>21</sup> Then, the average power π<sub>average</sub>(Ξ) over *M* candidate models is evaluated as

$$
\pi_{\text{average}}(\Xi) = \sum_{m=1}^{M} w_m \times \pi_m(\Xi)
$$

where w_m is the weight associated to each model m used in the computation of the optimality criteria (here 1/M,
with M = 4).

3.3 Results

Figure 2 displays, for each model, the efficiencies with respect to each criterion, obtained with the different combinations of the second and third samples out of four measurement times.

3.3.1 D- and D<sub>S</sub>-optimal designs for each model

Different models lead to different optimal designs. Table 3 reports the D- and D<sub>S</sub>-efficiencies, respectively, of the D-optimal design Ξ<sub>D,m</sub> and D<sub>S</sub>-optimal design Ξ<sub>D<sub>S</sub>,m</sub>, when the model used for data analysis is model *m* or another one. For instance, the design (0, 4, 5, 12) is D-optimal for M3 but underperforms for M4 with a low D-efficiency of only 0.646, and the design (0, 1, 11, 12) is D<sub>S</sub>-optimal for M2 but conduces to a low D<sub>S</sub>-efficiency of only 0.389 for M3.

In this work, the D<sub>S</sub>-optimality focuses on the estimation precision of the treatment effect β. We notice that using this criterion could lead to different optimal measurement times than those obtained with the D-criterion, as in the case of M2 or M4 (Figure 2). For M2, the allocations of informative measurement times leading to a D- and D<sub>S</sub>-efficiency above 0.8 are very similar, and the respective optimal designs are ξ<sub>D,2</sub> = (0, 1, 8, 12) vs. ξ<sub>D<sub>S</sub>,2</sub> = (0, 1, 11, 12). However, we notice quite different efficient designs between D- and D<sub>S</sub>-criterion for M4 (ξ<sub>D,4</sub> = (0, 6, 11, 12) vs. ξ<sub>D<sub>S</sub>,4</sub> = (0, 10, 11, 12)). This could be explained by the shape of the model curve (Figure 1, the second measurement time should be more tardive to better distinguish the two treatment groups) and by the higher number of parameters of this model compared to others (the importance of estimating β is thus relatively reduced when using D-optimality).

3.3.2 DD<sub>S</sub>-optimal design for each model

For each model, Figure 3 shows allocated measurement times as well as D- and D<sub>S</sub>-efficiencies of the DD<sub>S</sub>-optimal designs over α<sub>m</sub> possible values varying between 0 and 1. The third measurement time tends to 11 with increasing α<sub>m</sub> for all models. These results emphasize the necessity to get tardive measurements to better distinguish the responses between the control and the treated group, and thus to better estimate the treatment effect β. The product of D- and D<sub>S</sub>-efficiencies is, respectively, maximal when α<sub>1</sub> ∈ [0.2, 1], α<sub>2</sub> ∈ [0.3, 1], α<sub>3</sub> ∈ [0.25, 1] and α<sub>4</sub> ∈ [0.6, 1]. For these values of α<sub>m</sub>, D- and D<sub>S</sub>-efficiencies are over 0.8 excepted for M4. With M4, D- and D<sub>S</sub>-efficiencies are over 0.8 when 0.5 ≤ α<sub>4</sub> ≤ 0.55. In order to conserve the same range of importance for estimating the precision of parameters between the four models, we choose 0.5 as the same α<sub>m</sub> value with
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Heatmaps of design efficiencies, for each candidate model (M1 in red, M2 in blue, M3 in green, M4 in purple, by column) according to different considered optimality criterion (D-, D<sub>S</sub>- or DD<sub>S</sub>-optimality, by row). The fifth column (in grey) represents the robust criteria CD-, CD<sub>S</sub>- or CDD<sub>S</sub>-efficiencies. For each heatmap, the second measurement time is in abscissa and the third in ordinate. The filled black square shows the optimal design (i.e. efficiency = 1), which is also specified below each heatmap.</figcaption></figure>

each model. Thus, DD<sub>S</sub>-optimal design for each model is computed with $\alpha_m = 0.5$, with this value satisfying our two conditions.

We notice that the DD<sub>S</sub>-criterion could lead to different optimal measurement times than those obtained with the D- or D<sub>S</sub>-criterion (Figure 2), as in the case of M4 ($\xi_{D,4} = (0, 6, 11, 12)$, $\xi_{D_{S,4}} = (0, 10, 11, 12)$, $\xi_{DD_{S,4}} = (0, 9, 11, 12)$). These three designs lead to different predicted relative standard errors (RSE) of $\beta$, which for 100 patients, are 60%, 50%, 52%, respectively, illustrating the impact of using D<sub>S</sub>- or DD<sub>S</sub>-criterion on power to detect a treatment effect. For this model, we also studied the influence of the number of sampling times. We find that the DD<sub>S</sub>-optimal design with five samples is $\xi_{DD_S} = (0, 6, 10, 11, 12)$, leading to an RSE of $\beta$ of 50%, whereas the design with all 13 points lead to an RSE of 40%. For M1, M2, and M3, the DD<sub>S</sub>-efficient designs coincide with the D<sub>S</sub>-efficient designs. In the same way as for previous criteria (D- and D<sub>S</sub>-efficiencies), different models lead to different DD<sub>S</sub>-efficient designs. However, all the DD<sub>S</sub>-optimal designs include the 11th month, which again emphasizes the importance of tardive measurement times ($t_3$ and $t_4$) in estimating the treatment effect. Table 3 reports the DD<sub>S</sub>-efficiencies of the DD<sub>S</sub>-optimal design $\Xi_{DD_{S,m}}$ when the true model is *m* or another one. For instance, the design (0, 1, 11, 12) is DD<sub>S</sub>-optimal for M2 but underperforms for M3 with a low DD<sub>S</sub>-efficiency of 0.611. Therefore, this design provides poor precision on parameters of M3, especially $\beta$. With 100 individuals, the predicted RSE of $\beta$ is 65% with this design vs. 41% with the DD<sub>S</sub>-optimal for M3 ($\xi_{DD_{S,3}} = (0, 4, 11, 12)$).

### 3.3.3 Robust optimal designs

The CD-, CD<sub>S</sub>- and CDD<sub>S</sub>-efficiencies of every possible designs as well as the corresponding optimal design are reported in the last column of Figure 2. First, we can note that efficient robust designs are nearly the same using CD<sub>S</sub>- or CDD<sub>S</sub>-optimality.

Assuming model uncertainty leads to a robust CD-optimal design ($\xi_{CD} = (0, 5, 11, 12)$) which is different from the four D-optimal designs for each model. The most CD-efficient designs are closer to the D-efficient designs for the quadratic model (M3) than for other models. Table 3 presents D-efficiencies for each model. The robust CD-optimal design conducts to a D-efficiency above 0.8 regardless of the selected model. A model misspecification can lead to loss of efficiencies of up to 0.354.
<table><caption>**Table 3.** D-, Ds- and DDs-efficiencies for different optimal designs, accounting or not for model uncertainty.</caption><thead><tr><th colspan="5">Design Ξ<br>Ξ = {N = 100, ξ}</th></tr><tr><th></th><th>E<sub>D, 1</sub></th><th>E<sub>D, 2</sub></th><th>E<sub>D, 3</sub></th><th>E<sub>D, 4</sub></th></tr></thead><tbody><tr><td>ξ<sub>ES</sub> = (0, 4, 8, 12)</td><td>0.908</td><td>0.926</td><td>0.975</td><td>0.869</td></tr><tr><td>ξ<sub>D, 1</sub> = (0, 2, 11, 12)</td><td>1</td><td>0.898</td><td>0.812</td><td>0.706</td></tr><tr><td>ξ<sub>D, 2</sub> = (0, 1, 8, 12)</td><td>0.932</td><td>1</td><td>0.875</td><td>0.787</td></tr><tr><td>ξ<sub>D, 3</sub> = (0, 4, 5, 12)</td><td>0.916</td><td>0.839</td><td>1</td><td>0.646</td></tr><tr><td>ξ<sub>D, 4</sub> = (0, 6, 11, 12)</td><td>0.829</td><td>0.802</td><td>0.960</td><td>1</td></tr><tr><td>ξ<sub>CD</sub> = (0, 5, 11, 12)</td><td>0.860</td><td>0.805</td><td>0.994</td><td>0.958</td></tr><tr><td>ξ<sub>CDs</sub> = (0, 4, 11, 12)</td><td>0.910</td><td>0.828</td><td>0.977</td><td>0.882</td></tr><tr><td>ξ<sub>CDDs</sub> = (0, 4, 11, 12)</td><td>0.910</td><td>0.828</td><td>0.977</td><td>0.882</td></tr><tr><td></td><td>E<sub>Ds, 1</sub></td><td>E<sub>Ds, 2</sub></td><td>E<sub>Ds, 3</sub></td><td>E<sub>Ds, 4</sub></td></tr><tr><td>ξ<sub>ES</sub> = (0, 4, 8, 12)</td><td>0.841</td><td>0.646</td><td>0.835</td><td>0.393</td></tr><tr><td>ξ<sub>Ds, 1</sub> = (0, 2, 11, 12)</td><td>1</td><td>0.806</td><td>0.434</td><td>0.619</td></tr><tr><td>ξ<sub>Ds, 2</sub> = (0, 1, 11, 12)</td><td>0.751</td><td>1</td><td>0.389</td><td>0.627</td></tr><tr><td>ξ<sub>Ds, 3</sub> = (0, 5, 11, 12)</td><td>0.774</td><td>0.600</td><td>1</td><td>0.667</td></tr><tr><td>ξ<sub>Ds, 4</sub> = (0, 10, 11, 12)</td><td>0.529</td><td>0.519</td><td>0.611</td><td>1</td></tr><tr><td>ξ<sub>CD</sub> = (0, 5, 11, 12)</td><td>0.774</td><td>0.600</td><td>1</td><td>0.667</td></tr><tr><td>ξ<sub>CDs</sub> = (0, 4, 11, 12)</td><td>0.865</td><td>0.637</td><td>0.906</td><td>0.655</td></tr><tr><td>ξ<sub>CDDs</sub> = (0, 4, 11, 12)</td><td>0.865</td><td>0.637</td><td>0.906</td><td>0.655</td></tr><tr><td></td><td>E<sub>DDs, 1</sub></td><td>E<sub>DDs, 2</sub></td><td>E<sub>DDs, 3</sub></td><td>E<sub>DDs, 4</sub></td></tr><tr><td>ξ<sub>ES</sub> = (0, 4, 8, 12)</td><td>0.882</td><td>0.754</td><td>0.924</td><td>0.723</td></tr><tr><td>ξ<sub>DDs, 1</sub> = (0, 2, 11, 12)</td><td>1</td><td>0.863</td><td>0.644</td><td>0.765</td></tr><tr><td>ξ<sub>DDs, 2</sub> = (0, 1, 11, 12)</td><td>0.859</td><td>1</td><td>0.611</td><td>0.755</td></tr><tr><td>ξ<sub>DDs, 3</sub> = (0, 5, 11, 12)</td><td>0.827</td><td>0.722</td><td>1</td><td>0.947</td></tr><tr><td>ξ<sub>DDs, 4</sub> = (0, 9, 11, 12)</td><td>0.683</td><td>0.672</td><td>0.761</td><td>1</td></tr><tr><td>ξ<sub>CD</sub> = (0, 5, 11, 12)</td><td>0.827</td><td>0.722</td><td>1</td><td>0.947</td></tr><tr><td>ξ<sub>CDs</sub> = (0, 4, 11, 12)</td><td>0.893</td><td>0.751</td><td>0.954</td><td>0.894</td></tr><tr><td>ξ<sub>CDDs</sub> = (0, 4, 11, 12)</td><td>0.893</td><td>0.751</td><td>0.954</td><td>0.894</td></tr></tbody></table>

Note: D-, Ds- and DDs-efficiencies (respectively E<sub>D,m</sub>, E<sub>D,m</sub> and E<sub>DDs,m</sub>) and their corresponding compound efficiencies (E<sub>CD</sub>, E<sub>CDs</sub> and E<sub>CDDs</sub>) computed for each optimal design when ignoring model uncertainty (ξ<sub>D,m</sub>, ξ<sub>D,m</sub>, ξ<sub>DDs,m</sub>) or accounting for model uncertainty (ξ<sub>CD</sub>, ξ<sub>CDs</sub> and ξ<sub>CDDs</sub>) with each model m (m = 1, ..., 4). The efficiencies are highlighted in lightgrey if between 0.5 and 0.8 and in darkgrey if below 0.5.

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Results of the DDs-optimal designs as function of the weight $a_m$, which quantifies the balance between the D- and the Ds-criteria, varying between 0 and 1. Results for each candidate model are represented by column (M1 in orange, M2 in blue, M3 in green, M4 in purple). First row, location of DDs-optimal measurement times (in months): the symbol x represents the xth measurement in each design. Second row, corresponding D- (solid line) and Ds-efficiency (dotted line), with respect to the D- and Ds-optimal design, respectively.</figcaption></figure>
In the same way as using CD-optimality, assuming model uncertainty leads to CD<sub>S</sub>-optimal design close to D<sub>S</sub>-optimal design for model M3: ξ<sub>CD<sub>S</sub></sub> = (0, 4, 11, 12) and ξ<sub>D<sub>S</sub>, 3</sub> = (0, 5, 11, 12). Furthermore, CD<sub>S</sub>-efficient designs include mostly more tardive measurement times than the CD-efficient designs (Figure 2). The robust CD<sub>S</sub>-optimal design also conducts to a D-efficiency above 0.8 regardless of the candidate model and reduces in average the maximal loss of D<sub>S</sub>-efficiency as compared to the D<sub>S</sub>-optimal designs obtained separately for each model (Table 3, 0.363 vs. respectively 0.566, 0.619, 0.4 and 0.481 with ξ<sub>D<sub>S</sub>, 1</sub>, ξ<sub>D<sub>S</sub>, 2</sub>, ξ<sub>D<sub>S</sub>, 3</sub> and ξ<sub>D<sub>S</sub>, 4</sub>). The CD<sub>S</sub>-optimal design could be particularly appropriate to precisely estimate β, especially if the final model is M3 (D<sub>S</sub>-efficiency of 0.906).

Assuming model uncertainty leads to a robust CDD<sub>S</sub>-optimal design (ξ<sub>CDD<sub>S</sub></sub> = (0, 4, 11, 12)) which is different from the four local DD<sub>S</sub>-optimal designs for each model. The robust CDD<sub>S</sub>-optimal design is the same as the CD<sub>S</sub>- and thus conducts to a D-efficiency above 0.8 regardless of the candidate model. It also reduces the maximal loss of DD<sub>S</sub>-efficiency as compared to the DD<sub>S</sub>-optimal design obtained separately for each model (Table 3, 0.249 vs. respectively 0.356, 0.389, 0.278 and 0.328 with ξ<sub>DD<sub>S</sub>, 1</sub>, ξ<sub>DD<sub>S</sub>, 2</sub>, ξ<sub>DD<sub>S</sub>, 3</sub> and ξ<sub>DD<sub>S</sub>, 4</sub>), and shows DD<sub>S</sub>-efficiencies above 0.8 for M1, M3 and M4.

3.3.4 Expected average power using FIM

Expected power over number of subjects N, under each model and in average, are presented with different designs in Figure 4: the CDD<sub>S</sub>-optimal or robust design (ξ<sub>CDD<sub>S</sub></sub> = (0, 4, 11, 12)), the DD<sub>S</sub>-optimal design for M1 (ξ<sub>DD<sub>S</sub>, 1</sub> = (0, 2, 11, 12)) and the non-optimized equispaced design (ξ<sub>ES</sub> = (0, 4, 8, 12)). The exponential model (M4) presents the worst expected power among the four candidate models, especially with ξ<sub>ES</sub>, mainly because this model involves one more estimated parameter than the three other models. Among these three designs, the ξ<sub>CDD<sub>S</sub></sub> presents the closest curves of power under the four models. Moreover, to reach a π<sub>average</sub> of 0.9, the predicted number of subjects needed (NSN) is 274 with ξ<sub>CDD<sub>S</sub></sub> vs. 320 with ξ<sub>DD<sub>S</sub>, 1</sub> and 358 with ξ<sub>ES</sub>. In our example, the CDD<sub>S</sub>-criterion is able to propose design giving decent performances in terms of predicted power whatever the model chosen among the candidates.

4 Evaluation by simulation

In this section, we perform simulations to assess the relevance of the design strategy presented above. For this purpose, with simulated datasets of 274 subjects (i.e. the NSN to reach an average power of 0.9 with the robust design ξ<sub>CDD<sub>S</sub></sub>), we (1) compare performances of ξ<sub>CDD<sub>S</sub></sub> vs. a locally optimal design for the linear model M1 ξ<sub>DD<sub>S</sub>, 1</sub> and a non-optimized equispaced design ξ<sub>ES</sub> in terms of bias and precision of estimates and (2) evaluate the type I error and power as well as the adequacy between predictions and simulation results with the robust design.

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Predicted power for each model and in average with different designs, for varying total number of subjects from 50 to 450. From the left to the right, with the robust design ξ<sub>CDD<sub>S</sub></sub>, the design optimized for model M1 ξ<sub>DD<sub>S</sub>, 1</sub> and the non-optimized design ξ<sub>ES</sub>, the expected power from FIM for varying number of subjects is displayed in orange for M1, blue for M2, green for M3 and purple for M4. The average power over the four models is in dashed line. The dotted line crossing the x-axis shows the number of individuals needed to obtain an average power of 0.9.</figcaption></figure>
## 4.1 CTS

### 4.1.1 Designs and models

Repeated binary response trials of 274 individuals (137 per treatment group) are simulated under each model *m* (M1 to M4) with their respective population parameters values indicated in Table 2.

For the first objective (to compare different designs performances), with $\xi_{CDD_s}$, $\xi_{DD_{s,1}}$ and $\xi_{ES}$, $K_{1,m} = 500$ datasets are simulated with each model under $H_1$ (treatment effect $\beta = 5$) and each one of these three designs.

To achieve the second objective, with the robust design $\xi_{CDD_s}$, $K_{0,m} = 500$ additional datasets are simulated with each model under $H_0$ (without treatment effect, $\beta = 0$), in order to evaluate the type I error.

### 4.1.2 Evaluation

For each dataset, population parameters are estimated by the SAEM algorithm³² implemented in MONOLIX 2016R1 (www.lixoft.eu), a software devoted to maximum likelihood estimation of parameters in NLMEMs.

First, from the analysis of the $K_{1,m}$ simulated datasets, we calculate the inverse of the variance-covariance matrix computed from 500 estimated parameter vectors, for each of the three designs. We defined the observed relative D-efficiency between designs as the ratio of determinants of these matrices, normalized by the number of parameters. We also compare the distribution of the relative estimation error of parameters for each model (associated with the relative bias values) between the three designs.

Then, with the design $\xi_{CDD_s}$, the adequacy between the FIM predictions and the CTS observations is evaluated (in terms of parameter precision and power). The empirical standard errors ($SE_{CTS}$) of each parameter of each model are defined as the standard deviation of all estimated values. The empirical RSEs are then calculated as the ratio of the $SE_{CTS}$ to the parameter simulated values. The 95% confidence intervals for each empirical RSE are given as $[\sqrt{\frac{q_1 \times RSE^2}{K_{1,m-1}}}, \sqrt{\frac{q_2 \times RSE^2}{K_{1,m-1}}}]$, with $q_1$ and $q_2$, respectively, the 2.5% and 97.5% quantiles of the $\chi^2$ distribution with $K_{1,m-1} = 499$ degrees of freedom. To evaluate the relevance of MC/HMC-based approach to evaluate the expected FIM (see Section 2.1.3.) and the parameter estimation method, predicted RSEs obtained using the expected FIM are compared with empirical RSE and relative root mean square error (RRMSE) obtained with CTS. Thus, the type I error and the power are evaluated. From the $K_{0,m}$ datasets simulated with $\beta = 0$, the type I error of the Wald test on the treatment effect is evaluated as the proportion of trials for which $H_0$ is rejected and compared with 0.05. The 95% prediction interval (95% PI) for proportion $\pi$ of $H_0$ rejection (type I error or power) is computed by binomial test.³³ As explained in Section 3.2.1., we can evaluate the uncertainty in the computation of the D-criterion and the expected power derived from FIM, according to the MC and HMC samples. Under $H_1$, we can thus give a confidence interval as 2.5th and 97.5th percentiles of bootstrap obtained with our FIM computation approach, around the FIM expected power value. The observed proportion of significant Wald tests is compared with the FIM expected power value and its (respectively prediction or confidence) intervals, by binomial test or by bootstrap.

## 4.2 Results

### 4.2.1 Comparison of robust design with non-optimized and non-robust design

Informativeness of the three designs is investigated here. For each model, ratios of observed normalized determinant of the variance-covariance matrix of these designs are computed with the robust design $\xi_{CDD_s}$ as reference.

As presented in Table 4, overall, the estimations of parameters are more precise with the robust design than with the equispaced design for all the models (observed D-efficiency with respect to $\xi_{CDD_s}$ is below 1). With the optimal design for M1 $\xi_{DD_{s,1}}$, the parameters of M1 and M2 are slightly more precise than with the robust design, but this design induces important loss of precision for the parameters of M3 and M4 compared to the robust design.

The three designs show acceptable and similar distribution of relative estimation errors (REE, Figure 5) under M1, M2 and M3. REE under M4 are larger than other models with the three designs. The robust design performs better on M4 than the other designs, especially for the parameters $\mu_2$ and $\omega_2$ which present larger range of whiskers. Moreover, the maximal relative bias with $\xi_{CDD_s}$ is 15% but is over 20% with $\xi_{DD_{s,1}}$ (37% for $\mu_2$ and 35% for $\omega_2$ on M4) and $\xi_{ES}$ (47% for $\mu_2$ and 35% for $\omega_2$ on M4).
<table><caption>**Table 4.** Observed D-efficiency computed from clinical trial simulations for each model, with respect to the robust design.</caption><thead><tr><th></th><th>M1</th><th>M2</th><th>M3</th><th>M4</th></tr></thead><tbody><tr><td>ξ<sub>CDDs</sub>: reference</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>ξ<sub>DDs,1</sub></td><td>1.11</td><td>1.06</td><td>0.79</td><td>0.74</td></tr><tr><td>ξ<sub>ES</sub></td><td>0.98</td><td>0.94</td><td>0.98</td><td>0.84</td></tr></tbody></table>

Note: For each model M1 to M4, the inverse of the variance-covariance matrix of parameter estimates is calculated for each of the three designs: the DDs-optimal design for M1 (ξ<sub>DDs,1</sub>), the non-optimized equispaced design (ξ<sub>ES</sub>) and the robust design (ξ<sub>CDDs</sub>). The observed D-efficiency between designs is defined as the ratio of determinants of these matrices, normalized by the number of parameters.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Boxplot of the relative estimation error (REE, in %, whiskers from 5th to 95th percentile) for each parameter of the four models: M1 in orange, M2 in blue, M3 in green and M4 in purple. From the top to the bottom with 500 datasets and 274 subjects by dataset, the robust design ξ<sub>CDDs</sub>, the optimal design of M1 ξ<sub>DDs,1</sub> and the equispaced design ξ<sub>ES</sub>. The relative bias (RB, in %) of each parameter is positioned under each corresponding boxplot, the values in bold font highlight the RB over 20%. μ are the fixed effects, β the covariate effect and ω the standard deviation of random effects.</figcaption></figure>
### 4.2.2 Adequation between FIM predictions and CTS results with the robust design

As shown in Figure 6 for all models, the predicted RSE by the FIM computation ($RSE_{FIM}$) are in the same range as empirical RSE from CTS ($RSE_{CTS}$) especially for the fixed effect parameters ($\mu$) and for standard deviations of random effects parameters ($\omega$), even when the model is more complex (M4, excepted $\mu_3$). RSE of the covariate effect $\beta$ are, however, slightly under-predicted. Nevertheless, the maximal differences between $RSE_{FIM}$ and $RSE_{CTS}$ are at most 4%. Moreover, the population parameters of all candidate models could be estimated with reasonable estimation error for the proposed robust design (RSE around 30% for the fixed effects and 50% for the random effects) except M4 which contains six parameters: $\theta_2$ is quite difficult to be estimated with such a sparse design.

Table 5 displays the type I error of the Wald test on the treatment effect in each candidate model, using the estimated or empirical SE. For the Wald tests using estimated SE, it can be seen that the type I error lies in the 95% PI of the nominal level [0.033,0.073] for all models. Table 5 shows also the predicted power from FIM vs. observed power from CTS under $H_1$. Two different kinds of intervals are built around the FIM predicted power: as under $H_0$ using the binomial test or as percentiles of bootstrap performed with FIM evaluation. We can note that the interval built by bootstrap is smaller than the binomial test one. The choice of MC and HMC samples seems therefore appropriate for power predictions (see convergence plots in Supplementary Material, Figure S2). Nevertheless, for each model, the observed power is slightly above the predicted power. These results are partly due to the positive bias on $\beta$ (Figure 5) and developed further in the discussion.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Comparison of predicted and observed estimation errors on all parameters of the four candidate models with 274 subjects and the robust design $\xi_{CDDs}$. Relative standard error RSE (%) for the candidate models M1 to M4 are predicted using the expected FIM (full bar). Empirical RSE (hatched bar) are obtained by clinical trial simulations (CTS) with 500 datasets. $\mu$ are the fixed effects, $\beta$ the covariate effect and $\omega$ the standard deviation of random effects.</figcaption></figure>

<table><caption>**Table 5.** Wald tests predictions and simulations on the covariate effect, under $H_0$ and under $H_1$.</caption><thead><tr><th rowspan="3"></th><th colspan="2">Under H<sub>0</sub></th><th colspan="3">Under H<sub>1</sub></th></tr><tr><th rowspan="2">Nominal type I error</th><th rowspan="2">Observed type I error</th><th colspan="2">FIM predicted power</th><th rowspan="2">Observed power</th></tr><tr><th>Binomial P<sub>95</sub></th><th>Bootstrap CI<sub>95</sub></th></tr></thead><tbody><tr><td>M1</td><td></td><td>0.048</td><td>0.967</td><td>[0.949, 0.982]</td><td>[0.963, 0.971]</td><td>0.988</td></tr><tr><td>M2</td><td>0.05</td><td>0.060</td><td>0.889</td><td>[0.859, 0.916]</td><td>[0.873, 0.902]</td><td>0.988</td></tr><tr><td>M3</td><td>[0.033,0.073]</td><td>0.068</td><td>0.973</td><td>[0.959, 0.986]</td><td>[0.969, 0.976]</td><td>0.996</td></tr><tr><td>M4</td><td></td><td>0.036</td><td>0.769</td><td>[0.731, 0.806]</td><td>[0.741, 0.792]</td><td>0.860</td></tr></tbody></table>

Note: Under $H_0$, the 95% prediction interval calculated by binomial test is [0.033,0.073] for 500 simulations. For each model M1 to M4, the predicted power to detect a covariate effect $\beta$ under $H_1$ using the Fisher information matrix (FIM) is calculated with the design $\Xi = \{N = 274, \xi_{CDDs} = (0.4, 11, 12)\}$. Two intervals are presented around the predicted power: a prediction interval obtained by binomial test (PI) and an uncertainty interval (CI) using bootstrap from MC/HMC method for FIM evaluation. The observed power from clinical trial simulations (CTS) is given as the observed proportion of significant Wald tests performed on the covariate effect $\beta$ with estimated standard error from 500 simulated clinical trials with a type I error of 5%.
# 5 Application to a real study

## 5.1 Presentation of the study

We applied the proposed approach to a clinical study conducted by the National Institute on Aging³⁴ using the published model, i.e. not assuming model uncertainty. The study EPESE (Established Populations for Epidemiologic Studies of the Elderly) was a longitudinal study with 4162 persons aged 65 years or older involved. The outcome of interest was the instrumental activities of daily living (IADL). The IADL was measured once a year during three years. The evolution of IADL was compared between those who were cognitively impaired at baseline and those who were cognitively intact.

In Landerman et al.,³⁵ the IADL tasks (traveling, shopping, preparing meals, doing housework and managing finances without assistance) were dichotomized: individuals unable to perform four or five tasks were classified as disabled ($y_{ij}=1$) and those able to perform at least two tasks were not classified as disabled ($y_{ij}=0$). This dichotomous outcome was measured at 0, 1, 2 and 3 years. They studied the cognitive impairment at baseline on the outcome. The proposed model in Landerman et al.³⁵ to describe probability to be disabled at time $t_{ij}$ was

$$ \operatorname{logit}(p(y_{ij} = 1 | b_i)) = \theta_1 + \beta_1 \times 1_{BI} + (\theta_2 + \beta_2 \times 1_{BI})t_{ij} $$

where $1_{BI}$ is the cognitive impairment at baseline ($1_{BI} = 0$ if no impairment, $1_{BI} = 1$ if impairment). $\beta_1$ was the impairment effect on the intercept parameter $\theta_1$ and $\beta_2$ the impairment effect on the slope parameter $\theta_2$. $\theta_1$ follows a normal distribution: $\theta_1 = \mu_1 + b_1$, where $b_1 \sim N(0, \omega_1^2)$. No interindividual variability was considered on $\theta_2$ ($\theta_2 = \mu_2$).

## 5.2 Design evaluation and optimization

### 5.2.1 Methods

To evaluate the relevance of the FIM-based approach, using the EPESE design $\xi_{EPESE} = (0, 1, 2, 3)$ and the provided model described above, we compared the SE predicted using MC/HMC (see Section 2.1.3.) and the SE reported in Landerman et al.³⁵ after data fitting.

Then, based on the DD$_S$-criterion (with the subset of parameters of interests $\psi_S$ composed of $\mu_2$ and $\beta_2$), we propose an optimized design of four measuring times, assuming that $t_1$ and $t_4$ are fixed to 0 and 3 years (i.e. the end of the study), and that only two times $t_2$ and $t_3$ are optimized among 0.2 to 2.8 years by increment of 0.2 years. No repetition is considered, i.e. four measurement times must be different. Combinatorial optimization is performed for different optimality criteria. This corresponds to $\binom{14}{2} = 91$ possible elementary designs to be evaluated, and the optimal one is chosen according to each considered criterion.

Finally, we computed the power of the Wald test and the NSN is calculated (as explained in Section 3.2.3.) for the effect of cognitive impairment on the slope, $\beta_2$, for the slope parameter $\mu_2$.

### 5.2.2 Results

The comparison of the observed and the FIM computed SE of the parameters is reported in Table 6. SE of each fixed or random effect is close between FIM predictions and results reported in Landerman et al.³⁵

<table><caption>**Table 6.** Parameter estimation precision in EPESE analysis (Observed) in Landerman et al.³⁵ and FIM evaluation (Predicted).</caption><thead><tr><th></th><th>Parameter values</th><th>Observed SE</th><th>Predicted SE</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>-7.205</td><td>0.372</td><td>0.270</td></tr><tr><td>μ<sub>2</sub></td><td>1.677</td><td>0.107</td><td>0.074</td></tr><tr><td>β<sub>1</sub></td><td>5.604</td><td>0.384</td><td>0.367</td></tr><tr><td>β<sub>2</sub></td><td>0.034</td><td>0.140</td><td>0.157</td></tr><tr><td>ω<sup>2</sup><sub>1</sub></td><td>15.017</td><td>1.757</td><td>1.316</td></tr></tbody></table>

Note: $\mu$ are the fixed effects, $\omega$ the standard deviation of random effect and $\beta$ the effect of the cognitive impairment.
According to the D-criterion, the D-optimal design with the provided model is $\xi_D = (0, 1.2, 2.8, 3)$. The evolution of the probability to be disabled and the cognitive impairment effect on this evolution were particularly of interest in this study. Thus, we compute the $D_S$-criterion with $\mu_2$ and $\beta_2$ in the subset of parameters $\psi_S$ (see equation (5) in Table 1). The $D_S$-optimal design is $\xi_{D_S} = (0, 1.2, 2.8, 3)$, i.e. the same measuring times allocation of the D-optimal design. Moreover, using the $DD_S$-criterion with $\alpha$ between 0.4 (i.e. D-criterion) and 1 (i.e. $D_S$-criterion) leads to the same optimal design ($\xi_{DD_S} = (0, 1.2, 2.8, 3)$). As for the linear model M1 of Section 3, all the possible designs lead to D-efficiencies higher than 0.8. For instance, with the $\xi_{EPESE}$, the D-efficiency is 0.944.

The provided model in Landerman et al.³⁵ shows that the slope $\mu_2$ is non null. Moreover, only 31 individuals would have been needed with $\xi_{DD_S}$ (vs. 33 with $\xi_{EPESE}$) to reach a power of 0.8. However, the covariate $\beta_2$ is not significant. With the $\xi_{EPESE}$, the NSN to detect this covariate with a power of 0.8 would be 484,518 subjects. With the optimized design $\xi_{DD_S} = (0, 1.2, 2.8, 3)$, it would be 437,087 subjects. Recruiting this number of individuals, even with the optimized design, would be not feasible. Nevertheless, this result highlights that the $DD_S$-optimal design can reduce the NSN compared to a non-optimized design.

# 6 Discussion

In this paper, we propose a new methodological strategy based on MC/HMC and compound optimality theory to design longitudinal trials with discrete outcomes. Our approach accounts for model uncertainty when assuming a set of candidate models and ensures a compromise between the overall precision of estimation and the power of the Wald test to detect a covariate effect. To our knowledge, this is the first time that a robust optimal design approach, based on MC/HMC to compute the FIM without any linearization, is evaluated by extensive CTS. The relevance of using the FIM to efficiently predict the power of the Wald test has been evidenced before but only in the context of repeated continuous outcomes.²¹,³⁶ We have now confirmed the utility of this FIM-based method to predict average power across several candidate models, with an example of repeated discrete outcomes including two treatment groups. Although the considered example is quite theoretical here (Sections 3 and 4), the proposed method is illustrated with a published study analyzing real data (Section 5). In this study, the repeated binary outcome was the ability to undertake life tasks (based on the Instrumental Activities Daily Living tasks). The chosen settings can also easily be extended to clinical longitudinal studies with binary outcomes such as remission of dyskinesia in Parkinson's Disease. In such cases, it is not possible to obtain dense measuring times because medical tests to evaluate this scale are very time-consuming. In another context, a method to design longitudinal studies was proposed when dynamics is modelled by a binary Markov process, using examples of infection by *Streptococcus pneumoniae*.³⁷

To determine informative designs, we study different optimality criteria, according to different purposes of the study. The D-optimality is used to optimize the precision of the whole set of population parameters. The $D_S$-optimality allows to accommodate situations in which only a subset of the model parameters is of interest. It is particularly useful to minimize the standard error on the treatment effect, thus to maximize the power of the Wald test of the study and therefore to reduce the number of subjects needed. A parameter which is a primary end-point (i.e. treatment effect) or a key secondary end-point (i.e. time effect) of a study should be included in the subset of parameters of interest. This subset could also be defined based on the future use of the model. For example, in CTS, if outcomes are sensitive to some specific parameters, they should be included in the subset of parameters of interest. However, this criterion usually leads to designs that do not ensure the experimental identifiability.²² To overcome this issue, the $DD_S$-optimality allows the experimenter to find a compromise between the D- and the $D_S$-criterion,²⁰ by weighting each one according to the importance given to each model parameter. The weight for each criterion depends of the objective of the study, and its meaning is further discussed in Cook and Wong.³⁸ In spite of its usefulness, the $DD_S$-optimality was mostly used with continuous standard nonlinear regression models³⁹,⁴⁰ but has not been reported so far in the context of discrete NLMEMs, to our knowledge. Another novelty of this work is to extend these standard optimal criteria to account for model uncertainty, using their weighted product to propose compound CD-, $CD_S$- and $CDD_S$-criterion over a set of candidate models.

In the motivation example, we considered several candidate binary NLMEMs to describe the logit-probability of response over time: linear (M1), loglinear (M2), quadratic (M3) and exponential (M4). M4 contains one more parameter than other models and therefore provides information about the predicting abilities of the MC/HMC method with a more complex model. The non-optimized equispaced design in our example provides D-efficiencies > 0.8 under each candidate model. This design shows, however, very poor $D_S$-efficiency of 0.393 under M4 (Table 3). When optimizing designs, we find that accounting or not for uncertainty in models and
focusing or not on the estimation precision of the covariate effect may lead to different optimal locations of measurement times. Moreover, misspecification of models in case of local design optimization (i.e. assuming a given model) can lead to D-, D<sub>S</sub>- or DD<sub>S</sub>-efficiency, respectively, as low as 0.666, 0.389 or 0.611 for some models. The proposed compound CD-, CD<sub>S</sub>- or CDD<sub>S</sub>-optimal designs provides a better compromise for different candidate models and greatly reduces the loss of efficiency. These three compound optimal designs conduce all to D-efficiencies of at least 0.8 for each model. Furthermore, in terms of ability to detect the treatment effect, choosing the robust design ξ<sub>CDD<sub>S</sub></sub> instead of the design optimal for the linear model ξ<sub>DD<sub>S</sub>,1</sub> or either a non-optimized ξ<sub>ES</sub> reduces the required sample size of the study. The NSN for an average power of 0.9 is indeed 274 with ξ<sub>CDD<sub>S</sub></sub> vs. 320 with ξ<sub>DD<sub>S</sub>,1</sub> and 358 with ξ<sub>ES</sub>. The NSN is, however, close with either the CDD<sub>S</sub>- or the CD-design (274 vs. 280). This could be explained by the fact that these two designs only differ by their respective second time point (respectively, at fourth and fifth month). However, the gain in power could be larger in other examples, with more differences in allocation of sampling times between the CDD<sub>S</sub>- and the CD-design. The three optimal designs using compound criteria CD-, CD<sub>S</sub>- or CDD<sub>S</sub> are close, which is partly explained by the use of a combinatorial optimization, with only one possible measurement per month, limiting the total number of possible designs. Surprisingly, the second measurement time t<sub>2</sub> of the CD-design is latter than the one of the CD<sub>S</sub>- design, whereas the CD<sub>S</sub>-design maximizes the average power. It is explained by the loglinear model (M2) for which early measurement seems already informative to distinguish the logit probability of response between the two treatment groups. Overall, we notice that with this set of models, picking only tardive measurement could be very suboptimal to maximize the precision of the parameter β and thus the power to detect a treatment effect. The closeness of D- and D<sub>S</sub>-optimal measurement times under each model and averaging over the four models (using CD- and CD<sub>S</sub>-optimality) indicates that a precise estimation of all the parameters (and probably, even more the “trend” parameters μ<sub>2</sub> and ω<sub>2</sub> on which β acts) is necessary to reach a satisfactory power. In the optimization procedure, it would also be interesting to consider reducing the duration of the study.

We also evaluate our approach by simulating clinical trials under each candidate model. To assess the relevance of the robust design, we compare performances between three different designs: a non-optimized design with a measurement time every four months, a design optimal for a linear model (M1, maximizing the DDs-optimality) and a robust design accounting for model uncertainty (maximizing the compound DDs-optimality). In our example, choosing the robust design avoids important loss of information on parameters, especially when the right model was exponential (M4). Other scenarios could have been examined, e.g. simulating data under another model M5 and fit under each candidate model M1 to M4. We also could evaluate the power performances of model averaging<sup>23</sup> vs. model selection as in Buatois et al.,<sup>28</sup> but it is not the purpose of the present work. Furthermore, a method has been recently proposed to optimize designs when the analysis is done by model averaging.<sup>41</sup>

With the CDDs-design and 274 individuals, in spite of a slight under or over-prediction of RSE for some
parameters, the discrepancy is at most 4%; thus, the prediction using the FIM evaluated by MC/HMC is accurate.
This approach avoids extensive CTS and can efficiently help to easily detect non-informative designs in case of
large predicted RSE values. Using the predicted SE given by FIM, we predict the power of the Wald test to show
significant difference between treatment groups. Under H0, we note a good control of the type I error by CTS,
despite the asymptotic properties of the Wald test (explained in Dubois et al.⁴²) and our sparse measurement
design. Under H1, we found differences in predicted power between the different models: respectively 0.967, 0.889,
0.973 and 0.769 under M1, M2, M3 and M4. Whereas the SE for treatment effect (β) are not over-predicted,
powers are under-predicted under all the models M1 to M4 (especially M2 and M4): observed power are,
respectively, 0.988, 0.988, 0.996 and 0.86. However, we note that the power is never over-predicted. Moreover,
the FIM predicted Wald statistic is close to the median of the observed Wald statistic distribution, even with M2
and M4 (Supplementary Material, Figure S2). We investigated the possibility that the number of MC/HMC
samples would be insufficient under our design to adequately predict the precision on parameter β (and thus of
the power). However, the convergence plots of the determinant of the FIM and of the SE on β showed acceptable
results (Figure S3, Supplementary Material). Discrepancies between predicted and observed power are mainly
explained by the positive bias on β in CTS. If β would be estimated at its true simulated value (5), the estimated
power would be close to the predicted values (Supplementary Material, Table S1). Other settings of the SAEM
algorithm of MONOLIX were investigated, but none provided more satisfactory relative bias. To our knowledge,
an evaluation of different algorithms for repeated binary responses has never been published, as it has been done
for ordered categorical<sup>43</sup> or continuous data.<sup>44</sup> The analysis of this kind of data is challenging, and requires
informative design to avoid problem of bias. Further research works on the estimation and design method for
this kind of data are necessary. It could thus lead to fair measures of discrepancies between FIM prediction and
CTS observations, and also to a fair comparison of different designs in terms of average power.
We apply the described FIM based methods, using a published logistic model built on real data.³⁵ This model is describing the evolution of the disability probability over time, including parameters as the baseline (fixed and random effect), the time effect or the cognitive impairment effect. With the model and the design elements provided in the study, we find that the SE predicted by FIM are close to those provided in the study, showing that the FIM evaluation method by MC/HMC is relevant. The provided model is the only one mixed model which is reported. If other possible models would have been proposed, we would be able to compute compound criteria, with each model weight depending on its likelihood, in accordance with the model averaging principle.²³ Therefore, a robust design based on the CDD<sub>S</sub>-criterion would have been proposed. However, with the given model, we propose to optimize the choice of the measuring times according to the DD<sub>S</sub>-criterion. We also predict the number of subjects that would be needed to reach a power of 0.8 to detect the time effect and the cognitive impairment. We show that, instead of using a non-optimized equispaced design, optimizing the design using the DD<sub>S</sub>-criterion can reduce the number of subjects needed.

In this work, we perform combinatorial optimization by evaluating FIM for every possible design, an
alternative approach would consist in implementing an optimization algorithm to improve computing time as
suggested in Riviere et al.¹³ More efficient algorithms could also be explored, such as stochastic approximation,
annealing methods, randomized exchange algorithm,⁴⁵ multiplicative method⁴⁶ or particle swarm
optimization.⁴⁷,⁴⁸ It is also important to account for uncertainty in parameters in addition to uncertainty in
models,⁴⁹,⁵⁰ by evaluating the expectation of the FIM over the distribution of population parameters instead of
assuming known values of these parameters. We also plan to integrate these robust design methods (over the
parameter distribution and over the models) in a next version of PFIM, which is the software program for
designing longitudinal studies developed by our team (www.pfim.biostat.fr).⁵¹ Furthermore, robust design
methods studied in this paper could be combined with adaptive designs,⁵²,⁵³ by using the accumulated
information at each stage of the study to update knowledge of the candidate models and of the parameter
distributions to be taken into account in design optimization.

The proposed design strategy based on the expected FIM evaluated by MC/HMC and the compound optimality theory is a relevant approach, which enables, for the first time, optimization of informative sparse designs for longitudinal binary data, over several candidate models. This approach can also be applied to other type of discrete data models, such as multi-category or Poisson count data.⁵⁴ A robust design method and application using *MIXFIM* for Poisson models was studied in Loingeville et al.⁵⁰ but not for the comparison of two treatment groups, hence using CD-optimality. The approach described in this article also aims to propose design, providing good power to detect the treatment effect and acceptable precision of all parameters while accounting for model uncertainty.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this
article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article

ORCID iD

Jérémy Seurat https://orcid.org/0000-0002-8127-7751

Supplemental material

Supplemental material for this article is available online.

References

1. Lalonde RL, Kowalski KG, Hutmacher MM, et al. Model-based drug development. *Clin Pharmacol Ther* 2007; **82**: 21–32.

2. Mould DR and Upton RN. Basic concepts in population modeling, simulation, and model-based drug development. *CPT Pharmacometrics Syst Pharmacol* 2012; **1**: e6.
3. Mould DR and Upton RN. Basic concepts in population modeling, simulation, and model-based drug development – part 2: introduction to pharmacokinetic modeling methods. *CPT Pharmacomet Syst Pharmacol* 2013; **2**: e38.

4. Upton RN and Mould DR. Basic concepts in population modeling, simulation, and model-based drug development: part 3 – introduction to pharmacodynamic modeling methods. *CPT Pharmacomet Syst Pharmacol* 2014; **3**: e88.

5. Syrjänen RK, Kilpi TM, Kaijalainen TH, et al. Nasopharyngeal carriage of *Streptococcus pneumoniae* in Finnish children younger than 2 years old. *J Infect Dis* 2001; **184**: 451–459.

6. Mestre TA, Beaulieu-Boire I, Aquino CC, et al. What is a clinically important change in the Unified Dyskinesia Rating Scale in Parkinson's disease? *Parkinsonism Relat Disord* 2015; **21**: 1349–1354.

7. Atkinson A, Donev A and Tobias R. *Optimum experimental designs, with SAS*. Oxford, New York: Oxford University Press, 2007.

8. Pronzato L and Pázman A. *Design of experiments in nonlinear models: asymptotic normality, optimality criteria and small-sample properties*. New York: Springer Science & Business Media, 2013.

9. Fedorov V and Leonov S. *Optimal design for nonlinear response models*. Boca Raton: CRC Press, 2013, pp.20–29.

10. Mentré F, Mallet A and Baccar D. Optimal design in random-effects regression models. *Biometrika* 1997; **84**: 429–442.

11. Nyberg J, Bazzoli C, Ogungbenro K, et al. Methods and software tools for design evaluation in population pharmacokinetics-pharmacodynamics studies. *Br J Clin Pharmacol* 2015; **79**: 6–17.

12. Nguyen TT and Mentré F. Evaluation of the Fisher information matrix in nonlinear mixed effect models using adaptive Gaussian quadrature. *Comput Stat Data Anal* 2014; **80**: 57–69.

13. Riviere M-K, Ueckert S and Mentré F. An MCMC method for the evaluation of the Fisher information matrix for nonlinear mixed effect models. *Biostatistics* 2016; **17**: 737–750.

14. Mielke T. *Approximation of the Fisher information and design in nonlinear mixed effects models*. PhD Thesis. Otto-von-Guericke-Universität Magdeburg, Germany, 2012.

15. Waite TW and Woods DC. Designs for generalized linear models with random block effects via information matrix approximations. *Biometrika* 2015; **102**: 677–693.

16. Waite TW. *Design of experiments with mixed effects and discrete responses plus related topics*. PhD Thesis. University of Southampton, United Kingdom, 2012.

17. Ueckert S and Mentré F. A new method for evaluation of the Fisher information matrix for discrete mixed effect models using Monte Carlo sampling and adaptive Gaussian quadrature. *Comput Stat Data Anal* 2017; **111**: 203–219.

18. Riviere M-K and Mentré F. R package MIXFIM, version 1.0 http://mc-stan.org/ (accessed 17 August 2015).

19. Stan Development Team. Rstan: the R interface to Stan, version 2.12.0. http://mc-stan.org/ (accessed 19 January 2018).

20. Atkinson AC and Bogacka B. Compound D- and Ds-optimum designs for determining the order of a chemical reaction. *Technometrics* 1997; **39**: 347–356.

21. Retout S, Comets E, Samson A, et al. Design in nonlinear mixed effects models: optimization using the Fedorov-Wynn algorithm and power of the Wald test for binary covariates. *Stat Med* 2007; **26**: 5162–5179.

22. Sölkner J. Choice of optimality criteria for the design of crossbreeding experiments. *J Anim Sci* 1993; **71**: 2867–2873.

23. Buckland ST, Burnham KP and Augustin NH. Model selection: an integral part of inference. *Biometrics* 1997; **53**: 603–618.

24. Atkinson AC. DT-optimum designs for model discrimination and parameter estimation. *J Stat Plan Infer* 2008; **138**: 56–64.

25. Nguyen TT, Bénech H, Delaforge M, et al. Design optimisation for pharmacokinetic modeling of a cocktail of phenotyping drugs. *Pharm Stat* 2016; **15**: 165–177.

26. Schorning K, Bornkamp B, Bretz F, et al. Model selection versus model averaging in dose finding studies. *Stat Med* 2016; **35**: 4021–4040.

27. Aoki Y, Röshammar D, Hamrén B, et al. Model selection and averaging of nonlinear mixed-effect models for robust phase III dose selection. *J Pharmacokinet Pharmacodyn* 2017; **44**: 581–597.

28. Buatois S, Ueckert S, Frey N, et al. Comparison of model averaging and model selection in dose finding trials analyzed by nonlinear mixed effect models. *AAPS J* 2018; **20**: 56.

29. Lestini G, Ueckert S and Mentré F. Model-based optimal robust design in pharmacometrics. *PODE, Uppsala, Sweden*, 2016, www.maths.qmul.ac.uk/~bb/PODE/PODE2016_JLestini.pdf (accessed 30 May 2017).

30. Kiefer J and Wolfowitz J. The equivalence of two extremum problems. *Can J Math* 1960; **12**: 363–366.

31. Ogungbenro K and Aarons L. Population Fisher information matrix and optimal design of discrete data responses in population pharmacodynamic experiments. *J Pharmacokinet Pharmacodyn* 2011; **38**: 449–469.

32. Kuhn E and Lavielle M. Maximum likelihood estimation in nonlinear mixed effects models. *Comput Stat Data Anal* 2005; **49**: 1020–1038.

33. Clopper CJ and Pearson ES. The use of confidence or fiducial limits illustrated in the case of the binomial. *Biometrika* 1934; **26**: 404–413.

34. Blazer D, Burchett B, Service C, et al. The association of age and depression among the elderly: an epidemiologic exploration. *J Gerontol* 1991; **46**: M210–M215.

35. Landerman LR, Mustillo SA and Land KC. Modeling repeated measures of dichotomous data: testing whether the within-person trajectory of change varies across levels of between-person factors. *Soc Sci Res* 2011; **40**: 1456–1464.

36. Nguyen TT, Bazzoli C and Mentré F. Design evaluation and optimisation in crossover pharmacokinetic studies analysed by nonlinear mixed effects models. *Stat Med* 2012; **31**: 1043–1058.
37. Mehtälä J, Auranen K and Kulathinal S. Optimal designs for epidemiologic longitudinal studies with binary outcomes. *Stat Methods Med Res* 2015; **24**: 803–818.

38. Cook RD and Wong WK. On the equivalence of constrained and compound optimal designs. *J Am Stat Assoc* 1994; **89**: 687–692.

39. Yeatts SD, Gennings C and Crofton KM. Optimal design for the precise estimation of an interaction threshold: the impact of exposure to a mixture of 18 polyhalogenated aromatic hydrocarbons. *Risk Anal* 2012; **32**: 1784–1797.

40. Winkens B, Schouten HJA, van Breukelen GJP, et al. Optimal designs for clinical trials with second-order polynomial treatment effects. *Stat Methods Med Res* 2007; **16**: 523–537.

41. Alhorn K, Schorning K and Dette H. Optimal designs for frequentist model averaging. arXiv:180705234 [stat], http://arxiv.org/abs/1807.05234 (2018, accessed 1 August 2018).

42. Dubois A, Lavielle M, Gsteiger S, et al. Model-based analyses of bioequivalence crossover trials using the stochastic approximation expectation maximisation algorithm. *Stat Med* 2011; **30**: 2582–2600.

43. Savic RM, Mentré F and Lavielle M. Implementation and evaluation of the SAEM algorithm for longitudinal ordered categorical data with an illustration in pharmacokinetics-pharmacodynamics. *AAPS J* 2010; **13**: 44–53.

44. Plan EL, Maloney A, Mentré F, et al. Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models. *AAPS J* 2012; **14**: 420–432.

45. Harman R, Filová L and Richtárik P. A randomized exchange algorithm for computing optimal approximate designs of experiments. *Journal of the American Statistical Association*. In press. https://www.tandfonline.com/doi/full/10.1080/01621459.2018.1546588.

46. Yu Y. Monotonic convergence of a general algorithm for computing optimal designs. *Ann Statist* 2010; **38**: 1593–1606.

47. Wong WK, Chen R-B, Huang C-C, et al. A modified particle swarm optimization technique for finding optimal designs for mixture models. *PloS One* 2015; **10**: e0124720.

48. Kim S and Wong WK. Extended two-stage adaptive designs with three target responses for phase II clinical trials. *Stat Methods Med Res* 2018; **27**: 3628–3642.

49. Foo LK, McGree J, Eccleston J, et al. Comparison of robust criteria for D-optimal designs. *J Biopharm Stat* 2012; **22**: 1193–1205.

50. Loingeville F, Nguyen TT, Riviere M-K, et al. Robust designs in longitudinal studies accounting for parameter and model uncertainties - application to count data. *J Biopharm Stat*. In press. https://www.tandfonline.com/doi/abs/10.1080/10543406.2019.1607367?journalCode=lbps20.

51. Bazzoli C, Retout S and Mentré F. Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0. *Comput Methods Programs Biomed* 2010; **98**: 55–65.

52. Lestini G, Dumont C and Mentré F. Influence of the size of cohorts in adaptive design for nonlinear mixed effects models: an evaluation by simulation for a pharmacokinetic and pharmacodynamic model for a biomarker in oncology. *Pharm Res* 2015; **32**: 3159–3169.

53. Strömberg EA and Hooker AC. The effect of using a robust optimality criterion in model based adaptive optimization. *J Pharmacokinet Pharmacodyn* 2017; **44**: 317–324.

54. Fang X, Li J, Wong WK, et al. Detecting the violation of variance homogeneity in mixed models. *Stat Methods Med Res* 2016; **25**: 2506–2520.
Supplementary Material

<figure><img src="image_12.png" /><figcaption>**Figure S1.** Convergence plots of D-criterion with a non-optimized design. The D-criterion is plotted in black as a function of the number of Monte Carlo samples, with the equispaced design, for each model M1 to M4 by column. For the computation of the Fisher information matrix, two different algorithm settings are displayed by row (1000 Hamiltonian Monte-Carlo (HMC) samples with 1000 burn-in on the left or 200 HMC with 500 burn-in on the right). The red dotted lines represent the 2.5th and the 97.5th percentiles of the D-optimality values using bootstrap.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure S2.** Convergence plots of D-criterion and power with the robust design. For each model M1 to M4 by column, the D-criterion on the top part and the power on the bottom part, are plotted in black as a function of the number of Monte Carlo samples with the robust design. Different algorithm settings for the computation of the Fisher information matrix are represented (1000 Hamiltonian Monte-Carlo (HMC) samples and 1000 burn-in on the first row of each part or 200 HMC and 500 burn-in on the second row of each part). The red dotted lines represent the 2.5th and the 97.5th percentiles of the bootstrap values</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure S3.** Distribution of observed Wald test statistic with the robust design. For each model M1 to M4, the distribution of the 500 observed Wald test statistic W is represented as histograms. For each model, the solid line is the FIM predicted W and the dashed line is the median of observed W.</figcaption></figure>
<table><caption>**Table S1.** Observed power of the Wald test without bias on with the robust design.</caption><thead><tr><th></th><th>M1</th><th>M2</th><th>M3</th><th>M4</th></tr></thead><tbody><tr><td>FIM predicted power</td><td>0.968</td><td>0.889</td><td>0.973</td><td>0.77</td></tr><tr><td>CTS observed power</td><td>0.988</td><td>0.988</td><td>0.996</td><td>0.86</td></tr><tr><td>CTS observed power using simulated</td><td>0.976</td><td>0.914</td><td>0.992</td><td>0.746</td></tr></tbody></table>

FIM predicted power are obtained using the true value of (5) and the SE() calculated from FIM evaluation. CTS observed power are obtained using estimated and their SE() as given by the SAEM algorithm. CTS observed power using simulated are obtained using the true value of (5) and SE() as given by the SAEM algorithm.
## 2.2. Algorithmes pour l'optimisation de protocoles

### 2.2.1. Algorithmes multiplicatifs

#### 2.2.1.1. Résumé

Nous avons proposé d'étendre un algorithme dans le but d'optimiser des protocoles d'études longitudinales, en utilisant la matrice d'information de Fisher. Nous nous sommes intéressés aux algorithmes multiplicatifs, conçus pour l'optimisation de protocoles dans un espace de design discret. Les algorithmes multiplicatifs ont été principalement utilisés pour les régressions individuelles. Dans ce travail, ils ont été étendus aux modèles mixtes et ont été codés dans deux langages : R et C.

Contrairement à l'utilisation de l'algorithme de Fedorov-Wynn, l'utilisateur n'a pas besoin de définir un protocole pour l'initialisation de cet algorithme multiplicatif. Toutefois, il est nécessaire que la matrice de Fisher soit calculée pour chaque protocole élémentaire possible $\xi_q$, avec $q = (1, ..., Q)$. Les $\xi_q$ sont définis à partir de contraintes (par exemple, chaque $\xi_q$ peut être défini par 1 dose parmi 10 possibles et un vecteur de 5 temps de mesures parmi 8 possibles, soit un nombre total de protocoles élémentaires $Q = 10 \times 56 = 560$). A l'initialisation de l'algorithme, le même poids (c'est à dire le même nombre de sujets $N_q$) est attribué à chaque $\xi_q$. Ensuite, à partir du critère D, les poids optimaux de chaque $\xi_q$ sont calculés par l'algorithme suivant un processus itératif. La monotonicité et la convergence de l'algorithme ont été prouvées à partir d'une précédente démonstration de la convergence des algorithmes multiplicatifs (Yu 2010).

Deux exemples illustratifs d'optimisation de protocole d'études longitudinales ont été choisis afin d'appliquer l'algorithme développé. Le premier exemple traite d'une optimisation de temps de mesures et/ou de choix de doses pour une étude longitudinale sur la dégénérescence maculaire, en utilisant un modèle caractérisant l'évolution de l'acuité visuelle, définie comme une valeur continue (Rosenfeld et al. 2006; Buatois et al. 2018). La FIM a été évaluée par méthode FO. Deux cas d'optimisation de bras de doses parmi 7 possibles (et du nombre d'individus dans chaque bras) ont été réalisés, avec respectivement 26 et 4 temps de mesures fixés. L'optimisation des doses (et du nombre de sujets à attribuer à chaque dose) par l'algorithme a conduit à un gain de D-efficacité de respectivement 13 et 14 % par rapport aux protocoles de référence respectifs (non optimisés, avec patients équitablement répartis dans 4
bras de doses). Un exemple plus complexe a été étudié, avec recherche à la fois des doses et des temps de mesures les plus informatifs (1932 protocoles élémentaires possibles). Le protocole de population optimisé était composé de 4 groupes avec des doses et temps de mesures différents. Ce protocole optimisé conduisait à un gain d'efficacité de 26%. Cet exemple a permis aussi de comparer la vitesse des versions codées (en C et en R) de l'algorithme. L'algorithme en C était 14 fois plus rapide que l'algorithme en R pour cet exemple (31 s vs. 432 s, respectivement). Le deuxième exemple étudié, aussi utilisé pour appliquer la méthode de planification de protocoles robustes (section 2.1.), consiste en une optimisation de temps de mesures (2 parmi 4) pour des données binaires analysées par modèle logistique, avec 55 protocoles élémentaires possibles. L'algorithme a permis de trouver un protocole présentant un gain de D-efficacité de 18 % par rapport au protocole de référence, non optimisé (avec le même intervalle de temps entre deux mesures).

Ces travaux ont montré que l'utilisation d'un algorithme multiplicatif est une solution efficace pour l'optimisation de protocoles d'études longitudinales. Ils ont fait l'objet d'un article en révision dans *Computer Methods and Programs in Biomedicine*.

L'algorithme présenté a été intégré à la nouvelle version du programme PFIM (section 2.2.2.).

### 2.2.1.2. Article 2 (accepté)

Finding optimal design in nonlinear mixed effect models using multiplicative algorithms, publié
dans *Computer Methods and Programs in Biomedicine*, 2021

Jérémy Seurat*, Yuxin Tang*, France Mentré, Thu Thuy Nguyen

<footnote>* : co-premier auteur</footnote>
# Finding optimal design in nonlinear mixed effect models using multiplicative algorithms

Jérémy Seurat*, Yuxin Tang, France Mentré, Thu Thuy Nguyen

Université de Paris, INSERM, IAME, F-75006 Paris, France

## Abstract

**Background and objectives:** To optimize designs for longitudinal studies analyzed by nonlinear mixed effect models (NLMEMs), the Fisher information matrix (FIM) can be used. In this work, we focused on the multiplicative algorithms, previously applied in standard individual regression, to find optimal designs for NLMEMs. **Methods:** We extended multiplicative algorithms to mixed models and implemented the algorithm both in R and in C. Then, we applied the algorithm to find D-optimal designs in two longitudinal data examples, one with continuous and one with binary outcome. **Results:** For these examples, the algorithm performed faster in C than in R. Design optimization using the multiplicative algorithm led to designs with D-efficiency gains between 13% and 25% compared to non-optimized designs. **Conclusion:** We found that the multiplicative algorithm can be used efficiently to design longitudinal studies.

**Keywords:** nonlinear mixed effect model, Fisher information matrix, multiplicative algorithm, optimal design, D-optimality, longitudinal studies

The first two authors equally contributed.

<footnote>*Corresponding author
Email address: jeremy.seurat@inserm.fr (Thu Thuy Nguyen)</footnote>
# 1. Introduction

Nonlinear mixed effect models (NLMEMs) are frequently used in model-based drug development to analyze all longitudinal data obtained during clinical trials, notably pharmacokinetic/pharmacodynamic data [1, 2, 3, 4]. NLMEMs can compensate for the lack of individual information by borrowing the strength from the longitudinal data in the whole population and therefore allow for precise parameter estimation even with sparse designs, where few samples are collected from each subject, according to experimental constraints. Nevertheless, it is crucial to find informative experimental designs before modeling to get a precise estimation of parameters [5]. Indeed, the informativeness of a dataset for parameter estimation depends on the design elements as the number of subjects, the number and time of the measurements or the doses.

To evaluate and compare designs, the theory of optimum experimental design in classic nonlinear models was introduced [6, 7, 8], based on the expected Fisher information matrix (FIM), as an alternative to time-consuming clinical trial simulation. The inverse of the FIM, according to the Cramer-Rao inequality, is the lower bound of the variance covariance matrix of any unbiased estimators of the parameters. Then, the predicted standard errors for estimated parameters can be calculated, as the square roots of the diagonal elements of the inverse of the FIM.

Once the FIM is evaluated, criteria based on FIM can be computed to optimize designs. The D-optimality criterion is widely used for design optimization. It consists in maximizing the determinant of the FIM, normalized by the number of parameters to be estimated. Several algorithms have been developed for design optimization [6]. Particularly in PFIM (www.pfim.biostat.fr), the first R program dedicated to design evaluation and optimization of longitudinal studies [9], Simplex [10] and Fedorov-Wynn [11, 12, 13] algorithms were implemented. The Simplex algorithm optimized designs within a continuous design space while the Fedorov-Wynn algorithm works in a discrete design space. Each type of algorithm corresponds to a different optimization strategy. The Simplex