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

HAL
open science

# Stochastic Expectation Maximization Algorithm for Linear Mixed-Effects Model with Interactions in the Presence of Incomplete Data

Alandra Zakkour, Cyril Perret, Yousri Slaoui

► To cite this version:

Alandra Zakkour, Cyril Perret, Yousri Slaoui. Stochastic Expectation Maximization Algorithm for Linear Mixed-Effects Model with Interactions in the Presence of Incomplete Data. Entropy, 2023, 25 (3), pp.473. 10.3390/e25030473. hal-04026682

HAL Id: hal-04026682

https://hal.science/hal-04026682

Submitted on 13 Mar 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Stochastic Expectation Maximization Algorithm for Linear Mixed-Effects Model with Interactions in the Presence of Incomplete Data

Alandra Zakkour <fnref n="1" />,<fnref n="2" />, Cyril Perret <fnref n="1" />,<fnref n="2" /> and Yousri Slaoui <fnref n="1" />,*

**Abstract:** The purpose of this paper is to propose a new algorithm based on stochastic expectation maximization (SEM) to deal with the problem of unobserved values when multiple interactions in a linear mixed-effects model (LMEM) are present. We test the effectiveness of the proposed algorithm with the stochastic approximation expectation maximization (SAEM) and Monte Carlo Markov chain (MCMC) algorithms. This comparison is implemented to highlight the importance of including the maximum effects that can affect the model. The applications are made on both simulated psychological and real data. The findings demonstrate that our proposed SEM algorithm is highly preferable to the other competitor algorithms.

**Keywords:** linear mixed-effects model; interactions; missing data; censored data; EM algorithm; SEM algorithm

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

check for
updates

**Citation:** Zakkour, A.; Perret, C.;

Slaoui, Y. Stochastic Expectation
Maximization Algorithm for Linear
Mixed-Effects Model with
Interactions in the Presence of
Incomplete Data. *Entropy* **2023**, *25*,
473. https://doi.org/10.3390/
e25030473

Academic Editors: Sergio Curilef and
Francisco Calderón

Received: 9 February 2023

Revised: 27 February 2023

Accepted: 2 March 2023

Published: 8 March 2023

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

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

# 1. Introduction

The time between the presentation of a stimulus and a participant's motor response is the oldest and most widely used measure for exploring the functioning of the human mind. In 1869, ref. [1] theorized this duration, called reaction time (RT), as involving three sets of activities: perceptual mechanisms, cognitive processing and motor preparation. Based on the assumption that the first and last sets of processing can be considered as having virtually identical duration for the same task, any change in RTs between two experimental conditions is then interpreted as indicating a change in the duration of cognitive processing. RT is then considered by psychologists as a tool to explore cognitive processing mechanisms ([2]).

Psycholinguistics research on the cognitive mechanisms involved in language recognition or production frequently uses RT as a measure of behavior. Like all scientific disciplines, psycholinguistics relies on hypothesis testing to support its theoretical propositions. Since the early 2000s, researchers have taken up linear mixed-effects models (LMEMs). As described by [3], a LMEM allows for the proper consideration of one of the characteristics of psycholinguistic experiments: the presence of two random-effect variables. Indeed, the experimental structure of the experiments conducted in this field involves having a group of participants process a set of stimuli (see hereafter). Thus, the statistical analysis must allow the inclusion of these two random-effect variables, i.e., participants and items, in the structure of the model. The introduction of LMEM was thus an important methodological advance for psycholinguistics.

There is, however, one point that has not received much attention from psycholinguists. Experimentation with human beings is often subject to many vagaries. Imagine researchers whose goal is to understand how an adult retrieves the spelling of a word from their

<footnote n="1">¹ Laboratoire de Mathématiques et Applications, Université de Poitiers, 11 Boulevard Marie et Pierre Curie, 86962 Futuroscope Chasseneuil, CEDEX 9, 86073 Poitiers, France</footnote>

<footnote n="2">² CeRCA-CNRS UMR 7295, Université de Poitiers, 5 rue T. Lefebvre, MSHS, CEDEX 9, 86073 Poitiers, France</footnote>

<footnote>* Correspondence: yoursri.slaoui@math.univ-poitiers.fr</footnote>
memory ([4]). They ask a group of thirty participants to write down the names of 150 images presented on a computer screen. A device allows them to measure the time between the appearance of the image on the screen and the first writing gesture of the participant. These researchers can potentially collect 4500 RT values. However, between trials lost for technical reasons, spelling errors, and certain habits during data analysis (e.g., censoring of data greater than two standard deviations from the mean, right-censoring data), a non-negligible number of data are removed. For example, ref. [4] reported removing just over 20% of the data.

The issue of missing and censored data has received relatively little attention in psycholinguistics (see, however, [5,6]). This is especially true since the introduction of LMEMs. Psycholinguists assume that these models can be run on a sample of data with “holes”. The analysis strategy is then of the “keep the case empty” type, ignoring the bias that this introduces in the estimation of the parameters of the model and thus of the decisions taken. The objective of this work is thus to develop algorithms to manage the presence of missing and censored data in these psycholinguistic experiments. Two points are to be taken into account. On one hand, experimental designs in scientific psychology may involve assumptions about interactions between two or more fixed-effect variables. On the other hand, researchers suggest, for theoretical reasons, that all interactions between fixed-effect variables and random-effect variables made possible by the model should be included ([7]). The potential presence of these two types of interactions (fixed–fixed variables; random–fixed variables) are constraints on the development of the missing data procedure.

Let us first recall that the *LMEM* ([8]) is an extension of the simple linear model that allows both fixed and random effects to be represented. In 1861, the *LMEM* was introduced under the name of a one-way random-effects model ([9]), that is, a model with one random variable and without any fixed variable. From 1990 and onward, ref. [9] underlined that *LMEMs* became popular in many research applications including economics, sociology, insurance, agronomy, epidemiology, genetics, etc. They are used in longitudinal data analysis, multilevel modeling and small estimations. The analysis of this type of models is presented in [10]; for more literature about the methodology, theoretical results and software, see the books [11–14]. *LMEMs* are fitted and analyzed in R by using the package *lme4* or *lmertest* ([15]).

In [16], the author presented the fixed interactions between two factors in an *LMEM*, using the maximum product interaction F-test. Ref. [17] was interested in the sample size of an *LMEM*. He proposed a formula to estimate the sample size based on testing a mixed model that contained fixed interactions. In 2019, ref. [18] proposed an estimation method to recover the principal effects and interactions, because the existing method did not allow the integration of these effects in a mixed data frame. They approved that the proposed method gave optimal results to their applications’ conditions. On the same side, ref. [19] presented the estimation of the fixed interactions. These interactions are normally introduced by the product of the variables, but the algebraic transformations revealed that this technique did not produce a within-unit estimator. A Monte Carlo method confirmed that the FE (fixed-effects) estimator of (x, z) was biased, if one of these variables was correlated with another one. In order to present the interaction between x and z, it is possible to use the current syntax `x * z`. This consideration is called “double-demeaned”; it is less efficient than the standard FE and only works with *T* < 2. For the application in R, the Markov chain Monte Carlo method is applied using the package *mcmc*. Ref. [20] introduced the *MCMCpack* package that contains functions to perform Bayesian inference using a posterior simulation for a number of statistical models (https://CRAN.R-project.org/package=MCMCpack (accessed on 14 Jun 2011)), while [21] presented the package *MCMCglmm* (https://cran.r-project.org/web/packages/MCMCglmm/index.html (accessed on 2 February 2010)) for the MCMC method to fit the generalized linear mixed-models for multiple response vectors. This method was also used by [22] to identify unknown parameters in the biological field, such as detecting the concentration of target molecules, because of the importance of this
method for extracting the information. The MCMC method was applied in the work of [23]
to solve mechanical problems by using a Bayesian method.

Ref. [24] presented five types of LMEMs by introducing two random variables (see Appendix A.1), where the models were compared to show the importance of including the maximum effects that could affect the model.

One of the solutions that handles the missing data can be the expectation-maximization (EM) algorithm. It is a very useful algorithm for the estimation of the maximum likelihood function. This method can be a solution when the only data available do not allow the estimation of the parameters ([25]), and/or the expression of the likelihood is analytically impossible to maximize.

In other words, it aims to provide an estimator when the problem comes from the
presence of missing data. When the data are incomplete, the knowledge of these values
would make it possible to estimate the parameters. The EM algorithm takes its name from
the fact that at each iteration, it operates two distinct steps:

**Step E:** This step is called the expectation (E) step, where we are interested in finding the expected value of the unobserved or unknown variables given the observed data and the value of the parameters.

**Step M:** This step is called the maximization (M) step; in this step, we maximize the expected log-likelihood by using the estimation of the unknown data carried out in the previous step. These parameter estimates are then used to determine the distribution of the unknown variables in the next iteration.

At some points, the expectation or the maximization steps are impossible to apply directly ([26]); from there, the use of an extension form of EM is useful, such as MCEM, where the (E) step is replaced by a Monte Carlo simulation, or SAEM, where the (E) step is replaced by a stochastic approximation. SAEM was a solution of nonlinear mixed-effects models (NLMEM). Ref. [27] proposed a new methodology for maximum likelihood estimation in mixtures of nonlinear mixed-effects models. The resulting MSAEM (mixture SAEM) algorithm is now implemented in the Monolix software tool.

The aim of this paper was to perform the SEM algorithm, under an LMEM by including
two types of incomplete data (the censored and the MAR types) and by taking into consider-
ation for the first time the interactions, where our proposed model contains the interactions
between fixed variables and fixed-random variables. This document is organized as fol-
lows: In Section 2, we present three algorithms based on the expectation-maximization (EM)
method, the first one is called the SAEM algorithm, the second is our proposed SEM algorithm
and the third one is based on the MCMC method. We also present the incomplete data types,
divided into missing data and censored one. In Section 3, we define the proposed model
with some specific cases. In Section 4, we introduce a method to achieve the convergence.
In Section 5, we compare the results obtained from simulated psychological and real data.
In Section 6, we conclude the proposed study with some perspectives and future directions.

## 2. Methodology

### 2.1. *EM Algorithm*

As previously said, the EM algorithm is a widely used algorithm in the case of
incomplete data; in this situation, the maximum likelihood function is difficult or impossible
to use to estimate the parameter vector θ of the considered model. We formalize directly an
iteration from which we can understand clearly how this algorithm works:

- Let $y = (y_1, \dots, y_n)$ be the independent and identically distributed (i.i.d.) observations of likelihood $p(y|\theta)$.

- The maximization of $\log p(y|\theta)$ is impossible.

- We consider hidden data $z = (z_1, \dots, z_n)$ which make the maximization of the likelihood of the complete data possible when known.
- As we do not know these data z, we estimate the likelihood of the complete data by taking into account all the known information so the estimator is given as follows (E step):

$$Q(\theta|\theta_{k-1}) = \mathbb{E}[\log p(y, z; \theta|y; \theta_{k-1})], \quad (1)$$

where $\theta_{k-1}$ is the vector of the parameters at iteration $k-1$.

- Finally, we maximize this estimated likelihood to determine the new value of the parameter (M step). Thus, the transition from iteration $k-1$ to iteration $k$ in the algorithm consists in determining the parameters vector at iteration $k$, $\theta_k$:

$$\theta_k = \arg\max_{\theta} Q(\theta|\theta_{k-1}).$$

where $\theta_0$ is chosen arbitrarily. When one of these two steps are impossible to complete, we can consider an extension form of the EM algorithm such as the SAEM, SEM or MCEM algorithms. In the next subsection, we derive the SAEM algorithm into the formula of the EM algorithm.

### 2.1.1. SAEM Algorithm

The stochastic approximation expectation maximization (SAEM) algorithm was proposed by [26], in which the E step was replaced by a stochastic approximation. The stochastic approximation algorithm was first introduced by [28] and also used by [29], where the algorithm generates iterates of the form:

$$\theta_k = \theta_{k-1} - \gamma_k [h(\theta_{k-1}) + \epsilon_k],$$

where $\gamma_k$ is a sequence of positive step sizes, $h$ is a function of $\theta$, and $\epsilon_k$ is a constant such that $\epsilon_k = \mathbb{E}[h(\theta_k)]$. From this form, the SAEM algorithm is obtained where the E step of the EM algorithm is divided into two steps; consider the iteration $k$:

First, we sample a realization $z_k$ of the latent variable from the conditional distribution ($p_{\theta_{k-1}}(z|y)$) of $z$ given $y$, using the value of the parameter $\theta_{k-1}$ at iteration $k-1$.

Second, by using the realization $z_k$ from the first step, we update the value of $Q_k(\theta|\theta_{k-1})$ (see, (1)) through a stochastic approximation procedure.

Then, the algorithm continues as follows:

**Initialization step:** Initialize $\theta_0$ in a fixed compact set.

Then, for all $k \ge 1$, the $k$th iteration consists in three steps:

**Simulation step:** simulate $z_k$ from the conditional distribution $p_{\theta_{k-1}}(\cdot|y)$.

**Stochastic approximation step:** compute the quantity

$$Q_k(\theta) = Q_{k-1}(\theta) + \gamma_k \left( \frac{1}{m(k)} \sum_{j=1}^{m(k)} \log f(z_k(j); \theta) - Q_{k-1}(\theta) \right),$$

where $m_k$ is the number of simulations at each iteration.

**Maximization step:** update the parameter value according to $\theta_k = \arg\max_{\theta} Q_k(\theta)$.

The SAEM algorithm is more efficient compared to the MCEM, where, at each iteration, the simulation of the missing values is repeated and the ones obtained previously are not used. In the SAEM, all the simulated missing values contribute to the expectation step and that is the advantage of this algorithm where the maximization step is cheaper than the simulation one. Next, we consider the SEM algorithm to show how it differs from an SAEM.
### 2.1.2. SEM Algorithm

The SEM approach is a specific case of the SAEM algorithm. In the SEM algorithm, the step size is set to zero $γ_k = 0$ and the number of simulations by iteration is constant (usually, $m_k = 1$).

This method is an extension of the EM algorithm and is a solution when the maximum likelihood is impossible to complete. This algorithm was introduced by [30], where between the steps E and M, they used a stochastic step that consisted in generating a complete sample containing latent variables from the conditional density, based on the observed data. Starting with an initial position $θ_0$, an iteration k of the SEM algorithm, where, to each $θ_k$, a $θ_{k+1}$ is associated, is defined as follows:

**Step E:** Compute the conditional density $f(y|x; θ_k)$;

**Step S:** Draw $z_k$ from the conditional distribution, then obtain the complete sample $y_k = (x, z_k)$;

**Step M:** Update the parameters $θ_{k+1}$ by maximizing the likelihood function based on the complete vector $y_k$.

In the next subsection, we consider the MCEM algorithm to show how a Monte Carlo simulation is used in the EM algorithm.

### 2.1.3. MCEM Algorithm

In 1990, ref. [31] proposed to replace the expectation step to compute $Q(θ|θ_{k-1})$ by a Monte Carlo integration, hence the name MCEM algorithm (Monte Carlo EM). At iteration k, the E step is replaced by the following procedure:

**Simulation step (S-step):** generate $m(k)$ realizations $z_k(j)$ (with $j = 1, \dots, m(k)$)) of the missing data vector under the distribution function $p(z; θ_k)$

