{'fraction_non_alphanumeric': 0.04243800209570381, 'fraction_numerical': 0.015468290005488747, 'mean_word_length': 4.487063655030801, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mediation analysis with case-control sampling: Identification and estimation in the presence of a binary mediator

With reference to a stratified case-control procedure based on a binary variable of primary interest, we derive the expression of the distortion induced by the sampling design on the parameters of the logistic model of a secondary variable. This is particularly relevant when performing mediation analysis (possibly in a causal framework) with stratified case-control data in settings where both the outcome and the mediator are binary. Our identification result opens the way to M-estimation and Maximum Likelihood estimation. We then conduct a simulation study showing the gain in efficiency of the estimators of both the outcome and mediator model parameters w.r. to existing methods, based on weighting. As an illustrative example, we reanalyze a German case-control dataset in order to investigate whether the effect of reduced immunocompetency on listeriosis onset is mediated by the intake of gastric acid suppressors.


Introduction
Retrospective sampling schemes are common practice in studies when either the outcome of interest is rare or some covariates are difficult to measure. In this context, it is too expensive or too lengthy to random sample the whole population. Instead, more efficiently, random samples are extracted from a partition of the population according to the outcome and possibly other stratifying factors. When the outcome is binary, this procedure is known as (stratified) case-control sampling design.
Though the implications of the retrospective sampling are widely understood when the only interest is the relationship between the covariates and the outcome [5], less is known when the relationship between some of the covariates is also the object of inference. Typically, this situation arises in mediation analysis, where the aim is to decompose, in a counterfactual framework, the total effect of a treatment on the outcome into a direct and an indirect effect, the second one due to a mediator. A mediator is a variable that is a response of the treatment and in turn influences the outcome [28,39]. In this situation, two equations are of interest, the first linking the response to the treatment and the mediator, the second linking the mediator to the treatment. If the probability of a unit to be in the sample depends on the outcome (and possibly some other stratifying factors), the relationship between the treatment and the mediator is bound to be distorted and standard mediation methods fail.
To tackle this problem, VanderWeele and Vanstelaandt (2010) [43] propose to estimate the mediator-treatment equation on the sub-population of controls. Under the rare outcome assumption, this can be considered a random sample from the whole population; see also VanderWeele and Tchetgen Tchetgen (2016) [41]. Another strategy, also mentioned in VanderWeele and Vanstelaandt (2010) [43], is based on weighting, a proposal originally contained in some papers by Mansky and coauthors [21,22] in the context of choice-based samples, which are the analogous of stratified case-control samples within the econometric literature; see also VanderLaan (2008) [38].
With reference to a binary mediator and a binary outcome, we here propose a parametric approach to mediation analysis when the sampling scheme depends on the outcome, and possibly other stratifying covariates. We build on results already existing in the literature on parametric mediation analysis for a binary mediator and a binary outcome [35,9] to derive the explicit expression of the distortion induced by the sampling design. We exploit the knowledge provided by Directed Acyclic Graph (DAGs) [19], where an additional binary node is representing the sampling scheme. In this regard, we are in line with   [7]. We then propose two estimating procedures, one based on Maximum Likelihood (ML) and the other one on M-estimation, and compare them with the weighting method.
There are a number of additional situations where this investigation can be of interest. The first one concerns informative missingness in logistic regression. It turns out that our work extends the derivations in Wang et al. (2107) [47] to the continuous treatment context, also proposing two likelihood-based estimating procedures that can be easily implemented with standard statistical software. Another situation of interest arises in case control studies with secondary or additional outcomes [33], which is particularly common in genetic epidemiology, where information on a secondary phenotype is also of interest [46]. More generally, properly recovering the distortion induced on associations by conditioning on a collider node is fundamental for many other purposes. In this sense, this work extends to a parametric framework the results on odds-ratios of Nguyen et al. (2019) [26].
The reminder of the paper is organized as follows. In Section 2, we set notation and introduce the target regression models, showing that adjustments are needed to identify them in the presence of outcome-dependent sampling schemes. Our identification strategy is presented in Section 3, while the specifics of the two estimation methods deriving from it (M-estimation and ML) are described in Section 4. In Section 5, we apply the proposed methodology to data coming from a German case-control study on listeriosis. Our aim is to decompose, in a counterfactual framework, the effect of reduced immunocompetency on listeriosis into a direct and an indirect effect, with the latter mediated by gastric acid suppressor intake. Section 6 reports evidence from a simulation study conducted to compare M-and ML estimators to weighting estimators, while some concluding remarks are given in Section 7.

Background
Let A denote a generic treatment/exposure, while M and Y represent the binary mediator and the binary outcome, respectively. Like in other related approaches (see e.g. Didelez et al., 2010 [7]), outcome-dependent sampling is theoretically represented by a selection indicator variable, W . Data are available just on the selected units, i.e. those W = 1. A straightforward extension, that is not here considered, is when data can be seen as a random sample extracted from the population with W = 1. The binary variable W is influenced by Y via a known probabilistic mechanism, that can be either marginal or conditional on the strata formed by a set of categorical background covariates, B, that possibly affect A, M and/or Y as well. When B is empty, this corresponds to the unconditional Case-Control (CC) design. Otherwise, the Stratified Case-Control (SCC) design arises. We also indicate by S and Z the sets of covariates (of any nature) that influence M and/or Y , respectively, but not the selection mechanism W . Notice that S and Z do not necessarily have to be disjoint sets, i.e., variables in S might also be in Z and vice versa. The overall framework is represented in the two DAGs in Figure 1, where the W node is squared in order to pinpoint its nature of selection node.
In what follows, we account for the presence of B and thus refer to the SCC setting, with the CC one being a special case thereof. In order to ease the exposition, without loss of generality we think of B, S and Z as of singleton variables, with B taking values in {1, . . . , N b } to represent the N b different strata formed by background variables. Also, if ambiguities do not occur we occasionally use P (· | k) as a shorthand for P (· | K = k) . The logistic regression models assumed to hold in the population are and where β and δ are coefficient vectors, while x y and x m contain the values of the explanatory variables (including the constant term) as well as stratumspecific indicator variables I{B = b} for b = 2, . . . , N b (the usual corner-point parametrization taking the first stratum as reference is adopted). Clearly, such a general formulation allows for the presence, in both models, of firstand higher-order (e.g., quadratic) effects as well as of interaction terms. The impact of SCC sampling can be expressed in the language of DAGs and conditional independence [6] by stating that Y and M are not conditionally independent of W given the respective parent node sets, i.e., As a consequence, which shows that the models postulated in the right-hand sides of (1) and (2) are not identified in the sub-population of selected units, and that adjusting approaches are needed. However, by inspection of the DAG it is possible to see that and therefore Y can be seen as a selection bias breaking variable, a term introduced elsewhere in the literature [12]. This property, together with the explicit expressions of the distortion induced by the sampling design on the parameters of (1) and (2), forms the basis of our proposal. Notice that the minimum conditioning set for independence of M and W to hold is (B, Y ). However, since we are interested in the parametric formulation of the model for M, all variables in the conditioning set of (3) are relevant.

Identification
In order to identify the coefficients of model (1) from SCC data, it is wellknown that simple coefficient adjustments involving the conditional population prevalences π b = P (Y = 1 | B = b) suffice [30,10,4]. In detail, the model expression holding in the sub-population of selected units is given by where β ⋆ is equal to β except for modifications concerning the intercept (β INT ) and the coefficient sub be the proportion of cases in the SCC sample for stratum b, the replacing elements in β ⋆ are where Each of the correction terms above is known whenever π b is. Indeed, p b /(1 − p b ) is the stratum-specific CC ratio, which is fixed by design. We here assume π b (b = 1, . . . , N b ) to be known from external auxiliary information sources.
We turn now to the model for P (M = 1 | a, b, s, z, y, W = 1). From (3), we know that P (M = 1 | a, b, s, z, y, W = 1) = P (M = 1 | a, b, s, z, y). After some derivations, see Appendix A, it is possible to see that where is a correction term depending on y and the coefficient vector β. The parametric expression of it is (13) of Appendix A.

Estimation
Starting from the identification result in (4), inference on β ⋆ can be conducted from SCC data in a standard ML framework [30,2]. In parallel to this coefficient adjusting approach, weighting methods have also been proposed within the econometric literature dealing with choice-based samples [21,22], which can essentially be thought of as the analogous of SCC samples [5]. These methods consist in fitting the population model (1) to the SCC sample, assigning cases and controls in each stratum b a weight equal to π b /p b and (1 − π b )/(1 − p b ), respectively. The resulting estimators are consistent but, although possibly less prone to model misspecification [49], they were found to be typically less efficient than the ones obtained with parameter corrections [18]. Importantly, the weighting approach was also extended to variables different from the original outcome (sometimes termed secondary or additional outcomes [33]); see for example the more general work by VanderLaan (2008) [38]. Such an extension is of particular relevance in a mediation framework, where the population parameters of the mediator model also need to be recovered from SCC data in order to achieve effect decomposition [43]. In our setting, this would correspond to fitting the population model (2) to the SCC sample with the same weighting scheme as above.
Here, we introduce an alternative framework where estimation of the δ vector from SCC data exploits the identification result in (6). Since β is also involved in that equation, our framework is designed for estimating the combined parameter vector θ = (β ⊤ , δ ⊤ ) ⊤ altogether. This also fits with the mediation setting, where typically both β and δ have to be estimated at the same time. Two estimation strategies are discussed: in Section 4.1 we present an M-estimation approach, while in Section 4.2 the above mentioned ML framework is extended to our bivariate context.

M-estimation
In a formal M-estimation setting [16], the M- where d θ is the dimension of θ and s y (·) and s m (·) are the observed score functions associated to the models in (4) and (6), respectively; see Appendix B for their expressions. In practice, a two-step procedure is required. First, theβ M estimate is obtained by fitting model (4) and implementing the corrections for β INT and β b . Then,δ M is computed by fitting model (6), where o(y, x y ;β M ) is included as an offset term. Notice thatδ M is called a partial M-estimator [36], since it requires theβ M estimate of the first step to be plugged-in in place of the unknown β vector [32]. In this case, knowledge of the adjusting terms log(k 1 ) and log(k −1 ) is essential, since the estimated offset o(y, The variance-covariance matrix of theθ M estimator is given by and Ψ i (θ) is the score random vector of a generic unit i. The finite-sample estimate of V (θ M ) can be computed aŝ where A(y, θ) and B(y, θ) are the sample equivalent of A(Y i , θ) and B(Y i , θ) [36]. Their expressions are also reported in Appendix B.

Maximum Likelihood estimation
Since we are dealing with a SCC sample, the likelihood of the observed data corresponds to the conditional density of the (M, A, S, Z) random vector given (Y, B) and W = 1. As not only the β parameters of (1) are the object of inference, but also the δ parameters of (2), ML theory developed by Prentice and Pyke (1979) [30] should be extended. With reference to discrete choice models, Imbens (1992) [17] provides more general results that can be applied to this context. Specifically, given n independent sample units indexed by i = 1, . . . , n, it follows from Imbens (1992) that ML estimation of θ can be performed by maximizing An additional factorization of (9) leads to where , respectively. While the latter is related to the parameter vector θ via (6), the former involves marginalization over the binary mediator M. The corresponding model on the logistic scale can be written as where the first term is (4) evaluated at M = 0 while is a correction term which depends on x y , x m as well as on both parameter vectors. See Appendix C for the parametric expression of (11). Within this framework, the ML estimate of θ,θ ML , can be obtained by maximizing the log-likelihood in (10) via the usual iterative methods. Since these methods are likely to suffer from local maxima problems, it is advisable to set several starting points fluctuating around a sensible deterministic choice; see Section 5.2 for an account of the approach undertaken within our specific application. To further enhance the performance of optimization algorithms, it is typically useful to provide the expression of the log-likelihood gradient s(θ) = ∂ℓ(θ)/∂θ, which is reported in the first part of Appendix D.
In line with standard theory, the estimated variance-covariance matrix ofθ ML is given by Alternatively, sandwich estimation [34,48] can be performed via and s i (θ) is the individual contribution of the i-th unit to s(θ).

The Listeriosis dataset
To illustrate the proposed approach, we reconsider the dataset analyzed by Preußel et al. [31], where risk factors related to listeriosis in Germany are investigated. Listeriosis is an infection primarily contracted through the intake of food contaminated with Listeria monocytogenes bacteria. It mainly affects older adults, individuals with weakened immune systems and pregnant women [13]. Consequences might be quite severe, ranging from the development of life-threatening conditions for the fetus -in the case of pregnancyto severe illness or death in the other cases.
The study by Preußel et al. [31] focusses on sporadic non-pregnancy associated listeriosis. Many food-related risk factors, like consumption of cold cooked sausages and pre-sliced cheese, are considered in combination with personal characteristics as well as other risk factors such as Reduced Immunocompetency (RI) or the intake of Gastric Acid Suppressors (GASs). The latter act as effect modifiers and can also be thought of as risk factors themselves when focus goes beyond food-related factors [13,15,3,25]. In this application, we build upon this framework and extend it to a mediation setting in order to investigate whether -and to what extent -the additional vulnerability to listeriosis due to RI goes through the intake of GASs, which are commonly prescribed drugs for patients with RI [31,1]. In line with the original approach of Preußel et al. [31], the exposure (RI) is coded as a categorical variable with three mutually exclusive levels: 0=none, 1=RI due to immunocompromising diseases without immunosuppressive therapies (e.g., diabetes, autoimmune disorder), 2=RI due to immunosuppressive therapies (e.g., chemotherapy, radiation therapy). The mediator (GAS intake) is binary (0=no, 1=yes). RI and GAS intake are considered present if occurring within three months from listeriosis onset (cases) or interview (controls).
The dataset in Preußel et al. [31] contains 109 sporadic (i.e., not due to an outbreak) cases ascertained among over 40 men and non-pregnant women from all German Federal states (except Bremen, which accounts for around 0.8% of the German population) in the time span from March 2012 until December 2013. In addition, 1982 controls were enrolled from the population of subjects with no history of listeriosis which are accessible by telephone in Germany. In the original study, controls were sampled to get an equal number of individuals in the three age classes 40-65 years, 66-75 years and 76+ years, according to the age distribution of cases in the 2004-2011 years. However, the age distribution of the 109 cases in the dataset differs from this one. While this fact might be due to a number of reasons, we argue that such a discrepancy does not allow one to treat the original dataset as a proper SCC one. Indeed, in an SCC design involving age as a stratifying factor, a fixed number of controls within the same age class for each case are sampled, resulting in exactly the same age distribution for cases and controls. Conversely, in an unconditional CC design the age distribution of the controls should be close to the age distribution of the population (differing by sampling error only). Since the dataset has neither of these features, we obviate by randomly selecting a fixed number of controls (8 times the number of cases) in a way such that the age distribution of the German population in 2012 (according to the official webpage of the Federal Statistical Office of Germany) is matched, thereby generating an unconditional CC sample. The prevalence of sporadic non-pregnancy associated listeriosis is calculated assuming that 7 days is the average length of infection. The obtained rate is 1.70 × 10 −7 , which results in an adjusting correction term log(k 1 ) = 13.511.
In accordance with the data generating process postulated in Preußel et al. [31], we initially include age, educational level and gender as covariates both for the outcome and mediation model (Figure 2a). However, subsequent analyses lead to keeping only age (coded with the three classes introduced above) in the outcome model, and no covariates in the mediator model; see Figure 2b. With this regard, it is worth to underline that model selection is conducted separately for the two models, within an M-estimation framework. In detail, selection procedures for the mediator model are implemented  Table 1 reports results for two pairs of logistic models fitted with M-estimation, ML and weighting. Standard errors are obtained with a sandwich approach not only for M-and ML estimation (see Section 4), but also for weighting. ML estimation is performed through direct maximization of the log-likelihood function by means of the maxLik function in R [14]. Multiple starting values (10) are obtained by perturbing theθ M vector obtained from the M-estimates with a random normal deviation with null mean and standard deviation equal to 0.5. All attempts converge to the same solution. However, notice that unstable solutions might occur when starting values are determined totally at random.

Results
In the first pair of models (top part of Table 1), the outcome equation contains interaction effects between age classes and RI, whereas in the second pair (bottom part of the table) it does not. Conversely, the mediator equation always includes RI effects only, meaning that its selection strategy is not influenced by whether or not interactions are included in the outcome model (parameter estimates are in practice not sensitive to this change, too). The choice between the outcome model with and without interactions is not straightforward. Indeed, only the coefficient for the I{RI = 2} · I{AGE = 76+} product is significant. However, an overall ANOVA test comparing them (within the M-estimation framework) provides evidence in favor of the extended one (χ 2 deviance statistics equal to 11.706 with 4 degrees of freedom, p-value 0.020). Evidence is mixed also with regard to the Akaike and Bayesian Information Criteria, with the former favoring the interaction model (523.65 versus 527. 35 for the no interaction model), and the latter favoring the no interaction model (556.11 versus 571.58 for the interaction model). For this reason, both pairs are kept for reference. In order to ease interpretation, the outcome regression coefficients of AGE and RI are also linearly combined to obtain contrasts, on the log odds-ratio scale, between any pair of levels. These are reported, together with their standard errors, in Table 2, both for the interaction (top part) and no interaction (bottom part) setting. Clearly, in the former case the estimated contrasts of each factor are conditional on the level of the other. Tables 1 and 2 show that ML and M-estimation provide very similar results for both model pairs. In particular, the estimated coefficient for the effect of GAS intake on listeriosis development is always positive (very close to 0.99, p-value 0.002). This confirms that GAS intake remains positively associated with listeriosis, even when controlling for RI [31]. Similar conclusions can be drawn for RI effects due to both diseases (RI=1) and therapies (RI=2) on GAS intake, with all coefficients being around 0.84 (p-values lower than 0.002). In light of these results, we can conclude that there is no evidence of a difference between diseases and therapies with respect to their effect on GAS intake.
With regard to effects of age on the outcome, the interaction model shows an increasing pattern for the no RI (RI=0) group, although the shift from the 40 − 65 to the 66 − 75 age class is barely significant for M-estimation (p-value 0.061), and that from the 66 − 75 to the 76+ class is not significant (p-values greater than 0.25). For RI=1, a rather relevant positive difference is estimated for the contrast involving the second and the first age class. Conversely, the estimated difference between the third and the first class is smaller (and barely significant, p-values around 0.010). As a consequence, the risk for listeriosis in the second class is higher than that of the third one, with the difference being significant (p-value 0.030 for M-estimation and 0.042 for ML). For the RI=2 group, no significant differences among age classes are present. In the no interaction model, a unique pattern is estimated for all RI groups which is similar to the one described for the RI=1 group in the model with interactions, although with reduced magnitudes.
As for RI effects on the outcome, in the no interaction model we notice a significantly increased vulnerability to listeriosis when moving from RI=0 to both RI=1 and RI=2 (p-values lower than 0.001). The difference between the latter levels is significant, too, with immunosuppressive therapies further enhancing the risk for listeriosis with respect to immunocompromising diseases (p-values 0.002). Effect magnitudes are in line with the findings in the original study; see Preußel et al. [31]. In the model with interactions, such a pattern of effect directions is preserved within each age class. However, the shifts from RI=0 to RI=1 and from RI=1 to RI=2 are significant only for the second and the first class, respectively (p-values lower than 0.002). Overall, we can conclude that: i) only therapies have a significant effect in the 40 − 65 class, ii) differences between therapies and diseases are essentially null in the 66 − 75 class, and iii) RI effects considerably lessen in the 76+ class.
In terms of effect direction and statistical significance, results from the weighting approach are not extremely dissimilar to those previously reported. However, some noteworthy discrepancies emerge which concern GAS intake effects in the outcome model and RI effects in the mediator model. For the former, we observe a sensible increase in the coefficients with respect to M-and ML estimation. Moreover, it is important to remark that these estimates are somewhat less robust to model specification, with two rather distant values in the models with and without AGE-RI interactions (1.267 and 1.169, respectively). As for RI effects on GAS intake, the comparison with the other two approaches yields smaller coefficients, with no significant differences between any pair of levels.

