Delicatessen: M-Estimation in Python

M-estimation is a general large-sample statistical framework [1]. Widespread application of M-estimation is hindered by tedious derivative and matrix calculations. To address this barrier, we developed delicatessen, an open-source Python library that flexibly automates the mathematics of M-estimation. Delicatessen supports both built-in and custom, user-specified estimating equations. To help contextualize the value of M-estimation and delicatessen, we showcase several applications common to statistical analyses in life sciences research: regression with outliers, dose-response curves, and standardization.

To begin, we briefly review M-estimation; see Stefanski and Boos (2002) for a thorough introduction [1]. An M-estimator, $\hat{\bm{\theta}}$, is the solution for $\bm{\theta}$ to the vector equation $\sum_{i=1}^{n}\psi(\mathbf{Z}_{i};\hat{\bm{\theta}})=0$, where $\bm{\theta}=\{\theta_{1},\theta_{2},...,\theta_{v}\}$ is the parameter vector, $\psi(.)$ is a ($v\times 1$)-dimension estimating equation, and $\mathbf{Z}_{i}$ are the observed data for $n$ independent and identically distributed units. As an example, the M-estimator for the mean ($\mu$) of a variable ($Y_{i}$) is $\sum_{i=1}^{n}\left(Y_{i}-\mu\right)=0$, which is equivalent to the usual mean estimator, $\hat{\mu}=n^{-1}\sum_{i=1}^{n}Y_{i}$. To find $\hat{\bm{\theta}}$, root-finding procedures can be used, which iteratively search for the value(s) of $\hat{\bm{\theta}}$ where $\sum_{i=1}^{n}\psi(\mathbf{Z}_{i};\hat{\bm{\theta}})=0$. The variance of $\hat{\bm{\theta}}$ can be estimated using the empirical sandwich variance estimator (see Appendix). Key advantages of the sandwich estimator are reduced computational burden for variance estimation relative to common alternatives (e.g., bootstrapping [2, 3], Monte Carlo methods [4], etc.), automation of the delta method, and valid variance estimation for parameters that depend on other estimated parameters [1].

M-estimators are implemented in delicatessen via the MEstimator class. Input to MEstimator includes the estimating equation(s) and starting values for the root-finding algorithm. Root-finding algorithms implemented in SciPy [5], as well as custom root-finding algorithms, are supported. After finding $\hat{\bm{\theta}}$, the ‘bread’ of the sandwich estimator is calculated by numerically approximating the partial derivatives via the central difference formula and the ‘filling’ of the sandwich is calculated using NumPy [6]. Finally, the sandwich covariance matrix is constructed. While other M-estimation implementations exist [7, 8, 9], delicatessen provides greater support for both user-specified and pre-built estimating equations (Appendix Table 1).

To demonstrate application of delicatessen, we provide three examples common to life sciences research. First, consider linear regression with outliers [10, 11]. A common approach to handling outliers is to exclude them. However, exclusion ignores all information contributed by outliers, and should only be done when outliers are unambiguously a result of experimental error [11]. Yet, including outliers with simple linear regression can lead to unstable or unreliable estimates. Robust regression has been proposed as an alternative, whereby outliers contribute but their influence is curtailed [12, 13]. The estimating equations for the intercept and slope with robust linear regression are

$$\psi(Y_{i},X_{i},\bm{\alpha})=\begin{bmatrix}f_{k}(Y_{i}-(1,X_{i})\bm{\alpha}^{T})1\\ f_{k}(Y_{i}-(1,X_{i})\bm{\alpha}^{T})X_{i}\end{bmatrix}$$

where $Y_{i}$ is the independent variable, $X_{i}$ is the dependent variable(s), $\bm{\alpha}=(\alpha_{0},\alpha_{1})$ is the parameter vector for the regression model, and $f_{k}(\cdot)$ truncates or bounds the residuals at $-k,k$. To illustrate, 15 observations were simulated following Altman & Krzywinski (2016) [11] and all models were fit using delicatessen with built-in estimating equations. As a reference, a linear model was fit to the simulated data (Figure 1a, original reference slope: $0.72$). An outlying value of $Y_{i}$ was then generated by replacing $Y_{i}$ with $Y_{i}+3$ for the observation with the smallest $X_{i}$ value. When fitting a linear model with this outlier, the slope is noticeably decreased ($0.52$). With robust regression and $k=1.345$ [14], the influence of the outlier was curtailed and the estimated slope ($0.62$) was closer to the reference.