**Monte Carlo integration:** compute the approximation of $Q(θ|θ_{k-1})$ according to:

$$Q_k(\theta) = \frac{1}{m(k)} \sum_{j=1}^{m(k)} \log f(z_k(j); \theta).$$

The **maximization step** remains unchanged.

As said before, these algorithms can be used in such a way as to handle the unobserved values of the data that are defined in this next subsection.

## 2.2. Incomplete Data

We can talk about incomplete data when the values in our response vector are not all observed for many reasons; we can also say it is a question without an answer.

These types of data are a common problem recognized by statisticians. Ref. [32] presented many solutions to handle this problem based on simple and multiple imputations. Ref. [33] proposed to handle the incomplete data in a hierarchical model by using the SEM algorithm.

In this work, we are interested in data that contain two types of incomplete values: censored data denoted by $Y_c$ and missed data $Y_{miss}$. In the next subsection, the censored data are considered.

### 2.2.1. Censoring

Censored data are any data for which we do not know the exact event time. There are two types of censored data: right-censored and left-censored.

Right-censored data: these are data where the event has not yet been achieved when the study is finished.

Left-censored data: these are data when the event is achieved before the study starts.

We consider in this paper both censoring types of (left- and right-censoring) data with four percentages (0%, 5%, 10% and 20%). Let t be the censoring level and T the set of censored indices that can be written as $T = \{(i_1, \dots, i_{p+k}) \in I, Y_{i_1, \dots, i_{p+k}} < t \text{ or } Y_{i_1, \dots, i_{p+k}} > t\}$
$t$}, and let $Y = (Y_o, Y_c)$, with $Y_o = \{Y_{i_1, \dots, i_{p+k}, j'} (i_1, \dots, i_{p+k}) \notin T\}$ the vector of observed response variables and $Y_c = \{Y_{i_1, \dots, i_{p+k}, j'} (i_1, \dots, i_{p+k}) \in T\}$ defining the vector of censored response variables. The second type is the missing data, defined in the next subsection.

### 2.2.2. Missing Data

In general, missing values are produced when a value in the data is not represented for a given variable, for many reasons that can be linked to the objective of the study (for example, participants do not answer the questions). This appears in many research studies, particularly when collecting data and where participants are studied over a period of time.

At first, studies were developed assuming no missing values. In the late 1980s, with the advancement in technology, that problem attracted the attention of many researchers who wished to study several techniques for handling it. Depending on the reasons for their absence, these values could be divided into three types: missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR).

In Appendix A.2, we define the types of missing value and the methods to generate these types in a simulation study.

In this paper, we are interested in imputing the MAR type of missing data with four percentages (0%, 5%, 10% and 20%), which are crossed with four other percentages of censoring, resulting in 16 cases of incomplete data. For the missing data, we denote by $Y_o$ the vector of observed response variable and $Y_{miss}$ the missed one, so $Y = (Y_o, Y_{miss})$. In the next section, we present the main results of this work.

## 3. Main Results

### 3.1. The Proposed Model

Motivated by psychological data, in which we typically have two random variables, participants designed by $S$ and items designed by $I$, we consider a model that contains $p$ covariates $\{x_1; x_2; \dots; x_p\}$ and two random variables. The predicted variable, also named the variable of interest, is denoted by $Y$. Then, we consider the following linear mixed model:

$$Y = \beta_{0,s,i} + \beta_{1,s,i}x_1 + \beta_{2,s,i}x_2 + \dots + \beta_{p,s,i}x_p + r_{j,i,s} \quad (2)$$

with $r$ the residual of the model that follows a normal distribution. We give:

$$\left\{ \begin{array}{l} \beta_{0,s,i} = \delta_0 + S_{0,s} + I_{0,i} \\ \beta_{1,s,i} = \delta_1 + S_{1,s} + I_{1,i} \\ \beta_{2,s,i} = \delta_2 + S_{2,s} + I_{2,i} \\ \vdots \\ \beta_{p,s,i} = \delta_p + S_{p,s} + I_{p,i}. \end{array} \right.$$

By a replacement of these considerations in Equation (2), our model can be rewritten as follows:

$$\begin{aligned} Y = & \delta_0 + \delta_1 x_1 + \delta_2 x_2 + \dots + \delta_p x_p \\ & + S_{0,s} + S_{1,s} x_1 + S_{2,s} x_2 + \dots + S_{p,s} x_p \\ & + I_{0,i} + I_{1,i} x_1 + I_{2,i} x_2 + \dots + I_{p,i} x_p \\ & + r. \end{aligned}$$

By taking

$$X = [1 \ x_1 \ x_2 \ \dots \ x_{p+k}] = [1 \ x_j], \delta = (\delta_0 \ \delta_j)^T$$

and denoting by $T$ the transpose matrix,

$$Z = [z_j] = [x_j], u = [u_j] = [S_{j,l,s} + I_{j,l,i}]^T;$$
and $e_{jl} = [e_j]$.

This leads us to rewrite our model as follows:

$$Y = X\delta + \sum_{k=1}^{K} z_k u_k + e = X\delta + Zu + e,$$

where $E(Y) = X\delta$ and $Var(Y) = V = ZGZ^T + R$, with $u \sim N(0, G)$, $e \sim N(0, R)$ and $Cov(u, e^T) = 0$.

Now, by including all the interactions, the general model can be written as:

$$Y = F + I_{FF} + R + I_{FR} + e,$$

where we take:

$$F = \underbrace{\delta_0 + \delta_1 x_1 + \dots + \delta_p x_p}_{\text{fixed variables}},$$
$$I_{FF} = \underbrace{\delta_{p+1} x_{p+1} + \dots + \delta_{p+k} x_{p+k}}_{\text{interaction between fixed-fixed variables}},$$
$$R = \underbrace{(S_{1,s} x_1 + I_{1,i} x_i) + \dots + (S_{p,s} x_p + I_{p,i} x_i)}_{\text{random variables}},$$
$$I_{FR} = \underbrace{(S_{p+1,s} x_{p+1} + I_{p+1,i} x_{p+1}) + \dots + (S_{p+k,s} x_{p+k} + I_{p+k,i} x_{p+k})}_{\text{interaction between fixed-random variables}},$$
$$\text{and } e = \underbrace{(S_{0,s} + I_{0,i}) + r}_{\text{residuals}}.$$

We set $k$, the number of interactions, as $k = 2^p - p - 1$ and let $j, l \in [1, p + k]$.

Our random participants variable is given as follows:

$$S_{j,l,s} \sim N(0, \Sigma_S),$$

and the variance-covariance matrix of the participants is written as

$$\Sigma_S = \begin{cases} \sigma_{j,S}^2 & ; \text{if } j = l \\ \rho_S \sigma_{j,S} \sigma_{l,S} & ; \text{otherwise.} \end{cases}$$

We consider that the items random variable also follows a normal distribution:

$$I_{j,l,I} \sim N(0, \Sigma_I);$$

then, the variance-covariance matrix of the items is given by:

$$\Sigma_I = \begin{cases} \sigma_{j,I}^2 & ; \text{if } j = l \\ \rho_I \sigma_{j,I} \sigma_{l,I} & ; \text{otherwise.} \end{cases}$$

For $j, l \in [1, p + k]$, the variance-covariance matrix $G$ is equal to:

$$G = \Sigma_u = \begin{cases} \sigma_{j,S}^2 + \sigma_{j,I}^2 & ; \text{if } j = l \\ \rho_S \sigma_{j,S} \sigma_{l,S} + \rho_I \sigma_{j,I} \sigma_{l,I} & ; \text{otherwise.} \end{cases}$$
The goal is to predict our model by determining the parameters δ and u; we propose to use Henderson’s linear system ([34]) based on the maximum likelihood function. This approach is presented in Appendix A.3. The matrix R has the form:

$$R = \begin{cases} I & ; \text{if } j = l \le n - p - k \\ \sigma_{e_j}^2 & ; \text{if } j = l > n - p - k \\ 0 & ; \text{otherwise,} \end{cases}$$

where n represents the total number of observations. Therefore, with Z = (z₁, ..., zₚ, k), the variance of Y is equal to:

$$V = \left[ \begin{array}{l} \left( \sum_{j=1}^{p+k} z_j \Sigma_{u_1, u_j} \right) z_1 + \left( \sum_{j=1}^{p+k} z_j \Sigma_{u_2, u_j} \right) z_2 + \dots + \left( \sum_{j=1}^{p+k} z_j \Sigma_{u_p, u_j} \right) z_p + \dots \\ + \left( \sum_{j=1}^{p+k} z_j \Sigma_{u_{p+k}, u_j} \right) z_{p+k} \end{array} \right] + R.$$

The linear system equations of Henderson applied in this study is also presented in Appendix A.4. In the next subsection, we present some specific cases.

## 3.2. Specific Cases

In this section, we show how the model can be simplified if we take the following two cases: first, only the fixed-fixed interaction part is presented and second, there is no interaction.

### 3.2.1. Case 1: Fixed-Fixed Interaction

In this case, we consider the interactions between the fixed variables. Then, in the random part, the u vector is equal to 0:

$$ \begin{aligned} u &= \begin{bmatrix} u_1 & u_2 & \dots & u_p & u_{p+1} & \dots & u_{p+k} \end{bmatrix}^T \\ &= \begin{bmatrix} 0 & 0 & \dots & 0 & 0 & \dots & 0 \end{bmatrix}^T; \end{aligned} $$

therefore, the V matrix is equal to:

$$V = R = \begin{cases} I & ; \text{if } j = l \le n - p - k \\ \sigma_{e_j}^2 & ; \text{if } j = l > n - p - k \\ 0 & ; \text{otherwise.} \end{cases}$$

For more mathematical development, see Appendix A.4.

### 3.2.2. Case 2: No Interactions

In this case, there are no interactions between the variables; our model is simplified, where the random part with the fixed-fixed interaction is ignored:

$$X = [1 \ x_j]; \delta = (\delta_0 \ \delta_j)^T; Z = [z_j]; u = [u_j] = [0 \ 0 \ \dots \ 0]^T; e_{jl} = [e_j].$$

The V matrix is reduced to:

$$V = R = \begin{cases} I & ; \text{if } j = l \le n - p \\ \sigma_{e_j}^2 & ; \text{if } j = l > n - p \\ 0 & ; \text{otherwise.} \end{cases}$$

See the simplified system in Appendix A.4. Next, we present how the SEM algorithm can be applied in this work.
### 3.3. SEM Algorithm

In this paper, we propose to handle missing data and censored problems in the presence of interactions by using the SEM algorithm, proposed by [33]. The stochastic expectation maximization (SEM) algorithm is a method used to estimate the parameters when it is complicated to use the EM algorithm; it is a particular case of multiple imputations (MI).

To understand MI ([35]), we need to know the idea behind a simple imputation (SI) in which the nonobserved data {$Y_{no}$} (missed or censored) are replaced by a value and then the parameters are estimated using known methods to maximize the likelihood function. To obtain a robust estimate, the simple imputation can be repeated several times with different values of {$Y_{no}$} and the results can be combined. This method is called multiple imputations.

We applied the algorithm after crossing the MAR values with the censored ones (see Algorithm 1). In the SEM algorithm ([30]), the values of {$Y_{no}$} are drawn from the conditional distribution of the nonobserved data given the observed ones using the current values of the parameters. We generated samples from $P_{\Theta}(Y_{no}|Y_o)$ where this distribution was calculated using Gibbs's sampling; the procedure is underlined in Algorithm 1 (the SEM method using Gibbs's sampling was developed in [33]). For more information about the expectation maximization algorithm, see [36]. In the next subsection, we present the extension form of the SEM that leads us to the SAEM algorithm.

