Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
75,437 characters · 11 sections · 55 citation commands
Penalized Likelihood Inference with Survey Data
Survey data are widely used in many disciplines of social sciences. The statistical methodology for survey samples has been well-developed and culminated in a large body of literature Cameron-Trivedi(2009), Wooldridge(2010),Fuller(2011),Thompson(2012).
Despite the rapid development of machine learning/high-dimensional econometrics and the increasing availability of big datasets in recent years, the research on how to adapt/apply these high-dimensional statistical methods to survey data has been lagging. This paper aims to fill this gap in the literature by providing extensions of Lasso inference methods to survey environment. Our hope is to enrich the toolbox of practitioners who want to apply the high-dimensional regression methods to survey data.
Most prediction-oriented methods, including the Lasso, trade off bias and variance, and consequently, deliver a biased estimate that is not suitable for making inference on the model coefficients. Several post-Lasso-selection inference methods that mitigate this shortcoming have been proposed in the literature. Among others, Zhang-Zhang(2014) and Javanmard-Montanari(2014) propose a debiased Lasso (DB) method which is based on one-step iteration of the initial Lasso estimator. Belloni-Chernozhukov-Wei(2016) propose double selection and $C(\alpha)$-type methods in a generalized linear model (GLM) that satisfies sparsity assumptions. The latter is based on an estimating equation orthogonalized against the nuisance parameter “score" function.
Lee-Sun-Sun-Taylor(2016) propose a selective inference (SI) method for the parameters in a linear model selected by the Lasso. The method is extended to a homoskedastic GLM by Taylor-Tibshirani(2018). In SI, the target parameters are determined from the data as opposed to being fixed before the the selection events. This feature makes the post-selection method conceptually different from the $C(\alpha)$ and DB methods, where the target parameters are the population parameters.
This paper presents two rather straightforward results. We first extend the $C(\alpha)$, DB and SI methods to a GLM estimated by the Lasso to accommodate survey weights and/or heteroskedasticity. The survey framework we adopt is similar to that of Wooldridge(2001). Accounting for survey weights naturally leads to conditional heteroskedasticity which, in turn, brings about an extra challenge because the active and inactive constraints of the Karush-Kuhn-Tucker condition for the Lasso problem are no longer asymptotically independent, and conditioning only on the active constraints as considered by Taylor-Tibshirani(2018) for a homoskedastic GLM may lead to invalid inference.
Second, we establish the asymptotic validity of the above three methods for inference on nonlinear parameter functions such as the average marginal effects (AMEs) in a survey logit model. There exist very few studies on the application of Lasso methods to survey data. Mcconville-etal(2017) consider a survey-weighted linear Lasso regression and develop a finite population asymptotic theory for Lasso estimators with a fixed number of regressors. In contrast, we consider a survey-weighted GLM and establish the asymptotic validity inference procedures under the usual infinite population framework, see e.g. Wooldridge(2001),Wooldridge(2010) and Cameron-Trivedi(2009) for the latter. Additionally, we allow for a growing number of covariates in the survey extensions of the debiased Lasso and $C(\alpha)$ methods.
The paper is organized as follows. Section (ref) lays out the model framework. We propose extensions of the selective inference, debiased Lasso and $C(\alpha)$/orthogonalization methods in Section (ref). Section (ref) applies the proposed methods to inference on AMEs in a survey logit model. Section (ref) provides a simulation evidence on the properties of the proposed methods and Section (ref) presents an empirical application to Canadian Internet Use Survey 2020 data. We conclude in Section (ref).
\paragraph*{Notations and terminology} Let $1(\cdot)$ denote the indicator function, and $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the smallest and the largest eigenvalue of a symmetric matrix $A$, respectively. For a $k\times 1$ vector $a=(a_1,\dots, a_k)'$, we define $\Vert a\Vert_0\equiv \mathrm{supp}(a)$ (the number of nonzero components of the vector $a$) and $\Vert a\Vert_1\equiv \sum_{i=1}^k\vert a_j\vert$. For a real matrix $A=(a_{ij})$, let $\Vert A\Vert_{\infty}\equiv \max_{i,j}\left|a_{ij}\right|$, and $\Vert A\Vert=\sqrt{\mathrm{tr}(A'A)}$ and $\Vert A\Vert_2=\sqrt{\lambda_{\max}(A'A)}$ denote its Frobenius and spectral norms, respectively. The sub-Gaussian norm of a random variable $X$ is defined as
A random variable $X$ is called sub-Gaussian if $\Vert X\Vert_{\psi_2}\leq C<\infty$ for a constant $C>0$. A random vector $X\in \mathbb{R}^p$ is called sub-Gaussian if the one-dimensional marginals $X'b$ are sub-gaussian random variables for all $b\in\mathbb{R}^p$. The sub-Gaussian norm for the random vector is defined as $\Vert X\Vert_{\psi_2}\equiv \sup_{\Vert b\Vert=1}\Vert X'b\Vert_{\psi_2}$. The sub-exponential norm of a random variable $X\in \mathbb{R}$ is
Moreover, let $\bm{1}_{m}=(1,\dots, 1)'$ and $0_{m}=(0,\dots, 0)'$ denote the $m\times 1$ vector of ones and zeros, respectively, and $e_{jm}$ denote the $m\times 1$ unit vector whose $j$-th element is $1$ and the remaining elements are $0$. Let $F(x; \mu, \sigma^2, a,b)$ denote the CDF of a ${N}(\mu,\sigma^2)$ random variable truncated on the interval $[a,b]$, that is,
where $\Phi(\cdot)$ is the CDF of a $N(0,1)$ random variable. Also, let $\Lambda(z)\equiv \exp(z)/(1+\exp(z))$ denote the CDF of logistic distribution. We abbreviate central limit theorem and continuous mapping theorem as CLT and CMT, respectively.
We consider a GLM that specifies the conditional density of a scalar outcome variable $y_{i}$ given a $(p+1)\times 1$ vector of covariates $x_{i}$ which includes a constant as
where $\theta_0$ is the true value of the parameter vector $\theta\in\mathbb{R}^{p+1}$, and $a(\cdot)$ and $c(\cdot)$ are known functions. To each vector of observations $(y_i, x_i')', i=1,\dots, n,$ there corresponds a positive, bounded survey weight denoted as $w_{i}, i=1,\dots, n$.\footnote{In our framework, $\{(y_i, x_i', w_i)'\}_{i=1}^n$ actually forms a triangular array $\{\{(y_{ni}, x_{ni}', w_{ni})'\}_{i=1}^n$. We drop the index $n$ for notational simplicity.} Let $g(y, x'\theta)\equiv -\log f(y, x'\theta)$ and define the weighted log-likelihood function as follows:
As is well known, the weighted likelihood framework is commonly used in survey data analysis Manski-Lerman(1977),Cameron-Trivedi(2009),Wooldridge(2010), and accommodates, among others, the following stratification schemes.
Wooldridge(2001) established the asymptotic properties of $M$-estimator under the above two sampling schemes. We use the same sampling schemes to establish the asymptotic validity of the Lasso-based inference methods described below.
The score function, the sample information and negative Hessian matrices corresponding to (ref) are defined as
Moreover, we define $H(\theta_0)\equiv \operatorname{E}[\hat{H}(\theta_0)]$ and $I(\theta_0)\equiv \operatorname{E}[\hat{I}(\theta_0)]$.
Let us partition $x_i=(1, \tilde{x}_i')'\in\mathbb{R}^{p+1}$ and $\theta=(\alpha, \beta')'\in\mathbb{R}^{p+1}$, where $\tilde{x}_i=(\tilde{x}_{i1},\dots, \tilde{x}_{ip})'\in\mathbb{R}^p$, $\alpha\in\mathbb{R}$ and $\beta\in\mathbb{R}^p$ so that $x_i'\theta=\alpha+\tilde{x}_i'\beta$.
In this paper, the variable selection, estimation and inference are performed using a survey-weighted Lasso where the negative of the weighted log-likelihood function (ref) is minimized subject to $\ell_1$ penalty on the slope parameters:
where $\lambda\geq 0$ is a tuning parameter. Note here that, as it is standard in the Lasso literature, only the “slope" parameters in $\beta=(\beta_1,\dots, \beta_p)'$ are penalized. The $j$-th elements of $\theta$ and $\theta_0$ are denoted as $\theta_{(j)}$ and $\theta_{0(j)}$, respectively.
Hereafter, $M\subseteq\{1,\dots, p+1\}$ denotes the subset of regressors that includes the constant term and non-constant regressors with a vector of (non-zero) Lasso estimates $\hat{\beta}_M\in\mathbb{R}^{|M|-1}$ and $\hat{\mathtt{s}}_M\equiv \mathrm{sign}(\hat{\beta}_M)\in\{-1,1\}^{\vert M\vert-1}$. Also, let ${\beta}_M\in\mathbb{R}^{|M|-1}$ be the subvector of $\beta$ corresponding to $M$, ${\theta}_M=({\alpha}, {\beta}_M')'$ and $\hat{\theta}_M=(\hat{\alpha}, \hat{\beta}_M')'$.
The Lasso solution in (ref), with $\lambda$ fixed, returns a random subset of regressors $\hat{M}\subset\{1,\dots, p+1\}$. Since the intercept $\alpha$ is not penalized, $\hat{M}$ always includes the constant term. The target parameter vector in the selective inference considered in Section (ref) is ${\theta}_{M0}=({\alpha}_0, {\beta}_{M0}')'$, the true value of $\theta_M$ in the selected model $\hat{M}=M$. In contrast, the DB and $C(\alpha)$ in Sections (ref)--(ref) target the entire vector $\theta_0$.
Let $\mathtt{s}_M\equiv \mathrm{sign}({\beta}_{M0})\in\{-1,1\}^{\vert M\vert-1}$, and denote by $m_0\equiv \Vert\theta_0\Vert_0$ the number of nonzero elements of $\theta_0$. Moreover, let $\theta_{-M}\in\mathbb{R}^{p+1-|M|}$ be the subvector of parameters other than $\theta_M$ and $L_M(\theta_M)$ be the weighted log-likelihood function for the selected model with parameters $\theta_M$. It is clear that $L_M(\theta_M)$ can be obtained by evaluating $L(\theta)$ at $\theta^{*}$ whose non-zero elements are $\theta_M$ and remaining $p+1-|M|$ elements are 0. We partition (ref)-(ref) as follows:
where $\hat{H}_M(\theta_M)\in\mathbb{R}^{|M|\times |M|}$ and $S_M(\theta_M)\in\mathbb{R}^{|M|}$ denote the negative Hessian matrix and score functions corresponding to $\theta_M$, respectively. With the partitioning above, Lee-Sun-Sun-Taylor(2016) and Taylor-Tibshirani(2018) show that the event $\{\hat{M}=M, \hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M\}$ holds if and only if there exist random vectors $\hat{\theta}_M\in\mathbb{R}^{|M|}$ and $\mathtt{u}\in\mathbb{R}^{p+1-|M|}$ in the Karush-Kuhn-Tucker condition for the problem (ref) such that
We establish the asymptotic validity of the three inference methods under the following assumptions imposed directly on the loss function $g(y, t)$ which are similar to the assumptions employed in vandeGeer-etal(2014) and Xia-Nan-Li(2021).
The boundedness of the variables stated in Assumption (ref)(ref) is employed frequently in the literature, see Negahban-etal(2012), vandeGeer-etal(2014) and Xia-Nan-Li(2021). To deal with survey samples, we relax the i.i.d. assumption used in these papers, and although the proofs of validity of the inference procedures considered below are quite standard, much of the effort of the proof goes into verifying that the same results that hold in an i.i.d. setup carries over to independent non-identically distributed (i.n.i.d.) samples.
The weight $w_i$ is deterministic but we require that it is bounded from above and below away from $0$, so it rules out strata that become asymptotically degenerate. Moreover, the weights do not need to sum to 1. In the R package glmnet, the weights are rescaled to sum to $n$.
Assumption (ref)(ref) is a mild condition that ensures nonsingularity of the Hessian and information matrices in the case of slowly diverging number of covariates considered below. Assumption (ref)(ref) is standard and requires the convexity and boundedness of the first two derivatives and Lipschitz continuity of the second derivative of $g(y, t)$ with respect to $t$ uniformly in a neighborhood of $x_i'\theta_0$ (see vandeGeer-etal(2014) and Xia-Nan-Li(2021)).
The condition (ref) is essentially the the Quadratic Margin Condition needed for the consistency of the Lasso Buhlmann-vandeGeer(2011) and see also Negahban-etal(2012) for a related (stochastic) Restricted Strong Convexity condition. A sufficient condition for (ref) is that $\ddot{g}(y, x'\theta)$ is bounded away from zero locally around $x_i'\theta_0$ for all $i=1,\dots,n$.
In this section, we extend the selective inference argument of Taylor-Tibshirani(2018) for a homoskedastic GLM to a GLM with survey weights and/or heteroskedasticity. In the SI, the target parameters are the coefficients selected by the Lasso. As a result, they are random before the selection, but not so conditional on the Lassso selection events. This feature distinguishes the selective inference method from the $C(\alpha)$ and debiased Lasso inference where the target parameters are the population parameters. See Lee-Sun-Sun-Taylor(2016) for further discussions about the difference between the SI and other inference methods. As in Taylor-Tibshirani(2018), we fix $\lambda>0$ and consider the following one-step estimator
where $S_M(\hat{\theta}_M)=(0, \lambda\, \mathtt{s}_M')'$, from which we obtain the one-step estimator of $\beta_{M0}$:
The SI is based on the asymptotic distribution of $ \tilde{\beta}_M$ conditional on the selection event $\hat{M}=M$ and $\hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M$.\footnote{Lee-Sun-Sun-Taylor(2016) also propose a test statistic which is conditional on $\hat{M}=M$ only by taking the union of the events characterized by polyhedral constraints over all possible combinations of the signs of the selected coefficients.} From (ref) and (ref), it follows that
hence
As argued by Taylor-Tibshirani(2018) (see Equation (21) therein), in a homoskedastic GLM, the random quantities appearing in the active and inactive constraints (ref) and (ref) are asymptotically independent (after suitable normalizations). However, this no longer holds in our setup because the covariance matrix of the limiting Gaussian random variables is not block-diagonal in the presence of survey weights and heteroskedasticity. This entails conditioning not only on the active constraints, but also on the inactive constraints. In light of this, we next derive an affine constraint corresponding to (ref). Let
By the fact that $\Vert \texttt{u}\Vert_{\infty}<1$, and after some algebra, we can express the inactive constraints in (ref) as follows:
where the inequalities hold element-wise. The Lasso selection events in (ref), (ref) and (ref) can be rewritten in a compact form as
where
Assuming that the number of non-constant regressors, $p$, is fixed, one can establish the asymptotic normality of the one-step estimators before the Lasso selection (see Section (ref)):
where $\Sigma$ is a $p\times p$ asymptotic covariance matrix. The estimator of $\Sigma$ is
where
From (ref), we have the distributional approximation for $Z$
The latter combined with the affine constraints $\{AZ\leq b\}$ in (ref), is now amenable to application of Lemma (ref) below, which summarizes two key results of Lee-Sun-Sun-Taylor(2016) (Lemma 5.1 and Theorem 5.2). To describe the lemma, we define the following quantities for a general $k\times 1$ random vector $Z$, and $A\in\mathbb{R}^{k\times k}$, $b\in\mathbb{R}^k$ and $\eta\in\mathbb{R}^k$:
where $(Ac)_j$ denotes the $j$-th element of $Ac$. Lee-Sun-Sun-Taylor(2016) show the following the result.
In our setup, $b$ defined in (ref) is random whereas Lemma (ref) assumes constant $b$. In addition, we have an approximate normality in (ref) instead of the exact normality assumed in Lemma (ref). These lead to an asymptotic version of (ref), namely, as $n\to\infty$
(ref) can be established using the results of Markovic-etal(2017). Although we do not directly use (ref), it provides the basis of the inference procedures described below.
Suppose we wish to make inference on the $j$-th element of ${\beta}_{M0}$, $j=1,\dots, |M|-1$ (conditional on the Lasso selection event $\hat{M}=M$ and $\hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M$). Let $Z, A$ and $b$ be as in (ref), $\mu$ be as in (ref), and set $\eta= e_{jp}\in\mathbb{R}^{p}$, $c=c(Z, \eta)$ and $r=r(Z, \hat{\Sigma}, e_{jp})$ in (ref), where $\hat{\Sigma}$ is defined in (ref). Fix $\zeta\in(0,1)$. The SI confidence interval (CI) of level $1-\zeta$ is of the form $\mathrm{CI}_{\hat{M}j}\equiv [\tilde{q}_l, \tilde{q}_u]$, where $\tilde{q}_l$ and $\tilde{q}_u$ are the solutions to the following equations
The asymptotic validity of the above CI is established in the following proposition.
See Appendix (ref) for a proof. The assumption of fixed $p$ is commonly used in the literature on SI Lee-Sun-Sun-Taylor(2016),Tian-Taylor(2017),Taylor-Tibshirani(2018),KKK(2022). Taylor-Tibshirani(2018) provide a heuristic argument for the validity of the selective inference in a homoskedastic GLM. Proposition (ref) extends their argument to i.n.i.d. and possibly heteroskedastic survey samples. The asymptotic validity of the SI procedures typically entails showing that CLTs that hold before selection extend to selective inference under suitable assumptions Tian-Taylor(2017),KKK(2022). We establish the asymptotic validity of the selective inference procedure by verifying the conditions given in KKK(2022). \paragraph*{Inference on a nonlinear parameter function.} Next we consider inference on a scalar nonlinear parameter function $\rho_M(\theta_{M0})$ (which may depend on $n$) in the selected model with coefficients $\beta_M$ on the active variables. Such results are especially useful in the context of logit and probit models because the AMEs are often the objects of interest therein. Analogously to (ref), consider the one-step estimator
Standard arguments yield the distributional approximation
Again, an approach similar to those applied to the elements of $\beta$ allows us to define the augmented variables:
where $Z$, $A$ and $b$ are as defined in (ref), and
Then, for $\zeta\in(0,1)$, the level $1-\zeta$ CI for $\rho_M(\theta_{M0})$ can be constructed as in (ref) and (ref) by replacing $A$, $Z$, $b$ and $e_{jp}$ by $A_\rho$, $Z_\rho$, $b_\rho$ and $e_{j(p+1)}$, respectively, and letting $r_\rho=r(Z_\rho, \hat{\Sigma}_\rho, e_{j(p+1)})$ in (ref). We can also infer the parameter $\rho_M(\theta_{M0})$ by conditioning on the sign of the estimated parameter $\rho_{\hat{M}}(\hat{\theta}_{\hat{M}})$ in addition to the event $\{AZ\leq b\}$ considered previously in (ref). To this end, let $\mathtt{s}_M^\rho\equiv \mathrm{sign}(\rho_M(\theta_{M0}))$ and $\hat{\mathtt{s}}_M^\rho\equiv \mathrm{sign}(\rho_M(\hat{\theta}_M))$ and redefine
and keep $Z_\rho$ and $\hat{\Sigma}_\rho$ as defined in (ref) and (ref). Then, we can rewrite the event $\{\mathtt{s}_M^\rho=\hat{\mathtt{s}}_M^\rho\}=\{\mathtt{s}_M^\rho= \mathrm{sign}(\rho_M(\hat{\theta}_M))\}$ as
Therefore, the event $\{\hat{M}=M, \hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M, \hat{\mathtt{s}}_{\hat{M}}^\rho=\mathtt{s}_M^\rho\}$ is equivalent to the affine constraint $A_\rho Z_\rho\leq b_\rho$. We proceed similarly to the subvector case considered previously to obtain the SI CI for $\rho_M(\theta_{M0})$. Let $r_\rho=r(Z_\rho, \hat{\Sigma}_\rho, e_{j(p+1)})$ as in (ref), and fix $\zeta\in(0,1)$. The SI CI of level $1-\zeta$ for $\rho_M(\theta_{M0})$ is given by $\mathrm{CI}_{\hat{M}}^\rho\equiv [\tilde{q}_l^\rho, \tilde{q}_u^\rho]$, where $\tilde{q}_l^\rho$ and $\tilde{q}_u^\rho$ are the solutions respectively to the following equations
We summarize the asymptotic validity of the above CI in the next corollary which follows from the arguments similar to the proof of Proposition (ref).
The debiased Lasso method of Zhang-Zhang(2014) and Javanmard-Montanari(2014) is based on the one-step estimator constructed from the initial Lasso estimator $\hat{\theta}$:
This particular variant of the debiased Lasso that employs the standard Hessian is proposed by Xia-Nan-Li(2021) for a homoskedastic GLM. Similarly, we use $\hat{I}(\hat{\theta})$ to estimate the asymptotic variance of $n^{1/2}S(\theta_0)$ and $I(\theta_0)$. To show the consistency of $\hat{H}(\hat{\theta})$ and $\hat{I}(\hat{\theta})$, we first extend Corollary 5.50 of Vershynin(2010) to random matrices i.n.i.d. rows with non-identical second moment matrices in the following lemma.
See Appendix (ref) for a proof. Relative to Theorem 5.39 and Corollary 5.50 of Vershynin(2010), the invertibility of $\bar{\Sigma}_n$ is required in Lemma (ref), but the rows of the matrix $A$ can be heterogeneous with non-identical second moment matrices $\Sigma_i, i=1,\dots, n$. Lemma (ref) together with Lemma S2 of Xia-Nan-Li(2021) yields the following result.
The proof is provided in Appendix (ref) which essentially verifies that the argument of Xia-Nan-Li(2021) goes through with i.n.i.d. data. For inference on a $r\times 1$ vector nonlinear parameter function $\rho(\theta)$ (which may depend on $n$), we define a debiased Lasso (one-step estimator) as
We establish the asymptotic validity of Wald-type inference based on the debiased Lasso estimator above in the proposition below.
The proof is given in Appendix (ref). $\lambda=C\sqrt{\frac{\log p}{n}}$ is a standard assumption in the literature Buhlmann-vandeGeer(2011), Negahban-etal(2012),vandeGeer-etal(2014), HTW(2015). The assumptions imposed on the number of covariates $p$, and the model sparsity $m_0$ are the same as those in Xia-Nan-Li(2021). In particular, while the condition $m_0 \log p \sqrt{\frac{p}{n}}\to 0$ is stronger than the condition $m_0 \frac{\log p}{\sqrt{n}}\to 0$ assumed by vandeGeer-etal(2014), no assumption is imposed directly on the sparsity of the inverse Hessian (and information matrix) i.e. $\max_{j}m_j=o(n/\log p)$, where $m_j\equiv |\{k\neq j: (H(\theta_0)^{-1})_{jk}\neq 0\}|$ is the number of non-zero elements of the $j$-th row of $H(\theta_0)^{-1}$, as in vandeGeer-etal(2014). As noted by Xia-Nan-Li(2021), the condition $p^2/n\to 0$ is weaker than the condition $\max_{j}m_j=o(n/\log p)$, when $m_j$ is of the order $p$.
The assumption of locally Lipschitz Jacobian $\dot{\rho}(\theta)$ is slightly stronger than the usual continuous differentiability assumption required for testing nonlinear hypotheses (see e.g. Section 9 of Newey-McFadden(1994) and Hansen(2022b), Hansen(2022)). Under this assumption, an error term $n^{1/2}(\dot{\rho}(\hat{\theta})-\dot{\rho}(\bar{\theta}))'(\hat{\theta}-\theta_0)$, where $\bar{\theta}$ is a mean-value between $\hat{\theta}$ and ${\theta}_0$, that results from the estimation of $\theta_0$ and $\rho(\theta_0)$ becomes negligible.
Using Proposition (ref), we obtain confidence intervals for the elements of $\theta_0$ as well as the vector nonlinear parameter function $\rho(\theta_0)$. One can also consider a plug-in estimator $\rho(\tilde{\theta})$, where $\tilde{\theta}$ is the one-step estimator defined in (ref). This estimator is asymptotically equivalent to the one-step estimator $\tilde{\rho}$ in (ref). The proof is actually similar to that of Proposition (ref), thus is omitted. In addition, multi-step estimators of $\theta_0$ and $\rho(\theta_0)$ can also be considered.
Belloni-Chernozhukov-Wei(2016) develop subvector inference procedure in a high-dimensional GLM that satisfies sparsity assumptions. They construct an estimating equation orthogonalized against the direction of the nuisance parameter estimation which also underlies the Neyman(1959)'s $C(\alpha)$ test. Here, we consider a survey version of the $C(\alpha)$-type statistic for the $r\times 1$ nonlinear parameter functon $\rho(\theta)$ defined as
where $\tilde{\theta}^{*}$ is an auxiliary estimate that satisfies $\rho(\tilde{\theta}^{*})=\rho_0$. This test statistic is proposed, in a regular likelihood context, by Smith(1987c) and studied further by Dufour-Trognon-Tuvaandorj(2016) among others.
The proof is given in Appendix (ref). In general, determining an auxiliary estimator that satisfies the constraint $\rho(\tilde{\theta}^{*})=\rho_0$ may be difficult. However, as we show in the next section, when testing a restriction on the AME of a binary regressor in the logit model, such an estimator can be readily obtained.
This section applies the results established in the previous sections to inference on the logit model estimated by the Lasso from survey data. The standard logit specification for a dependent variable $y_i, i=1,\dots,n,$ is
where $x_i=(1, \tilde{x}_i')'\in\mathbb{R}^{p+1}$, $\tilde{x}_i=(\tilde{x}_{i1},\dots, \tilde{x}_{ip})'\in\mathbb{R}^p$, and $\theta=(\alpha, \beta')'\in\mathbb{R}^{p+1}$, $\alpha\in\mathbb{R}$, $\beta\in\mathbb{R}^p$. Given the survey weights $\{w_i\}_{i=1}^n$ on the observations $\{(y_i, x_i')'\}_{i=1}^n,$ the weighted log-likelihood function is
The score function, the sample information and negative Hessian functions are given by
In the context of the logit model, a key parameter of interest is the AME which is a nonlinear function of the model parameters. As such, this section focuses on the inference on AMEs. The marginal effect (ME) of a binary regressor $\tilde{x}_{ij}, j=1,\dots, p, i=1,\dots, n,$ with a coefficient $\theta_{(j)}$ is calculated by the change in $P[y_i=1\vert x_{i}]$ when the regressor $\tilde{x}_{ij}$ is switched from 0 to 1 holding all other variables constant:
The AME of the $j$-th regressor is defined as
where $\theta_0$ denotes the true value of $\theta$ and the expectation is taken with respect to the distribution of the regressors. Let us first consider the debiased Lasso inference for the AMEs. A natural estimator of $\text{AME}_j(\theta_0)$ is
where $\hat{\theta}=(\hat{\alpha}, \hat{\beta}')'$ is an estimator of $\theta_0$ e.g. the survey-weighted Lasso estimator. In the current context, the one-step estimator defined in (ref) specializes to
where
To obtain a confidence interval for $\text{AME}_j(\theta_0), j=2,\dots, p+1,$ we can then use
Next, we turn to the SI. Let $\text{AME}_M(\theta_M)=[\text{AME}_{M2}(\theta_M),\dots, \text{AME}_{MM}(\theta_M)]'\in\mathbb{R}^{|M|-1}$ denote the AMEs for the active variables selected by the survey-weighted Lasso with coefficients $\beta_M$. Then, from (ref) the SI for the AMEs in the selected model is based on the one-step estimator
Finally, we consider the $C(\alpha)$ statistic. Let $C_\alpha(\mathrm{AME}_{0j})$ denote the $C(\alpha)$ statistic for testing $H_0:\mathrm{AME}_{j}=\mathrm{AME}_{0j}$. To obtain an auxiliary estimate that satisfies $\mathrm{AME}_{j}(\tilde{\theta}^{*})=\mathrm{AME}_{0j}$, we only need to solve for a scalar ${\theta}_{(j)}$ in the following equation:
Testing the zero restriction $H_0:\mathrm{AME}_{j}=0$ is particularly simple. First, note that $\mathrm{AME}_{j}=0$ if $\theta_{(j)}=0$. Furthermore, the Jacobian used in the $C_\alpha(\mathrm{AME}_{0j})$ statistic is
Let $\tilde{\theta}^{*}$ denote the estimator when the $j$-th element of the Lasso estimator $\hat{\theta}$ is replaced by $0$. Then, we have
where $C_\alpha(\theta_{0(j)})$ is the $C(\alpha)$ statistic for testing the coefficient $H_0:\theta_{(j)}=0$. We summarize this simple observation in the following lemma.
This section presents a simulation evidence on the performance of the proposed procedures. We consider a logit model where the regressors and the dependent variables are generated as follows:
where $\theta_0=(1,1,1, 0_{1\times (p-2)})'$, $\tilde{x}_{ij}\sim i.i.d.\, \mathrm{Bernoulli}(\mathrm{prob})$, $j=1,\dots, p,$ $i=1,\dots, N$, $x_i=(1, \tilde{x}_i')'$ and $\pi_i=x_i'\theta_0$. We set the size of the population equal to $N=10,000$. Two sampling schemes are considered: standard stratified sampling and exogenous stratification with $\mathrm{prob}=0.5$ and $\mathrm{prob}=0.4$, respectively. For each scheme, we create 4 strata and consider two cases: $(n_s, n)\in\{(50, 200), (100, 400)\}$, where $n_s$ observations are drawn from each stratum with replacement yielding a stratified sample of size $n$. In the standard stratified sampling with $\mathrm{prob}=0.5$, the population is stratified into 4 strata of sizes $N_1=1000, N_2=2000, N_3=3000$ and $N_4=4000$, respectively. As a result, the weights on the observations are $w_i=0.1, 0.2, 0.3, 0.4$ corresponding to the four strata. In the exogenous stratification with $\mathrm{prob}=0.4$, the population is stratified according to the values of the first two non-constant regressors: $(\tilde{x}_{i1}, \tilde{x}_{i2}) \in\{(0,0),(0,1), (1,0), (1,1)\}$. The weights on the observations in the above four strata are $w_i=0.36,0.24, 0.24, 0.16$, respectively. To assess the effect of the dimension of the regressors, the values for $p$ are set such that $\frac{p}{n}\in\{0.01, 0.025, 0.05, 0.1, 0.25, 0.5\}$ for each $n\in\{200,400\}$. The true value of AME corresponding to the coefficient $\theta_{(2)}=\beta_1$ is $0.11$. The empirical size of the tests is examined by testing the following two restrictions separately:
To test the hypothesis on AME, we implement the two SI approaches, labeled as SI and SI2, with or without conditioning on the sign of the estimated AME, respectively, described in Section (ref). For the auxiliary estimate $\tilde{\theta}^{*}$ in the $C(\alpha)$ statistic, we used the one-step iteration of $(1,\hat{\theta}_{(-2)}')'$, where $1$ corresponds to the tested value and $\hat{\theta}_{(-2)}$ is the (weighted) logistic Lasso estimate of $\theta_{(-2)}$, the model coefficients other than $\theta_{(2)}$. Moreover, whenever the sample Hessian evaluated at $\tilde{\theta}^{*}$ in the $C(\alpha)$ statistic is found to be singular, we used the Moore-Penrose inverse. There was no such issue in the other test statistics. The model (ref) is fit using the R package glmnet. For the tuning parameter $\lambda$, we use the default value of the package which is chosen by 10-fold cross validation with loss function “auc” (area under the ROC curve).
Tables (ref) and (ref) report the empirical sizes of the tests under standard stratified sampling and exogenous stratification, respectively. The results under both sampling schemes are qualitatively similar. All tests show reasonable size control when the number of regressors is moderate i.e. $p/n=0.01, 0.025, 0.05, 0.1, 0.25$ for both hypotheses. We can also see that the SI tests tend to underreject in most cases while the $C(\alpha)$ test does so when $p/n=0.5$. The size distortions of the SI method could potentially be alleviated by considering an appropriate form of bootstrap. When $p/n=0.5$, that is, the number of covariates is large relative to the sample size, all tests tend to underreject. This may be attributed to the conditions imposed on the growth rate of $p$ relative to the degree of sparsity, the tuning parameter and the sample size which are needed for the asymptotic validity of the DB and $C(\alpha)$ tests given in Propositions (ref) and (ref).
Moreover, when $p/n=0.5$, the $C(\alpha)$ test exhibits a substantial size distortion, while the DB and SI tests show somewhat better performance despite the fact that, in this case, the number of covariates are too high relative to the sample size for our results to hold. It is also clear that the rejection rate of the survey $t$-test, denoted as $t_{\mathrm{svy}}$, deteriorates as the ratio $p/n$ grows, which is expected as the test is not robust to increasing number of covariates.
This section applies the proposed methods to Canadian Internet Use Survey (CIUS) 2020 data, and examines what demographic factors affect a person's access to a government program or service.\footnote{Available at \url{https://www150.statcan.gc.ca/n1/daily-quotidien/210622/dq210622b-eng.htm}} The dependent variable is a binary variable where respondents answered 1) yes; 2) no; 3) not stated to the question “During the past 12 months, what activities did you perform on the Internet to interact with the government in Canada? Was it: Accessed an account for a government program or service?" The covariates in this analysis are income, education, employment status, aboriginal identity, visible minority status, immigration status, gender, type of household, language spoken at home, and province. All have two or more categories. There are $n = 17,031$ observations in the survey.
The collection of CIUS 2020 is based on a stratified design employing probability sampling; the stratification is done at the province/census metropolitan area (CMA) and census agglomeration (CA) level where each of the ten Canadian provinces were divided into strata/geographic areas.\footnote{There are 151 strata with the largest stratum, Toronto, having 2,235,145 private dwellings and the smallest stratum, Elliot Lake, having 6,259 private dwellings as of 2016, see \url{https://www12.statcan.gc.ca/census-recensement/2016/dp-pd/hlt-fst/pd-pl/Table.cfm?Lang=Eng&T=201&SR=1&S=3&O=D&RPP=9999&PR=0}} Each record on a sampling frame used in CIUS 2020 is a group of one or several telephone numbers associated with the same address from the Census and various administrative sources with Statistics Canada's dwelling frame. The records\textemdash the groups of telephone numbers\textemdash were sampled independently without replacement from each stratum.
The initial weight on observations is the inverse of an adjusted version of the probability of selection equal to the number of records sampled in the stratum divided by the number of records in the stratum from the survey frame. The final person weight $w_i$ is an adjusted version of the initial weight that takes into account the household size and survey-response among others.\footnote{Further details of the weighting procedure can be found in Section 10 of Microdata user Guide, CIUS 2020 at \url{https://www23.statcan.gc.ca/imdb/p2SV.pl?Function=getSurvey$\&$SDDS=4432$\#$a2}}
The base categories are omitted in each model as the comparison category for the logit model. The representative individual in the base category has the following characteristics \textendash\, male, non-aboriginal, neither English nor French (e.g. English and non-official language) speaker, not employed, some post-secondary education, not a visible minority, family household with children under 18, income of $\$44,120$-$\$75,321$, landed immigrant (recent immigrant), and from the province Alberta.
Table (ref) reports the inference results for the logit coefficients. The survey logit Lasso selects French, Employed, High school or less, University degree, Visible Minority, Family household with no children under 18, and Single. The magnitude of the Lasso estimates are in line with the survey logit estimates, and the signs of the estimates also appear reasonable. All inference methods indicate that the coefficients on Employed, University degree and Visible minority are highly significant. The variable French is selected by the Lasso but the inference results show that its coefficient is far from being significant. It is interesting to note that although New Brunswick (NB) is not selected by the Lasso, its debiased Lasso estimate $-0.32$ is almost identical to the survey GLM estimate $-0.33$ and highly significant. Table (ref) displays the inference results for the AMEs. The employed are about $8$-$11$ percentage points more likely to use the government online service than those not employed. Moreover, the use of government online service in NB appears to be 6-7 percentage points lower than the level of Alberta (AB).
The debiased Lasso and $C(\alpha)$ test results in Table (ref) show that the family household without children under age 18 is less likely to use the government service than those with children under age 18. Moreover, low educational attainment and high income negatively affect the likelihood of an individual using the government online services. The variables with the largest (in absolute value) AMEs on whether a person uses government online services are whether or not a person is employed, whether or not a person is single, and if their educational attainment was a High school or less or a University degree.
This paper has provided two main results. First, we have extended Lasso inference methods to a GLM with survey weights and/or heteroskedasticity, and established their asymptotic validity. Second, we have considered inference on nonlinear parameter functions. The proposed extended inference methods were applied to the logit model and remain reliable when $p/n$ increases as illustrated in a simulation study with standard stratified sampling and exogenous stratification. An empirical illustration based on the CIUS 2020 data also confirms the relevance of the proposed approach.