A Multistage Binomial Model with Measurement Errors: Application to Protein Viability Prediction

1. Introduction

Binomial data are routine in quantitative sciences, naturally arising for instance in biology (e.g. presence or absence of a particular disease, a mutation leading to a viable or an unviable protein, survive or not survive the great migration in the Serengeti), economics (e.g. purchase or not purchase a product), material science (chiral or achiral compound, presence or absence of defects), and behavioral science (e.g. voting, school dropout, divorce). Binomial regression is the modern approach to analyzing such dichotomous outcomes, and most researchers turn to the logistic model [1,2,3].

In the Generalized Linear Model (GLM) framework [4], the logistic regression model for a binary outcome Y assumes a log-linear odds ratio for the success event $Y=1$ versus the failure $Y=0$, i.e. the expectation $\mu$ (also the success probability) of Y satisfies

$$\begin{array}{c}\hfill log\left({\displaystyle \frac{\mu }{1-\mu }}\right)={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_1{X}_p\end{array}$$

(1)

where $X_j$ ($j=1,\cdots ,p$) are explanatory variables for Y, and ${\beta }_j$ ($j=0,1,\cdots ,p$) are regression coefficients. Stated, otherwise, the effects of the predictors of the odds ratio ${\displaystyle \frac{\mu }{1-\mu }}$ of the event $Y=1$ are multiplicative. The logistic assumption (1) is mathematically convenient, allowing an unrestricted parameterization of the binomial GLM: the coefficients ${\beta }_j$ (which measure changes in the log-odds ratio) take values on the entire real line.

Despite the popularity of the standard logistic regression (multiplicative odds) for binary outcomes, the underlying assumption of this model does not cover phenomena resulting from many independent or successive discrete events. Indeed, complex systems often progress through a series of rate-limiting steps which can be subject to modifications by environmental exposures [5]. There can thus be a large number of routes to the failure of such systems, but only the final state is often observable. Indeed, many physical processes can be described by a number of specific events occurring randomly or in some predetermined order. In other words, the occurrence of a particular event often depends on the prior occurrence of many other (partially) related events.

Examples of such events are common along a manufacturing chain: the absence of defects on an end-line product obtained from assembling q different pieces depends on the occurrence of defects in each assembly line producing a piece. If the lines are independent, then the probability $\mu$ of a defect-free final item is the product of the no-defect probabilities ${\mu }_j$ in the q lines ($j=1,\cdots ,q$):

$$\begin{array}{c}\hfill \mu ={\mu }_1\times {\mu }_2\times \cdots \times {\mu }_q.\end{array}$$

(2)

Of course, the design of products is generally optimized so as to make errors additive and the effect of a single defect as insignificant as possible. However, if only the presence or absence of defects is of interest, then the multiplicative absolute risk model (2) reflects the process generating the binary outcome.

Processes with multiplicative risks are in fact common in biology, since many biological processes involve discrete events, at least at some levels of description and observation [6,7]. Consider, for instance, a mutated protein resulting from a single amino acid substitution, with functional viability depending upon two factors: the protein’s ability to fold into a native conformation (related to folding stability), and its ability to bind other subunits or ligands (related to binding stability). If such a protein does not fold properly (e.g. because of a highly destabilizing mutation), it is likely not going to bind (very low probability) since most proteins only bind their targets when folded [8]. The viability of the mutated protein is thus proportional to the probabilities of folding (${\mu }_1$) and binding (${\mu }_2$): $\mu \propto {\mu }_1\times {\mu }_2$ so that any predictor value leading to a small probability begets an overall prediction near zero (i.e. a mutation that renders a protein able to fold, but not bind, or vice versa, should be inviable). In this specific case, the multiplicative odds model (1) appears particularly inappropriate. Indeed, for $p=2$ with $X_1$ and $X_2$ representing folding and binding stabilities, respectively, if we fix folding stability to a value $x_1$, then it is always possible to find some value $x_2$ of binding stability so that the viability gets arbitrarily close to one (including an interaction term ${\beta }_{12}{X}_1{X}_2$ into (1) does not help because this is simply equivalent to adding ${\beta }_{12}{x}_1$ to ${\beta }_2$). This contradicts the fact that when folding stability is disrupted, the functional viability of the mutated protein cannot get arbitrarily close to one, regardless of the binding stability.

More generally, a multiplicative risk model arises naturally for any all-or-none type event. For instance, the presence or absence of a disease may result from the total absence of one or more protective events (e.g. diet deficiencies in required dietary constituents leading to a resultant disease) [9,10]. Conversely the presence or absence of a particular disease may result from one or more adverse events (e.g. encounters with carcinogenic agents leading to a specific cancer). Multiplicative risks are indeed widely used, for instance, in multi-stage cancer models [11,12,13,14,15,16,17], and for incidence or hazard rate models in survival analysis [18,19,20]. In social epidemiology, many studies reported evidence for multiplicative relations between the effects of risk factors and disease incidence or mortality [21,22,23,24,25].

However, to the best of our knowledge, there is no general binary data model with multiplicative absolute risks. To overcome interpretation difficulties related to the non-collapsibility of odds ratio [26,27,28,29,30], some works considered log-binomial models where, as in a Poisson GLM, covariates have multiplicative effects on the relative risk [31,32,33,34]. In addition, the Zero-Inflated Bernoulli (ZIB) model, a special case of the zero-inflated binomial model [35], is related to the multiplicative risk model (2); as we discuss in Section 2, the ZIB model can be described as a special case of our proposal. In spite of some recent developments in relative risk modeling where the probability of the target event is modeled as a ratio to the probability under a reference condition [36,37], these models not only lack the mathematical convenience of the logistic model (sample dependent constraints required on the parameter space to ensure that predicted probabilities are restricted to $[0,1]$, or restricted ranges of covariates), but also have their own set of interpretation issues [38,39,40,41].

We propose in this work an extension of the Traditional Binomial (TB) model to handle binomial data resulting from data generation processes with known multiplicative risk components. The proposed model assumes that in addition to an observed binomial outcome, we are aware of some unmeasurable prior events that determine the final outcome (with paths to success or failure involving sequential, or non-sequential steps, or a combination of both), and each prior event is known to be related to a subset of the observable explanatory variables. The proposal results from combining equation (2) with a standard binomial model inverse link function, i.e. each factor ${\mu }_j$ is given by ${\mu }_j=h\left({\eta }_j\right)$ with $h:\mathbb{R}\to [0,1]$ and ${\eta }_j\in \mathbb{R}$ a linear predictor obtained from explanatory variables related to ${\mu }_j$. Most popular binomial model link functions include the logit link function (1) with inverse the cumulative distribution function (cdf) of the logistic distribution,

$$\begin{array}{c}\hfill h\left({\eta }_j\right)={\left[1+exp\left(-{\eta }_j\right)\right]}^{-1},\end{array}$$

(3)

and the probit link with $h\left({\eta }_j\right)=\Phi \left({\eta }_j\right)$ where $\Phi$ is the cdf of the standard normal distribution. The resulting model is termed Multistage Binomial (MSB) model. It is worthwhile emphasizing that the MSB model (with multiple stages) is not a generic model for binomial data. Its use should be prescribed by some theoretical or empirical knowledge of the phenomenon being modeled. In case of only a (theoretical) suspicion of multiple stages, the MSB model should be tested against the TB model.

In addition to providing a general framework for modeling multiplicative risks in binomial data, nesting the TB and ZIB models, our proposed MSB model allows to control for two important sources of issues in statistical analysis: omitted variables and measurement errors in predictors. Indeed, researchers can seldom be aware of, let alone measure, all of the variables that determine a target outcome, and the impact of left out or omitted variables can be dramatic for interpretations in TB models [42,43,44,45,46,47,48,49,50]. Our formulation of the MSB model allows to control for omitted variables in some variable omission scenarios, i.e., when the explanatory variables which determine some factors ${\mu }_j$ in (2) are omitted. Moreover, observed predictors are often error-prone. To account for errors in observed explanatory variables, the proposed MSB model also allows for heteroscedastic additive Berkson measurement errors [51] when the variability of each numerical predictor around its observed values is known.

For inference in the MSB model, we developed a Penalized Maximum Likelihood (PML) estimation procedure. The PML method has two main advantages over a simple Maximum Likelihood (ML) approach: reduce small sample bias and circumvent estimation issues related to quasi or complete separation. On the one hand, although ML estimates in TB regression are asymptotically unbiased to first order, it is well known that the ML estimators are biased to first order in small sample ($n<200$ [52,53]). This occurs especially when the marginal success probability of the response approaches zero or one [54]. The ML estimators of the MSB model parameters are also expected to be biased to the first order in small samples because this is a general feature of binomial regression models [55].

On the other hand, complete separation occurs when one or more of the explanatory variables can perfectly predict some binary outcomes [56]. Quasi separation occurs when for any values of an explanatory variable (or a linear predictor) $x_i$ less than a threshold, say $x_i<t$, the binary outcome is always $y_i=0$ (or always $y_i=1$), but when $x_i\ge t$, $y_i$ might equal zero or one. The occurrence of (quasi) separation is generally related to small sample, sparse data, rare outcome ($E\left[Y_i\right]$ is very small), rare exposures (and more generally predictors with very small variability), large number of dichotomous risk factors [57], and highly correlated predictors (multicollinearity) [58], or covariates with strong effects [59].

Since the MSB model can be viewed as a combination of many TB regressions across stages of a complex process, the occurrence of separation becomes more likely not only as the number p of predictors in any stage increases, but also as the number q of stages increases, for not sufficiently large samples. As a consequence, in the MSB model framework, we anticipate an increased likelihood of unrealistically large ML estimates or related large (or undefined) standard errors as p or q increases, and as the sample size n or the maximum success probabilities decrease. Quasi or complete separation implies that maximum likelihood estimates approach infinity or zero [60,61]. In practice, separation generally results in a flat log-likelihood surface (likely with an asymptote) which leads to implausible estimates, and no or very large standard error estimates [62].

In this paper, we investigate some theoretical and numerical properties of MSB models. We first study the identifiability of the proposed model. Then, using simulation experiments, we study finite-sample properties of the PML estimates of the MSB model parameters. Specifically, we i) compare, for demonstration purposes, the Two-Stage Binomial (TSB) model to TB and ZIB models considering a simulated biological case study on the viability of a mutated protein; ii) assess the ability of PML estimators of MSB to recover the population values of regression parameters and predict event success probabilities in finite samples; iii) and evaluate finite-sample properties of a non-standard Likelihood Ratio (LR) test to motivate the choice of the MSB model over alternatives such as the TB or ZIB model.

The organization of the remainder of the paper is as follows. We specify the MSB model, discuss its identifiability and describe the PML inference in Section 2. Then, we present simulation results in Section 3. We discuss our results and finish with concluding remarks in Section 4.

2. The Multistage Binomial Model

The MSB model is built as an extension of the standard binomial GLM [4]. In this section, we specify and describe the main features of the MSB model, starting with the simpler model in the absence of measurement errors in observed covariates, and then introducing measurement errors. We then establish and discuss conditions for the identifiability of the model. Finally, we describe the estimation of model parameters using Jeffreys invariant prior as likelihood penalty, and an LR test to the detect omission of important process stages or predictors.

2.1. Measurement Error-Free Model

Let us consider a design with n units indexed i, each having $n_i$ (known positive integer) independent Bernoulli trials. The binomial outcome $Y_i$ is assumed to be generated through a multistage process with q known steps, each related to a $p_j$-vector of covariates ${\mathit{W}}_{ij}$ ($j=1,2,\cdots ,q$). The multiplicative risk binomial model is given for the ith outcome $Y_i$ as:

$$\begin{array}{cc}\hfill \left(Y_i\right|{\mathit{W}}_i={\mathit{w}}_i)& \stackrel{ind}{\sim }\mathcal{B}in(n_i,{\xi }_i),\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill {\xi }_i& ={\lambda }_i\left[{\alpha }_i+(1-{\alpha }_i)\prod _{j=1}^q{\xi }_{ij}\right],\hfill \end{array}$$

(4b)

$$\begin{array}{cc}\hfill {\xi }_{ij}& =h\left({\zeta }_{ij}\right),\hfill \end{array}$$

(4c)

$$\begin{array}{cc}\hfill {\zeta }_{ij}& ={\mathit{w}}_{ij}^{\top }{\mathit{\beta }}_j,\hfill \end{array}$$

(4d)

where ${\mathit{W}}_i^{\top }=({\mathit{W}}_{i1}^{\top },\cdots ,{\mathit{W}}_{iq}^{\top })$ is the vector of covariates; $\mathcal{B}in(n_i,{\xi }_i)$ denotes the binomial distribution with $n_i$ trials and success probability ${\xi }_i$; ${\lambda }_i\in (0,1]$ is the upper limit of the success probability (${\xi }_i$) of each Bernoulli trial given ${\mathit{W}}_i$ and ${\alpha }_i\in [0,1)$ is such that ${\alpha }_i{\lambda }_i$ is the lower limit of ${\xi }_i$; h is the inverse link function for the binomial outcome (h is assumed to be a continuous cdf, monotone, increasing and twice differentiable); and ${\mathit{\beta }}_j$ is a $p_j$-vector of real regression parameters (possibly including an intercept ${\beta }_{j0}$).

Model (4) describes a binomial process where the outcome results from q latent stages where each stage is represented by an independent Bernoulli experiment, and the final outcome corresponds to all experiments being simultaneously successful. As indicated above, the success probability ${\xi }_i$ is bounded as ${\xi }_i\in [{\alpha }_i{\lambda }_i,{\lambda }_i]$. These bounds represent theoretical limits that may arise from some unknown stages, or some known stages with no measured explanatory variable. Indeed, setting $n_i=1$ without loss of generality, and starting from multiplicative risks of the form (2) with $q+1$ stages (${\xi }_i={\prod }_{j=0}^q{\xi }_{ij}$) among which the first ($j=0$) is omitted, the conditional success probability of the outcome given the measured covariates ${\mathit{w}}_i$ is then

$$\begin{array}{cc}\hfill P(Y_i=1|{\mathit{W}}_i={\mathit{w}}_i)& ={\int }_{{\mathbb{R}}^{p_0}}h\left({\mathit{w}}_{i0}^{\top }{\mathit{\beta }}_0\right)\times \left(\prod _{j=1}^q{\xi }_{ij}\right)f_{{\mathit{w}}_0}\left({\mathit{w}}_{i0}\right)d{\mathit{w}}_{i0}\hfill \\ \hfill & =\lambda \times \prod _{j=1}^q{\xi }_{ij}\hfill \end{array}$$

(5)

where $\lambda ={\int }_{{\mathbb{R}}^{p_0}}h\left({\mathit{w}}_{i0}^{\top }{\mathit{\beta }}_0\right)f_{{\mathit{w}}_0}\left({\mathit{w}}_{i0}\right)d{\mathit{w}}_{i0}$ is the success probability at the omitted stage, averaged over of the distribution of the related $p_0$-vector ${\mathit{W}}_{i0}$ of explanatory variables, which has joint density function $f_{{\mathit{w}}_0}$ (integration will be replaced by summation where $f_{{\mathit{w}}_0}$ is a probability mass function (pmf), i.e., if ${\mathit{W}}_{i0}$ is discrete).

Stated in other words, the omission of explanatory variables related to a stage in the multiplicative-risk model (2) leads to (4b) with ${\alpha }_i=0$. It is worth noticing that while this rationale applies to unobserved confounders that vary across experimental or observational units, it remains valid when the vector of predictors ${\mathit{W}}_0$ of a stage is (willfully or not) held constant across all investigated/sampled units: in this case, $\lambda$ is not an average of success probabilities across the possible values of ${\mathit{W}}_0$, but rather the conditional success probability of this stage given the constant value ${\mathit{w}}_0$ of ${\mathit{W}}_0$: $\lambda =h\left({\mathit{w}}_0^{\top }{\mathit{\beta }}_0\right)$.

This development indicates that $\lambda$ can arise from one or both of two potential sources: uncontrolled confounders (marginalization), or a limit of the data generation mechanism in the specific investigated experimental or observational environment/condition (conditioning). In the general formulation (4) of the MSB model, we allow $\lambda$ to vary across observations to potentially account for conditioning when the study is replicated across some discrete experimental/environmental settings: ${\lambda }_i$ can capture the maximum success probability for all units i from the same setting. For the sake of symmetry, ${\alpha }_i\in (0,1)$ is allowed in model (4): since, one can equivalently model the failure of an event to occur instead of its success, the minimum success probability can also be greater than zero. By the same reasoning, this rationale can be extended to the omission and marginalization over many stages, or conditioning on one value of the vector of predictors at each of many stages.