**Algorithm 1 SEM algorithm:** N is the number of iterations of the SEM algorithm, M is the burn-in level, $Y_1, \dots, Y_n$ is the response vector, G is the number of iterations of Gibbs sampling, $Y[u_i]$ is the response vector with respect to the ith random variable, and fixed is the summation of the fixed effects with the fixed interaction part.

**Input:** N, M, $Y_1, \dots, Y_n$, and G.

1: Random initialization of $\Theta$

2: **for** j = 1, ..., N **do**

3: **for** g = 1, ..., G **do**

4: draw $e^{(g)}$ from $\mathcal{N}(0, \sigma_e^2)$, lower = $\min(e) - \text{fixed} - u^{(g-1)}$, upper = $\max(e) - \text{fixed} - u^{(g-1)}$)

5: draw $u_1^{(g)}$ from $\mathcal{N}(0, \Sigma_{u_1})$, lower = $\max(\min(Y[u_1]) - u_2^{(g-1)} - e^{(g)})$, upper = $\min(\max(Y[u_1]) - u_2^{(g-1)} - e^{(g)})$)

6: draw $u_2^{(g)}$ from $\mathcal{N}(0, \Sigma_{u_2})$, lower = $\max(\min(Y[u_2]) - u_1^{(g)} - e^{(g)})$, upper = $\min(\max(Y[u_2]) - u_1^{(g)} - e^{(g)})$)

7: **end for**

**output:** $\mathbf{Y}_{no}$ obtained from sampled $(u^{(G)}, e^{(G)})$

8: $\Theta_{j+1} = \operatorname*{argmax}_{\Theta} L(\Theta|\mathbf{Y}_o, \mathbf{Y}_{no})$

9: **end for**

10: $\hat{\Theta} = (N - M)^{-1} \sum_{j=M+1}^{N} \Theta_j$

**output:** $\hat{\Theta}$

### 3.4. SAEM Algorithm

In the SEM algorithm, $\Theta_{j+1}$ (the parameter at state j + 1) depends only on $\Theta_j$ and $Y_j$. By taking T the operator of the EM algorithm and M, the operator which associates $\Theta_{j+1}$ with $Y_j$ via step M, the updated parameter $\Theta_{j+1}$ can be written as follows (see [30]):

$$ \Theta_{j+1} = T(\Theta_j) + V(\Theta_j, Y_j), $$

where $V(\Theta_j, Y_j) = M(Y_j) - T(\Theta_j)$.
The SAEM algorithm (see also Algorithm 2) is an extension form of the SEM algorithm. Starting with an initial position $Θ_0$, the iteration of the SAEM algorithm, to which for each $Θ_j$, a $Θ_{j+1}$ is associated, is defined by the equation:

$$Θ_{j+1} = T(Θ_j) + γ_j V(Θ_j, Y_j),$$

where $γ_j$ is the step size, which is a decreasing sequence of positive real numbers starting with $γ_0 = 1$, $γ_j \mapsto 0$ when $j \mapsto \infty$. The theoretical and practical performance values of the SAEM algorithm are strongly dependent on the speed of convergence towards 0 and the regular decreasing of the sequence $γ_j$, which explains the importance of choosing the step size.

Ref. [37] performed many tests to see which sequence was more efficient, and they obtained two types of sequences that were significant, one of a slow mode:

$$γ_j = \cos(\frac{j \times \pi}{2N}),$$

the other type was of a linear mode:

$$γ_j = -\frac{(j+1)}{N}.$$

Ref. [38] chose the decreasing sequence from their previous experiments [39]; they took that sequence as

$$γ_j = \cos(N \times \alpha) 1_{\{0 \le j \le 20\}} + \frac{c}{\sqrt{j}} 1_{\{j > 20\}}, \text{ where } \cos(20 \times \alpha) = \frac{c}{\sqrt{20}} = 0.3.$$

That sequence was also chosen by [40,41] in a nonlinear mixed model as 1 at the first iteration and $(j - N)^{-1}$ during the last iterations.

**Algorithm 2** SAEM algorithm: N is the number of iterations of the SEM algorithm, M is the burn-in level, $Y_1, \dots, Y_n$ is the response vector, G is the number of iterations of Gibbs's sampling, $γ_j$ is the decreasing sequence, $Y[u_i]$ is the response vector with respect to the ith random variable, and fixed is the summation of the fixed effects with the fixed interaction part.

**Input:** N, M, $Y_1, \dots, Y_n$, and G.

1: Random initialization of $Θ$

2: **for** j = 1, ..., N **do**

3: **for** g = 1, ..., G **do**

4: draw $e^{(g)}$ from $N(0, \sigma_e^2)$, lower = $\min(e) - \text{fixed} - u^{(g-1)}$, upper = $\max(e) - \text{fixed} - u^{(g-1)}$)

5: draw $u_1^{(g)}$ from $N(0, \Sigma_{u_1})$, lower = $\max(\min(Y[u_1]) - u_2^{(g-1)} - e^{(g)}$, upper = $\min(\max(Y[u_1]) - u_2^{(g-1)} - e^{(g)})$)

6: draw $u_2^{(g)}$ from $N(0, \Sigma_{u_2})$, lower = $\max(\min(Y[u_2]) - u_1^{(g)} - e^{(g)}$, upper = $\min(\max(Y[u_2]) - u_1^{(g)} - e^{(g)})$)

7: **end for**

**output:** $Y_{no}$ obtained from sampled $(u^{(G)}, e^{(G)})$

8: $Θ_{j+1} = \operatorname{argmax}_{Θ} γ_j \times L(Θ|Y_0, Y_{no})$

9: **end for**

10: $\hat{Θ} = (N - M)^{-1} \sum_{j=M+1}^{N} Θ_j$

**output:** $\hat{Θ}$

When applying these two algorithms, we were confronted with a problem of convergence caused by the presence of the interaction between the fixed and random variables.
Therefore, in the next section, we present how to handle this problem by using the Hamiltonian Monte Carlo algorithm.

## 4. Convergence of Parameters

While introducing the fixed-random effect in the application section, we faced a convergence problem in the SEM and SAEM algorithms when using the standard form of Monte Carlo (MC). In order to solve this problem, we propose to use the hybrid MC algorithm ([42]), which uses the symmetric Metropolis-Hastings algorithm to accept or reject a proposal based on the Hamiltonian function, hence the name of the algorithm, the Hamiltonian Monte Carlo algorithm (HMC). In the next subsection, we present how to implement the HMC algorithm.

### Implementation

The Hamiltonian function $H(p, q)$ is written in term of the joint density $\pi(p, q)$:

$$H(p, q) = -\log \pi(p, q),$$

which decomposes into two terms:

$$\begin{aligned} H(p, q) &= -\log \pi(p|q) - \log \pi(q) \\ &= K(p, q) + V(p, q). \end{aligned}$$

The first term corresponds to the density of the target distribution, the second one is determined by the target distribution when $K(p, q)$ is unconstrained and must be specified by the implementation.

The acceptance probability of moving from state $(p_0, q_0)$ to $(p_i, q_i)$ is determined using the Metropolis-Hastings algorithm ([43,44]), with $Q(p_i, q_i|p_0, q_0)$ being the probability density defining each proposal:

$$a(p_i, q_i|p_0, q_0) = \min \left( 1, \frac{Q(p_0, q_0|p_i, q_i) \pi(p_i, q_i)}{Q(p_i, q_i|p_0, q_0) \pi(p_0, q_0)} \right).$$

Referring to the symmetric Metropolis ([45]), the acceptance probability is simplified to the form:

$$\begin{aligned} a(p_i, q_i|p_0, q_0) &= \min \left( 1, \frac{\pi(p_i, q_i)}{\pi(p_0, q_0)} \right) \\ &= \min \left( 1, \frac{\exp(-H(p_i, q_i))}{\exp(-H(p_0, q_0))} \right) \\ &= \min(1, \exp(-H(p_i, q_i) + H(p_0, q_0))). \end{aligned}$$

In this study, we considered that we had evidence of the convergence towards stationarity due to the Markov Chain (see, [45,46]).

In the next section, we present some numerical applied results based on a psychological simulation and then real data.

## 5. Numerical Experiment

This section aims to compare the SEM algorithm proposed in this paper with other methods in the presence of incomplete data in a complete LMEM. The method was first applied to simulated data, then to real ones.