Causal interpretation of effects
Causal mediation analysis requires a number of assumptions to hold [43]. These include consistency and composition [40,42], as well as more contextspecific assumptions concerning the absence of unobserved confounding for the following relationships: i) exposure-outcome, ii) mediator-outcome, iii) exposure-mediator. Also, iv) no mediator-outcome confounding should be generated by variables causally affected by the treatment. In this case, many effects of interest are not identified unless strong parametric assumptions are made, even when these variables are measured [43,29]. While consistency and composition are technical statements typically tenable in many applied settings, assumptions i)-iv) warrant some discussion about the introduced covariate sets. As mentioned in Section 5.1, for our application we directly rely on the data generating process of Preußel et al. [31], from which it seems reasonable to assume that the included covariates are sufficient for the assumptions i)-iii) to hold. Also, none of these variables is causally affected by the exposure, so that assumption iv) seems realistic, too.
When the outcome is binary, VanderWeele and coauthors [43,37] introduce the notion of Natural Direct Effect (NDE) and Natural Indirect Effect (NIE) on the odds-ratio scale. In particular, when also the mediator is binary, Valeri and VanderWeele [37] derive parametric expressions under the rare outcome assumption. Though this is also our context, we rely on the formulae in Doretti et al. (2022) [9], that extend them to settings where the exposure (and the mediator) interact with the covariates (also outside the rare outcome case). When the binary mediator does not interact with other variables, it can be shown that exposure effects in the outcome model (that is, RI effects in Table 2) are good approximations of log odds-ratio NDEs. We have also computed NIEs on the same scale: for our data, these depend very little on the covariate patterns (with differences at the sixth digit), even in the model with interactions.
With regard to M-and ML estimation, the NIEs related to the RI = 2 vs RI = 1 contrast are always very close to 0 and not significant. Clearly, this is due to the almost null difference between the effects of these two exposure levels in the mediator model (see the corresponding regression coefficients in Table 1). Conversely, when the reference level of RI is 0, NIEs lie around 0.068, with Delta-method [27] standard errors close to 0.042 (p-value 0.105). An exception is represented by the RI = 2 vs RI = 0 contrast in the setting with interactions, where the NIE has the same magnitude but the standard error is 0.037 (p-value 0.066). As for weighting, the NDE and NIE values are slightly different, reflecting the fact that the estimated regression coefficient differ from those of the other approaches; see the related discussion in Section 5.2. However, the interpretation of results is essentially the same. In Table 3, we report the proportion of RI effect mediated by GAS intake (for the 1 vs 0 and the 2 vs 0 contrasts) on the log odds-ratio scale. On such a scale, total causal effects are given by the sum of NDEs and NIEs.
This proportion is small but non-negligible, since it ranges from 2.4% to 14% across RI contrasts and age classes. Nevertheless, it is worth to remark that the significance of estimated NIEs is weak, and that total causal effects are non-significant whenever the corresponding NDEs also are.