Remark 1.

When model (4) has only one ($q=1$) or two stages ($q=2$), we obtain the following two important special cases.

• When $q=1$, ${\alpha }_i\to 0$, and ${\lambda }_i\to 1$, the MSB model is reduced to the TB model.
• For binary outcomes ($n_i=1$ for $i=1,2,\cdots ,n$), when $q=2$, ${\alpha }_i\to 0$, and ${\lambda }_i\to 1$, the MSB model is reduced to the ZIB model.
  The ZIB model has been proposed by Hall [35] to face situations where a binary outcome contains an excess of zeros as compared to the expected zeros under the TB model. The ZIB model assumes that the binary outcome is the product of two Bernoulli processes; a zero from either process (or both) results in a zero. An example of process generating such data in health sciences was investigated by Diop et al. [63]: a study population consisting of a mixture of susceptible subjects who can experience the binary outcome of interest, and cured subjects (cure fraction) who cannot experience the binary outcome. Complications arise from not knowing the immunity (cured or susceptible) status of each subject, and the ZIB model is used to jointly estimate the cure fraction on the one hand, and the logistic regression for the susceptible subjects on the other hand [64,65,66]. Using the MSB model framework ($q=2$, $n_i=1$, ${\alpha }_i=0$, and ${\lambda }_i=1$) to describe the ZIB model, the success probability of the binary outcome $Y_i$ is $E\left[Y_i\right]=h\left({\eta }_{i1}\right)h\left({\eta }_{i2}\right)$ and the probability that a subject i is cured is $1-h\left({\eta }_{i1}\right)$. Thus, in this case, the main difference between the MSB model and the ZIB model is that the MSB model allows the minimum probability ${\alpha }_i{\lambda }_i$ to be greater than zero, and the maximum success probability ${\lambda }_i$ to be less than one.
  When the susceptibility probability $h\left({\eta }_{i1}\right)$ is constant across subjects in the ZIB model, the latter can alternatively be considered as a one stage model ($q=1$) with ${\alpha }_i=0$ and $\lambda =h\left({\eta }_{i1}\right)$. In this case, the success probability of the binary outcome $Y_i$ is $E\left[Y_i\right]=\lambda h\left({\eta }_{i1}\right)$ and the probability that a subject i is cured is $1-\lambda$.

In complement to Remark 1, it is worthwhile emphasizing that the MSB model is not a zero inflation model, nor a mixture model for count data. These models only share some special cases with the MSB model when each and every count is binary. For general binomial data, the zero-inflated binomial model [67] is not a special case of the MSB model, even when $q=2$.

Remark 2.

An ovulation-time model was developed by Garel et al. [68] to investigate how the probability of ovulation and body growth co-evolve over time in moose. The special case of the probit link MSB model when $q=2$, ${\alpha }_i\to 0$ and ${\lambda }_i\to 1$, is closely related to the ovulation-time model. In the MSB model framework, the first stage in their model represents susceptibility to ovulation at a time, and the second stage is the occurrence of ovulation (this also points to a close relation with the ZIB model). However, unlike the fixed-scale link used in proper binomial models, the ovulation-time model employs a covariate-dependent probit link scale to capture spatiotemporal dynamics in the occurrence of ovulation.

2.2. Multistage Binomial Model Specification

In practice, the covariates ${\mathit{W}}_i$ in model (4) can be observed with measurement errors. We consider a situation where ${\mathit{W}}_i$ is error-prone: instead of ${\mathit{W}}_i$, we only have access to an observable ${\mathit{X}}_i^{\top }=({\mathit{X}}_{i1}^{\top },\cdots ,{\mathit{X}}_{iq}^{\top })$ such that ${\mathit{X}}_{ij}$ satisfies:

$$\begin{array}{c}\hfill {\mathit{X}}_{ij}={\mathit{W}}_{ij}+{\mathit{e}}_{ij}\end{array}$$

(6)

where ${\mathit{e}}_{ij}={\left(e_{ij1},\cdots ,e_{ij{p}_j}\right)}^{\top }$ is a measurement error vector independently following an unspecified zero-mean distribution,

$$\begin{array}{c}\hfill {\mathit{e}}_{ij}\stackrel{ind}{\sim }ME\left({\mathit{S}}_{ij}\right),\end{array}$$

(7)

which has covariance matrix ${\mathit{S}}_{ij}$ (a known semi-positive definite matrix), and a probability density function (pdf) $f_{\mathit{e}}(·|{\mathit{S}}_{ij})$. The measurement error $e_{ijk}$ in (6) is known as an additive Berkson error [51]. This error model is appropriate when a subset of the available explanatory variables are obtained from an estimation process such that the observed ${\mathit{X}}_{ij}$ can be regarded as an estimate of the truth ${\mathit{W}}_{ij}$, and ${\mathit{S}}_{ij}$ is an estimate of the covariance matrix of ${\mathit{X}}_{ij}$.

For inference purposes here, we are interested in the expected value of the outcome $Y_i$ given an observed covariate vector ${\mathit{x}}_i$: $E\left(Y_i\right|{\mathit{X}}_i={\mathit{x}}_i)=n_i{\mu }_i$ where $n_i$ is the number of trials, and ${\mu }_i$ is the success probability given ${\mathit{X}}_i={\mathit{x}}_i$. The Multistage Binomial (MSB) model is defined as:

$$\begin{array}{cc}\hfill \left(Y_i\right|{\mathit{X}}_i={\mathit{x}}_i)& \stackrel{ind}{\sim }\mathcal{B}in(n_i,{\mu }_i),\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill \mathrm{with}{\mu }_i& ={\lambda }_i\left[{\alpha }_i+(1-{\alpha }_i)\prod _{j=1}^q{\mu }_{ij}\right]\hfill \end{array}$$

(8b)

where ${\lambda }_i$ and ${\alpha }_i$ are as defined in (4), and ${\mu }_{ij}=\overline{h}\left({\mathit{x}}_{ij}|{\mathit{\beta }}_j,{\mathit{S}}_{ij}\right)$ is obtained for $j=1,2,\cdots ,q$ by integrating out the measurement error vector ${\mathit{e}}_{ij}$ as:

$$\begin{array}{cc}\hfill \overline{h}\left({\mathit{x}}_{ij}|{\mathit{\beta }}_j,{\mathit{S}}_{ij}\right)& ={\int }_{{\mathbb{R}}^{p_j}}h\left({\eta }_{ij}-{\mathit{e}}_{ij}^{\top }{\mathit{\beta }}_j\right)f_{\mathit{e}}\left({\mathit{e}}_{ij}\right|{\mathit{S}}_{ij})d{\mathit{e}}_{ij}\hfill \end{array}$$

(9)

with ${\eta }_{ij}={\mathit{x}}_{ij}^{\top }{\mathit{\beta }}_j$ (summation replaces integration where $f_{\mathit{e}}$ is a pmf, i.e., if $ME\left({\mathit{S}}_{ij}\right)$ is discrete). From the conditional expectation $E\left(Y_i\right|{\mathit{X}}_i)=n_i{\mu }_i$, the conditional variance of the binomial outcome $Y_i$ is

$$\begin{array}{cc}\hfill Var\left(Y_i\right|{\mathit{X}}_i)& =n_i{\mu }_i\left(1-{\mu }_i\right).\hfill \end{array}$$

(10)

When there is no measurement error (${\mathit{S}}_{ij}=\mathit{0}$), we have ${\mu }_{ij}=h\left({\eta }_{ij}\right)$ and ${\eta }_{ij}={\zeta }_{ij}$ from ${\mathit{e}}_{ij}=\mathit{0}$, hence ${\mu }_{ij}={\xi }_{ij}$. It follows that the measurement error-free model (4) is a special case of the MSB model (8) when ${\mathit{S}}_{ij}=\mathit{0}$ for all stages $j=1,\cdots ,q$; and across all units $i=1,\cdots ,n$.

A design for the general MSB model (8) includes the response vector $\mathit{y}={(y_1,\cdots ,y_n)}^{\top }$ of n independent outcomes, the associated $n\times p$ (with $p={\sum }_{j=1}^q{p}_j$) design matrix $\mathbf{X}={({\mathit{X}}_1,\cdots ,{\mathit{X}}_n)}^{\top }$, and the array $\mathbf{S}$ of $nq$ known $p_j\times p_j$ covariance matrices ${\mathit{S}}_{ij}$ of measurement errors. Because the MSB model allows the success probability bounds ${\lambda }_i$ and ${\alpha }_i{\lambda }_i$ to vary across experimental or observational units, the complete design includes two additional design matrices $\mathbf{U}={({\mathit{u}}_1,\cdots ,{\mathit{u}}_n)}^{\top }$ and $\mathbf{V}={({\mathit{v}}_1,\cdots ,{\mathit{v}}_n)}^{\top }$ where ${\mathit{u}}_i={(u_{i1},\cdots ,u_{id})}^{\top }$ is a d-vector and ${\mathit{v}}_i^{\top }=(v_{i1}^{\top },\cdots ,v_{ir}^{\top })$ is an r-vector such that ${\alpha }_i=h\left({\mathit{u}}_i^{\top }\mathit{\gamma }\right)$ and ${\lambda }_i=h\left({\mathit{v}}_i^{\top }\mathit{\delta }\right)$, and $\mathit{\gamma }\in {\mathbb{R}}^d$ and $\mathit{\delta }\in {\mathbb{R}}^r$ are parameters that determine the success probability bounds. In particular, each of ${\alpha }_i$ and ${\lambda }_i$ may have a fixed value across units, i.e., we can have $\mathbf{U}={\mathit{J}}_n$ (n-vector of all ones) and $\mathbf{V}={\mathit{J}}_n$ so that ${\alpha }_i={\alpha }_0$ and ${\lambda }_i={\lambda }_0$ with ${\alpha }_0=h\left({\gamma }_0\right)$, ${\lambda }_0=h\left({\delta }_0\right)$, ${\gamma }_0,{\delta }_0\in \mathbb{R}$ ($d=r=1$). In practice, we may also allow empty $\mathit{\delta }$ and $\mathit{\gamma }$ ($d=r=0$), i.e., ${\alpha }_i$ and ${\lambda }_i$ have fixed and known values, e.g. ${\alpha }_i=0$ and ${\lambda }_i=1$ for all $i=1,2,\cdots ,n$.

As motivated in Section 2.1, the additional covariates making ${\alpha }_i$ and/or ${\lambda }_i$ vary across units should be discrete, i.e. each ${\mathbf{U}}_i$ and ${\mathbf{V}}_i$ should be binary/contrast matrices (obtained by an appropriate coding of some categorical predictors) for the MSB model (8). As such, categorical predictors describing site conditions should be included into the design matrices $\mathbf{U}$ and $\mathbf{V}$ to define multiplicative intercepts for groups. Any stage of the considered multistage process can also include categorical variables. In the presence of categorical predictors, the MSB model becomes a multistage analogue of the analysis of covariance model for a binary response: categorical predictors define groups with different intercepts in the study population [69].

2.3. Conjugate Distribution for Measurement Errors

The possibly multivariable integrals (9) in the MSB model specification can make model estimation a difficult task requiring integral approximation methods such as Gaussian quadrature [70,71], or EM-type algorithm [72,73] which can be computationally demanding. To circumvent such difficulties, we consider an implicit conjugate distribution for measurement errors, such that each integral (9) is tractable. As per Proposition 1 below (see a proof in Appendix A.1), under this conjugate distribution of measurement errors, the integral (9) is given by the inverse link function h.

Proposition 1

(Conjugate measurement error distribution). Assume that the probability distribution corresponding to the link function h in the MSB model (8) has a finite variance ${\sigma }_o^2<\infty$. Then, there exists a family of distributions $ME\left({\mathit{S}}_{ij}\right)$ such that the function $\overline{h}$ defined in (9) satisfies

$$\begin{array}{cc}\hfill \overline{h}\left({\mathit{x}}_{ij}|{\mathit{\beta }}_j,{\mathit{S}}_{ij}\right)& =h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\hfill \end{array}$$

(11)

where ${\eta }_{ij}={\mathit{x}}_{ij}^{\top }{\mathit{\beta }}_j$ and ${\overline{\sigma }}_{ij}^2={\mathit{\beta }}_j^{\top }{\mathit{S}}_{ij}{\mathit{\beta }}_j$.

Remark 3.

The conjugate ME distribution is known for popular link functions.

• For the logit link, ${\sigma }_o^2={\displaystyle \frac{{\pi }^2}{3}}$ and the conjugate distribution $ME\left({\mathit{S}}_{ij}\right)$ of measurement errors ${\mathit{e}}_{ij}$ is the distribution such that the linear combination ${\mathit{e}}_{ij}^{\top }{\mathit{\beta }}_j$ follows a logistic Bridge distribution (which is the conjugate distribution of random-effects in marginalized mixed effects logistic models [74]).
• For the probit link, ${\sigma }_o^2=1$ and the conjugate measurement error distribution is the normal distribution with null mean vector and covariance matrix ${\mathit{S}}_{ij}$. This follows by Proposition 4 in [75].

From here on and afterwards, we restrict attention to the measurement error distribution $ME\left({\mathit{S}}_{ij}\right)$ implied by (11), unless otherwise specified. Under the conjugate $ME\left({\mathit{S}}_{ij}\right)$ assumption, the success probability ${\mu }_i$ in (8b) satisfies:

$$\begin{array}{cc}\hfill {\mu }_i& ={\lambda }_i\left[{\alpha }_i+(1-{\alpha }_i)\prod _{j=1}^{q}h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\right]\hfill \end{array}$$

(12)

where ${\lambda }_i=h\left({\mathit{u}}_i^{\top }\mathit{\delta }\right)$, ${\alpha }_i=h\left({\mathit{v}}_i^{\top }\mathit{\gamma }\right)$, ${\eta }_{ij}={\mathit{x}}_{ij}^{\top }{\mathit{\beta }}_j$, and ${\overline{\sigma }}_{ij}^2={\mathit{\beta }}_j^{\top }{\mathit{S}}_{ij}{\mathit{\beta }}_j$.

2.4. Regularity Conditions and Identifiability

The MSB model (8) is parameterized by the m-vector of all unknown parameters $\mathit{\theta }={\left({\mathit{\beta }}^{\top },{\mathit{\gamma }}^{\top },{\mathit{\delta }}^{\top }\right)}^{\top }$ with ${\mathit{\beta }}^{\top }=\left({\mathit{\beta }}_1^{\top },\cdots ,{\mathit{\beta }}_q^{\top }\right)$, $m=p+d+r$. If ${\alpha }_i$ or ${\lambda }_i$ have known values (e.g. ${\alpha }_i=0$ and ${\lambda }_i=1$), the corresponding sub-vectors ($\mathit{\delta }$ and $\mathit{\gamma }$) are empty ($d=0$ and $r=0$, and $\mathit{\theta }=\mathit{\beta }$). Classical inference on $\mathit{\theta }$ is achieved by the ML approach. For inference purposes, we present here the log-likelihood function of the MSB model. Recall that the pmf of a binomial distribution with $n_i$ trials and success probability ${\mu }_i$ has the form [76]

$$\begin{array}{c}\hfill P(Y_i=y_i)=\left({\displaystyle \genfrac{}{}{0pt}{}{n_i}{y_i}}\right){\mu }_i^{y_i}{\left(1-{\mu }_i\right)}^{n_i-y_i}.\end{array}$$

(13)

By the pmf (13), the log-likelihood of a parameter vector $\mathit{\theta }$ given a sample of n independent units $\left(y_i,{\mathbf{X}}_i,{\mathbf{U}}_i,{\mathbf{V}}_i,{\mathbf{S}}_i\right)$ reads

$$\begin{array}{c}\hfill {\ell }_n\left(\mathit{\theta }\right|\mathit{y})=\sum _{i=1}^n\left[log\left({\displaystyle \genfrac{}{}{0pt}{}{n_i}{y_i}}\right)+y_{i}log\left({\mu }_i\right)+(n_i-y_i)log\left(1-{\mu }_i\right)\right]\end{array}$$

(14)

where ${\mu }_i$ is as defined in Equation (12).