### 5.1. Simulated Data

The simulation data were created from an experimentally obtained database. In order to explore the cognitive process of retrieval in memory of the spelling of French words, 30 participants had to handwrite the label of 150 drawings of objects, constituting a database
of 4500 RTs ([4]). After removing the RTs corresponding to technical errors, spelling errors, and censoring data greater than 2.5 standard deviations from the mean, 3434 values remained. Ref. [4] performed a linear regression analysis involving nine fixed-effects variables. In the present work, we retained two of them. The first was an estimate of the age at which the image name was acquired in childhood (age of acquisition, AoA hereafter). This factor is one of the most important predictors of picture-naming RTs ([47]. The second fixed-effects variable corresponded to the number of letters in a picture label (Lett_cat). We performed a median split to obtain a categorical variable. Thus, image labels with a number of letters lower than the sample average were considered as short words and those with a higher number of letters were considered as long words. The interaction between these two fixed-effects variables was included in the model. Finally, there were two random effect variables: items, i.e., the 150 image labels produced by the 30 participants, and the participants themselves. The fixed-random interaction between AoA and the participants was also included in the LMEM.

Therefore, the model could be written as follows:

$$RTs = b_0 + b_1 \times AoA + b_2 \times Lett_cat + b_3 \times AoA \times Lett_cat + \text{participants} + \text{items} + AoA \times \text{participants} + e.$$

In this section, we illustrate the proposed algorithm by comparing the results to the SAEM and MCMC algorithms using different proportions of incomplete data. The comparison was performed by computing the mean absolute error (MAE), which measured the errors between the reference and the predicted vectors (with $\Theta$ as the reference value and $\hat{\Theta}$ as the predicted one); the standard equation of the MAE is given by:

$$MAE = \frac{1}{n} \sum_{i=1}^{n} |\hat{\Theta}_{i} - \Theta_{i}|.$$

We also computed the linear correlation between the parameters ($LCor = \frac{\text{Cov}(\hat{\Theta}_{i}, \Theta_{i})}{\sigma(\hat{\Theta}_{i})\sigma(\Theta_{i})}$) and the Spearman rank correlation ($RCor = 1 - \frac{6\sum d^2}{n(n^2-1)}$, where $d^2$ is the square of the difference in the ranks of the two coordinates).

In this numerical section, we provide the 16 cases to see the performance of the methods in the absence of one of these two considered types of unobserved values (MAR type/censored data) or in the presence of a low or high percentage. In addition, we used the simulated data of 4500 observations after testing our proposed methods on a larger number of observations. We obtained the same results as those for 4500 values.

After the implementation of our algorithm, we checked the convergence by illustrating the values of the estimated parameters obtained at each iteration; this is presented in Figure 1.

In Figure 2, we plotted the normalized residuals of each case. We observed a similar dispersion between the standard simulated data and the SEM algorithm, while we noticed a small but visible scattering by implementing the SAEM algorithm, and finally, we observed that when we used the MCMC method, the residuals were dispersed far from the normal line.

Our results are visualized quantitatively in Table 1. Based on the MAE value, we can observe that the SEM gave the lowest values in all cases except for the case of 5% MAR and 20% censoring, where the SAEM (8.108) beat the SEM (8.118) by a very narrow margin. By observing the full vectors, the SEM gave the best parameter values in two cases: (0%, 10%) and (20%, 5%). In the other cases, we observed that with a high percentage of MAR (20%), $b_2$ and $b_3$ in the SAEM method gave the closest values while for the other parameters, it was the SEM algorithm. In all cases, we see that the MCMC was not a good choice either for the parameters values or for the other measurements. Moreover, we can see that for LCor (respectively, Rcor), the SEM and the SAEM algorithms gave approximately the same results 0.999 (respectively, 1) in all cases. This showed that these two methods had a small
difference except for the case of (20%, 5%) where the RCor of the SAEM algorithm decreased from 1 to 0.9285. The MCMC had an RCor ranging between a maximum value of 1 and a minimum value of 0.738; that maximum value was obtained with 5% of MAR and the absence of censored data.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Value of the parameter estimates at each the SEM iteration with a burn-in period M = 10, we separate the parameters according to their values in (a), we plot b0, b1, sig1 and sig2 and in (b), we plot b2, b3, sig3 and sig4.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Normalized residual presented; (a)-standard simulated data, (b)-SEM algorithm, (c)-SAEM algorithm, (d)-MCMC algorithm.</figcaption></figure>
<table><caption>**Table 1.** Comparison of the SEM algorithm with the SAEM and MCMC algorithms in the complete model with the different percentages of incomplete data where the number of iterations was 50. sig₁, sig₂, sig₃ and sig₄ represent the standard deviations of the participants, items, the IFR part and the error, respectively. The percentages on the left present the MAR, and on the right are the percentages of censored data.</caption><thead><tr><th rowspan="2">Θ̂<sub>i</sub></th><th colspan="3">(0%, 0%)</th><th colspan="3">(0%, 5%)</th><th colspan="3">(0%, 10%)</th><th colspan="3">(0%, 20%)</th></tr><tr><th>Reference</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>1045</td><td>1044</td><td>1045</td><td>1068</td><td>1044</td><td>1043</td><td>1010</td><td>1047</td><td>1047</td><td>1047</td><td>1047</td><td>822.8</td></tr><tr><td>b<sub>1</sub></td><td>93.16</td><td>93.54</td><td>88.51</td><td>48.82</td><td>93.91</td><td>89.37</td><td>41.25</td><td>93.51</td><td>88.84</td><td>41.25</td><td>88.84</td><td>77.96</td></tr><tr><td>b<sub>2</sub></td><td>-52.28</td><td>-51.30</td><td>-42.33</td><td>-36.75</td><td>-52.50</td><td>-47.30</td><td>51.49</td><td>-30.27</td><td>-51.66</td><td>-30.27</td><td>-51.66</td><td>18.26</td></tr><tr><td>b<sub>3</sub></td><td>2.80</td><td>3.17</td><td>-0.18</td><td>6.92</td><td>4.43</td><td>5.52</td><td>-14.44</td><td>-1.79</td><td>5.02</td><td>-14.44</td><td>5.02</td><td>-18.20</td></tr><tr><td>sig<sub>1</sub></td><td>14.57</td><td>14.65</td><td>14.43</td><td>5.29</td><td>14.65</td><td>12.54</td><td>12.56</td><td>14.70</td><td>12.17</td><td>12.56</td><td>12.17</td><td>11.30</td></tr><tr><td>sig<sub>2</sub></td><td>230.8</td><td>231.7</td><td>225.9</td><td>350.5</td><td>231.1</td><td>223.0</td><td>430.5</td><td>233.2</td><td>219.0</td><td>430.5</td><td>219.0</td><td>538.7</td></tr><tr><td>sig<sub>3</sub></td><td>9.79</td><td>9.81</td><td>9.84</td><td>2.80</td><td>9.80</td><td>8.13</td><td>8.61</td><td>9.81</td><td>9.12</td><td>8.61</td><td>9.12</td><td>8.94</td></tr><tr><td>sig<sub>4</sub></td><td>282.0</td><td>282.0</td><td>275.5</td><td>283.1</td><td>282.0</td><td>269.4</td><td>277.8</td><td>282.0</td><td>253.5</td><td>269.4</td><td>253.5</td><td>253.9</td></tr><tr><td>MAE</td><td>0</td><td>0.421</td><td>3.637</td><td>28.07</td><td>0.505</td><td>4.628</td><td>51.84</td><td>3.935</td><td>6.575</td><td>4.628</td><td>6.575</td><td>83.64</td></tr><tr><td>LCor</td><td>1</td><td>0.999</td><td>0.999</td><td>0.992</td><td>0.999</td><td>0.999</td><td>0.971</td><td>0.999</td><td>0.999</td><td>0.971</td><td>0.999</td><td>0.915</td></tr><tr><td>RCor</td><td>1</td><td>1</td><td>1</td><td>0.904</td><td>1</td><td>1</td><td>0.738</td><td>1</td><td>1</td><td>0.738</td><td>1</td><td>0.833</td></tr></tbody></table><table><thead><tr><th rowspan="2">Θ̂<sub>i</sub></th><th colspan="3">(5%, 0%)</th><th colspan="3">(5%, 5%)</th><th colspan="3">(5%, 10%)</th><th colspan="3">(5%, 20%)</th></tr><tr><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>1044</td><td>1038</td><td>1004</td><td>1069</td><td>1037</td><td>1019</td><td>1069</td><td>1036</td><td>963</td><td>1078</td><td>1038</td><td>780</td></tr><tr><td>b<sub>1</sub></td><td>93.30</td><td>94.33</td><td>93.79</td><td>87.57</td><td>90.04</td><td>52.82</td><td>88.07</td><td>90.90</td><td>45.29</td><td>86.52</td><td>90.75</td><td>80.88</td></tr><tr><td>b<sub>2</sub></td><td>-46.69</td><td>-53.57</td><td>-49.05</td><td>-49.96</td><td>-42.94</td><td>-28.23</td><td>-52.68</td><td>-47.94</td><td>55.58</td><td>-32.18</td><td>-50.21</td><td>26.67</td></tr><tr><td>b<sub>3</sub></td><td>1.04</td><td>4.58</td><td>6.55</td><td>8.27</td><td>1.27</td><td>8.63</td><td>9.70</td><td>6.75</td><td>-11.90</td><td>3.82</td><td>5.61</td><td>-15.79</td></tr><tr><td>sig<sub>1</sub></td><td>15.19</td><td>13.65</td><td>16.59</td><td>15.25</td><td>17.02</td><td>8.34</td><td>15.36</td><td>13.93</td><td>16.19</td><td>15.79</td><td>12.80</td><td>4.10</td></tr><tr><td>sig<sub>2</sub></td><td>231.0</td><td>223.7</td><td>213.6</td><td>229.8</td><td>218.7</td><td>332.4</td><td>228.9</td><td>216.0</td><td>412.7</td><td>232.4</td><td>211.9</td><td>518.8</td></tr><tr><td>sig<sub>3</sub></td><td>9.94</td><td>7.93</td><td>11.07</td><td>9.97</td><td>8.63</td><td>5.05</td><td>9.95</td><td>7.52</td><td>10.61</td><td>10.11</td><td>8.42</td><td>6.56</td></tr><tr><td>sig<sub>4</sub></td><td>282.4</td><td>281.4</td><td>347.6</td><td>282.6</td><td>275.0</td><td>349.4</td><td>282.6</td><td>268.7</td><td>342.2</td><td>282.6</td><td>253.2</td><td>317.6</td></tr><tr><td>MAE</td><td>1.146</td><td>2.646</td><td>16.78</td><td>4.970</td><td>5.516</td><td>34.43</td><td>5.089</td><td>6.287</td><td>62.03</td><td>8.118</td><td>8.108</td><td>88.97</td></tr><tr><td>LCor</td><td>0.999</td><td>0.999</td><td>0.996</td><td>0.999</td><td>0.999</td><td>0.991</td><td>0.999</td><td>0.999</td><td>0.971</td><td>0.999</td><td>0.999</td><td>0.911</td></tr><tr><td>RCor</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0.928</td><td>1</td><td>1</td><td>0.738</td><td>1</td><td>1</td><td>0.809</td></tr></tbody></table><table><thead><tr><th rowspan="2">Θ̂<sub>i</sub></th><th colspan="3">(10%, 0%)</th><th colspan="3">(10%, 5%)</th><th colspan="3">(10%, 10%)</th><th colspan="3">(10%, 20%)</th></tr><tr><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>1070</td><td>1037</td><td>981</td><td>1069</td><td>1035</td><td>999</td><td>1070</td><td>1034</td><td>938</td><td>1080</td><td>1037</td><td>759</td></tr><tr><td>b<sub>1</sub></td><td>86.59</td><td>94.05</td><td>90.52</td><td>87.58</td><td>90.27</td><td>50.902</td><td>88.10</td><td>90.80</td><td>46.43</td><td>86.29</td><td>90.33</td><td>81.31</td></tr><tr><td>b<sub>2</sub></td><td>-45.96</td><td>-50.70</td><td>-62.01</td><td>-45.00</td><td>-39.50</td><td>-44.14</td><td>-47.79</td><td>-44.57</td><td>45.21</td><td>-28.34</td><td>-48.01</td><td>23.38</td></tr><tr><td>b<sub>3</sub></td><td>6.79</td><td>5.26</td><td>12.25</td><td>6.69</td><td>1.46</td><td>14.45</td><td>8.13</td><td>6.78</td><td>-8.83</td><td>2.53</td><td>5.81</td><td>-16.26</td></tr><tr><td>sig<sub>1</sub></td><td>14.80</td><td>24.43</td><td>4.36</td><td>14.97</td><td>27.35</td><td>25.00</td><td>15.15</td><td>22.09</td><td>30.51</td><td>15.93</td><td>21.39</td><td>18.64</td></tr><tr><td>sig<sub>2</sub></td><td>228.7</td><td>220.1</td><td>205.1</td><td>230.0</td><td>215.3</td><td>323.1</td><td>229.1</td><td>212.5</td><td>401.1</td><td>232.5</td><td>208.7</td><td>506.5</td></tr><tr><td>sig<sub>3</sub></td><td>9.79</td><td>11.46</td><td>5.94</td><td>9.89</td><td>12.10</td><td>6.336</td><td>9.88</td><td>10.16</td><td>13.87</td><td>10.11</td><td>10.50</td><td>8.73</td></tr><tr><td>sig<sub>4</sub></td><td>282.4</td><td>279.6</td><td>398.0</td><td>282.5</td><td>273.2</td><td>389.6</td><td>282.5</td><td>267.2</td><td>379.8</td><td>282.5</td><td>251.8</td><td>354.3</td></tr><tr><td>MAE</td><td>5.608</td><td>4.699</td><td>30.12</td><td>5.387</td><td>8.279</td><td>40.15</td><td>5.424</td><td>8.191</td><td>68.85</td><td>8.747</td><td>9.636</td><td>93.20</td></tr><tr><td>LCor</td><td>0.999</td><td>0.999</td><td>0.990</td><td>0.999</td><td>0.999</td><td>0.987</td><td>0.999</td><td>0.999</td><td>0.969</td><td>0.999</td><td>0.999</td><td>0.908</td></tr><tr><td>RCor</td><td>1</td><td>1</td><td>0.904</td><td>1</td><td>1</td><td>0.976</td><td>1</td><td>1</td><td>0.833</td><td>1</td><td>1</td><td>0.833</td></tr></tbody></table>
<table><caption>**Table 1.** Cont.</caption><thead><tr><th rowspan="2">Θ̂<sub>i</sub></th><th colspan="3">(20%, 0%)</th><th colspan="3">(20%, 5%)</th><th colspan="3">(20%, 10%)</th><th colspan="3">(20%, 20%)</th></tr><tr><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th><th>SEM</th><th>SAEM</th><th>MCMC</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td><b>1045</b></td><td>1035</td><td>896</td><td><b>1045</b></td><td>1044</td><td>908</td><td><b>1045</b></td><td>1036</td><td>857</td><td><b>1051</b></td><td><b>1039</b></td><td>687</td></tr><tr><td>b<sub>1</sub></td><td><b>93.34</b></td><td>94.05</td><td>100.8</td><td><b>93.70</b></td><td>80.78</td><td>64.43</td><td><b>94.05</b></td><td>90.14</td><td>57.35</td><td><b>92.78</b></td><td>89.41</td><td>91.59</td></tr><tr><td>b<sub>2</sub></td><td>-35.09</td><td><b>-38.36</b></td><td>2.54</td><td><b>-34.54</b></td><td>-128.87</td><td>21.66</td><td>-36.11</td><td><b>-34.79</b></td><td>97.19</td><td>-19.91</td><td><b>-39.24</b></td><td>85.55</td></tr><tr><td>b<sub>3</sub></td><td>-1.82</td><td><b>1.90</b></td><td>-16.60</td><td><b>-1.55</b></td><td>39.61</td><td>-13.36</td><td>-0.36</td><td><b>4.00</b></td><td>-31.09</td><td>-4.73</td><td><b>3.67</b></td><td>-42.98</td></tr><tr><td>sig<sub>1</sub></td><td><b>12.12</b></td><td>26.88</td><td>29.49</td><td><b>12.44</b></td><td>11.70</td><td>39.89</td><td><b>12.62</b></td><td>25.58</td><td>33.25</td><td><b>13.62</b></td><td>26.53</td><td>22.49</td></tr><tr><td>sig<sub>2</sub></td><td><b>230.9</b></td><td>211.4</td><td>204.7</td><td><b>231.8</b></td><td>118.2</td><td>315.9</td><td><b>231.3</b></td><td>204.4</td><td>387.9</td><td><b>233.1</b></td><td>200.9</td><td>485.3</td></tr><tr><td>sig<sub>3</sub></td><td>8.84</td><td>11.33</td><td><b>8.95</b></td><td><b>8.98</b></td><td>5.15</td><td>11.48</td><td><b>9.03</b></td><td>10.91</td><td>10.93</td><td><b>9.41</b></td><td>11.76</td><td>6.227</td></tr><tr><td>sig<sub>4</sub></td><td><b>282.3</b></td><td>276.2</td><td>465.9</td><td><b>282.3</b></td><td>265.3</td><td>453.6</td><td><b>282.3</b></td><td>264.0</td><td>442.2</td><td><b>282.4</b></td><td>249.2</td><td>411.6</td></tr><tr><td>MAE</td><td><b>3.261</b></td><td>8.002</td><td>57.07</td><td><b>3.374</b></td><td>32.83</td><td>67.38</td><td><b>2.994</b></td><td>10.905</td><td>93.01</td><td><b>6.343</b></td><td>12.47</td><td>117.2</td></tr><tr><td>LCor</td><td><b>0.999</b></td><td><b>0.999</b></td><td>0.969</td><td><b>0.999</b></td><td>0.991</td><td>0.969</td><td><b>0.999</b></td><td><b>0.999</b></td><td>0.944</td><td><b>0.999</b></td><td>0.998</td><td>0.873</td></tr><tr><td>RCor</td><td><b>1</b></td><td><b>1</b></td><td>0.976</td><td><b>1</b></td><td>0.9285</td><td>0.928</td><td><b>1</b></td><td><b>1</b></td><td>0.761</td><td><b>1</b></td><td><b>1</b></td><td>0.833</td></tr></tbody></table>

We marked in bold for each case the closest value to the reference among the three considered methods (the reference vector here was the simulated one that did not contain any missed value).

## 5.2. Real Data

In this section, we applied the three algorithms to the database used to create the simulated data. Of the 4500 RTs recorded during the experiment (30 participants producing the names of 150 items), 115 were removed by right censoring (2.56% of MNAR). Ref. [4] also removed 951 values, creating 21.13% of missing MCAR/MAR data. By applying the three proposed algorithms with some missing data to solve the problem in the presence of fixed and fixed-random interactions, we obtained the results presented in Table 2. The results were compared with respect to the initial vector where the missed values were treated by keeping the cases empty (KE). Interestingly, this is how missing and censored data are handled in psycholinguistic studies.

<table><caption>**Table 2.** Comparison of the initial parameters in the KE method versus the parameters estimated after the application of the SEM, SAEM and MCMC algorithms. The initial parameters presented here are the ones obtained from the model that contained 1066 NA (NA represent the incomplete data; here, we had 115 censored data and 951 missing data.</caption><thead><tr><th>Θ̂</th><th>KE</th><th>SEM</th><th>SAEM</th><th>MCMC</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>852.5</td><td>875.7</td><td>892.334</td><td>1298</td></tr><tr><td>b<sub>1</sub></td><td>146.6</td><td>134.2</td><td>120.4</td><td>-161.4</td></tr><tr><td>b<sub>2</sub></td><td>158.6</td><td>126.1</td><td>79.63</td><td>-335.3</td></tr><tr><td>b<sub>3</sub></td><td>-58.12</td><td>-42.08</td><td>-26.72</td><td>146.6</td></tr><tr><td>sig<sub>1</sub></td><td>184.3</td><td>196.0</td><td>203.5</td><td>216.9</td></tr><tr><td>sig<sub>2</sub></td><td>133.0</td><td>127.8</td><td>94.87</td><td>221.05</td></tr><tr><td>sig<sub>3</sub></td><td>37.05</td><td>28.67</td><td>32.84</td><td>66.13</td></tr><tr><td>sig<sub>4</sub></td><td>279.5</td><td>302.5</td><td>262.9</td><td>516.6</td></tr></tbody></table>

## 6. Conclusions

In this study, we proposed an algorithm based on the stochastic expectation maximization (SEM) to estimate the parameters under a linear mixed-effects model (LMEM) that contained fixed-fixed and fixed-random interactions in the presence of different percentages of incomplete data.

The simulated and real data showed that the proposed SEM algorithm gave the best results compared to other competitors, where the bias of the parameters was smaller than the bias in the SAEM and MCMC algorithms.

The simulation results were obtained by using statistical software R (see Appendix A.5 for some source code).
We plan to extend this approach by considering a more complicated level of interaction such as random–random effects and the interactions between more than two variables. Another direction could be to consider an extension of the present work by considering a generalized model under the logistic regression by taking one for the NA values and zero for the observed ones.

**Author Contributions:** Conceptualization, C.P. and Y.S.; methodology, Y.S.; software, A.Z.; validation, A.Z., C.P. and Y.S.; formal analysis, A.Z.; investigation, C.P. and Y.S resources, C.P.; data curation, C.P.; writing—original draft preparation, A.Z.; writing—review and editing, Y.S.; visualization, C.P. and Y.S.; supervision, C.P. and Y.S.; project administration, C.P and Y.S.; funding acquisition, Y.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Data is unavailable due to privacy.

**Acknowledgments:** We would like to thank the Editor-in-Chief of Entropy, the associate Editor and two anonymous reviewers for their helpful comments, which helped us to focus on improving the original version of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A

In the following section, we present how [24] defined the models by including the maximum effects. It is interesting to show the types of missing data with the implementation methods. Then, the approach of Henderson is introduced and our model matrices to estimate the parameters are defined. Finally, we present some source code of the SEM algorithm.

### Appendix A.1. Models Extension

Ref. [24] compared five types of models in order to prove the importance of taking the maximum effects. They presented the two random variables subject, designed by S, and items, designed by I. By taking i as the items' index, s as the subject's index and c as the observation, the first model was the maximum model that contained the fixed and random intercepts with the interactions, including also ρ_S and ρ_I:

$$ 
\begin{aligned}
Y_{c,i,s} &= \beta_0 + S_{0,s} + I_{0,i} + (\beta_1 + S_{1,s} + I_{1,i})X_c + \epsilon_{c,i,s}, \\
\begin{pmatrix} S_{0,s} \\ S_{1,s} \end{pmatrix} &\sim \mathcal{N} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00}^2 & \rho_S \tau_{00} \tau_{11} \\ \rho_S \tau_{00} \tau_{11} & \tau_{11}^2 \end{pmatrix} \right), \\
\begin{pmatrix} I_{0,i} \\ I_{1,i} \end{pmatrix} &\sim \mathcal{N} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} w_{00}^2 & \rho_I w_{00} w_{11} \\ \rho_I w_{00} w_{11} & w_{11}^2 \end{pmatrix} \right), \\
\epsilon_{c,i,s} &\sim (0, \sigma^2).
\end{aligned}
$$