As the second example, consider estimating the dose-response curve of herbicides on ryegrass root length [15, 16]. Since there is a natural lower bound of zero on root length, a three-parameter log-logistic (3PL) model was used. Parameters of the 3PL model are the halfway maximum effective concentration ($\gamma_{1}$), steepness of the dose-response curve ($\gamma_{2}$), and upper limit on the response ($\gamma_{3}$). Let $\psi_{PL}(D_{i},R_{i};\bm{\gamma})$ denote the estimating equations (see Appendix), where $D_{i}$ is the dose and $R_{i}$ is the response. Besides the halfway maximum effective concentration, other concentration levels may also be of interest. Let $\delta_{a}$ denote the $a$% effective concentration, where $0<a<100$. Here, estimation of the 20% effective concentration was of interest, with the estimating equation $\psi_{EC}(\delta_{20},\bm{\gamma})$. As $\delta_{a}$ is a transformation of the other parameters, this estimating equation does not depend on the observations but is instead a delta-method transformation of the parameters. Here, the estimating equations can be ‘stacked’ together

$$\psi(R_{i},D_{i},\bm{\theta})=\begin{bmatrix}\psi_{PL}(R_{i},D_{i};\bm{\gamma})\\ \psi_{EC}(\delta_{a},\bm{\gamma})\end{bmatrix}$$

where $\bm{\theta}=(\bm{\gamma},\delta_{a})$. The variance for $\delta_{a}$ is automatically estimated through the sandwich. The sandwich variance can then be used to construct Wald-type 95% confidence intervals (CI). The estimated dose-response curve using delicatessen and the built-in estimating equations is shown in Figure 1b. The estimated 20% effective concentration was 1.86 (95% CI: 1.58, 2.14). This example highlights how M-estimation can be used to apply the delta method, with delicatessen automating the process.

Finally, consider the problem of standardization. Often the available data is not a random sample of the population. Consider the biomarker data from a convenience sample of HIV patients ($n=57$) in Kamat et al. (2012) [17]. When comparing the prevalence of recent drug use (either cocaine or opiates) to a more generalizable cohort,[18] drug use was notably higher in Kamat et al.’s sample (70% versus 8%). As Kamat et al. reported differential biomarker expression by cocaine use, differential patterns in drug use between data sets indicates that the summary statistics on biomarker expression may not be generalizable. To account for the discrepancy in drug use, we standardized the mean for select biomarkers to the distribution of drug use in the secondary data source using inverse odds weights [19]. Inverse odds weights were estimated using a logistic regression model [20], with the estimating equations

$$\psi_{w}(S_{i},X_{i},\bm{\beta})=\begin{bmatrix}\left(S_{i}-\text{expit}((1,X_{i})\bm{\beta}^{T})\right)1\\ \left(S_{i}-\text{expit}((1,X_{i})\bm{\beta}^{T})\right)X_{i}\end{bmatrix}$$

where $X_{i}$ indicates drug use for individual $i$, $S_{i}$ indicates if the individual was in the Kamat et al. study ($S_{i}=1$) or in the cohort ($S_{i}=0$), and $\text{expit}(.)$ is the inverse logit. The estimating equation for the weighted mean is

$$\psi_{m}(B_{i},S_{i},\mu_{m},\bm{\beta})=I(S_{i}=1)\times w_{i}(X_{i};\bm{\beta})\times(\text{log}(B_{i}^{m})-\mu_{m})$$

where $\text{log}(B_{i}^{m})$ is the log-transformed biomarker $m$, $\mu_{m}$ is the mean for biomarker $m$, and

$$w_{i}(X_{i};\bm{\beta})=\left(1-\text{expit}(X_{i}^{T}\bm{\beta})\right)/\text{expit}(X_{i}^{T}\bm{\beta})$$

is the inverse odds weight. While $\bm{\beta}$ is not of primary interest, the estimates of $\mu_{m}$ depend on $\bm{\beta}$ through the inverse odds weights. This dependence means that the estimated variance of $\bm{\beta}$ should carry forward into the estimated variance of $\mu_{m}$. Ignoring this dependence puts us in danger of underestimating the uncertainty in the means for biomarker expression. M-estimators address this issue via the sandwich, where the uncertainty of parameters are propagated through the stacked estimating equations. Therefore, estimating equations for the logistic model and for the biomarker means were stacked together and $\bm{\theta}=(\bm{\beta},\mu_{1},\mu_{2},...,\mu_{m})$. The stacked estimating equations were implemented in delicatessen by a combination of built-in and user-specified estimating equations. Results for select biomarkers are presented in Figure 1c. Notably, the means for log-transformed sIL-2R and IL-12 decreased after standardization. While these results were standardized by drug use, definitions varied between studies and other important differences between the populations the samples were drawn from likely exist. Therefore, these results should only be viewed as an illustration of how M-estimation can be used.

To summarize, M-estimation is a flexible framework. To automate the mathematics of M-estimators, we developed delicatessen, which supports both pre-built and user-created estimating equations. Key features of delicatessen were highlighted through examples in life science research. Further examples can be found on our website.