Let us denote ${\mathit{\theta }}_{\circ }={\left({\mathit{\beta }}_{\circ }^{\top },{\mathit{\gamma }}_{\circ }^{\top },{\mathit{\delta }}_{\circ }^{\top }\right)}^{\top }$ as the vector of true model parameter. We require the following regularity conditions to ensure the identifiability of the MSB model, i.e., the ability to infer the parameter vector ${\mathit{\theta }}_{\circ }$ from observed data.

C1-

All predictors in the MSB model are bounded, i.e., there exist compact sets ${\mathcal{C}}_j\in {\mathbb{R}}^{p_j}$ (for $j=1,2,\cdots ,q$), ${\mathcal{C}}_u\in {\mathbb{R}}^d$ and ${\mathcal{C}}_v\in {\mathbb{R}}^r$ such that ${\mathbf{X}}_{ij}\in {\mathcal{C}}_j$, ${\mathbf{U}}_i\in {\mathcal{C}}_u$, and ${\mathbf{V}}_i\in {\mathcal{C}}_v$ for every $i=1,2,\cdots ,n$. For every $i=1,2,\cdots ,n$, $Var\left[X_{ijk}\right]>0$ for $k=1,2,\cdots ,p_j$, $Var\left[U_{ik}\right]>0$ for $k=1,2,\cdots ,d$, and $Var\left[V_{ik}\right]>0$ for $k=1,2,\cdots ,r$. For every $i=1,2,\cdots ,n$, $j=1,2,\cdots ,q$, and $k=1,2,\cdots ,p_j$, the $X_{ijk}$ are linearly independent. For every $i=1,2,\cdots ,n$ and $k=1,2,\cdots ,d$, and $k=1,2,\cdots ,d$ the $U_{ik}$ are linearly independent, and for every $i=1,2,\cdots ,n$ and $k=1,2,\cdots ,r$, the $V_{ik}$ are linearly independent.

C2-

The true parameter vectors ${\mathit{\beta }}_{j\circ }$, ${\mathit{\gamma }}_{\circ }$ and ${\mathit{\delta }}_{\circ }$ lie in the interior of known compact sets ${\mathcal{B}}_j\in {\mathbb{R}}^{p_j}$ for $j=1,2,\cdots ,q$, ${\mathcal{B}}_d\in {\mathbb{R}}^d$, and ${\mathcal{B}}_r\in {\mathbb{R}}^r$ respectively.

C3-

There exists, for each stage j, a continuous predictor $T_j$ which is in ${\mathit{X}}_j$ but not in any ${\mathit{X}}_k$ (for all $k\ne j$), nor in $\mathbf{U}$, neither in $\mathbf{V}$. Moreover, there exists at least one continuous predictor, either $T_u$ in $\mathbf{U}$ but not in $\mathbf{V}$ nor in ${\mathit{X}}_j$ (for all $j=1,2,\cdots ,q$), or $T_v$ in $\mathbf{V}$ but not in $\mathbf{U}$ nor in ${\mathit{X}}_j$ (for all $j=1,2,\cdots ,q$).

The conditions C1–C3 are extensions of the corresponding conditions in Diop et al. [63]. The assumptions C1 and C2 are classical requirements for identifiability in the TB model. With regard to condition C3, it ensures that the product ${\prod }_{j=1}^q{\mu }_{ij}$ in (8b) is not invariant under any permutation of the parameter vectors ${\mathit{\beta }}_j$ ($j=1,\cdots ,q$), $\mathit{\gamma }$ and $\mathit{\delta }$.

Remark 4.

As emphasized by Diop et al. [63], the condition C3 is not particularly stringent.

-

In fact, the motivation of the MSB model assumes that we have some prior knowledge of the steps of the multistage process that generates the binary outcome $Y_i$. This implies that, for each step that we are aware of, we know about at least one specificity that differentiates it from other known steps. This is reflected by condition C3. If our prior knowledge of the steps of the process is confused to the point that two stages share the exact same predictors or risk factors, then we should consider that the two stages are only one known (meaningful) stage.

-

If all stages of the modeled process are known and accounted for, the model approaches the limiting special case where ${\alpha }_i=0$ and ${\lambda }_i=1$. In this situation, the distinct continuous predictor in each of the q stages (condition C3 is only required for $q-1$ stages. That is, one of the q stages may not have a distinct continuous predictor.

-

A continuous predictor is required in the design matrix of either ${\alpha }_i$ ($T_u$ in $\mathbf{U}$) or λ ($T_v$ in $\mathbf{V}$) only if both ${\alpha }_i$ and ${\lambda }_i$ are being estimated. For instance, in the important special cases where either ${\alpha }_i=0$ or ${\lambda }_i=1$, no distinct continuous predictor is required in any of ${\alpha }_i$ and ${\lambda }_i$ design matrices, if present.

Based on the conditions C1–C3, Proposition 2 establishes identifiability of the MSB model (8) (see a proof in Appendix A.2).

Proposition 2

(Identifiability of an MSB model). Under the conditions C1–C3, the MSB model (8) with a logit or probit link function and the conjugate measurement error distribution is identifiable; i.e., ${\ell }_n\left(\mathit{\theta }\right|\mathit{y})={\ell }_n\left({\mathit{\theta }}^{*}\right|\mathit{y})$ almost surely implies $\mathit{\theta }={\mathit{\theta }}^{*}$.

Remark 5.

Proposition 2 considers the case of conjugate measurement error distribution with known error covariance matrices ${\mathit{S}}_{ij}$. Note that this includes the special case of no measurement error (${\mathit{S}}_{ij}=\mathit{0}$). Some extensions are also readily available.

• For a logit link function, Proposition 2 remains true if the conjugate measurement error distribution were replaced by the (multivariate) normal distribution with null mean vector. This follows from Theorem 1 of Shklyar [77], given that the error covariance matrices ${\mathit{S}}_{ij}$ are assumed known in the MSB model (8).
• In the case $q=1$, the model remains identifiable if all covariance matrices ${\mathit{S}}_{ij}$ are only known up to a unique scaling factor (see Theorem 2 of Kuchenhoff [78]). In particular, if the unique stage includes only one error-prone continuous predictor, then under conditions C1–C3, a homogeneous variance of the measurement errors is identifiable along with all other model parameters. Extension to $q>1$ requires further investigations.

It is important to observe that a null true model parameter ${\mathit{\theta }}_0=\mathit{0}$ is not generally identifiable (unless $q=1$, and ${\alpha }_i$ and ${\lambda }_i$ are known constants, which essentially corresponds to the standard binomial model). Indeed, the condition C3 implies that some continuous predictors in the design matrices must have non-zero regression slopes. That is, the vector ${\mathit{\beta }}_j$ of true regression coefficients in each stage j must have at least one non-zero slope.

Also note that along with C1 and C2, condition C3 is sufficient for identifiability of the MSB model (8), but may not be required. That is, identifiability may also be established under some alternative conditions replacing C3. These may be physical or biological constraints on the predictors in each stage of the model, or restrictions on model parameters. For instance, in their probit link ovulation-time model, Garel et al. [68] imposed constraints on the parameter space (some coefficients were set to zero) to ensure numerical identifiability.

2.5. Penalized Maximum Likelihood Estimation

The ML estimates of the MSB model (8) parameters can be obtained by solving the score equation ${\dot{\mathit{\ell }}}_n=\mathit{0}$. The score function ${\dot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})={\displaystyle \frac{\partial {\ell }_n\left(\mathit{\theta }\right|\mathit{y})}{\partial \mathit{\theta }}}$ of the MSB model is given in Appendix B. To reduce bias in MSB parameter estimates for small samples and address quasi or complete separation issues, we consider a likelihood penalization approach which is the common route to first order bias correction in the standard logistic regression model [55,60,79,80,81]. For general exponential family models under canonical parametrization, Firth showed that using Jeffreys invariant prior as likelihood penalty removes the first-order bias term in the asymptotic expansion of the ML estimator [55]. The Jeffreys penalty, which is the square root of the determinant of the Fisher information matrix [82], is asymptotically negligible. As a result, the maximum a posteriori estimator is equivalent to the ML estimator in large samples.

For Penalized Maximum Likelihood (PML) estimation of MSB parameters based on the Jeffreys invariant prior, the penalized log-likelihood function ${\tilde{\ell }}_n$ is given by

$$\begin{array}{c}\hfill {\tilde{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})={\ell }_n\left(\mathit{\theta }\right|\mathit{y})+{\displaystyle \frac{1}{2}}log\left|{\mathcal{E}}_n\left(\mathit{\theta }\right)\right|\end{array}$$

(15)