The second model differed by its correlation parameters, where S_{0,s} and S_{1,s} (respectively, I_{0,i} and I_{1,i}) were independent. ρ_S (respectively, ρ_I) was set to zero:

$$ 
\begin{aligned}
Y_{c,i,s} &= \beta_0 + S_{0,s} + I_{0,i} + (\beta_1 + S_{1,s} + I_{1,i})X_c + \epsilon_{c,i,s}, \quad (A1) \\
\begin{pmatrix} S_{0,s} \\ S_{1,s} \end{pmatrix} &\sim \mathcal{N} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00}^2 & 0 \\ 0 & \tau_{11}^2 \end{pmatrix} \right), \\
\begin{pmatrix} I_{0,i} \\ I_{1,i} \end{pmatrix} &\sim \mathcal{N} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} w_{00}^2 & 0 \\ 0 & w_{11}^2 \end{pmatrix} \right), \\
\epsilon_{c,i,s} &\sim (0, \sigma^2).
\end{aligned}
$$
In the third model, they reduced Equation (A1) by ignoring the items' random slopes (the interaction between the items' fixed and random effects), where $w_{11}$ was set to 0:

$$ 
\begin{aligned}
Y_{c,i,s} &= \beta_0 + S_{0,s} + I_{0,i} + (\beta_1 + S_{1,s})X_c + \epsilon_{c,i,s}, \\
\binom{S_{0,s}}{S_{1,s}} &\sim \mathcal{N}\left(\binom{0}{0}, \binom{\tau_{00}^2}{0} \binom{0}{\tau_{11}^2}\right), \\
I_{0,i} &\sim (0, w_{00}^2), \\
\epsilon_{c,i,s} &\sim (0, \sigma^2).
\end{aligned}
$$