Simulation study
In this section, we present a simulation study conducted to investigate the finite-sample properties of the M-and ML estimators introduced in Section 4. The performance of these estimators is compared to that of the weighting estimator, implemented for both the outcome and mediator model. In Section 6.1, the structure of the study is described, while results are reported in Section 6.2.

Results
The regression parameter estimatesβ andδ of the ML and M-estimator are quite similar and show that these estimators perform slightly better than the weighting estimator; see Table 4 for the setting with 5 controls per case (with 8 controls similar results are obtained). Importantly, these findings suggest that the estimators derived within our analytical framework generally overcome the weighting approach not only for the outcome model parameters (see Section 4), but also for the mediator model parameters. The ML estimator is the most efficient when the maximization is stable. For 21 out of 1000 repetitions of Scenario 2, the maximization was however unstable, resulting in one or more extreme estimated standard errors (up to 3.5 million times the size of the corresponding standard error of the weighting estimator). The M-estimator is more robust compared to the ML estimator when it comes to such extreme results and performs slightly better compared to the weighting estimator when it comes to precision; see Figure 3. The bias is very low for almost all parameters, with two exceptions. First, the I{RI = 1}·I{AGE = 76+} coefficient (scenario 2); the empirical coverage is however as desired (at least 95%). Second, the intercept from the outcome model of both scenarios and sample sizes is biased for the weighting estimator, often resulting in slightly too low empirical coverage (see the values in Table 4 for the setting with 5 controls per case, while results for the larger sample size are not reported).