where $log|·|$ returns the logarithm of the determinant of its matrix argument and ${\mathcal{E}}_n$ denotes the Fisher (expected) information matrix: ${\mathcal{E}}_n\left(\mathit{\theta }\right)=Var\left[{\displaystyle \frac{\partial {\ell }_n\left(\mathit{\theta }\right|\mathit{y})}{\partial \mathit{\theta }}}\right]$ [83]. The PML estimators of the MSB model parameters are obtained by solving the penalized score equation ${\dot{\tilde{\ell }}}_n=\mathit{0}$. The Fisher information ${\mathcal{E}}_n\left(\mathit{\theta }\right)$ of the MSB model as well as the penalized score ${\dot{\tilde{\ell }}}_n\left(\mathit{\theta }\right|\mathit{y})={\displaystyle \frac{\partial \tilde{\ell }\left(\mathit{\theta }\right|\mathit{y})}{\partial \mathit{\theta }}}$ are given in Appendix B. For inference in the MSB model, we implemented both the ML and PML procedures in the open source library msbreg (available at https://github.com/Chenangnon/msbreg) in R freeware [84].

2.6. Testing for Unit Asymptote

One important feature of the proposed MSB model framework is the ability to infer the absence or not of unknown stages that significantly affect the binomial outcome. Indeed, when there is no controlled experimental condition making the maximum success probability ${\lambda }_i$ vary across experimental units, the MSB model includes a multiplicative intercept $\lambda$ such that ${\lambda }_i=\lambda$ for all units $i=1,2,\cdots ,n$ (i.e. the design matrix for ${\lambda }_i$ is a vector of all ones: $\mathbf{V}={\mathit{J}}_n$). In this situation, the question arises of whether confounding variables significantly affect the binomial outcome at unknown stages.

The absence or presence of unknown stages can be inferred via a statistical test with the null hypothesis $H_0:\lambda =1$ (unit asymptote). For instance, the nesting relationships between TB and MSB, or ZIB and MSB models allow to always check whether a dataset at hand supports $\lambda =1$ or not. For $q=1$, we can test the TB ($\lambda =1$) model against the more general MSB model ($\lambda <1$). Likewise, when $q=2$, we can test the ZIB model (with subject-dependent susceptibility probability and $\lambda =1$) against the general MSB model ($\lambda <1$). More generally ($q\ge 1$), testing for a unit asymptote will help achieve a parsimonious fit: when $H_0$ is rejected, it motivates identifying additional stages (predictors) for the binomial outcome.

Testing the null hypothesis of a unit asymptote $H_0$ can be achieved, for instance, via an LR test assuming a large sample $chi$-$square$ approximation to the sampling distribution of the LR statistic under $H_0$. The LR statistic is given by

$$\begin{array}{c}\hfill R_{\ell }\left(\hat{\mathit{\theta }},{\hat{\mathit{\theta }}}_0|\mathit{y}\right)=2\left[{\ell }_n\left(\hat{\mathit{\theta }}\right|\mathit{y})-{\ell }_n\left({\hat{\mathit{\theta }}}_0\right|\mathit{y})\right]\end{array}$$

(15)

where $\hat{\mathit{\theta }}$ is the full ML estimate and ${\hat{\mathit{\theta }}}_0$ is the ML estimate under the restriction $H_0$. Since $H_0$ lies on the boundary of the parameter space ($\lambda =1$), the naive large sample ${\chi }_1^2$ approximation to the null distribution of $R_{\ell }$ lacks power [85,86]. Indeed, the appropriate large sample null distribution is the so-called $chi$-$bar$-$square$ distribution ${\overline{\chi }}_{(0,1)}^2=0.5·{\chi }_0^2+0.5·{\chi }_1^2$, i.e. the mixture of ${\chi }_0^2$ (point mass at zero) and ${\chi }_1^2$ distributions, each with weight 0.5 [86,87,88,89]. In practice, this simply implies multiplying the p-value obtained from the ${\chi }_1^2$ distribution by 0.5.

Under PML inference, the unit asymptote null hypothesis $H_0$ can be tested using the corresponding penalized LR statistic $R_{\ell }\left(\tilde{\mathit{\theta }},{\tilde{\mathit{\theta }}}_0|\mathit{y}\right)$ where $\tilde{\mathit{\theta }}$ and ${\tilde{\mathit{\theta }}}_0$ are respectively the full and restricted PML estimates of $\mathit{\theta }$. Note that because the Jeffreys prior considered for PML inference is asymptotically negligible, the LR statistic remains valid for testing $H_0$ in large samples under either ML or PML estimation [90]. In small samples ($n<100$), one can consider a bootstrap approximation of the null distribution of the LR statistic under $H_0$ (see e.g. [91]).

3. Simulation Studies

In this section, we showcase the value of the proposed MSB model, starting with simulated case studies, and then exploring important statistical aspects including bias and accuracy, sample size determination, confidence interval and significance of model parameters. For demonstration purposes, we consider a simulated biological case study where the interest is in the viability of a mutant protein resulting from a single amino acid substitution. The outcome variable Y, which is binary (viable/unviable), depends upon the protein’s ability to fold into a native conformation, and its ability to bind other subunits or ligands. Free energy changes ($\Delta \Delta G_{fold}$, denoted $X_1$ and $\Delta \Delta G_{bind}$, denoted $X_2$) indicate the extent to which a mutation destabilizes protein folding and binding, and we model how they affect the outcome Y.

3.1. Simulated Examples on Mutated Protein Viability

3.1.1. Example 1: Simple Logistic Regression Analysis

For simplicity in the presentation, we start with a model having only one observed predictor: free energy of folding ($X_1$). We use one simulated dataset (sample size $n=300$) generated from the logit link model:

$$\begin{array}{cc}\hfill Y& \stackrel{ind}{\sim }\mathcal{B}er\left(\mu \right)\hfill \\ \hfill \mathrm{where}\mu & =\lambda ·h\left({\beta }_0+{\beta }_1{X}_1\right),\hfill \end{array}$$

(16)

h is the inverse logit link function defined in (3), ${\beta }_0=3$, ${\beta }_1=-2$, $X_1$ is generated from $\mathcal{U}nif(-1,4)$ (i.e. a uniform distribution on $[-1,4]$), and $\lambda =0.75$. The parameter $\lambda <1$ means that there are unmeasured variables (e.g. $X_2$) that affect the protein viability curve.

The scatter plot of the simulated data and the underlying true viability curve are displayed on Figure 1. The traditional approach in the face of the dataset is to fit a Standard Logistic Regression (SLR) model to estimate the parameters ${\beta }_0$ and ${\beta }_1$, to obtain the viability curve as a function of the predictor $\Delta \Delta G_{fold}$. This fit results in the biased estimates ${\hat{\beta }}_0=0.93$ and ${\hat{\beta }}_1=-1.14$ (Table 1). The corresponding fitted viability curve (Figure 1) shows large discrepancies with the true viability curve.

The SLR model fit wrongly suggests for instance that the protein viability can approach one as the folding stability $\Delta \Delta G_{fold}$ decreases below $X_1=-1$$kcal/mol$ (whereas the true viability plateaus at $\lambda =0.75$). The fit also incorrectly indicates that varying the folding stability $\Delta \Delta G_{fold}$ from $-1$ to 0 $kcal/mol$ will result into important changes in the protein viability whereas the true viability is nearly constant in this region. Clearly, mis-specifying the SLR model for the Multistage Logistic Regression (MLR) model can lead to biased estimates of regression parameters and probability curve, and ultimately to incorrect biological interpretations.

From both our understanding of protein viability and the scatter in Figure 1, we know or at least suspect that the expected protein viability does not approach one as $X_1$ decreases continuously below $X_1=-1$$kcal/mol$. The MLR model approach can thus yield more realistic estimates. We note that the MLR regression parameter estimates shown in Table 1 are closer to the true parameter values, ${\hat{\beta }}_0=2.95$ and ${\hat{\beta }}_1=-1.94$. Along with the estimate of the asymptote $\hat{\lambda }=0.76$, these estimates produce a fitted curve showing only mild discrepancies with the true viability curve (Figure 1). We can measure the overall discrepancy between the fitted and true curves using the mean absolute deviation between the predicted viability and the true viability, $mad\left(\hat{\mu }\right)$, and their squared Pearson correlation, $R^2\left(\hat{\mu }\right)$. We find $mad\left(\hat{\mu }\right)=0.01$ and $R^2\left(\hat{\mu }\right)=1.00$ for the MLR model fit, whereas the SLR model fit gives $mad\left(\hat{\mu }\right)=0.06$ and $R^2\left(\hat{\mu }\right)=0.95$. This indicates that the MLR model predicted viability is on average 1% below or above the true viability, versus 6% for the SLR model predictions. By allowing one more parameter ($\lambda$), the MLR model provides a better fit to the data.

In practice, the true viability curve is unknown, hence $mad\left(\hat{\mu }\right)$ or $R^2\left(\hat{\mu }\right)$ could not be computed. However, the Akaike’s Information Criterion (AIC) is smaller for the MLR model fit (AIC = 276.30) as compared to the SLR model fit (AIC = 281.43), indicating that the MLR model fit is better. The estimate of the asymptote ($\hat{\lambda }=0.76$) has the approximate 95% confidence interval $C{I}_{95\%}=[0.63,0.86]$. The LR test comparing the fits of the MLR model and the SLR model (which is nested into the MLR model) indicates with 95% confidence that the improvement achieved by estimating one additional parameter in the MLR model is significant (${\overline{\chi }}_{(0,1)}^2=7.16$, $P=0.0037$) and is thus worth the lost of one degree of freedom. Along with the AIC statistic and the 95% confidence interval for $\lambda$, this test provides evidence for $\lambda <1$, suggesting there are important but missing stages or covariates. Indeed, one advantage of the MLR model framework over the SLR model is the ability to infer or rule out the absence of important stages or covariates along the process generating the data.

Since the results discussed above are based on a unique dataset simulated from a MLR setting, they do not indicate that the MLR fit is always better than the SLR fit, even when the data comes from a MLR process. They nevertheless suggest that one always ponder the SLR model assumption which requires the success probability of the response to arbitrarily approach one under some conditions only dependent upon measured covariates. If the physical or biological process under study may not agree with such an assumption, it is then worthwhile comparing the SLR model and MLR fits and selecting the best model using an objective statistical procedure (e.g. AIC, LR test).

In order to make comparisons robust to randomness in a single dataset, we computed the relative bias, defined as

$$\begin{array}{cc}\hfill bias\left({\hat{\beta }}_j\right)& =100\times {\displaystyle \frac{E\left[{\hat{\beta }}_j\right]-{\beta }_j}{|{\beta }_j|}}\hfill \end{array}$$

and the accuracy, measured by the relative root mean square error,

$$\begin{array}{cc}\hfill rmse\left({\hat{\beta }}_j\right)& =100\times {\displaystyle \frac{\sqrt{E\left[{({\hat{\beta }}_1-{\beta }_1)}^2\right]}}{|{\beta }_j|}}.\hfill \end{array}$$

The expectations in the $bias$ and $rmse$ formulae are approximated using averages over Monte Carlo replicates of the dateset. To ensure that the standard error of the approximated $bias$ and $rmse$ are less than $0.5\%$, we estimated the minimum number of Monte Carlo repetitions to be $95\times (100-95)/(0.5^2)=1900$ [93]. We thus used 2000 dateset replicates.

As it is expected under a MLR data model, the MLR estimates are less biased and more accurate than logistic estimates, resulting in smaller discrepancies between true and fitted viability curves for the MLR fits. Indeed, replicating the working example 2000 times, we found $bias\left({\hat{\beta }}_0\right)=3\%$, $rmse\left({\hat{\beta }}_0\right)=30\%$, $bias\left({\hat{\beta }}_1\right)=-2\%$, and $rmse\left({\hat{\beta }}_1\right)=21\%$ for the MLR regression coefficient estimators. The corresponding performance measures for the SLR model estimators are $bias\left({\hat{\beta }}_0\right)=-71\%$, $rmse\left({\hat{\beta }}_0\right)=71\%$, $bias\left({\hat{\beta }}_1\right)=4\%$, and $rmse\left({\hat{\beta }}_1\right)=44.28\%$, showing a switch of the directions of the biases and larger absolute value of $bias$ and $rrmse$. In addition, the expected $mad$ and $R^2$ are $mad\left(\hat{\mu }\right)=0.03$ and $R^2\left(\hat{\mu }\right)=0.99$ for the MLR fits, versus $mad\left(\hat{\mu }\right)=0.06$ and $R^2\left(\hat{\mu }\right)=0.94$ for the SLR fits.

3.1.2. Example 2: Two-Stage Logistic Regression Analysis

The regression analysis in Example 1 is based on only one predictor. We consider here a model for protein viability with free energy of folding ($X_1$) and binding ($X_2$). In this case, the MSB model does not reduce to any other known binomial model. For demonstration, we again use one simulated dataset ($n=300$) generated from the MLR model:

$$\begin{array}{cc}\hfill Y& \stackrel{ind}{\sim }\mathcal{B}er\left(\mu \right)\hfill \\ \hfill \mathrm{where}\mu & =\lambda ·h\left({\beta }_{01}+{\beta }_1{X}_1\right)·h\left({\beta }_{02}+{\beta }_2{X}_2\right),\hfill \end{array}$$

${\beta }_{01}={\beta }_{02}=3$, ${\beta }_1=-2$, ${\beta }_2=-1$, $X_1\sim \mathcal{U}nif(-1,4)$ and $X_2\sim \mathcal{U}nif(-3,7)$ are independent, and $\lambda =0.80$. Here, $\lambda <1$ stands for further (unmeasured) variables (e.g. residue sites) that affect the viability of the mutated protein.

Fitting a traditional logistic regression model (i.e. two predictors with additive-effects) gives ${\hat{\beta }}_0=0.68$ ($C{I}_{95\%}=[0.25,1.12]$) and ${\hat{\beta }}_1=-0.97$ ($C{I}_{95\%}=[-1.23,-0.70]$) and ${\hat{\beta }}_2=-0.36$ ($C{I}_{95\%}=[-0.48,-0.25]$) with AIC $=259.80$. This fit results in predicted success probabilities $\hat{\mu }$ with mean absolute deviation from the truth $mad\left(\hat{\mu }\right)=0.07$ and predictive power $R^2\left(\hat{\mu }\right)=0.89$. Under a multiplicative-effect assumption based on the biology of protein viability, the fitted Zero-Inflated Logistic (ZIL) and MLR models are shown in Table 2. As expected, the ZIL and MLR models provide better fits to the data (higher deviance reduction ratio $R_{dev}^2$, lower AIC) as compared to the traditional logistic fit. It appears that the MLR model is fitter than the ZIL model based on its lower AIC and higher $R^2$. The $95\%$ confidence interval for $\lambda$ ($C{I}_{95\%}=[0.63,0.91]$) as well as the LR test for the ZIL fit ($\lambda =1$) versus the MLR fit ($0<\lambda <1$) also point to the choice of the MLR model (${\overline{\chi }}_{(0,1)}^2=14.06$, $P<0.001$).

Beyond this one dataset simulation, the MLR model estimates are less biased and more accurate as compared to the ZIL estimates, as expected for data generated from a MLR process. Indeed, from Table 3–which shows the mean estimates, relative biases, root mean square errors of parameter estimates, and predictive power (computed from ZIL and MLR model fits on 2000 replicates of the MLR data)– we notice that the MLR model estimates have much less bias than ZIL estimates, and slightly lower $rmse$.

3.2. General Performance in Synthetic Data

In the examples discussed in Section 3.1, we have considered $\lambda \in \left\{0.75,0.80\right\}$, intercept equal to 3 and slopes in $\left\{-1,-2\right\}$, set the sample size to $n=300$, and used a uniform distribution for the predictors. Our goal here is to assess the behavior of the PML estimators of the MSB model when these parameters vary. We considered regression intercepts in $\left\{3,8\right\}$ and slopes in $\left\{-1,-2,-3\right\}$, and the predictor distributions ($PD$) include normal, and log-normal in addition to uniform. We generally investigated sample sizes $n\in [100,2000]$ and asymptote values $\lambda \in [0.5,1]$, although we occasionally used n up to $n=10000$ and $\lambda$ down to $\lambda =0.25$ to illustrate some specific trends.

We summarize the key results on the performance of PML inference in simulated MSB data. We consider aspects such as the bias and accuracy of the PML estimators, their implication for sample size determination under an MSB design, the ability of the LR test procedure to detect the absence of unknown stages ($\lambda <1$), and the coverage and significance probabilities of 95% confidence interval for a regression slope. Detailed simulation design and results are available in Appendix C.

3.2.1. Bias and Accuracy

Figure 2 and Figure 3 show, respectively, the relative $bias$ and the $rmse$ of the PML estimates of $\lambda$, ${\beta }_0$ and ${\beta }_1$ as functions of the true asymptote $\lambda$, the sample size n, and the $PD$ using data from the one stage model (17) with ${\beta }_0=3$ and ${\beta }_1=-2$. It appears that for the estimates $\hat{\lambda }$, ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$, both the relative $bias$ and $rmse$ decrease and approach zero as n increases

• Bias in $\hat{\lambda }$. At $\lambda =1$, the relative $bias$ of $\hat{\lambda }$ is close to zero for all considered $PD$ and sample sizes $n\in [100,2000]$. As $\lambda$ decreases, the estimate $\hat{\lambda }$ exhibits a slight $bias$, reaching a maximum of about 10% in absolute value at $n=100$ when $\lambda =0.5$ (Figure 2). We notice a pattern in the $bias$: under a uniform $PD$ for instance, the $bias$ is negative at $\lambda =0.9$ and positive for $\lambda \in \left\{0.5,0.75\right\}$. Additional simulations with $\lambda$ values down to 0.25 (see Appendix C) reveal that under a uniform or log-normal $PD$, the $bias$ in $\hat{\lambda }$ is negative (under-estimation) for high $\lambda$ values ($0.75<\lambda <1$), essentially zero for medium $\lambda$ values ($0.5\le \lambda \le 0.75$), and positive (over-estimation) for low $\lambda$ values ($0.25\le \lambda <0.5$). For the normal $PD$, the pattern persists, but the $bias$ is negative for $0.8<\lambda <1$.
• Accuracy of $\hat{\lambda }$. At $\lambda =1$, the relative $rmse$ of $\hat{\lambda }$ is close to zero for all $PD$ and sample sizes. The accuracy also decreases with $\lambda$ (Figure 3). Indeed, under the uniform $PD$ with $n=100$, the $rmse$ of $\hat{\lambda }$ is about 6% for $\lambda =0.9$, but reaches about 21% for $\lambda =0.5$. We observe a similar trend under the log-normal $PD$, with a $rmse$ around 5% at $\lambda =0.9$, and 19% at $\lambda =0.5$ ($n=100$). The estimate $\hat{\lambda }$ is comparatively less accurate under a normal $PD$, with a $rmse$ around 10% at $\lambda =0.9$, and 25% at $\lambda =0.5$ ($n=100$).
• Bias and accuracy of ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$. At $\lambda =1$ (all $PD$), the $bias$ of regression parameter estimates ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$ is also close to zero (Figure 2), and the $rmse$ (which is then essentially the standard deviation for each estimator) is about $15\%$ at $n=300$, $10\%$ at $n=500$, and $5\%$ at $n=1500$ (Figure 3). As $\lambda$ decreases, the relative $rmse$ increases, reaching $32\%$ for ${\hat{\beta }}_0$ and $22\%$ for ${\hat{\beta }}_1$ when $\lambda =0.5$ and $n=500$. This indicates that the lower $\lambda$, the larger the sample size required to achieve the same level of accuracy for regression coefficients.
• Effect of predictor distribution. The distribution of the predictor does not affect the accuracy of the estimates ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$ for $\lambda$ close to one (Figure 3). But as $\lambda$ decreases, the relative $rmse$ is ceteris paribus minimum under uniform $PD$ and maximum under log-normal $PD$ for both ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$. At $\lambda =0.5$ for instance, $rmse\left({\hat{\beta }}_1\right)=15\%$ at $n=1000$ under a uniform predictor. In contrast, $rmse\left({\hat{\beta }}_1\right)=18\%$ at $n=2000$ under a log-normal predictor.
• Impact of regression coefficients. The trends in the relative $bias$ and $rmse$ shown for ${\beta }_0=3$ and ${\beta }_1=-2$ (Figure 2 and Figure 3) are similar across combinations of ${\beta }_0\in \left\{3,8\right\}$ and ${\beta }_1\in \left\{-1,-2,-3\right\}$. Specifically, the $rrmse$ of the estimates ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$ increase with $\lambda$, and for $\lambda <1$, the $rrmse$ is ceteris paribus minimum under uniform $PD$ and maximum under log-normal $PD$ (see Appendix C).
• Beta regression on the accuracy of ${\hat{\beta }}_1$. To summarize the general trends for the slope estimate ${\hat{\beta }}_1$, we fitted a beta regression to the values of the relative $rmse$ of ${\hat{\beta }}_1$ as a function of the simulation factors. Table 4 presents the results of the fit which explain about 87% of the variability in $rmse\left({\hat{\beta }}_1\right)$. From Table 4, we observe that the leading factors in the variations of $rmse\left({\hat{\beta }}_1\right)$ are $\lambda$, $PD$ and ${\beta }_0$, with significant interactions. For instance, decreasing $\lambda$ by 0.1 results on average in a 6.5% increase in $rmse\left({\hat{\beta }}_1\right)$ under a uniform $PD$ and a small sample ($n<200$). Switching then to a log-normal $PD$, we have on average a 16% increase in $rmse\left({\hat{\beta }}_1\right)$ for small $\lambda$ values, and about a 26% increase for $\lambda$ values close to one. Under scenarios with ${\beta }_0=8$, we observe an average 7% decrease in $rmse\left({\hat{\beta }}_1\right)$ as compared with ${\beta }_0=3$.

The presence of a significant negative interaction between n and $\lambda$ (Table 4) indicates that larger sample sizes are generally required when $\lambda <1$ as compared to when $\lambda =1$. Indeed, the $rmse$ generally decreases with n, but the interactive effect of $\lambda$ with the sample size n is also negative and relatively large (twice the main effect of n) as compared to the marginal effect of n (Table 4). For instance, increasing n by 100 would result in a 5.5% average decrease in $rmse\left({\hat{\beta }}_1\right)$ if $\lambda =0.5$, versus a 8.3% average decrease if $\lambda =1$.

3.2.2. Sample Size Determination

The trends in the accuracy of PML estimates indicate that the sample size in an MSB design depends on $\lambda$. To illustrate the importance of such trends for sample size calculation when considering an MSB framework for data analysis, we run additional simulations with sample sizes up to $n=10000$ under a uniform $PD$ with true slope ${\beta }_1=-2$ and intercept ${\beta }_0=3$. Figure 4 shows for different $\lambda$ values, the sample size n required to achieve a specified root mean square error for the estimate ${\hat{\beta }}_1$.

It appears that for $\lambda =1$, a sample size of about $n=393$ is required to achieve $rmse\left({\hat{\beta }}_1\right)=10\%$. Reaching the same level of accuracy with $\lambda =0.9$ requires a sample of size $n=853$. The required sample size becomes $n=1313$ for $\lambda =0.75$, and reaches $n=2417$ when $\lambda =0.5$. Hence, in this particular simulation setting, halving $\lambda$ demands about six times more samples to achieve the same 10% level of accuracy. To achieve $rmse\left({\hat{\beta }}_1\right)=5\%$, $n=1425$ samples are required when $\lambda =1$, $n=3609$ when $\lambda =0.9$, $n=5508$ when $\lambda =0.75$ and $n=9949$ when $\lambda =0.5$. Hence for a 5% level of accuracy, halving $\lambda$ demands about seven times more samples.

We observe similar trends under normal and log-normal $PD$ (see details in Appendix C). For instance, under the normal $PD$, to achieve $rrmse\left({\hat{\beta }}_1\right)=10\%$, about $n=480$ samples are required when $\lambda =1$. With $\lambda =0.5$, the required sample size becomes $n=3819$, that is, almost eight fold increase of the sample size requirement. Under the log-normal $PD$, about $n=497$ samples are required to achieve $rrmse\left({\hat{\beta }}_1\right)=10\%$ when $\lambda =1$. With $\lambda =0.5$, the required sample size becomes $n=5779$, that is, 11.6 fold increase of the sample size requirement.

3.2.3. Testing for Unit Asymptote

Figure 5 shows the rejection rate of the LR test of the null hypothesis ($H_0:\lambda =1$) at the significance level of 5% for true $\lambda$ values in $[0.5,1]$ and regression slopes ${\beta }_1=-3,-1$. It appears that when the TB model is the truth ($\lambda =1$) and the sample size is large ($n=2000$), the rejection rate of the LR test is close to the nominal level of 5% (5.5% for uniform $PD$, 4.5% for normal $PD$ and 4.9% for log-normal $PD$). This shows that the LR test detects the TB model when it is the correct model and sample size is large. Interestingly, this remains true in small sample with rejection rates under 6.7% (largest value observed under log-normal $PD$ with ${\beta }_1=-3$).

For true $\lambda \le 0.9$, we observed on Figure 5 that the power of the LR test is 1 under uniform $PD$ in samples of size $n\ge 1000$. In smaller samples, the power decreases with decreasing sample size, but increases as the true $\lambda$ decreases, reaching 59% at $\lambda =0.5$. In other words, the LR test becomes conservative in small samples. We also notice that the rejection rates are lower under normal and log-normal $PD$ as compared to the uniform $PD$. Finally, the rejection rates are lower for ${\beta }_1=-3$ as compared to ${\beta }_1=-1$ when the $PD$ deviates from uniform.

3.2.4. Statistical Inference on Regression Coefficients

We present the finite sample properties of inference based on the large sample normal approximation of the distributions of PML estimates under the MSB model. Figure 6 shows the Coverage Probability (CP) of approximate 95% confidence interval using the inverse negative Fisher information matrix as covariance matrix for PML estimates. We observed that the CP is close to the nominal level (95%) and stays above 90% when either $\lambda =1$, or the sample size is large, or the predictor distribution is uniform or normal. For a log-normal predictor however, the CP decreases with $\lambda$, reaching 88% at $n=1000$ for $\lambda =0.5$.

Figure 7 displays the Significance Probability (SP) of the approximate 95% confidence interval. It appears that, irrespective of the $PD$ and $\lambda$ values, the significance of the slope ${\beta }_1$ is close to 1 for any $n\ge 200$. The intercept coefficient ${\beta }_0$ is more sensitive to both $\lambda$ and sample size, the SP decreasing and reaching 0.22 under normal $PD$ at $n=100$.

4. Discussion

In this paper, we propose a binomial regression framework for modeling binary data arising from processes with multiplicative steps. The use of the model is prescribed when modeling a complex binary process with no data on the detailed discrete stages of the phenomenon. As soon as we have data on each specific stage, a multivariate binomial GLM [95,96] becomes more appropriate for the data analysis. We derived conditions for the identifiability of model parameters. Not surprisingly, it is sufficient that each stage has a continuous predictor only present at that stage to ensure identifiability. The maximum likelihood estimates of the MSB model parameters are asymptotically unbiased, but subject to bias in small samples. To reduce the small sample while addressing separation issues, we developed a penalized maximum likelihood estimation procedure based on Jeffreys invariant prior. We have demonstrated the use context and usefulness of the modeling framework using synthetic data on mutant protein viability (see Figure 1).

Interestingly, both the standard binomial GLM and the ZIB model are special instances of the MSB model. Indeed, though it is not usually described as based on a multiplicative risk process, the ZIB model is based on (2), but limited to $q=2$. It is apparent from (2) that the overall success probability $\mu$ of the binary outcome can be very small. It follows that data generated from an MSB model may have a larger occurrence of zeros as compared to the corresponding standard binomial model. This observation immediately points to the ZIB model which has become popular for modeling binary data with an excess of zeros relative to the expectation under the standard binomial model [63,64,65,66]. The ZIB model stems from a mixture of two distinct subpopulations where one subpopulation is subject to a risk factor (group of “susceptible” units) whereas the other subpopulation is not (group of “immune” units). This can be mathematically described as a sequence of two Bernoulli processes. A first Bernoulli process assigns each individual to one of the two subpopulations. Then, the susceptible subpopulation, when exposed, is subject to a risk through a second Bernoulli process, while the immune subpopulation is not subject to any risk, even if exposed. The first stage of the process is considered as a source of nuisance and is not of interest to the data analysis [63]. It is here important to notice that the MSB model (8) is not a zero-inflation model. The model simply restricts the success probability of the binomial outcome $Y_i$ to the subset $({\alpha }_i{\lambda }_i,{\lambda }_i)$ of $(0,1)$, equally allowing for inflation of both zeros and ones.

The asymptote parameter $\lambda$ included in the MSB model (maximum success probability for the binary outcome) can capture unknown and missing stages. In particular, when only one stage is known for the process under investigation, $\lambda$ provides an estimate of the importance of unknown or unaccounted stages or components in the determination of the binomial outcome. In health data context for instance, $1-\lambda$ can represent a cure fraction [63], or a maximum survival probability. The standard binomial model would over estimate survival or risk if we omit $\lambda$, leading to poor prediction (Figure 1), and possibly to poor policy design. We investigated the finite sample properties of a non-standard likelihood ratio test for detecting $\lambda =1$. The test have shown rejection rates (of the assumption $\lambda =1$) close to the nominal level when $\lambda =1$ (Figure 5). The test is generally conservative and will point to $\lambda <1$ only when the latter is strongly supported by the data. Hence, when studying a phenomenon without a strong justification that the process is additive, it is worthwhile testing the MSB model against the standard binomial model to identify which model is supported by the available data. Rejecting $\lambda =1$ would suggest the absence of some important stages or covariates.

We investigated the finite-sample properties of the penalized maximum likelihood estimates via simulations. It turns out that $\lambda$ affects the estimation of regression coefficients. In particular, lower $\lambda$ values are associated with higher mean square errors in regression coefficient estimates. This implies that larger samples are required for lower $\lambda$ values (Figure 4). Intuitively, $\lambda <1$ means more failures in the binary response as compared to $\lambda =1$. However, unlike the failures under $\lambda =1$ which provide as much information on regression parameters as successes, the additional failures due to $\lambda <1$ do not provide any information on regression parameters. Stated otherwise, a lower $\lambda$ means that each data point is less informative as compared to when $\lambda =1$. Hence the effective size of the sample contributing to the estimation of regression parameters is only a small proportion of n when $\lambda$ is very low. We also established that larger sample sizes are required under log-normal distribution as compared to uniform or normal distribution of predictors.

We finally point to a few open questions that deserve attention. First, a rigorous theory of the existence, consistency and asymptotic normality of the (penalized) maximum likelihood estimators of MSB model parameters should be established. A related aspect of interest is the identifiability and inference when predictor dimension is much larger than the sample size. Second, the bias in parameters estimates relates to the particular parameterization of the MSB model. Nevertheless, first order bias corrected estimates can perform substantially better (in terms of accuracy, i.e. mean square error) than crude maximum likelihood estimates [97]. Although the use of the Jeffreys invariant prior as likelihood penalty ensures finite estimates and reduced bias, it does not completely remove first order bias from the maximum likelihood estimates of the MSB model parameters. Further bias reduction is thus possible using recent general purpose methods developed for nonlinear models (see e.g [98]). Third, diagnostic statistics such as influence measures [99] should be developed for real data analysis under the MSB model framework. Fourth, the power of the test for unit asymptote may be improved by implementing the Bartlett correction [100] for the likelihood ratio statistic. Moreover, the score test, an alternative to likelihood ratio test, may offer a faster procedure for model selection since it only requires fitting the null model (i.e. under $H_0$: $\lambda =1$) [101]. Last, but not least, we assumed that measurement error standard deviations are known in this work. Future research should investigate to what extent the proposed conjugate measurement error distribution is robust to distribution mis-specification.

Appendix A.1. Proof of Proposition 1

Proof. 

Each factor ${\xi }_{ij}$ in (4b) can be considered as the success probability of the independent Bernoulli variable $Y_{ij}$ associated with stage j given a covariate vector ${\mathit{w}}_{ij}$, i.e.,

$$\begin{array}{cc}\hfill \left(Y_{ij}\right|{\mathit{W}}_{ij}={\mathit{w}}_{ij})& \stackrel{ind}{\sim }\mathcal{B}er\left({\xi }_{ij}\right).\hfill \end{array}$$

(A1)

From (A1), the Bernoulli variable $Y_{ij}$ can further be represented by truncating a continuous distribution, e.g. truncated normal, logistic or scale mixture of skew normal variables [102]. In this framework, $Y_{ij}$ has the stochastic representation

$$\begin{array}{c}\left(Y_{ij}\right|Z_{ij}=z_{ij},W_{ij}=w_{ij})\stackrel{d}{=}{I}_{(0,\infty )}({\zeta }_{ij}+z_{ij})\\ \hspace{1em}\left(Z_{ij}\right|W_{ij}=w_{ij})\sim H(0,{\sigma }_o^2)\end{array}$$

(A2)

where $\stackrel{d}{=}$ means “equal in distribution to”, and H denotes the distribution with cdf h. For the logit link for instance, $H(\eta ,{\sigma }_o^2)$ is the logistic distribution with mean $\eta$ and variance ${\sigma }_o^2={\pi }^2/3$; and in the case of the probit link, $H(\eta ,{\sigma }_o^2)$ is the normal distribution with mean $\eta$ and variance ${\sigma }_o^2=1$. From Equation (6), we have ${\mathit{x}}_{ij}={\mathit{w}}_{ij}+{\mathit{e}}_{ij}$. It follows that ${\zeta }_{ij}+z_{ij}$ is a function of ${\mathit{e}}_{ij}$ as ${\zeta }_{ij}+z_{ij}={\eta }_{ij}+\left(z_{ij}-{\mathit{e}}_{ij}^{\top }{\mathit{\beta }}_j\right)$. It then appears that the integral in (9) can be simplified by choosing the distribution $ME\left({\mathit{S}}_{ij}\right)$ of the measurement error ${\mathit{e}}_{ij}$ such that the distribution of $z_{ij}-{\mathit{e}}_{ij}^{\top }{\mathit{\beta }}_j$ is in the same family as $H(0,{\sigma }_o^2)$. Indeed, since $Z_{ij}$ and ${\mathit{e}}_{ij}$ are independent, $Z_{ij}-{\mathit{e}}_{ij}^{\top }{\mathit{\beta }}_j$ has variance $s_{ij}^2={\sigma }_o^2+{\overline{\sigma }}_{ij}^2$ with ${\overline{\sigma }}_{ij}^2={\mathit{\beta }}_j^{\top }{\mathit{S}}_{ij}{\mathit{\beta }}_j$, and since both $Z_{ij}$ and ${\mathit{e}}_{ij}$ have null mean, the representation (A2) becomes

$$\begin{array}{cc}\hfill \left(Y_{ij}\right|{\overline{T}}_{ij}={\overline{t}}_{ij})& \stackrel{d}{=}{I}_{(0,\infty )}({\eta }_{ij}+{\overline{t}}_{ij})\hfill \\ \hfill {\overline{T}}_{ij}& \sim H\left(0,{\overline{\sigma }}_{ij}^2\right).\hfill \end{array}$$

In light of (A2), the distribution H in the stochastic representation must have variance ${\sigma }_o^2$ in order to write e.g. ${\xi }_{ij}=h\left({\zeta }_{ij}\right)$. From ${\eta }_{ij}+{\overline{t}}_{ij}>0\iff {\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}+{\displaystyle \frac{{\overline{t}}_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}>0$, we get by setting $T_{ij}={\displaystyle \frac{{\overline{T}}_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}$:

$$\begin{array}{cc}\hfill \left(Y_{ij}\right|T_{ij}=t_{ij})& \stackrel{d}{=}{I}_{(0,\infty )}\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}+t_{ij}\right)\hfill \\ \hfill T_{ij}& \sim H\left(0,{\sigma }_o^2\right).\hfill \end{array}$$

(A3)

It follows that the expectation of $Y_{ij}|{\mathit{x}}_{ij}$ is given by

$$\begin{array}{c}\hfill {\mu }_{ij}=h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right).\end{array}$$