The fourth model kept the random items' slope and excluded that of the subject, where $\tau_{11} = 0$:

$$ 
\begin{aligned}
Y_{c,i,s} &= \beta_0 + S_{0,s} + I_{0,i} + (\beta_1 + I_{1,i})X_c + \epsilon_{c,i,s}, \\
S_{0,s} &\sim (0, \tau_{00}^2), \\
\binom{I_{0,i}}{I_{1,i}} &\sim \mathcal{N}\left(\binom{0}{0}, \binom{w_{00}^2}{0} \binom{0}{w_{11}^2}\right), \\
\epsilon_{c,i,s} &\sim (0, \sigma^2).
\end{aligned}
$$

The fifth model excluded the two random slopes (for items and subject), where $\tau_{11} = w_{00} = 0$:

$$ 
\begin{aligned}
Y_{c,i,s} &= \beta_0 + S_{0,s} + I_{0,i} + \beta_1 X_c + \epsilon_{c,i,s}, \\
S_{0,s} &\sim (0, \tau_{00}^2), \\
I_{0,i} &\sim (0, w_{00}^2), \\
\epsilon_{c,i,s} &\sim (0, \sigma^2).
\end{aligned}
$$

In the next section, we present the different types of missing data and how to generate them using different methods.

## Appendix A.2. Missing Data Types and Imputations

Let $M$ be the missing data indicator that defines what is known and what is missing. The response variable is denoted by $Y$ in the complete data, where $Y = \{Y_o, Y_{no}\}$ with $Y_o$ is the observed part and $Y_{no}$ is the nonobserved one. Ref. [48] denoted the missing part by $Y_{miss}$. Consider a dataset of $p$ variables and $n$ subjects, with $Y = \{Y_1, \dots, Y_p\}$. $Y$ should look like a matrix of $n$ rows and $p$ columns. In the literature on missing data, see also [49,50].

We show in this section the three types of missing values and how to generate them in a simulation study based on a rule in each case ([51]).

We start with the first type named MCAR.

The missing completely at random type (MCAR) is a special case of the MAR type that is presented later. In this type, the value is missed due to chance and the absence is not related to the subject, so the distribution of $M$ does not depend on $Y_o$ or on $Y_{miss}$ and is identical for all the observations:

$$ \mathbb{P}(M|Y) = \mathbb{P}(M), $$

The statistical advantage of MCAR data is that the analysis remains unbiased, none of the variables is affected more than another.

Suppose that each subject has a probability $\pi$ of being missed from a variable, the missing data rule is given by $\mathbb{P}(M = 1) = \pi$. From this rule, we can determine various properties associated with the MCAR data such as the expected percentage of missing values and the expected number of missing data patterns in the MCAR data.
If we are working with univariate patterns, let $n$ be the number of subjects in the data and $K$ the random variable indicating the number of subjects with missing values in the data that follows a binomial distribution given by $K \sim B(n, \pi)$, where $0 \le \pi \le 1$. Since $\mathbb{E}(K) = n\pi$ and $\mathbb{V}(K) = n\pi(1 - \pi)$, the expected percentage of missing values is:

$$ \begin{aligned} \mathbb{E}(\Pi) &= \frac{1}{n}\mathbb{E}(K) \\ &= \pi, \end{aligned} $$

where $\Pi = \frac{K}{n}$ is the random variable that defines the estimated percentage of missing values in a sample. The variance of this estimated percentage is given by:

$$ \begin{aligned} \mathbb{V}(\Pi) &= \frac{1}{n^2}\mathbb{V}(K) \\ &= \frac{\pi(1 - \pi)}{n}. \end{aligned} $$

In this case, we have two missing data patterns, pattern 1 includes subjects with complete data; pattern 2 includes subjects with missing data. Let $I_j$ be the indicator variable of the event, where pattern $j$ is present in at least one subject in the sample with $j \in \{1, 2\}$. The probability that pattern 1 is present in at least one subject is $\mathbb{P}(I_1 = 1) = \mathbb{E}(I_1) = 1 - \pi^n$. The probability that pattern 2 is present in at least one subject is $\mathbb{P}(I_2 = 1) = \mathbb{E}(I_2) = 1 - (1 - \pi)^n$. Let $D$ be the number of distinct missing data patterns, $D = \sum_{j=1}^2 I_j$; we can determine the expected value:

$$ \begin{aligned} \mathbb{E}(D) &= \sum_{j=1}^{2} \mathbb{E}(I_j) \\ &= (1 - \pi^n) + (1 - (1 - \pi)^n) \\ &= 2 - \pi^n - (1 - \pi)^n. \end{aligned} \qquad (A2) $$

Therefore, in MCAR data that contain more than two missing data patterns, the expected number of distinct pattern in a sample is given by:

$$ \mathbb{E}(D) = m - \sum_{j=1}^{m} (1 - \eta_j)^n, $$

with $\eta_1, \dots, \eta_m$ are the corresponding probabilities for patterns $1, \dots, m$, and depends only on the probability of missing values in each variable $i$ (i.e., $\pi_1, \dots, \pi_l$). $m = 2^l$ is the number of patterns and $l$ is the total number of variables. Generally, the most used method to implement the MCAR type is randomly deleting the desired percentage of missing values.

The second type is MAR.

Missing at random is the most classic case; Rubin and D.B (1976) defined missing data to be MAR if the distribution of missing data did not depend on $Y_{miss}$:

$$ \mathbb{P}(M|Y) = \mathbb{P}(M|Y_o). $$

MAR data occur when the absence is not random but can be explained by the variables for which complete information exists. More generally, we are in the MAR type if the cause of missing values is not related to the goal of the study. If $Y_1$ has missing values, then it is regressed on other variables $Y_2$ to $Y_l$. The missing values in $Y_1$ are then replaced by the obtained predictive values. Similarly, if $Y_2$ has missing values, then the $Y_1, Y_3$ to $Y_p$ variables are used in the prediction model as independent variables. Thus, the missing values are replaced with the predicted ones. The missing data rules for MAR data can be organized into several categories: single-cutoff method, multiple-cutoff method, percentile method and logistic regression method.
Here, we define these three rules.

**Single-cutoff method:** Consider $Y_1$ as the variable with missing value and $Y_2$ is the missing data predictor with a cutoff point $a$. If a subject has $Y_2 \ge a$, then its probability of being missing from $Y_1$ is $\pi_1$, and if $Y_2 < a$, then its probability of being missing from $Y_1$ is $\pi_2$. Let $U$ be an indicator that takes the value 1 ($U = 1$) when $Y_2 \ge a$ and $U = 0$ when $Y < a$. The missing data rule is given as:

$$ \mathbb{P}(M = 1|U = 1) = \pi_1; \mathbb{P}(M = 1|U = 0) = \pi_2, $$

where the parameters associated with this rule are $\pi_1$ and $\pi_2$, and $\pi_0$ is the probability that $Y_2$ is equal to or greater than the cutoff: $\mathbb{P}(Y_2 \ge a) = \mathbb{P}(U = 1) = \pi_0$. To find the expected percentage of missing values, we first calculate the unconditional probability of a subject being missing from $Y_1$:

$$ \begin{aligned} \pi_{miss} &= \mathbb{P}(M = 1) \\ &= \mathbb{P}(M = 1|U = 1)\mathbb{P}(U = 1) + \mathbb{P}(M = 1|U = 0)\mathbb{P}(U = 0) \\ &= \pi_1\pi_0 + \pi_2(1 - \pi_0). \end{aligned} $$

The expected percentage of missing values in a sample is given by:

$$ \begin{aligned} \mathbb{E}(\Pi) &= \pi_{miss} \\ &= \pi_1\pi_0 + \pi_2(1 - \pi_0), \end{aligned} $$

and the variance of this estimated percentage of missing values is written as:

$$ \mathbb{V}(\Pi) = \frac{\pi_{miss}(1 - \pi_{miss})}{n}. $$

**Multiple-cutoff method:** One advantage of this method is that it can be used to create a nonlinear relationship between the missing data indicator and the missing data predictor. To create a nonlinear relationship between the indicator and the predictor, we need to fix an upper cutoff and a lower cutoff. Suppose we fix the two cutoff points $a$ and $-a$, $U = 1$ when $Y_2 \ge a$ or $Y_2 \le -a$, and $U = 0$ when $-a < Y_2 < a$. The missing data rule can be written as:

$$ \mathbb{P}(M = 1|U = 1) = \pi_1; \mathbb{P}(M = 1|U = 0) = \pi_2, \qquad (A3) $$

and we can see that this missing data rule is the same as the missing data rule in the single-cutoff method. In the case of a linear relationship, we need to specify two or more cutoff points in the missing data predictor. Suppose we take the quartile points ($Q_1$, $Q_2$ and $Q_3$) as a cutoff, let $V$ be a discrete uniform random variable created based on the values of $Y_2$ and that takes the following values:

$$ V = \begin{cases} 1 & \text{if } Y_2 < Q_1 \\ 2 & \text{if } Q_1 \le Y_2 < Q_2 \\ 3 & \text{if } Q_2 \le Y_2 < Q_3 \\ 4 & \text{if } Y_2 \ge Q_3. \end{cases} \qquad (A4) $$

Then, for the linear relation case, the missing data rule is given by:

$$ \mathbb{P}(M = 1|V = 1) = \pi_1; \mathbb{P}(M = 1|V = 2) = \pi_2; $$

$$ \mathbb{P}(M = 1|V = 3) = \pi_3; \mathbb{P}(M = 1|V = 4) = \pi_4. $$
Let $n$ be the total number of subjects and $n_j$ be the number of subjects in each quartile group. The variance of the estimated $\pi_0$ is:

$$V(\Pi_0) = \frac{\pi_0(1 - \pi_0)}{n},$$

and the variance of the estimated $\pi_j$ where $j \in \{1, 2, 3, 4\}$ is defined as follows:

$$V(\Pi_j) = \frac{\pi_j(1 - \pi_j)}{n_j}.$$

We can calculate each subject's probability of being missing by calculating the marginal probability of $M = 1$:

$$\begin{aligned} \pi_{miss} &= P(M = 1) \\ &= \pi_1\pi_0 + \pi_2\pi_0 + \pi_3\pi_0 + \pi_4\pi_0, \end{aligned}$$

and the expected percentage of missing values can be written as:

$$\begin{aligned} \mathbb{E}(\Pi) &= \pi_{miss} \\ &= \pi_1\pi_0 + \pi_2\pi_0 + \pi_3\pi_0 + \pi_4\pi_0; \end{aligned}$$

thus, the variance of this estimated percentage is given by:

$$V(\Pi) = \frac{\pi_{miss}(1 - \pi_{miss})}{n}.$$

To implement the missing data rule in (A4), researchers typically delete $\pi_1$ subjects with $Y_2 < Q_1$, $\pi_2$ subjects with $Q_1 \le Y_2 < Q_2$, and so on.

**Percentile method:** This case is an extension of the multiple-cutoff method, where each subject has a probability of being missing that depends on its percentile rank in the missing data predictor.

If there is a direct relation between the missing data indicator and the predictor, then the missing data rule is defined by:

$$P(M = 1 | Y_2 = q_k) = \frac{k}{100},$$

where $q_k$ is the $Y_2$ value corresponding to its $k$th percentile. If there is an indirect relationship, then the missing data rule is defined as:

$$P(M = 1 | Y_2 = q_k) = 1 - \frac{k}{100}.$$

**Logistic regression method:** When generating MAR data by using the logistic regression method, the logistic regression model is considered as the missing data rule and the population regression coefficients of the model are the parameters associated with the missing data rule. Then, the logistic regression model for subject $i$ is written as:

$$\log \left( \frac{P(M_i = 1 | y_{2,i})}{1 - P(M_i = 1 | y_{2,i})} \right) = \beta_0 + \beta_1 y_{2,i},$$

where $y_{2,i}$ is the subject $i$'s value on $Y_2$. The probability of being missing for each subject is given by:

$$P(M_i = 1 | y_{2,i}) = \frac{1}{1 + e^{-\beta_0 - \beta_1 y_{2,i}}}.$$

Because the above function is continuous, it means the probability of being missing for $Y_1$ gradually increases or decreases as the value of $Y_2$ increases. With the logistic regression,
there is no simple formula for calculating the expected percentage of missing data, so we estimate the expected percentage of missing data by calculating the mean of the probabilities in a sample with a large sample size:

$$ \pi_{miss} = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{1 + e^{-\beta_0 - \beta_1 Y_{2,i}}} $$

The third and last type is the MNAR data:

In missing not at random data, the cause of missing data is related to the variable of interest. In this type, the probability that an observation is missing depends on the observed $Y_0$ and the missing data $Y_{miss}$. If $Y_1$ is neither MCAR nor MAR, then it is MNAR. Generating MNAR data is the same as generating MAR data; the only difference is that in MNAR data the probability of missing values depends on the variable's own value and not on the observed values of the other variables. Therefore, we can change the missing data predictor to the variable with missing values, then we use one of the methods presented above to generate MAR data. For example, to change rule (A3) that generates MAR data to one that generates MNAR data, we only need to replace the variable $Y_2$ by $Y_1$; therefore, the corresponding missing data rule for generating MNAR data is when $Y_1$ value is above a cutoff point a. $Y_1$ has a $\pi_1$ probability of being missing, otherwise, $Y_1$ has a $\pi_2$ probability of being missing.

Next, we define the linear system of Henderson to deal with the parameters' estimation.

## Appendix A.3. Henderson's Approach

In order to estimate the parameters $\hat{\delta}$ and $\hat{u}$, we need to maximize the log function $(\hat{\delta}, \hat{u}) = \arg\max_{\delta,u} \ln f(Y, u)$. We define $f(Y, u) = f(Y|u)f(u)$ with $Y|u \sim N(X\delta + Zu, R)$ and $u \sim N(0, G)$, so we can write:

$$ \begin{aligned} -2 \ln f(Y|u) &= N \ln(2\pi) + \ln |R| + (Y - X\delta - Zu)' R^{-1} (Y - X\delta - Zu), \\ -2 \ln f(u) &= q \ln(2\pi) + \ln |G| + u' G^{-1} u. \end{aligned} $$

Denote $l(\delta, u; Y) = \ln f(Y, u)$; to maximize this function, we need to cancel the first derivatives of $\ln f(Y, u)$ with respect to $\delta$ and $u$:

$$ \begin{aligned} \frac{\partial[-2l(\delta, u; y)]}{\partial\delta} &= -2X'R^{-1}(y - X\delta - Zu) = 0, \\ \frac{\partial[-2l(\delta, u; y)]}{\partial u} &= -2Z'R^{-1}(y - X\delta - Zu) + 2G^{-1}u = 0. \end{aligned} \quad (A5) $$

From the two equations in (A5), Henderson presented the linear mixed model under the system:

$$ \begin{bmatrix} X'R^{-1}X & X'R^{-1}Z \\ Z'R^{-1}X & Z'R^{-1}Z + G^{-1} \end{bmatrix} \begin{bmatrix} \hat{\delta} \\ \hat{u} \end{bmatrix} = \begin{bmatrix} X'R^{-1}y \\ Z'R^{-1}y \end{bmatrix} \quad (A6) $$

In the next part, we present how we applied Henderson's approach to our proposed model.

## Appendix A.4. Model Simplification

**Complete model:** We consider the following matrices, where $r$ determines the rows, $c$ determines the columns, and $l, h \in [1, p+k]$:

$$ (X'R^{-1}X)_{M \times N} = \begin{cases} (n-p-k) + \sigma_{e_1}^{-2} + \dots + \sigma_{e_p}^{-2} + \sigma_{e_{p+1}}^{-2} + \dots + \sigma_{e_{p+k}}^{-2} & ; \text{if } r = c = 1 \\ \sum_{j=1}^{n-p-k} x_{l,j} + \sum_{o=1}^{p+k} x_{l,n-p-k+o} \sigma_{e_o}^{-2} & ; \text{if } r = 1; c \neq 1 \\ \sum_{j=1}^{n-p-k} x_{l,j} x_{h,j} + \sum_{o=1}^{p+k} x_{l,n-p-k+o} x_{h,n-p-k+o} \sigma_{e_o}^{-2} & ; \text{otherwise.} \end{cases} $$
$$ (X'R^{-1}y)_{M \times 1} = \begin{cases} \sum_{j=1}^{n-p-k} y_j + \sum_{o=1}^{p+k} y_{n-p-k+o} \sigma_{e_o}^{-2} & ; \text{ if } r = 1 \\ \sum_{j=1}^{n-p-k} x_{l,j} y_j + \sum_{o=1}^{p+k} x_{l,n-p-k+o} y_{n-p-k+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (X'R^{-1}Z)_{M \times N} = \begin{cases} \sum_{j=1}^{n-p-k} z_{l,j} + \sum_{o=1}^{p+k} z_{l,n-p-k+o} \sigma_{e_o}^{-2} & ; \text{ if } r = 1 \\ \sum_{j=1}^{n-p-k} x_{l,j} z_{h,j} + \sum_{o=1}^{p+k} x_{l,n-p-k+o} z_{h,n-p-k+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (Z'R^{-1}Z)_{M \times N} = \begin{cases} \sum_{j=1}^{n-p-k} z_{l,j}^2 + \sum_{o=1}^{p+k} z_{l,n-p-k+o} \sigma_{e_o}^{-2} & ; \text{ if } r = c \\ \sum_{j=1}^{n-p-k} z_{l,j} z_{h,j} + \sum_{o=1}^{p+k} z_{l,n-p-k+o} z_{h,n-p-k+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (Z'R^{-1}y)_{M \times 1} = \left[ \sum_{j=1}^{n-p-k} z_{l,j} y_j + \sum_{o=1}^{p+k} z_{l,n-p-k+o} y_{n-p-k+o} \sigma_{e_o}^{-2} \right]. $$

**Case of no interactions:** In this case, all the matrices of the linear system are simplified to reach *p* instead of *p* + *k*. $\forall l, h \in [1, p]$:

$$ (X'R^{-1}X)_{M \times N} = \begin{cases} (n-p) + \sigma_{e_1}^{-2} + \dots + \sigma_{e_p}^{-2} & ; \text{ if } r = c = 1 \\ \sum_{j=1}^{n-p} x_{l,j} + \sum_{o=1}^{p} x_{l,n-p+o} \sigma_{e_o}^{-2} & ; \text{ if } r = 1, c \neq 1 \\ \sum_{j=1}^{n-p} x_{l,j} x_{h,l} + \sum_{o=1}^{p} x_{l,n-p+o} x_{h,n-p+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (X'R^{-1}y)_{M \times 1} = \begin{cases} \sum_{j=1}^{n-p} y_j + \sum_{o=1}^{p} y_{n-p+o} \sigma_{e_o}^{-2} & ; \text{ if } r = 1 \\ \sum_{j=1}^{n-p} x_{l,j} y_j + \sum_{o=1}^{p} x_{l,n-p+o} y_{n-p+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (X'R^{-1}Z)_{M \times N} = \begin{cases} \sum_{j=1}^{n-p} z_{l,j} + \sum_{o=1}^{p} z_{l,n-p+o} \sigma_{e_o}^{-2} & ; \text{ if } r = 1 \\ \sum_{j=1}^{n-p} x_{l,j} z_{h,j} + \sum_{o=1}^{p} x_{l,n-p+o} z_{h,n-p+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (Z'R^{-1}Z)_{M \times N} = \begin{cases} \sum_{j=1}^{n-p} z_{l,j}^2 + \sum_{o=1}^{p} z_{l,n-p+o} \sigma_{e_o}^{-2} & ; \text{ if } r = c \\ \sum_{j=1}^{n-p} z_{l,j} z_{h,j} + \sum_{o=1}^{p} z_{l,n-p+o} z_{h,n-p+o} \sigma_{e_o}^{-2} & ; \text{ otherwise.} \end{cases} $$

$$ (Z'R^{-1}y)_{M \times 1} = \left[ \sum_{j=1}^{n-p} z_{l,j} y_j + \sum_{o=1}^{p} z_{l,n-p+o} y_{n-p+o} \sigma_{e_o}^{-2} \right]. $$

## Appendix A.5. R Source Code

Finally, we give the source code in R for the SEM algorithm by considering our real data.

```R
library(lme4); # for lmer
library(msm); # for rnorm
# initial model
model_initial=lmer(RT~AoA*Lett_cat+(1+AoA|participants)+(1|items), data,
  control = lmerControl(optimizer="bobyqa", calc.derivs = FALSE));

# random initialization of parameters
b0=fixef(model_initial)[1];
b1=fixef(model_initial)[2];
b2=fixef(model_initial)[3];
b3=fixef(model_initial)[4];
```
```R
upp=min(max(!is.na(data[participants == participants[a],5]))-Rand_items[Subset1]
-Rand_parts.AoA[Subset1]-Err[Subset1]-fixed[Subset1]);
upp=ifelse(upp-low>0,upp,upp+abs(upp-low)+100);
newRand_parts[Subset1]=rtnorm(1,sd=sig_parts, lower=low,upper=upp);
}
}
p(Rand_parts,newRand_parts);
# Rand_items with constraints
for (a in 1:lengthfactor(items)){
Subset2=(items == items[a]) & censored;
newRand_items=Rand_items;
if (sum(Subset2)>0){
low=max(min(!is.na(data[items == items[a],5]))-Rand_parts[Subset2]
-Rand_parts.AoA[Subset2]-Err[Subset2]-fixed[Subset2]);
upp=min(max(!is.na(data[items == items[a],5]))-Rand_parts[Subset2]
-Rand_parts.AoA[Subset2]-Err[Subset2]-fixed[Subset2]);
upp=ifelse(upp-low>0,upp,upp+abs(upp-low)+100);
newRand_items[Subset2]=rtnorm(1,sd=sig_items, lower=low,upper=upp);
}
}
p(Rand_items,newRand_items);
# Rand_parts.AoA with constraints
for (a in 1:lengthfactor(participants)){
Subset1=(participants == participants[a]) & censored;
newRand_parts.AoA=Rand_parts.AoA;
if (sum(Subset1)>0) {
low=max(min(!is.na(data[participants == participants[a],5]))-Rand_items[Subset1]
-Rand_parts[Subset1]-Err[Subset1]-fixed[Subset1]);
upp=min(max(!is.na(data[participants == participants[a],5]))-Rand_items[Subset1]
-Rand_parts[Subset1]-Err[Subset1]-fixed[Subset1]);
upp=ifelse(upp-low>0,upp,upp+abs(upp-low)+100);
newRand_parts.AoA[Subset1]=rtnorm(1,sd=sig_parts.AoA, lower=low,upper=upp);
}
}
p(Rand_parts.AoA,newRand_parts.AoA);
# end of Gibbs sampling algorithm
}
# compute the new response vector
RT[censored]=round(fixed[censored]+Rand_parts[censored]+Rand_items[censored]
+Rand_parts.AoA[censored]+Err[censored],0);
# train the model with the complete data
model_final<-lmer(RT~AoA*Lett_cat+(1+AoA|participants)+(1|items), data,
control = lmerControl(optimizer="bobyqa",calc.derivs = FALSE));
# end of SEM algorithm
}
```

## References

1. Donders, F.C. On the speed of mental processes. *Acta Psychol.* **1969**, *30*, 412–431. [CrossRef] [PubMed]

2. Roelofs, A. One hundred fifty years after Donders: Insights from unpublished data, a replication, and modeling of his reaction times. *Acta Psychol.* **2018**, *191*, 228–233. [CrossRef] [PubMed]

3. Baayen, R.H.; Davidson, D.J.; Bates, D.M. Mixed-effects modeling with crossed random effects for subjects and items. *J. Mem. Lang.* **2008**, *29*, 390–412. [CrossRef]

4. Perret, C.; Solier, C. Application of a Bayesian approach for exploring the impact of syllable frequency in handwritten picture naming. *J. Cogn. Psychol.* **2022**, *34*, 622–635. [CrossRef]
5. Ratcliff, R. Methods for dealing with reaction time outliers. *Psychol. Bull.* **1993**, *114*, 510–532. [CrossRef]

6. Van Selst, M.; Jolicoeur, P. A solution to the effect of the size on outlier elimination. *Q. J. Exp. Psychol.* **1994**, *47*, 631–650. [CrossRef]

7. Barr, D.J.; Levy, R.; Scheepers, C.; Tily, H.J. Random effects structure for confirmatory hypothesis testing: Keep it maximal. *J. Mem. Lang.* **2013**, *68*, 255–278. [CrossRef]

8. Dempster, A.P.; Rubin, D.B.; Tsutakawa, R.K. Estimation in Covariance Components Models. *J. Amer. Statist. Assoc.* **1981**, *76*, 341–353. [CrossRef]

9. West, T.B.; Welch, B.K.; Galecki, T.A. *Linear Mixed Models, A Practical Guide Using Statistical Software*; Chapman & Hall: New York, NY, USA, 2006.

10. Jiang, J.; Wand, M.P.; Bhaskaran, A. Usable and precise asymptotics for generalized linear mixed model analysis and design. *J. R. Stat. Soc. Ser. B* **2022**, *84*, 55–82. [CrossRef]

11. Faraway, J.J. *Extending the Linear Model with R*, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2016.

12. Jiang, J. *Asymptotic Analysis of Mixed Effects Models*; CRC Press: Boca Raton, FL, USA, 2017.

13. McCulloch, C.E.; Searle, S.R.; Neuhaus, J.M. *Generalized, Linear, and Mixed Models*, 2nd ed.; John Wiley & Sons: New York, NY, USA, 2008.

14. Stroup, W.W. *Generalized Linear Mixed Models*; CRC Press: Boca Raton, FL, USA, 2013.

15. Bates, D.; Mächler, M.; Bolker, M.B.; Walker, C.S. Fitting Linear Mixed-Effects Models Using lme4. *J. Stat. Softw.* **2015**, *67*, 1–48. [CrossRef]

16. Boik, J.R. The Analysis of Two-Factor Interactions in Fixed Effects Linear Models. *J. Educ. Stat.* **1993**, *18*, 1–40. [CrossRef]

17. Yi, Q.; Panzarella, T. *Estimating Sample Size for Tests on Trends across Repeated Measurements with Missing Data Based on the Interaction Term in a Mixed Model*; Elsevier: Amsterdam, The Netherlands, 2002; Volume 23, pp. 481–96.

18. Robin, G.; Josse, J.; Moulines, E.; Tibshirani, R. Main effects and interactions in mixed and incomplete data frames. *J. Am. Statist. Assoc.* **2019**, *115*, 1292–1303. [CrossRef]

19. Giesselmann, M.; Schmidt-Catran, A.W. Interactions in fixed effects regression models. *Sociol. Methods. Res.* **2022**, *51*, 1100–1127. [CrossRef]

20. Martin, D.A.; Quinn, M.K.; Park, H.J. MCMCPack: Markov chain Monte Carlo in R. *J. Stat. Softw.* **2011**, *42*, 22. [CrossRef]

21. Hadfield, J.D. MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package. *J. Stat. Softw.* **2010**, *33*, 1–22. [CrossRef]

22. Khodadadian, A.; Parvizi, M.; Teshnehlab, M.; Heitzinger, C. Rational Design of Field-Effect Sensors Using Partial Differential Equations, Bayesian Inversion, and Artificial Neural Networks. *Sensors* **2022**, *22*, 4785. [CrossRef]

23. Noii, N.; Khodadadian, A.; Ulloa, J.; Aldakheel, F.; Wick, T.; Francois, S.; Wriggers, P. Bayesian Inversion with Open-Source Codes for Various One-Dimensional Model Problems in Computational Mechanics. *Arch. Comput. Methods Eng.* **2022**, *29*, 4285–4318. [CrossRef]

24. Matuschek, H.; Kliegl, R.; Vasishth, S.; Baayen, H.; Bates, D. Balancing Type I error and power in linear mixed models. *J. Mem. Lang.* **2017**, *94*, 305–315. [CrossRef]

25. Dempster, A.; Laird, N.; Rubin, D. Maximum-likelihood from incomplete data via the EM algorithm. *J. R. Statist. Soc. Ser. B* **1977**, *39*, 1–38.

26. Delyon, B.; Lavielle, M.; Moulines, E. Convergence of a stochastic approximation version of the EM algorithm. *Ann. Stat.* **1999**, *1*, 94–128. [CrossRef]

27. Lavielle, M.; Mbogning, C. An improved SAEM algorithm for maximum likelihood estimation in mixtures of non linear mixed effects models. *Statist. Comput.* **2014**, *24*, 693–707. [CrossRef]

28. Robbins, H.; Monro, S. A Stochastic Approximation Method. *Ann. Math. Stat.* **1951**, *22*, 400–407. [CrossRef]

29. Bhatnagar, S.; Prasad, H.; Prashanth, L. *Stochastic Recursive Algorithms for Optimization*; Springer: Heidelberg, Germay, 2013; Volume 434, pp. 17–28.

30. Celeux, G.; Diebolt, J. The SEM algorithm: A probabilistic teacher algorithm derived from the EM algorithm for the mixture problem. *Comput. Stat. Q.* **1985**, *2*, 73–82.

31. Wei, G.; Tanner, M. A Monte-Carlo implementation of the EM algorithm and the Poor's Man's data augmentation algorithm. *J. Am. Statist. Assoc.* **1990**, *85*, 699–704. [CrossRef]

32. Bennett, D.A. How can I deal with missing data in my study. *Aust. N. Z. J. Public Health* **2001**, *25*, 464–469. [CrossRef] [PubMed]

33. Slaoui, Y.; Nuel, G. Parameter Estimation in a Hierarchical Random Intercept Model with Censored Response: An Approach using a SEM Algorithm and Gibbs Sampling. *Sankhya Indian J. Statist.* **2014**, *76*, 210–233. [CrossRef]

34. Foulley, J.L. Mixed Model Methodology. Part I: Linear Mixed Models. In *Technical Report*; Université de Montpellier: Montpellier, France, 2015.

35. Grund, S.; Lüdtke, O.; Robitzsch, A. Multiple imputation of missing data in multilevel models with the R package mdmb: A flexible sequential modeling approach. *Behav. Res. Meth.* **2021**, *53*, 2631–2649. [CrossRef]

36. Santos, F. *L'algorithme EM: Une Courte Présentation*; CNRS, UMR, 5199: Talence, France, 2015.

37. Celeux, G.; Diebolt, J. Une Version de Type Recuit Simule de L'algorithme EM. Ph.D. Thesis, INRIA, Le Chesnay, France, 1989.

38. Celeux, G.; Chauveau, D.; Diebolt, J. On Stochastic Versions of the EM Algorithm. Ph.D. Thesis, INRIA, Le Chesnay, France, 1995.

39. Celeux, G.; Diebolt, J. A stochastic approximation type EM algorithm for the mixture problem. *Stochastics Int. J. Probab. Stoch. Process.* **1992**, *41*, 119–134. [CrossRef]
40. Richard, F.; Samson, A.; Cuénod, A.C. A SAEM algorithm for the estimation of template and deformation parameters in medical image sequences. *Statist. Comput.* **2009**, *19*, 465–478. [CrossRef]

41. Panhart, X.; Samson, A. Extension of the SAEM algorithm for nonlinear mixed models with 2 levels of random effects. *Biostatistics* **2008**, *10*, 121–135. [CrossRef]

42. Betancourt, M. A conceptual introduction to Hamiltonian Monte Carlo. *arXiv* **2017**, arXiv:1701.02434.

43. Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. Equation of State Calculations by Fast Computing Machines. *J. Chem. Phys.* **1953**, *21*, 1087–1092. [CrossRef]

44. Hastings, W.K. Monte Carlo sampling methods using Markov Chains and their applications. *Biometrika* **1970**, *57*, 97–109. [CrossRef]

45. Roberts, O. G.; Rosenthal, S. J. General state space Markov chains and MCMC algorithms. *Probab. Surv.* **2004**, *1*, 20–71. [CrossRef]

46. Andrieu, C.; De Freitas, N.; Doucet, A.; Jordan, M. An Introduction to MCMC for Machine Learning. *KAP* **2003**, *50*, 5–43.

47. Perret, C.; Bonin, P. Which variables should be controlled for to investigate picture naming in adults? A Bayesian meta-analysis. *Behav. Res.* Methods **2019**, *51*, 2533–2545. [CrossRef]

48. Schafer, L.J.; Graham, W.J. Missing Data: Our View of the State of the Art. *Psychol. Methods* **2002**, *7*, 147–177. [CrossRef]

49. Enders, C.K. *Applied Missing Data Analysis*; Guilford Press: New York, NY, USA, 2010.

50. Little, R.J.; Rubin, D.B. *Statistical Analysis with Missing Data*; John Wiley & Sons: Hoboken, NJ, USA, 2019; Volume 793.

51. Zhang, X. Tutorial: How to Generate Missing Data For Simulation Studies. 2021. Available online: [https://files.osf.io/v1/resources/rq6yb/providers/osfstorage/60e93adf600da501810a8ea4?format=pdf&action=download&direct&version=1](https://files.osf.io/v1/resources/rq6yb/providers/osfstorage/60e93adf600da501810a8ea4?format=pdf&action=download&direct&version=1) (accessed on 8 February 2023).

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.