Discussion
We consider Stratified Case-Control (SCC) designs defined on a certain variable (primary outcome) and we address the problem of modeling another binary variable (secondary outcome) from the resulting dataset. In detail, we assume that a logistic regression model for the secondary outcome holds in the target population, and we recover the parameters of such a model via a simple offset adjustment. Like in other frameworks, knowledge of the (conditional) prevalence(s) of the primary outcome in the population is required.
Our analytical solution allows one to perform mediation analysis from SCC data in settings where both the outcome and the mediator are binary. As shown in Section 4, the derived offset correction leads to M-and Maximum Likelihood (ML) estimation of the joint (i.e., outcome and mediator) model parameter vector, provided that the intercept and/or the coefficients of background variables in the outcome model are suitably adjusted. These estimation methods are opposed to weighting, where such adjustments are not needed. With this respect, a noticeable finding of the simulation study in Section 6 is that the efficiency gain of M-and ML estimation compared to weighting is not limited to the outcome model vector, but it extends to the mediator model parameters. Although we have explored a single scenario (close to that of the application described in Section 5), we reasonably believe that a similar tendency would be observed for alternative parameter configurations, too.
The proposed approach is exemplified via the analysis of a CC dataset gathered within a German study on listeriosis conducted by Preußel et al. [31]. Specifically, we take a novel perspective and focus on evaluating whether the effect of Reduced Immunocompetency (RI) on listeriosis development is mediated by the intake of Gastric Acid Suppressors (GASs). We try to answer this question also in causal terms by decomposing the total causal effect of RI on listeriosis into the Natural Direct Effect (NDE) and the Natural Indirect Effect (NIE), though we acknowledge that such a terminology might be somewhat questionable from a strict causal mediation standpoint. This is because RI and GAS intake do not indeed cause listeriosis; the ingestion of food contaminated with Listeria monocytogenes bacteria (predominantly) or other unknown vehicles (to a lesser extent) do. Nevertheless, since infection means are not of primary interest, the disentanglement of the RI effect on listeriosis in causal terms appears sensible, due to the clear additional vulnerability to listeriosis for people with RI as well as for people taking GASs, and to the fact that GAS intake is a direct consequence of RI. Results show that the NIE of RI is almost negligible.
It is important to stress that any conclusion from the reported analyses has to be drawn cautiously, for several reasons. First, in the original study by Preußel et al. [31], 732 patients with listeriosis were potentially eligible as cases. Of these, only 109 entered the study due to participation refusal, decease, inability to answer/be contacted and other unknown reasons. Also, like in many telephone surveys recruited controls were found to systematically differ from the target population with respect to socio-economic status. Although included and not included cases present not too dissimilar distributions for many variables, and socio-economic status is essentially accounted for via the (initial) inclusion of educational level as its proxy [31], we cannot exclude that some degree of non-informative missingness [24] might affect the results.
Furthermore, as discussed in Section 5.2.1, though the typical assumptions of causal mediation analysis i)-iv) seem to be tenable in our application, it is worth to remark that these requirements about confounding might be more difficult to justify in other contexts. To deal with violations of i)-iii), sensitivity analyses [20] as well as methods embedded in an instrumental variable framework [8,23] have been introduced which could be possibly adapted to the present setting. As for iv), alternative causal estimands identifiable also in the presence of exposure-induced mediator-outcome confounding (named interventional effects) were proposed [44,45].
It is well-known that in SCC designs the number of strata can in principle increase until the limit situation of containing a single case (sometimes termed exact matching). In these settings, traditional estimation might be problematic, with conditional likelihood methods typically preferable [11,5]. Extension to deal with these cases can be of interest.
A Proof of (6) and (7) By standard probability results it is possible to express the left-hand side of (6) as For the outcome model (1), the partition x y = (x ⊤ y 0 , x ⊤ y 1 ) ⊤ is introduced, where x y 0 (x y 1 ) denotes the sub-vector of covariate values not involving (involving) M. The coefficient vector partition β = (β ⊤ 0 , β ⊤ 1 ) ⊤ is defined accordingly. We assume our model to be hierachical, and therefore in the most general setting x y 1 = m · x y 0 . However, in many applications β 1 and x y 1 are likely to be smaller-dimensional than β 0 and x y 0 , with some elements of the latter vectors not having a counterpart in the former. In these cases, the additional vectorβ 1 has to be defined by suitably expanding β 1 with zeros, so that the sum of conformable vectors β + = β 0 +β 1 can be computed. It then follows that Then,