(A4)

□

Appendix A.2. Proof of Proposition 2

To prove Proposition 2, we make use of the following two Lemmas.

Lemma A1.

Consider the function ${\phi }_1:{\mathbb{R}}^p\to {\mathbb{R}}^p$ defined for $\mathbf{z}\in {\mathbb{R}}^p$ as

$$\begin{array}{c}\hfill {\phi }_1\left(\mathbf{z}\right)={\displaystyle \frac{\mathbf{z}}{\sqrt{1+{(\overline{\sigma }/{\sigma }_o)}^2}}}\end{array}$$

(A5)

where ${\overline{\sigma }}^2={\mathbf{z}}^{\top }\mathit{S}\mathbf{z}$, $\mathit{S}\in {\mathbb{R}}^{p\times p}$ is a semi-positive definite matrix, and ${\sigma }_o>0$ is a positive constant. The function ${\phi }_1$ is injective, and accordingly, for any $\mathbf{a}\in {\mathbb{R}}^p$, the equation ${\phi }_1\left(\mathbf{z}\right)={\phi }_1\left(\mathbf{a}\right)$ has a unique solution $\mathbf{z}=\mathbf{a}$.

Proof. 

Denoting ${\mathit{I}}_p$ the $p\times p$ identity matrix, and setting $\upsilon =\sqrt{1+{(\overline{\sigma }/{\sigma }_o)}^2}$, the Jacobian matrix $\dot{{\phi }_1}\left(\mathbf{z}\right)={\displaystyle \frac{\partial {\phi }_1\left(\mathbf{z}\right)}{\partial {\mathbf{z}}^{\top }}}$ of ${\phi }_1\left(\mathbf{z}\right)$ is given by $\dot{{\phi }_1}\left(\mathbf{z}\right)={\displaystyle \frac{1}{\upsilon }}{\mathit{I}}_p+\mathbf{z}{\displaystyle \frac{\partial }{\partial {\mathbf{z}}^{\top }}}\left({\displaystyle \frac{1}{\upsilon }}\right)$. Further calculus yields

$$\begin{array}{cc}\hfill \dot{{\phi }_1}\left(\mathbf{z}\right)& ={\displaystyle \frac{1}{\upsilon }}\left[{\mathit{I}}_p-{\displaystyle \frac{1}{{\sigma }_o^2{\upsilon }^2}}\mathbf{z}{\mathbf{z}}^{\top }\mathit{S}\right].\hfill \end{array}$$

By the matrix determinant lemma [103], the determinant of $\dot{{\phi }_1}$ is given by

$$\begin{array}{cc}\hfill |\dot{{\phi }_1}\left(\mathbf{z}\right)|& ={\upsilon }^{-p}\left(1-{\displaystyle \frac{1}{{\sigma }_o^2}}{\displaystyle \frac{{\mathbf{z}}^{\top }\mathit{S}\mathbf{z}}{\left(1+\frac{{\mathbf{z}}^{\top }\mathit{S}\mathbf{z}}{{\sigma }_o^2}\right)}}\right)\hfill \\ \hfill & ={\upsilon }^{-p}{\left(1+{\displaystyle \frac{{\mathbf{z}}^{\top }\mathit{S}\mathbf{z}}{{\sigma }_o^2}}\right)}^{-1}={\upsilon }^{-p-2}.\hfill \end{array}$$

Since $\mathit{S}$ is semi-positive definite, ${\mathbf{z}}^{\top }\mathit{S}\mathbf{z}\ge 0$ and $\upsilon >0$ for any $\mathbf{z}\in {\mathbb{R}}^p$. It follows that $|\dot{{\phi }_1}\left(\mathbf{z}\right)|\ne 0$ for all $\mathbf{z}\in {\mathbb{R}}^p$. □

Lemma A2.

Consider the function ${\phi }_2:{\mathbb{R}}^p\to {\mathbb{R}}^p$ defined for $\mathbf{z}\in {\mathbb{R}}^p$ and $\mathbf{z}\ne \mathit{0}$ as

$$\begin{array}{c}\hfill {\phi }_2\left(\mathbf{z}\right)=\left({\mathbf{z}}_0^{\top }\mathbf{z}\right)·\mathbf{z}\end{array}$$

(A6)

where ${\mathbf{z}}_0\in {\mathbb{R}}^p$ with linearly independent components. The function ${\phi }_2$ is injective, and accordingly, for any $\mathbf{z}\in {\mathbb{R}}^p$ and $\mathbf{z}\ne \mathit{0}$, the equation ${\phi }_2\left(\mathbf{a}\right)={\phi }_2\left(\mathbf{z}\right)$ has a unique solution $\mathbf{a}=\mathbf{z}$.

Proof. 

The Jacobian matrix of ${\phi }_2\left(\mathbf{z}\right)$ is given by $\dot{{\phi }_2}\left(\mathbf{z}\right)=\left({\mathbf{z}}_0^{\top }\mathbf{z}\right)·\left({\mathit{I}}_p+{\displaystyle \frac{\mathbf{z}{\mathbf{z}}_0^{\top }}{{\mathbf{z}}_0^{\top }\mathbf{z}}}\right)$. Note that the second term in the square bracket is well defined since ${\mathbf{z}}_0$ have linearly independent components and $\mathbf{z}\ne \mathit{0}$ imply that ${\mathbf{z}}_0^{\top }\mathbf{z}\ne 0$. The determinant of $\dot{{\phi }_2}$ is thus

$$\begin{array}{cc}\hfill |\dot{{\phi }_2}\left(\mathbf{z}\right)|& ={\left({\mathbf{z}}_0^{\top }\mathbf{z}\right)}^p\left[1+{\displaystyle \frac{{\mathbf{z}}_0^{\top }\mathbf{z}}{{\mathbf{z}}_0^{\top }\mathbf{z}}}\right]=2{\left({\mathbf{z}}_0^{\top }\mathbf{z}\right)}^p.\hfill \end{array}$$

It follows that $|\dot{{\phi }_2}\left(\mathbf{z}\right)|\ne 0$ for all non null $\mathbf{z}\in {\mathbb{R}}^p$. □

Proof of Proposition 2. 

Let us denote $(\Omega ,\mathcal{A},\mathbb{P})$ the probability space on which the random vector $\left(Y,{\mathit{X}}_1,\cdots ,{\mathit{X}}_q,\mathbf{u},\mathbf{v}\right)$ is defined. Suppose that $\ell \left(\mathit{\theta }\right)=\ell \left({\mathit{\theta }}^{*}\right)$ almost surely. Then assumptions C1 and C2 imply the existence of a positive constant $c_1$ such that for every ${\mathit{x}}_j\in {\mathbb{C}}_j$, $\mathbf{u}\in {\mathbb{C}}_u$, $\mathbf{v}\in {\mathbb{C}}_v$, and $\mathit{\theta }\in {\mathcal{B}}_p\times {\mathcal{B}}_d\times {\mathcal{B}}_r$, $c_1<\mathbb{P}(Y=1|{\mathit{x}}_1,\cdots ,{\mathit{x}}_q,\mathbf{U},\mathbf{V})<1-c_1$. It follows that there exists $\omega \in \Omega$, outside the region (of measure zero) where $\ell \left(\mathit{\theta }\right)\ne \ell \left({\mathit{\theta }}^{*}\right)$ such that $Y\left(\omega \right)=1$ given that ${\mathbf{X}}_j={\mathit{x}}_j$, $\mathbf{U}=\mathit{u}$ and $\mathbf{V}=\mathit{v}$. For such a $\omega$, the equality $\ell \left(\mathit{\theta }\right)=\ell \left({\mathit{\theta }}^{*}\right)$ reads $\mu ={\mu }^{*}$.

For simplicity in the presentation we define $\varrho =\lambda {\prod }_{j=1}^q{\mu }_j$ so that $\mu$ can be rewritten as $\mu =\alpha \lambda +(1-\alpha )\varrho$. We first show that the quantities $\alpha$, $\lambda$ and $\varrho$ are individually identifiable. In case either $\alpha$ and $\lambda$ are known constants in the MSB model, we only need $\varrho$ to be identifiable. Otherwise, from $\alpha =h\left({\mathit{u}}^{\top }\mathit{\gamma }\right)$ or $\lambda =h\left({\mathit{v}}^{\top }\mathit{\delta }\right)$, the identifiability of $\alpha$ or $\lambda$ is reduced to the identifiability problem in the standard binomial regression model. Specifically, assumptions C1 and C2 ensure that $\alpha ={\alpha }^{*}$ and $\lambda ={\lambda }^{*}$ respectively imply that $\mathit{\gamma }={\mathit{\gamma }}^{*}$ and $\mathit{\delta }={\mathit{\delta }}^{*}$ (see e.g. [104]). As for $\varrho$, $\varrho ={\varrho }^{*}$ reads

$$\begin{array}{cc}\hfill \lambda \prod _{j=1}^{q}h\left({\displaystyle \frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{\sqrt{1+{({\overline{\sigma }}_j/{\sigma }_o)}^2}}}\right)& ={\lambda }^{*}\prod _{j=1}^{q}h\left({\displaystyle \frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{\sqrt{1+{({\overline{\sigma }}_j^{*}/{\sigma }_o)}^2}}}\right)\hfill \end{array}$$

where ${\overline{\sigma }}_j=\sqrt{{\mathit{\beta }}_j^{\top }\mathit{S}{\mathit{\beta }}_j}$, ${\overline{\sigma }}_j^{*}=\sqrt{{\mathit{\beta }}_j^{*\top }\mathit{S}{\mathit{\beta }}_j^{*}}$. Set ${\upsilon }_j=\sqrt{1+{({\overline{\sigma }}_j/{\sigma }_o)}^2}$ and ${\upsilon }_j^{*}=\sqrt{1+{({\overline{\sigma }}_j^{*}/{\sigma }_o)}^2}$. By the positivity of all h values under condition C2, and given that $\lambda ={\lambda }^{*}$, we can write for any $j=1,2,\cdots ,q$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)}{h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}& ={\displaystyle \frac{{\prod }_{k\ne j}h\left(\frac{{\mathit{x}}_k^{\top }{\mathit{\beta }}_k^{*}}{{\upsilon }_k^{*}}\right)}{{\prod }_{k\ne j}h\left(\frac{{\mathit{x}}_k^{\top }{\mathit{\beta }}_k}{{\upsilon }_k}\right)}}.\hfill \end{array}$$

(A7)

Differentiating with respect to $T_j$, a continuous predictor in ${\mathit{x}}_j$ with regression coefficient ${\beta }_{jt}\ne 0$ (or ${\beta }_{jt}^{*}\ne 0$) under condition C3, we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{\frac{{\beta }_{jt}}{{\upsilon }_j}\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)-\frac{{\beta }_{jt}^{*}}{{\upsilon }_j^{*}}h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{{\left[h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)\right]}^2}}& =0\hfill \end{array}$$

(A8)

where $\dot{h}\left(x\right)={\displaystyle \frac{\partial h\left(x\right)}{\partial x}}$ (density function corresponding to the link h). The right hand side of the result is zero because the right hand side of (A7) does not depend on ${\mathit{x}}_j$. The equality (A8) implies

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\beta }_{jt}{\upsilon }_j^{*}}{{\beta }_{jt}^{*}{\upsilon }_j}}& ={\displaystyle \frac{h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}.\hfill \end{array}$$

(A9)

Further differentiation of both sides of Equation (A9) with respect to $T_j$ gives

$$\begin{array}{cc}\hfill 0·{\displaystyle \frac{\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)}}+{\displaystyle \frac{h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)}{h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}{\displaystyle \frac{\frac{{\beta }_{jt}^{*}}{{\upsilon }_j^{*}}\ddot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)-\frac{{\beta }_{jt}}{{\upsilon }_j}\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)\ddot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)}{{\left[\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\right]}^2}}& =0\hfill \end{array}$$

where the factor zero in the first term on the left hand side results from (A8), and the right hand side is zero because the left hand side of (A9) does not depend on ${\mathit{x}}_j$. The last equality implies

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\beta }_{jt}{\upsilon }_j^{*}}{{\beta }_{jt}^{*}{\upsilon }_j}}& ={\displaystyle \frac{\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\ddot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{\ddot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}.\hfill \end{array}$$

(A10)

Equations (A9) and (A10) gives the double equality

$$\begin{array}{c}\hfill {\displaystyle \frac{{\beta }_{jt}{\upsilon }_j^{*}}{{\beta }_{jt}^{*}{\upsilon }_j}}={\displaystyle \frac{h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)h\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}={\displaystyle \frac{\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\ddot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{\ddot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)\dot{h}\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}.\end{array}$$

which is central to our proof of ${\mathit{\beta }}_j={\mathit{\beta }}_j^{*}$. We next consider specific cases of link functions.

• Case of logit link: $h\left(x\right)={\left(1+e^{-x}\right)}^{-1}$
  Here, $\dot{h}\left(x\right)=h\left(x\right)[1-h\left(x\right)]$, hence ${\displaystyle \frac{h\left(x\right)}{\dot{h}\left(x\right)}}={\displaystyle \frac{1}{1-h\left(x\right)}}=1+e^x$. Equation (A9) becomes

  $$\begin{array}{cc}\hfill {\displaystyle \frac{{\beta }_{jt}{\upsilon }_j^{*}}{{\beta }_{jt}^{*}{\upsilon }_j}}& ={\displaystyle \frac{1+exp\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)}{1+exp\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}}.\hfill \end{array}$$

  Differentiating both sides of the last equality with respect to $T_j$ yields ${\left({\displaystyle \frac{{\mathit{\beta }}_j}{{\upsilon }_j}}-{\displaystyle \frac{{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}}\right)}^{\top }{\mathit{x}}_j=0$. Along with the condition C1, this implies that ${\displaystyle \frac{{\mathit{\beta }}_j}{{\upsilon }_j}}={\displaystyle \frac{{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}}$, and by Lemma A1, we get ${\mathit{\beta }}_j={\mathit{\beta }}_j^{*}$.
• Case of probit link: $h\left(x\right)={\int }_{-\infty }^x{\displaystyle \frac{1}{\sqrt{2\pi }}}{e}^{-t^2/2}dt$
  We have $\dot{h}\left(x\right)={\displaystyle \frac{1}{\sqrt{2\pi }}}{e}^{-x^2/2}$, and $\ddot{h}\left(x\right)=-x\dot{h}\left(x\right)$, hence ${\displaystyle \frac{\ddot{h}\left(x\right)}{\dot{h}\left(x\right)}}=-x$. Note that by the condition C1 and the assumption ${\beta }_{jt}\ne 0$, we have ${\mathit{x}}_j^{\top }{\mathit{\beta }}_j\ne 0$. Equation (A10) thus reads

  $$\begin{array}{cc}\hfill {\displaystyle \frac{{\beta }_{jt}{\upsilon }_j^{*}}{{\beta }_{jt}^{*}{\upsilon }_j}}& ={\displaystyle \frac{\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}\right)}{\left(\frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}\right)}},\mathrm{which}\mathrm{is}\mathrm{rearranged}\mathrm{as}{\displaystyle \frac{{\beta }_{jt}}{{\upsilon }_j}}{\displaystyle \frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j}{{\upsilon }_j}}={\displaystyle \frac{{\beta }_{jt}^{*}}{{\upsilon }_j^{*}}}{\displaystyle \frac{{\mathit{x}}_j^{\top }{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}}.\hfill \end{array}$$

  Then, because the assumption ${\beta }_{jt}\ne 0$ gives ${\mathit{\beta }}_j\ne \mathit{0}$ and condition C1 ensures ${\mathit{\beta }}_j^{\top }{\mathit{x}}_j\ne 0$, by Lemma A2, we have ${\displaystyle \frac{{\mathit{\beta }}_j}{{\upsilon }_j}}={\displaystyle \frac{{\mathit{\beta }}_j^{*}}{{\upsilon }_j^{*}}}$, which then gives ${\mathit{\beta }}_j={\mathit{\beta }}_j^{*}$ by Lemma A1. We have overall shown that under the logit or probit link function, $\varrho ={\varrho }^{*}$ implies that ${\mathit{\beta }}_j={\mathit{\beta }}_j^{*}$ for $j=1,2,\cdots ,q$.

It remains to show that $\mu ={\mu }^{*}$ implies that $\alpha ={\alpha }^{*}$, $\lambda ={\lambda }^{*}$ and $\varrho ={\varrho }^{*}$. In case both $\alpha$ and $\lambda$ are known constants in the MSB model, $\mu =\alpha \lambda +(1-\alpha )\varrho$ implies that $\mu$ as a function of $\varrho$ is invertible ($\varrho ={\displaystyle \frac{\mu -\alpha \lambda }{(1-\alpha )}}$ given that $\alpha \in (0,1)$ and $\mu \ge \alpha \lambda$). Hence $\mu ={\mu }^{*}$ implies $\varrho ={\varrho }^{*}$. If both $\alpha$ and $\lambda$ are not known constants, we have for all $j=1,2,\cdots ,q$, ${\displaystyle \frac{\partial }{\partial T_j}}{\displaystyle \frac{{\varrho }_0}{{\varrho }_0^{*}}}=0$ (where ${\varrho }_0={\prod }_{j=1}^q{\mu }_j$) which leads to ${\varrho }_0={\varrho }_0^{*}$ by the above argument used for ${\mathit{\beta }}_j$. In case exactly one of $\alpha$ and $\lambda$ is a known constant (say $\alpha$ is a known constant), the second (say $\lambda$) can be uniquely recovered from knowing $\mu$ and ${\varrho }_0$.

If none of $\alpha$ and $\lambda$ is a known constant, then the condition C3 ensures that there is a continuous predictor $T_{\alpha }$ in the linear predictor of $\alpha$ or $T_{\lambda }$ in the linear predictor of $\lambda$. Then, $\mu ={\mu }^{*}$ reads $\alpha \lambda +(1-\alpha )\lambda {\varrho }_0={\alpha }^{*}{\lambda }^{*}+(1-{\alpha }^{*}){\lambda }^{*}{\varrho }_0^{*}$. Suppose that we have $T_{\alpha }$, then we get ${\displaystyle \frac{\partial }{\partial T_{\lambda }}}{\displaystyle \frac{\lambda }{{\lambda }^{*}}}=0$ and $\delta ={\delta }^{*}$ follows. This results in $\lambda ={\lambda }^{*}$. Since ${\varrho }_0={\varrho }_0^{*}$ and $\lambda ={\lambda }^{*}$, it follows that $\alpha ={\alpha }^{*}$. If we have $T_{\lambda }$, then we get ${\displaystyle \frac{\partial }{\partial T_{\alpha }}}{\displaystyle \frac{\lambda }{{\lambda }^{*}}}=0$ giving $\alpha ={\alpha }^{*}$. From ${\varrho }_0={\varrho }_0^{*}$ and $\alpha ={\alpha }^{*}$, it follows that $\lambda ={\lambda }^{*}$. □

Appendix B.1. Maximum Likelihood

The Fisher (expected) information matrix required in the penalized log-likelihood function (15) is defined as the variance-covariance matrix of the score vector ${\dot{\ell }}_n={\displaystyle \frac{\partial \ell \left(\mathit{\theta }\right|\mathit{y})}{\partial \mathit{\theta }}}$ [83]:

$$\begin{array}{c}\hfill {\mathcal{E}}_n\left(\mathit{\theta }\right)=Var\left[{\dot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})\right].\end{array}$$

(A11)

From (14), the score vector is given by

$$\begin{array}{cc}\hfill {\dot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})& =\sum _{i=1}^n{\displaystyle \frac{y_i-n_i{\mu }_i}{{\mu }_i(1-{\mu }_i)}}{\dot{\mathit{\mu }}}_i\hfill \end{array}$$

(A12)

where ${\dot{\mathit{\mu }}}_i={\displaystyle \frac{\partial {\mu }_i}{\partial \mathit{\theta }}}$ is the vector of first order partial derivatives of ${\mu }_i$ with respect to the elements of the parameter vector $\mathit{\theta }={\left({\mathit{\beta }}_1^{\top },\cdots ,{\mathit{\beta }}_q^{\top },{\mathit{\gamma }}^{\top },{\mathit{\delta }}^{\top }\right)}^{\top }$. From (12), ${\dot{\mathit{\mu }}}_i$ is given by:

$$\begin{array}{cc}\hfill {\dot{\mathit{\mu }}}_i& ={\left({\displaystyle \frac{\partial {\mu }_i}{\partial {\mathit{\beta }}_1^{\top }}},\cdots ,{\displaystyle \frac{\partial {\mu }_i}{\partial {\mathit{\beta }}_q^{\top }}},{\displaystyle \frac{\partial {\mu }_i}{\partial {\mathit{\gamma }}^{\top }}},{\displaystyle \frac{\partial {\mu }_i}{\partial {\mathit{\delta }}^{\top }}}\right)}^{\top }\hfill \end{array}$$