B Details of M-estimation
In order to obtain the expressions of the observed score function vectors s y (·) and s m (·), it is worth to rely on matrix notation. Specifically, we index sample units by i = 1, . . . , n and introduce the sample vectors m = (m 1 , . . . , m n ) ⊤ and y = (y 1 , . . . , y n ) ⊤ , as well as the sample design matrices We also denote by d δ , d β , d β 0 and d β 1 the number of columns of X m , X y , X y 0 and X y 1 , respectively. Letting expit(·) = exp(·)/{1 + exp(·)}, for a given value θ = (β ⊤ , δ ⊤ ) ⊤ we have where β ⋆ is like in Section 3 and o is a column vector collecting the offset terms o(y i , x y 0 ,i ; β) of every sample unit i = 1, . . . , n. Notice that the vector p my = (p my,1 , . . . , p my,n ) ⊤ contains the p my,i probabilities (i = 1, . . . , n) already introduced in Section 4.2, while p y ⋆ = (p y ⋆ ,1 , . . . , p y ⋆ ,n ) ⊤ differs from p y ⋆ (m) = (p y ⋆ (m) ,1 , . . . , p y ⋆ (m) ,n ), which is the vector collecting the p y ⋆ (m) ,i probabilities (also used in 4.2). A compact form for ψ(θ) = (s y (θ) ⊤ , s m (θ) ⊤ ) ⊤ is given by The formula for the A(y, θ) matrix is where each generic sub-matrix is defined as H qα = ∂s q /∂α ⊤ and In the above, the D matrix is given by whereX y 0 is the matrix obtained by extracting the columns of X y 0 corresponding to the non-null elements ofβ 1 while v 0 and v 1 are two column vectors of length n, whose i-th element is given by .
Notice that in the formulae above · denotes element-wise product, be it between column vectors or between a matrix and a column vector. In the latter case, the result is a conformable matrix where every column is element-wise multiplied by the vector. Finally, B(y, θ) is given by with ψ i (θ) being the contribution of unit i to the observed score vector.