(A13)

where

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mu }_i}{\partial {\mathit{\beta }}_j}}& ={\lambda }_i(1-{\alpha }_i)\left[\prod _{k\ne j}h\left({\displaystyle \frac{{\eta }_{ik}}{\sqrt{1+{({\overline{\sigma }}_{ik}/{\sigma }_o)}^2}}}\right)\right]{\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_j}}(\mathrm{for}j=1,\cdots ,q),\hfill \end{array}$$

(A14a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mu }_i}{\partial \mathit{\gamma }}}& ={\lambda }_i\left[1-\prod _{j=1}^{q}h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\right]\dot{h}\left({\mathbf{U}}_i^{\top }\mathit{\gamma }\right){\mathbf{U}}_i,\hfill \end{array}$$

(A14b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mu }_i}{\partial \mathit{\delta }}}& =\left[{\alpha }_i+(1-{\alpha }_i)\prod _{j=1}^{q}h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\right]\dot{h}\left({\mathbf{V}}_i^{\top }\mathit{\delta }\right){\mathbf{V}}_i\hfill \end{array}$$

(A14c)

with $\dot{h}\left(x\right)={\displaystyle \frac{\partial h\left(x\right)}{\partial x}}$ (density function corresponding to the link h),

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_j}}& =\dot{h}\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right){\mathit{\tau }}_{ij},\mathrm{and}\hfill \end{array}$$

(A14d)

$$\begin{array}{cc}\hfill {\mathit{\tau }}_{ij}& ={\displaystyle \frac{1}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\left[{\mathit{x}}_{ij}-{\displaystyle \frac{{\eta }_{ij}}{{\sigma }_o^2+{\overline{\sigma }}_{ij}^2}}{\mathit{S}}_{ij}{\mathit{\beta }}_j\right].\hfill \end{array}$$

(A14e)

Note that in the absence of measurement errors (${\mathit{S}}_{ij}=\mathit{0}$), ${\mathit{\tau }}_{ij}={\mathit{x}}_{ij}$ so that ${\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_j}}=\dot{h}\left({\eta }_{ij}\right){\mathit{x}}_{ij}$. Since ${\dot{\mathit{\mu }}}_i$ does not depend on $y_i$, and $Var\left[{\displaystyle \frac{y_i-n_i{\mu }_i}{{\mu }_i(1-{\mu }_i)}}\right]={\displaystyle \frac{Var\left[y_i\right]}{{\mu }_i^2{(1-{\mu }_i)}^2}}={\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}$, we obtain the Fisher information matrix

$$\begin{array}{cc}\hfill {\mathcal{E}}_n\left(\mathit{\theta }\right)& =\sum _{i=1}^n{\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}{\dot{\mathit{\mu }}}_i{\dot{\mathit{\mu }}}_i^{\top }.\hfill \end{array}$$

(A15f)

When the ML estimator $\hat{\mathit{\theta }}={\left({\hat{\mathit{\beta }}}^{\top },{\hat{\mathit{\gamma }}}^{\top },{\hat{\mathit{\delta }}}^{\top }\right)}^{\top }$ of $\mathit{\theta }$ exists, large sample inference on $\hat{\mathit{\theta }}$ relies on the asymptotic covariance matrix as usually in GLMs. A candidate for the asymptotic covariance matrix is the inverse of the observed information matrix (evaluated at $\hat{\mathit{\theta }}$) [105], given for the MSB model (8) by ${\Sigma }_n={\left[{\mathcal{I}}_n\left(\hat{\mathit{\theta }}\right|\mathit{y})\right]}^{-1}$ where

$$\begin{array}{cc}\hfill {\mathcal{I}}_n\left(\mathit{\theta }\right|\mathit{y})& =-{\ddot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})\hfill \end{array}$$

with ${\ddot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})={\displaystyle \frac{\partial {\dot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})}{\partial {\mathit{\theta }}^{\top }}}$. By differentiating (A12) with respect to ${\mathit{\theta }}^{\top }$, we obtain

$$\begin{array}{cc}\hfill {\mathcal{I}}_n\left(\mathit{\theta }\right|\mathit{y})& =\sum _{i=1}^n\left[\left({\displaystyle \frac{y_i}{{\mu }_i^2}}+{\displaystyle \frac{n_i-y_i}{{(1-{\mu }_i)}^2}}\right){\dot{\mathit{\mu }}}_i{\dot{\mathit{\mu }}}_i^{\top }-{\displaystyle \frac{y_i-n_i{\mu }_i}{{\mu }_i(1-{\mu }_i)}}{\ddot{\mathit{\mu }}}_i\right]\hfill \end{array}$$

(A16)

where ${\ddot{\mathit{\mu }}}_i={\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\top }}}$ is obtained from (A13) as

$$\begin{array}{cc}\hfill {\ddot{\mathit{\mu }}}_i& =\left[\begin{array}{cccccc}{\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_1\partial {\mathit{\beta }}_1^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_1\partial {\mathit{\beta }}_2^{\top }}}& \cdots & {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_1\partial {\mathit{\beta }}_q^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_1\partial {\mathit{\gamma }}^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_1\partial {\mathit{\delta }}^{\top }}}\\ \\ {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_2\partial {\mathit{\beta }}_1^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_2\partial {\mathit{\beta }}_2^{\top }}}& \cdots & {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_2\partial {\mathit{\beta }}_q^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_2\partial {\mathit{\gamma }}^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_2\partial {\mathit{\delta }}^{\top }}}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_q\partial {\mathit{\beta }}_1^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_q\partial {\mathit{\beta }}_2^{\top }}}& \cdots & {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_q\partial {\mathit{\beta }}_q^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_q\partial {\mathit{\gamma }}^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_q\partial {\mathit{\delta }}^{\top }}}\\ \\ {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\beta }}_1^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\beta }}_2^{\top }}}& \cdots & {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\beta }}_q^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\gamma }}^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\delta }}^{\top }}}\\ \\ {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\delta }\partial {\mathit{\beta }}_1^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\delta }\partial {\mathit{\beta }}_2^{\top }}}& \cdots & {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\delta }\partial {\mathit{\beta }}_q^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\delta }\partial {\mathit{\gamma }}^{\top }}}& {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\delta }\partial {\mathit{\delta }}^{\top }}}\end{array}\right]\hfill \end{array}$$

(A17)

where for $j=1,\cdots ,q$ and on setting $\ddot{h}\left(x\right)={\displaystyle \frac{{\partial }^{2}h\left(x\right)}{\partial x^2}}$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_j{\mathit{\beta }}_j^{\top }}}& ={\lambda }_i(1-{\alpha }_i)\left[\prod _{k\ne j}h\left({\displaystyle \frac{{\eta }_{ik}}{\sqrt{1+{({\overline{\sigma }}_{ik}/{\sigma }_o)}^2}}}\right)\right]{\displaystyle \frac{{\partial }^2{\mu }_{ij}}{\partial {\mathit{\beta }}_j\partial {\mathit{\beta }}_j^{\top }}},\hfill \end{array}$$

(A18a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\gamma }}^{\top }}}& ={\lambda }_i\left[1-\prod _{j=1}^{q}h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\right]\ddot{h}\left({\mathbf{U}}_i^{\top }\mathit{\gamma }\right){\mathbf{U}}_i{\mathbf{U}}_i^{\top },\hfill \end{array}$$

(A18b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\delta }\partial {\mathit{\delta }}^{\top }}}& =\left[{\alpha }_i+(1-{\alpha }_i)\prod _{j=1}^{q}h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\right]\ddot{h}\left({\mathbf{V}}_i^{\top }\mathit{\delta }\right){\mathbf{V}}_i{\mathbf{V}}_i^{\top },\hfill \end{array}$$

(A18c)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_j\partial {\mathit{\beta }}_k^{\top }}}& ={\lambda }_i(1-{\alpha }_i)\left[\prod _{l\ne j,k}h\left({\displaystyle \frac{{\eta }_{il}}{\sqrt{1+{({\overline{\sigma }}_{il}/{\sigma }_o)}^2}}}\right)\right]{\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_j}}{\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_k^{\top }}}(\mathrm{for}k\ne j),\hfill \end{array}$$

(A18d)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_j\partial {\mathit{\gamma }}^{\top }}}& =-{\lambda }_i\left[\prod _{k\ne j}h\left({\displaystyle \frac{{\eta }_{ik}}{\sqrt{1+{({\overline{\sigma }}_{ik}/{\sigma }_o)}^2}}}\right)\right]\dot{h}\left({\mathbf{U}}_i^{\top }\mathit{\gamma }\right){\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_j}}{\mathbf{U}}_i^{\top },\hfill \end{array}$$

(A18e)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial {\mathit{\beta }}_j\partial {\mathit{\delta }}^{\top }}}& =(1-{\alpha }_i)\left[\prod _{k\ne j}h\left({\displaystyle \frac{{\eta }_{ik}}{\sqrt{1+{({\overline{\sigma }}_{ik}/{\sigma }_o)}^2}}}\right)\right]\dot{h}\left({\mathbf{V}}_i^{\top }\mathit{\delta }\right){\displaystyle \frac{\partial {\mu }_{ij}}{\partial {\mathit{\beta }}_j}}{\mathbf{V}}_i^{\top },\hfill \end{array}$$

(A18f)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_i}{\partial \mathit{\gamma }\partial {\mathit{\delta }}^{\top }}}& =\left[1-\prod _{j=1}^{q}h\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right)\right]\dot{h}\left({\mathbf{U}}_i^{\top }\mathit{\gamma }\right)\dot{h}\left({\mathbf{V}}_i^{\top }\mathit{\delta }\right){\mathbf{U}}_i{\mathbf{V}}_i^{\top },\hfill \end{array}$$

(A18g)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{\mu }_{ij}}{\partial {\mathit{\beta }}_j\partial {\mathit{\beta }}_j^{\top }}}& =\ddot{h}\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right){\mathit{\tau }}_{ij}{\mathit{\tau }}_{ij}^{\top }+\dot{h}\left({\displaystyle \frac{{\eta }_{ij}}{\sqrt{1+{({\overline{\sigma }}_{ij}/{\sigma }_o)}^2}}}\right){\dot{\mathit{\tau }}}_{ij},\hfill \end{array}$$

(A18h)

$$\begin{array}{cc}\hfill \mathrm{with}{\dot{\mathit{\tau }}}_{ij}& ={\displaystyle \frac{{\sigma }_o}{{\left({\sigma }_o^2+{\overline{\sigma }}_{ij}^2\right)}^{3/2}}}\left[{\displaystyle \frac{3{\eta }_{ij}}{{\sigma }_o^2+{\overline{\sigma }}_{ij}^2}}{\mathit{S}}_{ij}{\mathit{\beta }}_j{\mathit{\beta }}_j^{\top }{\mathit{S}}_{ij}-{\eta }_{ij}{\mathit{S}}_{ij}-{\mathit{x}}_{ij}{\mathit{\beta }}_j^{\top }{\mathit{S}}_{ij}-{\mathit{S}}_{ij}{\mathit{\beta }}_j{\mathit{x}}_{ij}^{\top }\right].\hfill \end{array}$$

(A18i)

An alternative to ${\mathcal{I}}_n\left(\mathit{\theta }\right|\mathit{y})$ used in ${\Sigma }_n={\left[{\mathcal{I}}_n\left(\hat{\mathit{\theta }}\right|\mathit{y})\right]}^{-1}$ is the Fisher information matrix ${\mathcal{E}}_n\left(\mathit{\theta }\right)$. Using ${\mathcal{E}}_n$, we have the potential asymptotic covariance matrix ${\Gamma }_n={\left[{\mathcal{E}}_n\left(\hat{\mathit{\theta }}\right)\right]}^{-1}$ as an alternative to ${\Sigma }_n$. For the purpose of constructing asymptotic confidence intervals for estimates, recent studies suggest that the Fisher information ${\mathcal{E}}_n\left(\hat{\mathit{\theta }}\right)$ performs at least as well as the observed information matrix ${\mathcal{I}}_n\left(\hat{\mathit{\theta }}\right|\mathit{y})$ [106,107,108]. Interestingly, ${\mathcal{E}}_n$ is also less expensive to compute as compared to ${\mathcal{I}}_n$. Indeed, Equation (A17) shows that even after computing the score (A12), a lot of additional computations are required to obtain ${\mathcal{I}}_n\left(\mathit{\theta }\right|\mathit{y})$. However, ${\mathcal{I}}_n\left(\hat{\mathit{\theta }}\right|\mathit{y})$ is usually available as a by-product from optimization routines in popular statistical softwares (see e.g. optim in R [84]).

Appendix B.2. Penalized Maximum Likelihood

For PML estimation based on ${\tilde{\ell }}_n$ (15), the penalized score ${\dot{\tilde{\ell }}}_n\left(\mathit{\theta }\right|\mathit{y})={\displaystyle \frac{\partial {\tilde{\ell }}_n\left(\mathit{\theta }\right)}{\partial \mathit{\theta }}}$ is given by

$$\begin{array}{c}\hfill {\dot{\tilde{\ell }}}_n\left(\mathit{\theta }\right|\mathit{y})={\dot{\ell }}_n\left(\mathit{\theta }\right|\mathit{y})+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial log|{\mathcal{E}}_n\left(\mathit{\theta }\right)|}{\partial \mathit{\theta }}}\end{array}$$

(A19)

where ${\dot{\ell }}_n$ is the score given in (A12), and by the chain rule, the double score penalty satisfies

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial log|{\mathcal{E}}_n\left(\mathit{\theta }\right)|}{\partial {\mathit{\theta }}_j}}& =\mathrm{Tr}\left[{\left[{\mathcal{E}}_n\left(\mathit{\theta }\right)\right]}^{-1}{\displaystyle \frac{\partial {\mathcal{E}}_n\left(\mathit{\theta }\right)}{\partial {\mathit{\theta }}_j}}\right].\hfill \end{array}$$

(A20)

The required derivatives of ${\mathcal{E}}_n\left(\mathit{\theta }\right)$ are obtained from ${\displaystyle \frac{\partial {\mathcal{E}}_n\left(\mathit{\theta }\right)}{\partial {\mathit{\theta }}_j}}={\sum }_{i=1}^n{\displaystyle \frac{\partial }{\partial {\mathit{\theta }}_j}}\left({\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}{\dot{\mathit{\mu }}}_i{\dot{\mathit{\mu }}}_i^{\top }\right)$ as

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{E}}_n\left(\mathit{\theta }\right)}{\partial {\mathit{\theta }}_j}}& =\sum _{i=1}^n\left({\displaystyle \frac{n_i(2{\mu }_i-1)}{{\mu }_i^2{(1-{\mu }_i)}^2}}{\dot{\mathit{\mu }}}_{i,j}\right){\dot{\mathit{\mu }}}_i{\dot{\mathit{\mu }}}_i^{\top }+{\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}\mathit{M}(i,j),\hfill \end{array}$$

where ${\dot{\mathit{\mu }}}_{i,j}$ is the jth element of ${\dot{\mathit{\mu }}}_i$ given by Equation (A13), and the matrix $\mathit{M}(i,j)={\displaystyle \frac{\partial }{\partial {\mathit{\theta }}_j}}\left({\dot{\mathit{\mu }}}_i{\dot{\mathit{\mu }}}_i^{\top }\right)$ has u,vth element $\mathit{M}{(i,j)}_{uv}={\displaystyle \frac{\partial }{\partial {\mathit{\theta }}_j}}\left({\dot{\mathit{\mu }}}_{i,u}{\dot{\mathit{\mu }}}_{i,v}\right)$, i.e.,

$$\begin{array}{cc}\hfill \mathit{M}(i,j)& ={\ddot{\mathit{\mu }}}_{i,·j}{\dot{\mathit{\mu }}}_i^{\top }+{\dot{\mathit{\mu }}}_i{\ddot{\mathit{\mu }}}_{i,·j}^{\top },\hfill \end{array}$$

with ${\ddot{\mathit{\mu }}}_{i,·j}$ the jth column of the hessian matrix ${\ddot{\mathit{\mu }}}_i$ given in Equation (A17). Note that Equation (A20) can be rewritten as

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial log|{\mathcal{E}}_n\left(\mathit{\theta }\right)|}{\partial {\mathit{\theta }}_j}}& ={\dot{\mathcal{E}}}_{n,j}^{\top }\left(\mathit{\theta }\right)\mathrm{vec}\left({\left[{\mathcal{E}}_n\left(\mathit{\theta }\right)\right]}^{-1}\right)\hfill \end{array}$$