C Parametric expression of (11)
The parametric expression of g(x y 0 , x m ; β, δ) can be derived from (12) and (13). After some simplifications, it follows that where the elements in the right-hand side are defined as in Appendix A. We then have with β ⋆ 0 being the modification of β 0 accounting for the correction for β INT and β b .

D Gradient function for ML estimation
The log-likelihood gradient s(θ) can be conveniently expressed in matrix form. To this end, we maintain notation as in the body of the paper and in previous appendices, and we introduce the matrices where In the above, ⊗ denotes Kronecker product, while ∂β INT g(x y 0 ,n , x m,n ; β, δ) G βm = ∂ ∂βm g(x y 0 ,1 , x m,1 ; β, δ) · · · ∂ ∂βm g(x y 0 ,n , x m,n ; β, δ) G δ INT = ∂ ∂δ INT g(x y 0 ,1 , x m,1 ; β, δ) · · · ∂ ∂δ INT g(x y 0 ,n , x m,n ; β, δ) are three row vectors collecting the derivatives of the g(x y 0 ,i , x m,i ; β, δ) terms (i = 1, . . . , n) with respect to β INT , β m and δ INT , where the latter two are the coefficient of M in the outcome model and the intercept of the mediation model, respectively. These derivatives can be computed by noticing that the argument of the logarithm can be written as (q 1 q 2 q 3 + q 4 )/(q 2 q 3 + q 4 ), with q 1 = exp(x ⊤ y 0 ,iβ 1 ), q 2 = exp(x ⊤ m,i δ), q 3 = 1 + exp(x ⊤ y 0 ,i β 0 ) and q 4 = 1 + exp(x ⊤ y 0 ,i β + ). Consequently, we have Finally, letting the log-likelihood gradient is given by   Table 2 continues.  Table 3: Proportion of RI effect mediated by GAS intake Figure 3: Box plots of the quotient between the estimated standard errors of the ML / M-estimator and the weighting estimator for the regression parameters β and δ (i.e. 1000 replications of 9 or 13 parameters). Note that we have changed the limits of the y-axis in order to make the figure more clear, the number of observations above the window is printed above each box.