where ${\dot{\mathcal{E}}}_{n,j}\left(\mathit{\theta }\right)=\mathrm{vec}\left({\displaystyle \frac{\partial {\mathcal{E}}_n\left(\mathit{\theta }\right)}{\partial {\mathit{\theta }}_j}}\right)={\displaystyle \frac{\partial \mathrm{vec}\left({\mathcal{E}}_n\left(\mathit{\theta }\right)\right)}{\partial {\mathit{\theta }}_j}}$ and $\mathrm{vec}$ is the usual operator which stacks the columns of its matrix argument. Denoting ⨂ the Kronecker (direct) product, and applying Equation (2) in [109], we get

$$\begin{array}{c}\hfill {\dot{\mathcal{E}}}_{n,j}\left(\mathit{\theta }\right)=\sum _{i=1}^n\left({\displaystyle \frac{n_i(2{\mu }_i-1)}{{\mu }_i^2{(1-{\mu }_i)}^2}}{\dot{\mathit{\mu }}}_{i,j}\right){\dot{\mathit{\mu }}}_i⨂{\dot{\mathit{\mu }}}_i+{\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}\left[{\dot{\mathit{\mu }}}_i⨂{\ddot{\mathit{\mu }}}_{i,·j}+{\ddot{\mathit{\mu }}}_{i,·j}⨂{\dot{\mathit{\mu }}}_i\right].\end{array}$$

Further using Equation (1) in [109] yields

$$\begin{array}{c}\hfill {\dot{\mathcal{E}}}_{n,j}\left(\mathit{\theta }\right)=\sum _{i=1}^n\left({\displaystyle \frac{n_i(2{\mu }_i-1)}{{\mu }_i^2{(1-{\mu }_i)}^2}}{\dot{\mathit{\mu }}}_{i,j}\right){\dot{\mathit{\mu }}}_i⨂{\dot{\mathit{\mu }}}_i+{\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}\left({\mathit{K}}_{m,m}+{\mathit{I}}_{m^2}\right)\left({\ddot{\mathit{\mu }}}_{i,·j}⨂{\dot{\mathit{\mu }}}_i\right)\end{array}$$

where ${\mathit{K}}_{m,m}$ is the commutation matrix [110] and m denotes the length of the parameter vector $\mathit{\theta }$. Collecting the column vectors ${\dot{\mathcal{E}}}_{n,j}$, we obtain

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial log|{\mathcal{E}}_n\left(\mathit{\theta }\right)|}{\partial \mathit{\theta }}}={\dot{\mathcal{E}}}_n^{\top }\left(\mathit{\theta }\right)\mathrm{vec}\left({\left[{\mathcal{E}}_n\left(\mathit{\theta }\right)\right]}^{-1}\right)\hfill \end{array}$$

(A21a)

where ${\dot{\mathcal{E}}}_n$ is a $m^2\times m$ jerk matrix defined as ${\dot{\mathcal{E}}}_n\left(\mathit{\theta }\right)={\displaystyle \frac{\partial }{\partial {\mathit{\theta }}^{\top }}}\mathrm{vec}\left({\mathcal{E}}_n\left(\mathit{\theta }\right)\right)$ and given by

$$\begin{array}{cc}\hfill {\dot{\mathcal{E}}}_n\left(\mathit{\theta }\right)& =\sum _{i=1}^n{\displaystyle \frac{n_i}{{\mu }_i(1-{\mu }_i)}}\left[{\displaystyle \frac{2{\mu }_i-1}{{\mu }_i(1-{\mu }_i)}}\left({\dot{\mathit{\mu }}}_i⨂{\dot{\mathit{\mu }}}_i\right){\dot{\mathit{\mu }}}_i^{\top }+\left({\mathit{K}}_{m,m}+{\mathit{I}}_{m^2}\right)\left({\ddot{\mathit{\mu }}}_i⨂{\dot{\mathit{\mu }}}_i\right)\right].\hfill \end{array}$$

(A21b)

Solving the penalized score equation ${\dot{\tilde{\ell }}}_n\left(\mathit{\theta }\right|\mathit{y})=\mathit{0}$ (with ${\dot{\tilde{\ell }}}_n$ given by Equation (A19)) produces the PML estimates $\tilde{\mathit{\theta }}={\left({\tilde{\mathit{\beta }}}^{\top },{\tilde{\mathit{\gamma }}}^{\top },{\tilde{\mathit{\delta }}}^{\top }\right)}^{\top }$ of $\mathit{\theta }$. For large sample inference on $\tilde{\mathit{\theta }}$, the Fisher information based asymptotic covariance matrix is given by ${\Gamma }_n={\left[{\mathcal{E}}_n\left(\tilde{\mathit{\theta }}\right)\right]}^{-1}$. Indeed, the penalized Fisher information matrix (covariance matrix of ${\dot{\tilde{\ell }}}_n$) is, by Equation (A11) (the same as for ${\dot{\ell }}_n$), given by Equation (A15) because ${\dot{\tilde{\ell }}}_n$ is equal to ${\dot{\ell }}_n$ up to an additive constant (independent of the response Y) which does not affect covariance matrix calculation.

Appendix C.1. Simulation Design

We consider the following extension of the MSB data generation process considered in Section 3:

$$\begin{array}{ccc}\hfill \hspace{1em}\hspace{1em} \hspace{1em}& \hfill Y_i& \stackrel{ind}{\sim }\mathcal{B}in\left(1,{\mu }_i\right),\hfill \end{array}$$

(A22a)

$$\begin{array}{cc}\hfill \mathrm{where}{\mu }_i& =\lambda \prod _{j=1}^{q}h\left({\eta }_{ij}\right),\hfill \end{array}$$

(A22b)

$$\begin{array}{ccc}\hfill \hspace{1em} \hspace{1em}\hspace{1em}& \hfill {\eta }_{ij}& ={\beta }_0+{\beta }_1{X}_{ij}.\hfill \end{array}$$

(A22c)

In (A22), $Y_i$ is the binary response for unit i, $\lambda \in (0,1]$ is an upper limit for the success probability of the binary response Y, h is the logit link function defined in (3), X denotes the explanatory variable following a predictor distribution denoted $PD$, and ${\beta }_0$ and ${\beta }_1$ are intercept and slopes.

Note that if $\lambda =1$ with $q=1$, then (A22) is a standard simple logistic data generation process. For $\lambda <1$ and $q=1$, (A22) represents a multistage process where only one stage is observed. We varied the upper limit $\lambda$ of the success probability of the response Y from 0.5 to 1, varying regression coefficients as ${\beta }_0=3,8$ and ${\beta }_1=-3,-2,-1,0$, and the number n of units from 100 to 2000.

Under model (A22), the true success probability curves are shown in Figure A1 for ${\beta }_0=3$, ${\beta }_1=-1$ and selected $\lambda$ values from 0.1 to 1. For a predictor X, we first considered a uniform distribution. For instance, with ${\beta }_0=3$ and ${\beta }_1=-1$, we simulated X values uniformly from the interval $[-2,8]$ (Figure A2), having expectation ${\mu }_x=3$ and standard deviation ${\sigma }_x=2.89$, and spanning 99% probability mass of the standard logistic distribution centered at 3.

It is well established that the distribution of the explanatory variable has a large impact on inference in logistic regression [111]. To investigate this effect in our MSB framework, we considered two alternative continuous probability distributions for X, including the normal (symmetric) and the lognormal (asymmetric) distributions. Normal covariate values were simulated with mean ${\mu }_w=3$ and standard deviation ${\sigma }_x=1.6$. We used a shifted lognormal distribution, obtained by subtracting two from log normal variates with mean ${\mu }_x=2.65$ and standard deviation ${\sigma }_x=1.41$ (log scale mean and standard are 1.1 and 0.5, respectively). For both normal and lognormal covariates, the chosen distribution parameters ensure that X values cover the interval $[-2,8]$ with probability 99% when ${\beta }_0=3$ and ${\beta }_1=-1$ (Figure A2). The lognormal distribution giving more weight to X values along the asymptote of the success probability curve ($X\in [-0.89,2.91]$; Figure A2) is expected to result in more precise (lower variance) estimates of the asymptote $\lambda$ (irrespective of the true $\lambda$ value), but less accurate (more biased/less precise) estimates of ${\beta }_0$ and ${\beta }_1$ when the true $\lambda$ value is low.

A summary of our simulation factors, including $\lambda$, ${\beta }_0$, ${\beta }_1$, $PD$, and sample size (n) is shown in Table A1. As performance measures, we computed the relative bias and root mean square error ($rmse$) of the estimates of model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$, and the $mad\left(\hat{\mu }\right)$ and $R^2\left(\hat{\mu }\right)$ of predicted success probabilities as defined in Section 3.1, based on 2000 independent replicates of each simulation setting.

Table A1. Levels of simulation factors

Simulation factor	Notation	Levels
Sample size	$n$	100, 200, 300, 500, 750, 1000, 1500, 2000
Upper limit of success probability	$\lambda$	0.5, 0.75, 0.9, 1
Intercept	${\beta }_0$	3, 8
Slope	${\beta }_1$	-3, -2, -1
Distribution of the predictor	$PD$	Uniform, Normal, Log-normal

Table A1. Levels of simulation factors

Simulation factor	Notation	Levels
Sample size	$n$	100, 200, 300, 500, 750, 1000, 1500, 2000
Upper limit of success probability	$\lambda$	0.5, 0.75, 0.9, 1
Intercept	${\beta }_0$	3, 8
Slope	${\beta }_1$	-3, -2, -1
Distribution of the predictor	$PD$	Uniform, Normal, Log-normal

Figure A1. Success probability from model (A22) as a function the explanatory variable X for selected asymptote values $\lambda$ and ${\beta }_0=3$ and ${\beta }_1=-1$.

Figure A1. Success probability from model (A22) as a function the explanatory variable X for selected asymptote values $\lambda$ and ${\beta }_0=3$ and ${\beta }_1=-1$.

Figure A2. Probability distributions considered for the covariate X in model (A22) with ${\beta }_0=3$ and ${\beta }_1=-1$. All density curves cover the region $w\in [-2,8]$ with probability 99% or more (100% for the uniform distribution). The normal density is above the uniform density curve over the region $X\in [0.84,5.16]$ where the normal has 90% probability against 43% for the uniform. The lognormal density is above the uniform density curve over the region $X\in [-0.89,2.91]$ where the lognormal has 81% probability against 38% for the uniform. For general ${\beta }_0$ and ${\beta }_1$ values, X is rescaled (linear transformation) so as to preserve the logit scale range of the linear predictor $\eta ={\beta }_0+{\beta }_{1}X$.

Figure A2. Probability distributions considered for the covariate X in model (A22) with ${\beta }_0=3$ and ${\beta }_1=-1$. All density curves cover the region $w\in [-2,8]$ with probability 99% or more (100% for the uniform distribution). The normal density is above the uniform density curve over the region $X\in [0.84,5.16]$ where the normal has 90% probability against 43% for the uniform. The lognormal density is above the uniform density curve over the region $X\in [-0.89,2.91]$ where the lognormal has 81% probability against 38% for the uniform. For general ${\beta }_0$ and ${\beta }_1$ values, X is rescaled (linear transformation) so as to preserve the logit scale range of the linear predictor $\eta ={\beta }_0+{\beta }_{1}X$.

Appendix C.2. Additional Simulation Results

Appendix C.2.1. Bias

Figure A3–Figure A14 present the relative bias in the estimates $\hat{\lambda }$, ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$ for model parameters $\lambda \in [0.25,1]$, ${\beta }_0\in \left\{3,8\right\}$ and ${\beta }_1\in \left\{-3,-2,-1\right\}$ and sample size $n\in [100,2000]$. The relative bias generally approaches zero as n increases. For $\hat{\lambda }$, the bias is generally negative for high true $\lambda$ values (close to one), and positive for low $\lambda$ values ($0.25\le \lambda <0.5$).

Figure A3. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A3. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A4. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A4. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A5. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A5. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A6. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A6. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A7. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A7. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A8. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A8. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A9. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A9. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A10. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A10. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A11. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A11. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A12. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A12. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A13. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A13. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A14. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A14. Relative bias of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Appendix C.2.2. Accuracy

Figure A15–Figure A26 present the relative root mean square error ($rrmse$) in the estimates $\hat{\lambda }$, ${\hat{\beta }}_0$ and ${\hat{\beta }}_1$ for model parameters $\lambda \in [0.25,1]$, ${\beta }_0\in \left\{3,8\right\}$ and ${\beta }_1\in \left\{-3,-2,-1\right\}$ and sample size $n\in [100,2000]$. Figure A27 and Figure A28 shows the $rrmse$ in ${\hat{\beta }}_1$ grouped for all parameter settings. The $rrmse$ generally approaches zero as n increases.

Figure A15. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A15. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A16. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A16. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A17. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A17. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A18. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A18. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A19. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A19. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A20. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A20. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A21. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A21. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \ge 0.8$.

Figure A22. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A22. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-1$, and $\lambda \le 0.75$.

Figure A23. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A23. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \ge 0.8$.

Figure A24. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A24. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-2$, and $\lambda \le 0.75$.

Figure A25. Relative root mean square error ($rmse$A22 ) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A25. Relative root mean square error ($rmse$A22 ) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \ge 0.8$.

Figure A26. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A26. Relative root mean square error ($rmse$) of Penalized Maximum Likelihood estimates of multistage binomial model parameters $\lambda$, ${\beta }_0$ and ${\beta }_1$ under the MSB model data generation process in Equation (A22) with ${\beta }_0=8$, ${\beta }_1=-3$, and $\lambda \le 0.75$.

Figure A27. Relative root mean square error ($rrmse$) of Penalized Maximum Likelihood estimate of the slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with $q=1$ (one predictor) for $\lambda \ge 0.8$. The figure columns show the investigated pairs of true ${\beta }_0$ and ${\beta }_1$ values, and rows correspond to investigated distributions of predictor.

Figure A27. Relative root mean square error ($rrmse$) of Penalized Maximum Likelihood estimate of the slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with $q=1$ (one predictor) for $\lambda \ge 0.8$. The figure columns show the investigated pairs of true ${\beta }_0$ and ${\beta }_1$ values, and rows correspond to investigated distributions of predictor.

Figure A28. Relative root mean square error ($rrmse$) of Penalized Maximum Likelihood estimate of the slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with $q=1$ (one predictor) for $\lambda \le 0.75$. The figure columns show the investigated pairs of true ${\beta }_0$ and ${\beta }_1$ values, and rows correspond to investigated distributions of predictor.

Figure A28. Relative root mean square error ($rrmse$) of Penalized Maximum Likelihood estimate of the slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with $q=1$ (one predictor) for $\lambda \le 0.75$. The figure columns show the investigated pairs of true ${\beta }_0$ and ${\beta }_1$ values, and rows correspond to investigated distributions of predictor.

Appendix C.2.3. Sample Size Determination

In complement to Figure 4, Figure A29 and Figure A30 present the sample size n required to achieve a specified relative root mean square error for the PML estimate ${\hat{\beta }}_1$ of the regression slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with ${\beta }_0=3$, ${\beta }_1=-2$ for normal and log-normal predictor distributions, respectively.

Figure A29. Sample size n required to achieve a specified relative root mean square error for the Penalized Maximum Likelihood estimate ${\hat{\beta }}_1$ of the regression slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with a normal predictor distribution, ${\beta }_0=3$, ${\beta }_1=-2$ and selected $\lambda$ values.

Figure A29. Sample size n required to achieve a specified relative root mean square error for the Penalized Maximum Likelihood estimate ${\hat{\beta }}_1$ of the regression slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with a normal predictor distribution, ${\beta }_0=3$, ${\beta }_1=-2$ and selected $\lambda$ values.

Figure A30. Sample size n required to achieve a specified relative root mean square error for the Penalized Maximum Likelihood estimate ${\hat{\beta }}_1$ of the regression slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with a log-normal predictor distribution, ${\beta }_0=3$, ${\beta }_1=-2$ and selected $\lambda$ values.

Figure A30. Sample size n required to achieve a specified relative root mean square error for the Penalized Maximum Likelihood estimate ${\hat{\beta }}_1$ of the regression slope ${\beta }_1$ under the MSB data generation process in Equation (A22) with a log-normal predictor distribution, ${\beta }_0=3$, ${\beta }_1=-2$ and selected $\lambda$ values.