The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
75,436 characters
Penalized Likelihood Inference with Survey Data
\title{{Penalized Likelihood Inference with Survey Data}\thanks{
The authors acknowledge financial support of the Catalyzing Interdisciplinary Research Clusters Initiative, York University. This research is part of the project ``Digital Currencies'', approved and funded by the office of VPRI and participating faculties of York University. The authors acknowledge access to the data provided by RDC and
Statistics Canada (Project 21-MAPA YRK-721).}}
\date{\today}
\author{
Joann Jasiak\thanks{York University, \texttt{[email removed].}}
\and
Purevdorj Tuvaandorj\thanks{York University, \texttt{[email removed].}}
}
\maketitle
\begin{abstract}
This paper extends three Lasso inferential methods, Debiased Lasso, $C(\alpha)$ and Selective Inference to a survey environment.
We establish the asymptotic validity of the inference procedures in generalized linear models with survey weights and/or heteroskedasticity.
Moreover, we generalize the methods to inference on nonlinear parameter functions e.g. the average marginal effect in survey logit models.
We illustrate
the effectiveness of the approach in simulated data and Canadian Internet Use Survey 2020 data.
\begin{description}
\item[Keywords:]Survey data, Survey weights, Lasso, Logit, Average Marginal Effect, Post-Selection Inference
\end{description}
\end{abstract}
\section{Introduction}
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 \citep[see e.g.][]{Cameron-Trivedi(2009), Wooldridge(2010),Fuller(2011),Thompson(2012)}.\par
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.\par
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, \cite{Zhang-Zhang(2014)} and \cite{Javanmard-Montanari(2014)} propose a
debiased Lasso (DB) method which is based on one-step iteration of the initial Lasso estimator. \cite{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.\par
\cite{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 \cite{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.\par
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 \cite{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 \cite{Taylor-Tibshirani(2018)}
for a homoskedastic GLM may lead to invalid inference.\par
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.\par
There exist very few studies on the application of Lasso methods to survey data. \cite{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. \cite{Wooldridge(2001),Wooldridge(2010)} and
\cite{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.\par
The paper is organized as follows. Section \ref{sec: model} lays out the model framework.
We propose extensions of the selective inference, debiased Lasso and $C(\alpha)$/orthogonalization methods in Section \ref{sec: infer}.
Section \ref{sec: logit} applies the proposed methods to inference on AMEs in a survey logit model.
Section \ref{sec: simul} provides a simulation evidence on the properties of the proposed methods and
Section \ref{sec: app} presents an empirical application to Canadian Internet Use Survey 2020 data. We conclude in Section \ref{sec: conc}.\par
\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
\begin{equation}
\Vert X\Vert_{\psi_2}\equiv \sup_{m\geq 1}m^{-1/2}(\operatorname{E}[\vert X\vert^m])^{1/m}.
\end{equation}
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
\begin{equation}
\Vert X\Vert_{\psi_1} \equiv \inf \{t>0: \operatorname{E}[\exp(|X|/t)]\leq 2\}.
\end{equation}
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$.\par
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,
\begin{equation*}
F(x; \mu, \sigma^2, a,b)\equiv \frac{\Phi((x-\mu)/\sigma)-\Phi((a-\mu)/\sigma)}
{\Phi((b-\mu)/\sigma)-\Phi((a-\mu)/\sigma)},
\end{equation*}
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.\par
We abbreviate central limit theorem and continuous mapping theorem as CLT and CMT, respectively.
\section{Model}\label{sec: model}
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
\begin{equation*}
f(y_i\vert x_i,\theta_0)=\exp(y_ix_i'\theta_0-a(x_i'\theta_0))c(y_i),\quad i=1,\dots, n,
\end{equation*}
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:
\begin{equation}\label{eq: WLL}
L(\theta)\equiv -n^{-1}\sum_{i=1}^nw_i g(y_{i}, x_{i}'\theta).
\end{equation}
As is well known, the weighted likelihood framework is commonly used in survey data analysis \citep{Manski-Lerman(1977),Cameron-Trivedi(2009),Wooldridge(2010)}, and accommodates, among others, the following stratification schemes.
\begin{example}[Standard stratified sampling]\label{ex: SS}
Let $\mathcal{Z}$ be the population for $z=(y,{x}')'$ which is assumed to be infinite (or contain a large number of units). $\mathcal{Z}$ is stratified into $J$, nonempty, mutually exclusive and exhaustive strata such that $\mathcal{Z}=\bigcup_{j=1}^J\mathcal{Z}_j$. $n_j$ observations $\{z_{ij}\}_{i=1}^{n_j}=\{(y_{ij}, x_{ij}')'\}_{i=1}^{n_j}$ are sampled randomly from each stratum $\mathcal{Z}_j, j=1,\dots, J$. The strata sample sizes, $n_j$s, are non-random, and the population frequencies
$q_j=P[z\in \mathcal{Z}_j]>0, j=1,\dots, J,$ are assumed to be known. The weights on the observations from the $j$-th stratum are given by
$w_{n_0+\dots +n_{j-1}+1}=\dots= w_{n_0+\dots +n_{j-1}+n_j}=q_j/(n_j/n), j=1,\dots, J,$ with $n_0=0$ and $n=\sum_{j=1}^Jn_j$.
Let us re-label the observations as $z_{ij}=z_{n_0+\dots +n_{j-1}+i}, i=1,\dots, n_j, j=1,\dots, J$.
The corresponding
likelihood function is then
\begin{equation}\label{eq: SS loglik}
L(\theta)=-\sum_{j=1}^Jq_j \left(n_j^{-1}\sum_{i=1}^{n_j}g(y_{ij}, x_{ij}'\theta)\right)= -n^{-1}\sum_{i=1}^nw_i g(y_{i}, x_{i}'\theta).
\end{equation}
\end{example}
\begin{example}[Exogenous stratification]\label{ex: ES}
Let $\mathcal{Z}=\mathcal{Y}\times \mathcal{X}$, where $\mathcal{Y}$ and $\mathcal{X}$ are the sample spaces for $y$ and $x$.
The population is stratified into $J$ strata according to a deterministic function of ${x}_i$:
$\mathcal{X}=\cup_{j=1}^J\mathcal{X}_j$, where $\mathcal{X}_j, j=1\dots, J,$ are mutually exclusive.
The population frequencies
$q_j=P[z\in \mathcal{Z}_j]=P[x\in\mathcal{X}_j]>0, j=1,\dots, J,$ are assumed to be known.
Given $n=\sum_{j=1}^Jn_j$ observations
$\{z_{ij}\}_{i=1,\dots, n_j,\, j=1,\dots, J}=\{(y_{ij}, x_{ij}')'\}_{i=1,\dots,n_j,\, j=1,\dots, J}$, where $\{z_{ij}\}_{i=1}^{n_j}=\{(y_{ij}, x_{ij}')'\}_{i=1}^{n_j}$ sampled randomly from each stratum $\mathcal{Z}_j, j=1,\dots, J$, the likelihood function can be formulated as in
\eqref{eq: SS loglik}.
\end{example}
\noindent \cite{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.\par
The score function, the sample information and negative Hessian matrices corresponding to \eqref{eq: WLL} are defined as
\begin{align}
S(\theta)&\equiv \frac{\partial L(\theta)}{\partial \theta}
=-n^{-1}\sum_{i=1}^nw_ix_{i}\dot{g}(y_i,x_i'\theta),\quad
\dot{g}(y, t) \equiv \frac{\partial g(y, t)}{\partial t},\label{eq: scoref}\\
\hat{I}(\theta)
&\equiv n^{-1}\sum_{i=1}^nw_i^2x_ix_i'\dot{g}(y_i,x_i'\theta)^2,\label{eq: info}\\
\hat{H}(\theta)
&\equiv -\frac{\partial^2 {L}(\theta)}{\partial \theta\partial \theta'}=n^{-1}\sum_{i=1}^nw_ix_ix_i'\ddot{g}(y_i,x_i'\theta),\quad \ddot{g}(y, t)\equiv \frac{\partial^2 g(y, t)}{\partial t^2}.\label{eq: hessian}
\end{align}
Moreover, we define $H(\theta_0)\equiv \operatorname{E}[\hat{H}(\theta_0)]$ and $I(\theta_0)\equiv \operatorname{E}[\hat{I}(\theta_0)]$.\par
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$.\par
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 \eqref{eq: WLL} is minimized subject to $\ell_1$ penalty on the slope parameters:
\begin{equation}\label{eq: logitlasso}
\min_{\theta=(\alpha, \beta')'\in\mathbb{R}^{p+1}} \left(-{L}(\theta)+\lambda \Vert\beta\Vert_1\right),
\end{equation}
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.\par
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')'$.\par
The Lasso solution in \eqref{eq: logitlasso}, 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{subsub: selective} 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{subsec: db Lasso}--\ref{subsec: Calpha} 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 \eqref{eq: scoref}-\eqref{eq: hessian} as follows:
\begin{align*}
S(\theta)&=[S_M(\theta)', S_{-M}(\theta)']',\ S_M(\theta)\in\mathbb{R}^{|M|},\ S_{-M}(\theta)\in\mathbb{R}^{p+1-|M|},\\
\hat{H}(\theta)
&=
\begin{bmatrix}
\hat{H}_{M}(\theta_M)& \hat{H}_{M(-M)}(\theta_M)\\
\hat{H}_{-MM}(\theta_M)& \hat{H}_{-M}(\theta_M)
\end{bmatrix},\\
\hat{I}(\theta)
&=
\begin{bmatrix}
\hat{I}_{M}(\theta_M)& \hat{I}_{M(-M)}(\theta_M)\\
\hat{I}_{-MM}(\theta_M)& \hat{I}_{-M}(\theta_M)
\end{bmatrix},
\end{align*}
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.\par
With the partitioning above,
\cite{Lee-Sun-Sun-Taylor(2016)} and \cite{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 \eqref{eq: logitlasso} such that
\begin{align}
\frac{\partial L_M(\hat{\theta}_M)}{\partial \theta_M}-(0, \lambda\, \mathtt{s}_M')'
&=S_M(\hat{\theta}_M)-(0, \lambda\, \mathtt{s}_M')'=0,\quad \mathtt{s}_M=\mathrm{sign}(\hat{\beta}_M)\in\{-1,1\}^{\vert M\vert-1},\label{eq:KKT1}\\
\frac{\partial L_M(\hat{\theta}_M)}{\partial \theta_{-M}}-\lambda \mathtt{u}
&=S_{-M}(\hat{\theta}_M)-\lambda \texttt{u}=0,\quad \mathtt{u}\in\mathbb{R}^{p+1-|M|},\quad \Vert \mathtt{u}\Vert_{\infty}<1.\label{eq:KKT2}
\end{align}
\section{Post-Lasso selection inference}\label{sec: infer}
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 \cite{vandeGeer-etal(2014)} and
\cite{Xia-Nan-Li(2021)}.
\begin{assumption}[Asymptotic validity]\label{A: AsyValid}
\leavevmode
\begin{enumerate}[label={(\alph*)}]
\item \label{AsyValid max}
$\{(y_i, x_i')'\}_{i=1}^n$ are independent with $\max_{1\leq i\leq n}a_i<C_u<\infty$ a.s. where
\begin{equation*}
a_i\in\{\Vert x_i\Vert_{\psi_2}, \Vert x_{i}\Vert_{\infty}, \Vert X\theta_0\Vert_{\infty}\}.
\end{equation*}
Moreover, $w_i$ is non-random with $0<C_l<w_i<C_u$ for all $n, i$.
\item \label{AsyValid eval}
For $A\in \{H(\theta_0), I(\theta_0), \operatorname{E}[n^{-1}X'X]\}$, there exist positive constants $\lambda_l$ and $\lambda_u$ such that
$0<\lambda_l\leq \lambda_{\min}(A)\leq \lambda_{\max}(A)\leq \lambda_{u}<\infty$.
\item \label{AsyValid Lip}
The function $g(y, t)\equiv a(t)-yt-\log c(y)$ is convex in $t\in\mathbb{R}$ for all $y$, and
twice differentiable with
$\dot{g}(y, t) \equiv {\partial g(y, t)}/{\partial t}$ and $\ddot{g}(y, t)\equiv {\partial^2 g(y, t)} / {\partial t^2} $ for all $(y,t)$.
There exist a positive definite matrix $H$ and $\eta>0$ such that $\lambda_{\min}(H)>\lambda_l>0$ and
\begin{equation}\label{A: Strong Con}
n^{-1}\sum_{i=1}^n\operatorname{E}[w_i(g(y_i,x_i'\theta)-g(y_i,x_i'\theta_0))]\geq \Vert H^{1/2}(\theta-\theta_0)\Vert^2
\end{equation}
for all $\Vert X(\theta-\theta_0)\Vert_{\infty}<\eta$. Furthermore, $\ddot{g}(y,t)$ is Lipschitz with some constant $L_0>0$:
\begin{equation}\label{A: dg lip}
\max_{t_0 \in \{x_i' \theta_0\}} \sup_{\max(|t-t_0|, |\tilde{t} - t_0|) \leq \eta} \sup_{y \in \mathcal{Y}}\frac{| \ddot{g}(y,t)-\ddot{g}(y,\tilde{t})|}{|t-\tilde{t}|}\leq L_0,
\end{equation}
and
\begin{align}
&\max_{t_0 \in \{x_i' \theta_0\}}\sup_{y \in \mathcal{Y}} |\dot{g}(y,t_0)| \leq C_u,\label{A: dg1}\\
&\max_{t_0 \in \{x_i' \theta_0\}} \sup_{|t-t_0|\leq \eta} \sup_{y \in \mathcal{Y}} |\ddot{g}(y,t)| \leq C_u.\label{A: dg2}
\end{align}
\end{enumerate}
\end{assumption}
The boundedness of the variables stated in Assumption \ref{A: AsyValid}\ref{AsyValid max} is employed frequently in the literature, see \cite{Negahban-etal(2012)}, \cite{vandeGeer-etal(2014)} and
\cite{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.\par
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 \texttt{R} package \texttt{glmnet}, the weights are rescaled to sum to $n$.\par
Assumption \ref{A: AsyValid}\ref{AsyValid eval} 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{A: AsyValid}\ref{AsyValid Lip} 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 \cite{vandeGeer-etal(2014)} and \cite{Xia-Nan-Li(2021)}).\par
The condition \eqref{A: Strong Con} is essentially the the Quadratic Margin Condition needed for the consistency of the Lasso \citep{Buhlmann-vandeGeer(2011)} and see also \cite{Negahban-etal(2012)} for a related (stochastic) Restricted Strong Convexity condition. A sufficient condition for \eqref{A: Strong Con} is that $\ddot{g}(y, x'\theta)$ is bounded away from zero locally around $x_i'\theta_0$ for all $i=1,\dots,n$.
\subsection{Selective inference}\label{subsub: selective}
In this section, we extend the selective inference argument of \cite{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 \cite{Lee-Sun-Sun-Taylor(2016)} for further discussions about the difference between the SI and other inference methods.\par
As in \cite{Taylor-Tibshirani(2018)}, we fix $\lambda>0$ and consider the following one-step estimator
\begin{equation}\label{eq: onestep}
\tilde{\theta}_M\equiv \hat{\theta}_M+\hat{H}_M(\hat{\theta}_M)^{-1}S_M(\hat{\theta}_M),
\end{equation}
where
$S_M(\hat{\theta}_M)=(0, \lambda\, \mathtt{s}_M')'$, from which we obtain the one-step estimator of $\beta_{M0}$:
\begin{equation}\label{eq: one-step2}
\tilde{\beta}_M=\hat{\beta}_M+[0_{|M|-1}, I_{|M|-1}]\hat{H}_M(\hat{\theta}_M)^{-1}S_M(\hat{\theta}_M).
\end{equation}
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{\cite{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 \eqref{eq:KKT1} and \eqref{eq: one-step2}, it follows that
\begin{equation*}
\mathtt{s}_M=\mathrm{sign}\left(\tilde{\beta}_M-[0_{|M|-1}, I_{|M|-1}]\hat{H}_M(\hat{\theta}_M)^{-1}(0, \lambda\, \mathtt{s}_M')'\right),
\end{equation*}
hence
\begin{equation}\label{eq: activecon}
\mathrm{diag}(\mathtt{s}_M)\left(\tilde{\beta}_M-[0_{|M|-1}, I_{|M|-1}]\hat{H}_M(\hat{\theta}_M)^{-1}(0, \lambda\, \mathtt{s}_M')'\right)\geq 0.
\end{equation}
As argued by \cite{Taylor-Tibshirani(2018)} (see Equation (21) therein),
in a homoskedastic GLM, the random quantities appearing in the active and inactive constraints \eqref{eq:KKT1} and \eqref{eq:KKT2}
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 \eqref{eq:KKT2}. Let
\begin{equation}
\tilde{S}_{-M}(\hat{\theta}_M)\equiv S_{-M}(\hat{\theta}_{M})-\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}
S_{M}(\hat{\theta}_{M}).
\end{equation}
By the fact that $\Vert \texttt{u}\Vert_{\infty}<1$, and after some algebra, we can express the inactive constraints in \eqref{eq:KKT2} as follows:
\begin{align}
\tilde{S}_{-M}(\hat{\theta}_{M})
&\leq \lambda(\bm{1}_{p+1-|M|}-\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}(0, \mathtt{s}_M')'),\label{eq: inactive con1}\\
-\tilde{S}_{-M}(\hat{\theta}_{M})
&\leq \lambda(\bm{1}_{p+1-|M|}+\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}(0, \mathtt{s}_M)'),\label{eq: inactive con2}
\end{align}
where the inequalities hold element-wise. The Lasso selection events in \eqref{eq: activecon}, \eqref{eq: inactive con1} and \eqref{eq: inactive con2} can be
rewritten in a compact form as
\begin{equation}\label{eq: AZb}
\{AZ\leq b\},
\end{equation}
where
\begin{align}
A
&\equiv
\begin{bmatrix}
-\mathrm{diag}(\mathtt{s}_M)&0_{(|M|-1)\times (p+1-|M|)}\\
0_{(p+1-|M|)\times (|M|-1)}& I_{p+1-|M|}\\
0_{(p+1-|M|)\times (|M|-1)}& -I_{p+1-|M|}
\end{bmatrix}\in\mathbb{R}^{(2p+1-|M|)\times p},\quad
Z
\equiv n^{1/2}
\begin{bmatrix}
\tilde{\beta}_M\\
\tilde{S}_{-M}(\hat{\theta}_{M})
\end{bmatrix}\in\mathbb{R}^{p},\notag
\\
b
&\equiv n^{1/2}
\begin{bmatrix}
-\mathrm{diag}(\mathtt{s}_M)[0_{|M|-1}, I_{|M|-1}]\hat{H}_M(\hat{\theta}_M)^{-1}(0, \lambda\, \mathtt{s}_M')'\\
\lambda(\bm{1}_{p+1-|M|}-\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}(0, \mathtt{s}_M')')\\
\lambda(\bm{1}_{p+1-|M|}+\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}(0, \mathtt{s}_M')')
\end{bmatrix}\in\mathbb{R}^{2p+1-|M|}.\label{def: Ab}
\end{align}
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{subsec: proof SI}):
\begin{equation}\label{eq: AN}
\begin{bmatrix}
n^{1/2}(\tilde{\beta}_M-\beta_{M0})\\
n^{1/2}\tilde{S}_{-M}(\hat{\theta}_{M})
\end{bmatrix}
\displaystyle \stackrel{d}{\longrightarrow} {N}(0, \Sigma),
\end{equation}
where $\Sigma$ is a $p\times p$ asymptotic covariance matrix. The estimator of $\Sigma$ is
\begin{align}
&\hat{\Sigma}
=
\begin{bmatrix}
\hat{\Sigma}_{\beta\beta}&\hat{\Sigma}_{\beta s}\\
\hat{\Sigma}_{\beta s}'&\hat{\Sigma}_{ss}
\end{bmatrix},\label{def: Sigma}
\end{align}
where
\begin{align}
&\hat{\Sigma}_{\beta\beta}\equiv [0_{|M|-1}, I_{|M|-1}]\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_M(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}[0_{|M|-1}, I_{|M|-1}]',\notag\\
&\hat{\Sigma}_{\beta s}\equiv [0_{|M|-1}, I_{|M|-1}]\left[\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_{M(-M)}(\hat{\theta}_M)-\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_{M}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}\hat{H}_{M(-M)}(\hat{\theta}_M)\right],\notag\\
&\hat{\Sigma}_{ss}\equiv[I_{p+1-|M|}, -\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}]
\begin{bmatrix}
\hat{I}_{-M}(\hat{\theta}_M)& \hat{I}_{-MM}(\hat{\theta}_M)\\
\hat{I}_{M(-M)}(\hat{\theta}_M)& \hat{I}_{M}(\hat{\theta}_M)
\end{bmatrix}\notag\\
&\qquad\enspace\
[I_{p+1-|M|}, -\hat{H}_{-MM}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}]'.
\end{align}
From \eqref{eq: AN}, we have the distributional approximation for $Z$
\begin{equation}\label{eq: approx N}
Z\displaystyle \stackrel{a}{\sim} N(\mu, \Sigma),\quad \mu
\equiv n^{1/2}[\beta_{M0}', 0_{p+1-|M|}']'.
\end{equation}
The latter combined with the affine constraints $\{AZ\leq b\}$ in \eqref{eq: AZb}, is now amenable to application of Lemma \ref{lem: Poly} below, which summarizes two key results of \cite{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$:
\begin{align}
c=c(\Sigma, \eta)
&\equiv \Sigma \eta(\eta'\Sigma\eta)^{-1},\quad r=r(Z, \Sigma, \eta)\equiv (I_k-c\eta')Z,\label{eq:cr}\\
\mathcal{V}^{-}(r)
&\equiv \max_{j: (Ac)_j<0}\frac{b_j-(Ar)_j}{(Ac)_j},\label{eq:Vminus}\\
\mathcal{V}^{+}(r)
&\equiv \min_{j: (Ac)_j>0}\frac{b_j-(Ar)_j}{(Ac)_j}, \label{eq:Vplus}\\
\mathcal{V}^{0}(r)
&\equiv \min_{j: (Ac)_j=0}b_j-(Ar)_j\label{eq:Vzero},
\end{align}
where $(Ac)_j$ denotes the $j$-th element of $Ac$.
\cite{Lee-Sun-Sun-Taylor(2016)} show the following the result.
\begin{lemma}[Polyhedral lemma and truncated Gaussian pivot \citep{Lee-Sun-Sun-Taylor(2016)}]\label{lem: Poly}
Let $Z\sim {N}(\mu, \Sigma)$ and $A\in\mathbb{R}^{k\times k}$, $b\in\mathbb{R}^k$
and $\eta\in\mathbb{R}^k$ be fixed quantities. If $c, r, \mathcal{V}^{-}(r), \mathcal{V}^{+}(r)$ and $\mathcal{V}^{0}(r)$ are defined as in \eqref{eq:cr}-\eqref{eq:Vzero}, then
\begin{enumerate}[label={(\alph*)}]
\item $\eta'z$ is independent of $\mathcal{V}^{-}(r)$, $\mathcal{V}^{+}(r)$ and $\mathcal{V}^{0}(r)$, and
the following events are equivalent:
\begin{equation}\label{eq: PolyConst}
\{AZ\leq b\}
=\{\mathcal{V}^{-}(r)\leq \eta'Z\leq
\mathcal{V}^{+}(r), \mathcal{V}^{0}(r)\geq 0\}.
\end{equation}
\item\label{lem: TGP} Furthermore,
\begin{equation}\label{eq: TNP}
F(\eta'Z; \eta'\mu, \eta'\Sigma\eta, \mathcal{V}^{-}(r),\mathcal{V}^{+}(r))\vert \{AZ\leq b\}\sim U(0,1).
\end{equation}
\end{enumerate}
\end{lemma}
\noindent In our setup, $b$ defined in \eqref{def: Ab}
is random whereas Lemma \ref{lem: Poly} assumes constant $b$. In addition, we have an approximate normality in \eqref{eq: approx N} instead of the exact normality assumed in Lemma \ref{lem: Poly}.
These lead to an asymptotic version of \eqref{eq: TNP}, namely, as $n\to\infty$
\begin{equation}\label{TM2}
F(e_{jp}'Z, e_{jp}'\mu, e_{jp}'\hat{\Sigma} e_{jp}, \mathcal{V}^{-}(r),\mathcal{V}^{+}(r))\vert \{AZ\leq b\}\displaystyle \stackrel{d}{\longrightarrow} U(0,1).
\end{equation}
\eqref{TM2} can be established using the results of \cite{Markovic-etal(2017)}. Although
we do not directly use \eqref{TM2}, it provides the basis of the inference procedures described below.
\par
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 \eqref{def: Ab}, $\mu$ be as in \eqref{eq: approx N},
and set $\eta= e_{jp}\in\mathbb{R}^{p}$, $c=c(Z, \eta)$ and $r=r(Z, \hat{\Sigma}, e_{jp})$ in \eqref{eq:cr}, where $\hat{\Sigma}$ is defined in \eqref{def: Sigma}.
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
\begin{align}
&F(n^{1/2}q, e_{jp}'Z, e_{jp}'\hat{\Sigma} e_{jp}, \mathcal{V}^{-}(r),\mathcal{V}^{+}(r))=\frac{\zeta}{2},\label{eq: SI CI1}\\
&F(n^{1/2}q, e_{jp}'Z, e_{jp}'\hat{\Sigma} e_{jp}, \mathcal{V}^{-}(r),\mathcal{V}^{+}(r))=1-\frac{\zeta}{2}. \label{eq: SI CI2}
\end{align}
The asymptotic validity of the above CI is established in the following proposition.
\begin{proposition}\label{prop1}
Let Assumption \ref{A: AsyValid} hold with $p$ fixed, $\lambda=Cn^{-1/2}$, where $C=O(1)$, and
$H(\theta_0)$ and $I(\theta_0)$ converge to nonsingular matrices.
Then, it holds that for $\zeta\in(0,1)$
$$\displaystyle \operatornamewithlimits{\lim\inf \ }_{n\to\infty}P[e_{j\hat{M}}'\beta_{\hat{M}0}\in \mathrm{CI}_{\hat{M}j}\vert \hat{M}=M, \hat{\mathtt{s}}_M=\mathtt{s}_M]=1-\zeta.$$
\end{proposition}
See Appendix \ref{subsec: proof SI} for a proof. The assumption of fixed $p$ is commonly used in the literature on SI
\citep[see e.g.][]{Lee-Sun-Sun-Taylor(2016),Tian-Taylor(2017),Taylor-Tibshirani(2018),KKK(2022)}.
\cite{Taylor-Tibshirani(2018)} provide a heuristic argument for the validity of the selective inference in a homoskedastic GLM.
Proposition \ref{prop1} 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 \citep{Tian-Taylor(2017),KKK(2022)}.
We establish the asymptotic validity of the selective inference procedure
by verifying the conditions given in \cite{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 \eqref{eq: one-step2}, consider the one-step estimator
\begin{equation}\label{eq: si one-step}
\tilde{\rho}_M
\equiv \rho(\hat{\theta}_M)
+\dot{\rho}_M(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}S_M(\hat{\theta}_M),\quad \dot{\rho}_M({\theta}_M)\equiv \frac{\partial \rho_M({\theta}_M)'}{\partial \theta_M}.
\end{equation}
Standard arguments yield the distributional approximation
\begin{equation}\label{eq: one-stepb}
n^{1/2}\tilde{\rho}_M
\displaystyle \stackrel{a}{\sim} N\left(n^{1/2}\rho_M(\theta_{M0}), \dot{\rho}_M({\theta}_{M0})'\,\Sigma\, \dot{\rho}_M({\theta}_{M0})\right).
\end{equation}
Again, an approach similar to those applied to the elements of $\beta$ allows us to define the augmented variables:
\begin{align}
A_\rho
&\equiv
\begin{bmatrix}
0& 0_{p}'\\
0_{2p+1-|M|}& A
\end{bmatrix}\in\mathbb{R}^{(2p+2-|M|)\times (p+1)},\quad
Z_{\rho}
\equiv
\begin{bmatrix}
n^{1/2}\tilde{\rho}_M\\
Z
\end{bmatrix}\in\mathbb{R}^{p+1},\label{def: Ab rho1}\\
b_\rho
&\equiv
\begin{bmatrix}
0\\
b
\end{bmatrix}\in\mathbb{R}^{2p+2-|M|},\quad
\hat{\Sigma}_\rho
\equiv
\begin{bmatrix}\
\hat{\Sigma}_{\rho\rho}& \hat{\Sigma}_{\rho\beta}& \hat{\Sigma}_{\rho s}\\
\hat{\Sigma}_{\beta\rho}& \hat{\Sigma}_{\beta\beta}& \hat{\Sigma}_{\beta s}\\
\hat{\Sigma}_{s\rho}& \hat{\Sigma}_{s\beta}& \hat{\Sigma}_{ss}
\end{bmatrix}\in\mathbb{R}^{(p+1)\times (p+1)},\label{def: Ab rho2}
\end{align}
where $Z$, $A$ and $b$ are as defined in \eqref{def: Ab}, and
\begin{align*}
\hat{\Sigma}_{\rho\rho}
&\equiv \dot{\rho}_M(\hat{\theta}_{M})'\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_M(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}\dot{\rho}_M(\hat{\theta}_{M}),\\
\hat{\Sigma}_{\rho\beta}
&\equiv \dot{\rho}_M(\hat{\theta}_{M})'\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_M(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}[0_{|M|-1}, I_{|M|-1}]',\\
\hat{\Sigma}_{\rho s}
&\equiv \dot{\rho}_M(\hat{\theta}_{M})'\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_{M(-M)}(\hat{\theta}_M)-\dot{\rho}_M(\hat{\theta}_{M})'\hat{H}_M(\hat{\theta}_M)^{-1}\hat{I}_{M}(\hat{\theta}_M)\hat{H}_M(\hat{\theta}_M)^{-1}\hat{H}_{M(-M)}(\hat{\theta}_M).
\end{align*}
Then,
for $\zeta\in(0,1)$, the level $1-\zeta$ CI for $\rho_M(\theta_{M0})$ can be constructed as in \eqref{eq: SI CI1} and \eqref{eq: SI CI2}
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 \eqref{eq:cr}.\par
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 \eqref{eq: AZb}. 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
\begin{align*}
A_\rho
\equiv
\begin{bmatrix}
-\mathtt{s}_M^\rho& 0_{p}'\\
0_{2p+1-|M|}& A
\end{bmatrix}\in\mathbb{R}^{(2p+2-|M|)\times (p+1)},\,
b_\rho
\equiv
\begin{bmatrix}
-\mathtt{s}_M^\rho\,\dot{\rho}_M(\hat{\theta}_{M})'\hat{H}_M(\hat{\theta}_M)^{-1}(0, \lambda\, \mathtt{s}_M')'\\
b
\end{bmatrix}\in\mathbb{R}^{2p+2-|M|},
\end{align*}
and keep $Z_\rho$ and $\hat{\Sigma}_\rho$ as defined
in \eqref{def: Ab rho1} and \eqref{def: Ab rho2}. 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
\begin{align}
\{\mathtt{s}_M^\rho=\mathrm{sign}(\rho_M(\hat{\theta}_M))\}
&=\{\mathtt{s}_M^\rho\,\rho_M(\hat{\theta}_M)>0\}\notag\\
&=\{\mathtt{s}_M^\rho\,(\tilde{\rho}_M-\dot{\rho}_M(\hat{\theta}_M)'\hat{H}_M(\hat{\theta}_M)^{-1}(0, \lambda\mathtt{s}_M')')>0\}\notag\\
&=\{-\mathtt{s}_M^\rho\,\tilde{\rho}_M<-\lambda\,\mathtt{s}_M^\rho\,\dot{\rho}_M(\hat{\theta}_M)'\hat{H}_M(\hat{\theta}_M)^{-1}(0, \mathtt{s}_M')'\}.
\end{align}
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 \eqref{eq:cr}, 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
\begin{align}
&F(n^{1/2}q, n^{1/2}\tilde{\rho}_M, e_{j(p+1)}'\hat{\Sigma}_\rho e_{j(p+1)}, \mathcal{V}^{-}(r_\rho),\mathcal{V}^{+}(r_\rho))=\frac{\zeta}{2},\label{eq: SI CI3}\\
&F(n^{1/2}q, n^{1/2}\tilde{\rho}_M, e_{j(p+1)}'\hat{\Sigma}_\rho e_{j(p+1)}, \mathcal{V}^{-}(r_\rho),\mathcal{V}^{+}(r_\rho))=1-\frac{\zeta}{2}.\label{eq: SI CI4}
\end{align}
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{prop1}.
\begin{corollary}\label{corr1}
Suppose that the conditions of Proposition \ref{prop1} hold, and
the scalar nonlinear parameter function $\rho_M(\theta_M)$ is continuously differentiable in a neighborhood of $\theta_{M0}$ with
$\dot{\rho}_M({\theta}_{M0})'\dot{\rho}_M({\theta}_{M0})>\lambda_l>0,$
Then, it holds that for $\zeta\in(0,1)$
$$\displaystyle \operatornamewithlimits{\lim\inf \ }_{n\to\infty}P[\rho_{\hat{M}}\in \mathrm{CI}_{\hat{M}}^\rho \vert \hat{M}=M, \hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M, \hat{\mathtt{s}}_{\hat{M}}^\rho=\mathtt{s}_M^\rho]=1-\zeta.$$
\end{corollary}
\subsection{Debiased Lasso inference}\label{subsec: db Lasso}
The debiased Lasso method of \cite{Zhang-Zhang(2014)} and \cite{Javanmard-Montanari(2014)} is based
on the one-step estimator constructed from the initial Lasso estimator $\hat{\theta}$:
\begin{equation}\label{eq: one-step0}
\tilde{\theta}
=\hat{\theta}+\hat{H}(\hat{\theta})^{-1}S(\hat{\theta}).
\end{equation}
This particular variant of the debiased Lasso that employs the standard Hessian is
proposed by \cite{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 \cite{Vershynin(2010)} to random matrices i.n.i.d.
rows with non-identical second moment matrices in the following lemma.
\begin{lemma}[Covariance matrix consistency for i.n.i.d. random vectors.]\label{lem: covm}
Let $A$ be an $n\times p$ matrix whose rows $A_i'$ are
independent sub-Gaussian random vectors in $\mathbb{R}^p$ with $\operatorname{E}[A_i]=\mu_i$, $\operatorname{E}[A_iA_i']=\Sigma_i$ and $0<\lambda_{l}< \lambda_{\min}(\bar{\Sigma}_n)<\infty,$ where $\bar{\Sigma}_n\equiv n^{-1}\sum_{i=1}\Sigma_i$. Then for every $t\geq 0$, with
probability at least $1-2\exp(-t^2)$ it holds that
\begin{equation}\label{eq: HD LLN}
\Vert n^{-1}A'A-\bar{\Sigma}_n\Vert_2\leq C_K\max(\delta, \delta^2)\Vert\bar{\Sigma}_n\Vert_2,\quad \delta\equiv c\left(\sqrt{\frac{p}{n}}+\frac{t}{\sqrt{n}}\right),
\end{equation}
where $c$ is an absolute constant and $C_K>0$ is a constant that depend only on the sub-Gaussian norm $K=\max_{i}\Vert A_i\Vert_{\psi_2}<\infty$ of the rows and $\lambda_l$.
\end{lemma}
See Appendix \ref{proof covm} for a proof. Relative to Theorem 5.39 and Corollary 5.50 of \cite{Vershynin(2010)}, the invertibility of $\bar{\Sigma}_n$ is required in Lemma \ref{lem: covm}, but the rows of the matrix $A$ can be heterogeneous with non-identical second moment matrices $\Sigma_i, i=1,\dots, n$. Lemma \ref{lem: covm} together
with Lemma S2 of \cite{Xia-Nan-Li(2021)} yields the following result.
\begin{lemma}[The rate of convergence of the Hessian and information matrices]\label{lem: HI rate}
Under Assumption \ref{A: AsyValid},
\begin{align}
&\Vert\hat{H}(\hat{\theta})-{H}(\theta_0)\Vert_2=O_p\left(
\sqrt{\frac{p}{n}}+m_0\lambda
\right),\label{eq: Hhat consistency}\\
&\Vert\hat{H}(\hat{\theta})^{-1}-H(\theta_0)^{-1}\Vert_2
=O_p\left(
\sqrt{\frac{p}{n}}+m_0\lambda
\right),\label{eq: Hhat inverse consistency}
\\
&\Vert\hat{I}(\hat{\theta})-{I}(\theta_0)\Vert_2=O_p\left(
\sqrt{\frac{p}{n}}+m_0\lambda
\right),\label{eq: Ihat consistency}\\
&\Vert\hat{I}(\hat{\theta})^{-1}-{I}(\theta_0)^{-1}\Vert_2=O_p\left(
\sqrt{\frac{p}{n}}+m_0\lambda
\right).\label{eq: Ihat inverse consistency}
\end{align}
\end{lemma}
The proof is provided in Appendix \ref{subsec: proof HI rate} which
essentially verifies that the argument of \cite{Xia-Nan-Li(2021)} goes through with i.n.i.d. data.\par
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
\begin{equation}\label{eq: db one-step}
\tilde{\rho}
\equiv \rho(\hat{\theta})
+\dot{\rho}(\hat{\theta})'\hat{H}(\hat{\theta})^{-1}S(\hat{\theta}),\quad \dot{\rho}({\theta})\equiv \frac{\partial \rho({\theta})'}{\partial \theta}.
\end{equation}
We establish the asymptotic validity of Wald-type inference based on the debiased Lasso estimator
above in the proposition below.
\begin{proposition}[Asymptotic validity of Survey Debiased Lasso test]\label{prop: DB}
Let Assumption \ref{A: AsyValid} hold and
assume that $\lambda=C\sqrt{\frac{\log p}{n}}$ with $C=O(1)$ and $p\geq 1$, $p^2/n\to 0$ and $m_0 \log p \sqrt{\frac{p}{n}}\to 0$ as $n\to \infty$.
If the $r\times 1$ function $\rho(\theta)$ is differentiable in
a neighborhood of $\theta_0$ with a locally Lipschitz Jacobian $\dot{\rho}(\theta)$ and
$\lambda_{\min}\left(\dot{\rho}({\theta}_0)'\dot{\rho}({\theta}_0)\right)>\lambda_l>0,$ where $r(<(p+1))$ is fixed, then
\begin{equation}\label{eq: db psi}
\left(\dot{\rho}(\hat{\theta})'
\hat{H}(\hat{\theta})^{-1}\hat{I}(\hat{\theta}) \hat{H}(\hat{\theta})^{-1}\dot{\rho}(\hat{\theta})\right)^{-1/2}n^{1/2}(\tilde{\rho}-\rho({\theta}_0))
\displaystyle \stackrel{d}{\longrightarrow} N(0, I_r).
\end{equation}
\end{proposition}
The proof is given in Appendix \ref{subsec: proof DB}. $\lambda=C\sqrt{\frac{\log p}{n}}$ is a standard assumption
in the literature \citep[see e.g.][]{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 \cite{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 \cite{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 \cite{vandeGeer-etal(2014)}. As noted by \cite{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$.\par
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 \cite{Newey-McFadden(1994)} and \cite{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.\par
Using Proposition \ref{prop: DB}, 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 \eqref{eq: one-step0}. This estimator is asymptotically equivalent to
the one-step estimator $\tilde{\rho}$ in \eqref{eq: db one-step}. The proof is actually similar to that of Proposition \ref{prop: DB}, thus is omitted. In addition, multi-step estimators
of $\theta_0$ and $\rho(\theta_0)$ can also be considered.
\subsection{$C(\alpha)$/Orthogonalization inference}\label{subsec: Calpha}
\cite{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 \cite{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
\begin{equation}
C_\alpha(\rho_0)
\equiv n\, S(\tilde{\theta}^{*})'\hat{H}(\tilde{\theta}^{*})^{-1}
\dot{\rho}(\tilde{\theta}^{*})
\left(\dot{\rho}(\tilde{\theta}^{*})'
\hat{H}(\tilde{\theta}^{*})^{-1}\hat{I}(\tilde{\theta}^{*})\hat{H}(\tilde{\theta}^{*})^{-1}
\dot{\rho}(\tilde{\theta}^{*})\right)^{-1}
\dot{\rho}(\tilde{\theta}^{*})'
\hat{H}(\tilde{\theta}^{*})^{-1}
S(\tilde{\theta}^{*}),
\end{equation}
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 \cite{Smith(1987c)} and
studied further by \cite{Dufour-Trognon-Tuvaandorj(2016)} among others.
\begin{proposition}[Asymptotic validity of Survey $C(\alpha)$ test]\label{prop: Ca}
Let Assumption \ref{A: AsyValid} hold and
assume that $\lambda=C \sqrt{\frac{\log p}{n}}$ with $C=O(1)$, $p^2/n\to 0$ and $m_0 \log p \sqrt{\frac{p}{n}}\to 0$ as $n\to \infty$.
Let $\tilde{\theta}^{*}$ be an auxiliary estimator that satisfies
$\Vert\tilde{\theta}^{*}-\theta_0\Vert^2=O_p(m_0\lambda^2)$ and $\rho(\tilde{\theta}^{*})=\rho_0$, where
the nonlinear parameter function $\rho(\theta)\in\mathbb{R}^r$ satisfies the conditions given in
Proposition \ref{prop: DB}.
Then, under $H_0:\rho(\theta_0)=\rho_0$
\begin{equation}
C_\alpha(\rho_0)\displaystyle \stackrel{d}{\longrightarrow} \chi^2_{r}.
\end{equation}
\end{proposition}
\noindent The proof is given in Appendix \ref{subsec: Proof Calpha}. 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.
\section{Survey logit}\label{sec: logit}
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
\begin{equation}
P[y_i=1\vert x_i]
=\Lambda(x_i'\theta),
\end{equation}
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
\begin{equation}\label{eq: WLL logit}
L(\theta) =n^{-1}\sum_{i=1}^{n} w_i(y_{i} x_i'\theta - \log(1 + \exp(x_{i}'\theta))).
\end{equation}
The score function, the sample information and negative Hessian functions are given by
\begin{align}
S(\theta)
&=\frac{\partial {L}(\theta)}{\partial \theta}=n^{-1}\sum_{i=1}^nw_ix_i(y_i-\Lambda(x_i'\theta)),\label{eq: score logit}\\
\hat{I}(\theta)
&=n^{-1}\sum_{i=1}^nw_i^2x_ix_i'\Lambda(x_i'\theta)(1-\Lambda(x_i'\theta)),\label{eq: info logit}\\
\hat{H}(\theta)
&=-\frac{\partial^2 {L}(\theta)}{\partial \theta\partial \theta'}=n^{-1}\sum_{i=1}^nw_i x_ix_i'\Lambda(x_i'\theta)(1-\Lambda(x_i'\theta)).\label{eq: hessian logit}
\end{align}
\subsection{Inference on average marginal effects}
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:
\begin{align}
\text{ME}_{ij}(\theta)\equiv\Lambda({x}_{i}'\theta)\vert_{\tilde{x}_{ij}=1}-\Lambda({x}_{i}'\theta)\vert_{\tilde{x}_{ij}=0}.
\end{align}
The AME of the $j$-th regressor is defined as
\begin{equation*}
\text{AME}_j=\text{AME}_j(\theta_0)
\equiv \operatorname{E}\left[\frac{1}{\sum_{i=1}^nw_i}\sum_{i=1}^nw_i\text{ME}_{ij}(\theta_0)\right],
\end{equation*}
where $\theta_0$ denotes the true value of $\theta$ and the expectation is taken with respect to the distribution of the regressors.\par
Let us first consider the debiased Lasso inference for the AMEs.
A natural estimator of $\text{AME}_j(\theta_0)$ is
\begin{equation*}
\widehat{\text{AME}}_j(\hat{\theta})
\equiv \frac{1}{\sum_{i=1}^nw_i}
\sum_{i=1}^nw_i\left(\Lambda({x}_{i}'\hat{\theta})\vert_{\tilde{x}_{ij}=1}-\Lambda({x}_{i}'\hat{\theta})\vert_{\tilde{x}_{ij}=0}\right),
\end{equation*}
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 \eqref{eq: one-step0} specializes to
\begin{equation*}
\widetilde{\text{AME}}_j
=\widehat{\text{AME}}_j(\hat{\theta})
+\frac{\partial \widehat{\text{AME}}_j(\hat{\theta})}{\partial \theta'}\hat{H}(\hat{\theta})^{-1}S(\hat{\theta}),
\end{equation*}
where
\begin{equation*}
\frac{\partial \widehat{\text{AME}}_j(\hat{\theta})}
{\partial\theta}
\equiv \frac{1}{\sum_{i=1}^nw_i}
\sum_{i=1}^nw_i\left\{\left[x_i\Lambda({x}_{i}'\hat{\theta})(1-\Lambda({x}_{i}'\hat{\theta}))\right]\vert_{\tilde{x}_{ij}=1}-\left[x_i\Lambda(x_{i}'\hat{\theta})(1-\Lambda({x}_{i}'\hat{\theta}))\right]\vert_{\tilde{x}_{ij}=0}\right\}.
\end{equation*}
To obtain a confidence interval for $\text{AME}_j(\theta_0), j=2,\dots, p+1,$ we can then use
\begin{align*}
\left(\frac{\partial \widehat{\mathrm{AME}}_j(\hat{\theta})}{\partial\theta'}\hat{H}(\hat{\theta})^{-1}\hat{I}(\hat{\theta}) \hat{H}(\hat{\theta})^{-1}\frac{\partial \widehat{\mathrm{AME}}_j(\hat{\theta})}{\partial\theta}\right)^{-1/2}n^{1/2}(\widetilde{\mathrm{AME}}_j-\mathrm{AME}_j)
&\displaystyle \stackrel{d}{\longrightarrow} N(0, 1).
\end{align*}
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 \eqref{eq: si one-step} the SI for the AMEs in the selected model is based on the one-step estimator
\begin{equation}\label{eq: SI AME one-step}
\widetilde{\text{AME}}_M
=\widehat{\text{AME}}_M(\hat{\theta}_M)
+\frac{\partial \widehat{\text{AME}}_M(\hat{\theta}_M)}{\partial \theta_M'}\hat{H}_M(\hat{\theta}_M)^{-1}S_M(\hat{\theta}_M).
\end{equation}
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:
\begin{equation*}
\frac{1}{\sum_{i=1}^nw_i}
\sum_{i=1}^nw_i\left(\Lambda\left({\theta}_{(j)}+{x}_{i(-j)}'\hat{\theta}_{(-j)}\right)
-\Lambda\left({x}_{i(-j)}'\hat{\theta}_{(-j)}\right)\right)=\mathrm{AME}_{0j}.
\end{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
\begin{align*}
\frac{\partial \widehat{\text{AME}}_j(\theta)}
{\partial\theta}
&\equiv \frac{1}{\sum_{i=1}^nw_i}
\sum_{i=1}^nw_i\left[0,\dots,\left(\Lambda({x}_{i}'\theta)(1-\Lambda({x}_{i}'\theta))\right)\vert_{\tilde{x}_{ij}=1}, \dots, 0\right]'\\
&=\frac{1}{\sum_{i=1}^nw_i}
\sum_{i=1}^nw_i\Lambda({x}_{i}'\theta)(1-\Lambda({x}_{i}'\theta))\vert_{\tilde{x}_{ij}=1}e_{j(p+1)}.
\end{align*}
Let $\tilde{\theta}^{*}$ denote the estimator when the $j$-th element of the Lasso estimator $\hat{\theta}$ is replaced by $0$.
Then, we have
\begin{align*}
&C_\alpha(\mathrm{AME}_{0j})\\
&= n\, S(\tilde{\theta}^{*})'\hat{H}(\tilde{\theta}^{*})^{-1}
e_{j(p+1)}
\left(e_{j(p+1)}'
\hat{H}(\tilde{\theta}^{*})^{-1}\hat{I}(\tilde{\theta}^{*})\hat{H}(\tilde{\theta}^{*})^{-1}
e_{j(p+1)}\right)^{-1}
e_{j(p+1)}'
\hat{H}(\tilde{\theta}^{*})^{-1}
S(\tilde{\theta}^{*})\\
&=C_\alpha(\theta_{0(j)}),
\end{align*}
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.
\begin{lemma}\label{lem: Ca equiv.}
The $C(\alpha)$ statistic for testing the coefficient $H_0:\theta_{(j)}=0$ on a binary regressor $\tilde{x}_{ij}$ based on the auxiliary estimator $\tilde{\theta}^{*}$ is equivalent
to the $C(\alpha)$ statistic for testing the corresponding AME $H_0:\mathrm{AME}_{j}=0$ based on $\tilde{\theta}^{*}$.
\end{lemma}
\section{Simulations}\label{sec: simul}
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:
\begin{equation}
{y}_i\sim \mathrm{Bernoulli}(\pi_i),
\end{equation}
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$.\par
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.\par
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:
\begin{equation}\label{eq: H0 sim}
H_0: \theta_{(2)}=1,\quad H_0:\mathrm{AME}_2=0.11.
\end{equation}
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{subsub: selective}. 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.\par
The model \eqref{eq: logitlasso} is fit using the \texttt{R} package \texttt{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
``\texttt{auc}'' (area under the ROC curve).\par
Tables \ref{tab:Level0} and \ref{tab:LevelES} 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{prop: DB} and \ref{prop: Ca}.\par
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.\par
\begin{table}[htbp]
\caption{Empirical rejection frequencies of the tests for $H_0: \theta_{(2)}=1$ and $H_0:\mathrm{AME}_2=0.11$ at $5\%$ level. Standard stratified sampling.}
\label{tab:Level0}
\begin{center}
\begin{tabular}{lcccccc}
\toprule
\multicolumn{1}{l}{Tests} & \multicolumn{1}{c}{$p=2$} &\multicolumn{1}{c}{$p=5$} & \multicolumn{1}{c}{$p=10$} & \multicolumn{1}{c}{$p=20$} & \multicolumn{1}{c}{$p=50$}& \multicolumn{1}{c}{$p=100$}\\
\midrule
\multicolumn{7}{c}{$H_0:\theta_{(2)}=1$, $n=200$}\\
\midrule
DB & 5.0 & 4.4 & 3.7 & 3.1 & 4.5 & 3.3 \\
$C(\alpha)$ & 5.5 & 4.1 & 3.1 & 2.8 & 4.2 & 0.3 \\
SI & 3.9 & 2.5 & 2.3 & 2.6 & 3.6 & 3.6 \\
$t_{\mathrm{svy}}$ & 6.2 & 6.4 & 8.0 & 8.7 & 36.0 & 94.9 \\
\midrule
\multicolumn{7}{c}{$H_0:\mathrm{AME}_{2}=0.11$, $n=200$}\\
\midrule
DB & 5.4 & 5.3 & 4.6 & 3.7 & 3.5 & 1.4 \\
$C(\alpha)$ & 6.1 & 6.4 & 4.9 & 5.4 & 4.2 & 1.3 \\
SI & 4.2 & 2.6 & 2.2 & 2.8 & 3.6 & 4.7 \\
SI2 & 4.2 & 2.6 & 2.4 & 2.8 & 3.5 & 4.3 \\
$t_{\mathrm{svy}}$ & 5.7 & 7.7 & 7.4 & 8.2 & 50.9 & 93.3 \\
\bottomrule
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{lcccccc}
\toprule
\multicolumn{1}{l}{Tests} & \multicolumn{1}{c}{$p=4$} & \multicolumn{1}{c}{$p=10$} & \multicolumn{1}{c}{$p=20$} & \multicolumn{1}{c}{$p=40$} & \multicolumn{1}{c}{$p=100$} & \multicolumn{1}{c}{$p=200$}\\
\midrule
\multicolumn{7}{c}{$H_0:\theta_{(2)}=1$, $n=400$}\\
\midrule
DB & 4.8 & 4.4 & 6.0 & 3.7 & 5.6 & 3.9 \\
$C(\alpha)$& 4.7 & 4.2 & 4.2 & 4.2 & 5.7 & 0.5 \\
SI & 5.5 & 3.6 & 3.7 & 3.0 & 3.9 & 2.9 \\
$t_{\mathrm{svy}}$ & 5.0 & 5.1 & 6.3 & 15.9 & 40.4 & 98.3 \\
\midrule
\multicolumn{7}{c}{$H_0:\mathrm{AME}_{2}=0.11$, $n=400$}\\
\midrule
DB & 4.5 & 4.9 & 5.8 & 5.0 & 4.6 & 3.3 \\
$C(\alpha)$ & 5.2 & 6.9 & 8.5 & 3.1 & 3.4 & 1.0 \\
SI & 5.1 & 4.4 & 4.6 & 2.8 & 3.9 & 3.2 \\
SI2 & 5.1 & 4.4 & 4.6 & 2.9 & 3.7 & 2.8 \\
$t_{\mathrm{svy}}$ & 5.3 & 6.9 & 9.1 & 10.8 & 46.8 & 93.7 \\
\bottomrule
\end{tabular}
\end{center}
\footnotesize{Notes: $n=200, 400$ and 1000 simulation replications. DB, $C(\alpha)$, SI and $t_{\mathrm{svy}}$ denote
the debiased Lasso, $C(\alpha)$, selective inference and standard survey-weighted $t$ tests respectively. For the restriction $H_0:\mathrm{AME}_2=0.11$, SI is conditional on the sign of
the estimated AME in addition to $\hat{M}=M, \hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M$ while SI2 is conditional on the latter only.}
\end{table}
\begin{table}[htbp]
\caption{Empirical rejection frequencies of the tests for $H_0: \theta_{(2)}=1$ and $H_0:\mathrm{AME}_2=0.11$ at $5\%$ level. Exogenous stratification.}
\label{tab:LevelES}
\begin{center}
\begin{tabular}{lcccccc}
\toprule
\multicolumn{1}{l}{Tests} & \multicolumn{1}{c}{$p=2$} &\multicolumn{1}{c}{$p=5$} & \multicolumn{1}{c}{$p=10$} & \multicolumn{1}{c}{$p=20$} & \multicolumn{1}{c}{$p=50$}& \multicolumn{1}{c}{$p=100$}\\
\midrule
\multicolumn{7}{c}{$H_0:\theta_{(2)}=1$, $n=200$}\\
\midrule
DB & 4.9 & 4.8 & 3.1 & 4.1 & 7.3 & 4.6 \\
$C(\alpha)$ & 6.4 & 5.3 & 4.0 & 4.0 & 8.4 & 1.8 \\
SI & 4.4 & 2.1 & 2.7 & 2.2 & 2.9 & 4.1 \\
$t_{\mathrm{svy}}$ & 5.1 & 5.1 & 6.6 & 6.1 & 31.8 & 95.1 \\
\midrule
\multicolumn{7}{c}{$H_0:\mathrm{AME}_{2}=0.11$, $n=200$}\\
\midrule
DB & 5.4 & 4.9 & 5.3 & 3.8 & 5.6 & 4.3 \\
$C(\alpha)$ & 6.3 & 5.5 & 3.9 & 4.5 & 6.3 & 1.7 \\
SI & 4.1 & 1.9 & 3.0 & 2.4 & 2.8 & 5.2 \\
SI2 & 4.1 & 2.0 & 3.1 & 2.6 & 2.7 & 5.0 \\
$t_{\mathrm{svy}}$& 5.9 & 6.1 & 8.8 & 7.6 & 43.4 & 93.6 \\
\bottomrule
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{lcccccc}
\toprule
\multicolumn{1}{l}{Tests} & \multicolumn{1}{c}{$p=4$} & \multicolumn{1}{c}{$p=10$} & \multicolumn{1}{c}{$p=20$} & \multicolumn{1}{c}{$p=40$} & \multicolumn{1}{c}{$p=100$} & \multicolumn{1}{c}{$p=200$}\\
\midrule
\multicolumn{7}{c}{$H_0:\theta_{(2)}=1$, $n=400$}\\
\midrule
DB & 5.7 & 4.9 & 8.4 & 4.7 & 7.1 & 4.0 \\
$C(\alpha)$ & 5.1 & 4.6 & 5.7 & 4.7 & 9.1 & 7.0 \\
SI & 7.0 & 5.7 & 4.1 & 4.2 & 3.1 & 2.9 \\
$t_{\mathrm{svy}}$ & 4.8 & 5.2 & 6.0 & 9.2 & 29.9 & 98.4 \\
\midrule
\multicolumn{7}{c}{$H_0:\mathrm{AME}_{2}=0.11$, $n=400$}\\
\midrule
DB & 4.7 & 5.7 & 4.3 & 5.1 & 4.8 & 4.6 \\
$C(\alpha)$ & 4.6 & 5.3 & 5.6 & 4.5 & 4.3 & 0.3 \\
SI & 3.8 & 4.2 & 3.3 & 3.0 & 3.3 & 5.1 \\
SI2 & 3.8 & 4.2 & 3.3 & 3.0 & 3.2 & 4.8 \\
$t_{\mathrm{svy}}$ & 6.0 & 6.9 & 5.3 & 8.9 & 43.9 & 91.0 \\
\bottomrule
\end{tabular}
\end{center}
\footnotesize{Notes: $n=200, 400$ and 1000 simulation replications. DB, $C(\alpha)$, SI and $t_{\mathrm{svy}}$ denote
the debiased Lasso, $C(\alpha)$, selective inference and standard survey-weighted $t$ tests respectively. For the restriction $H_0:\mathrm{AME}_2=0.11$, SI is conditional on the sign of
the estimated AME in addition to $\hat{M}=M, \hat{\mathtt{s}}_{\hat{M}}=\mathtt{s}_M$ while SI2 is conditional on the latter only.}
\end{table}
\section{Empirical application}\label{sec: app}
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) \emph{yes}; 2) \emph{no}; 3) \emph{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?"\par
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.\par
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. \par
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\, \emph{male, non-aboriginal, neither English nor French} (e.g. English and non-official language) speaker, \emph{not employed, some post-secondary education, not a visible minority, family household with children under 18,} income of $\$44,120$-$\$75,321$, \emph{landed immigrant} (recent immigrant), and from the province Alberta.\par
Table \ref{tab: coefs} reports the inference results for the logit coefficients.
The survey logit Lasso selects \emph{French, Employed, High school or less, University degree, Visible Minority, Family household with no children under 18}, and \emph{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 \emph{Employed, University degree} and \emph{Visible minority} are highly significant. The variable \emph{French} is selected by the Lasso but the inference results show that its coefficient is far from being significant.\par
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. \par
Table \ref{tab: AME} 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).
\par
The debiased Lasso and $C(\alpha)$ test results in Table \ref{tab: AME} 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 \emph{High school or less} or a \emph{University degree}.\par
\begin{table}[htbp]
\small
\caption{Point estimates of $\theta_0$ and test p-values for $H_0:\theta_{0(j)}=0$ in Lasso Logistic Regression for Government Online Service Access.}
\begin{center}
\label{tab: coefs}
\begin{tabular}{lLLLLLLLL}
\toprule
\multicolumn{1}{l}{} & \multicolumn{4}{c}{Estimator} & \multicolumn{4}{c}{p-values}\\
\cmidrule(lr){2-5}\cmidrule(lr){6-9}
\multicolumn{1}{l}{Variable} & \multicolumn{1}{l}{GLM} & \multicolumn{1}{l}{Lasso} & \multicolumn{1}{l}{DB} &
\multicolumn{1}{l}{SI} & \multicolumn{1}{l}{$t_{\mathrm{svy}}$} & \multicolumn{1}{l}{DB} & \multicolumn{1}{l}{$C(\alpha)$} &
\multicolumn{1}{l}{SI} \\
\midrule
Intercept & -0.53 & -0.51 & -0.52 & -0.37 & 0.23 & 0.18 & 0.18 & 0.00 \\
Female & -0.02 & 0.00 & -0.02 & - & 0.61 & 0.61 & 0.61 & - \\
Aboriginal & 0.11 & 0.00 & 0.10 & - & 0.43 & 0.45 & 0.45 & - \\
Aboriginal n.s. & 0.85 & 0.00 & 0.71 & - & 0.11 & 0.16 & 0.16 & - \\
English & 0.30 & 0.00 & 0.28 & - & 0.49 & 0.44 & 0.44 & - \\
French & -0.18 & -0.22 & -0.15 & -0.59 & 0.68 & 0.70 & 0.78 & 1.00 \\
English and French & 0.48 & 0.00 & 0.46 & - & 0.27 & 0.22 & 0.22 & - \\
Language n.s. & -0.15 & 0.00 & -0.14 & - & 0.80 & 0.79 & 0.79 & - \\
Employed & 0.36 & 0.30 & 0.36 & 0.34 & 0.00 & 0.00 & 0.00 & 0.00 \\
Employment n.s. & 0.35 & 0.00 & 0.32 & - & 0.29 & 0.33 & 0.33 & - \\
High school or less & -0.53 & -0.41 & -0.51 & -0.51 & 0.00 & 0.00 & 0.00 & 1.00 \\
University degree & 0.37 & 0.32 & 0.37 & 0.35 & 0.00 & 0.00 & 0.00 & 0.00 \\
Education n.s. & -0.99 & 0.00 & -0.78 & - & 0.01 & 0.06 & 0.06 & - \\
Visible minority & 0.27 & 0.17 & 0.27 & 0.25 & 0.00 & 0.00 & 0.00 & 0.00 \\
Visible minority n.s. & 0.40 & 0.00 & 0.33 & - & 0.33 & 0.41 & 0.41 & - \\
Family household w.o.c.u 18 & -0.29 & -0.08 & -0.28 & -0.28 & 0.00 & 0.00 & 0.00 & 1.00 \\
Single & -0.72 & -0.31 & -0.70 & -0.65 & 0.00 & 0.00 & 0.00 & 1.00 \\
Other household type & -0.07 & 0.00 & -0.07 & - & 0.62 & 0.64 & 0.64 & - \\
Family n.s. & -0.17 & 0.00 & -0.17 & - & 0.48 & 0.48 & 0.48 & - \\
\$44,119 and less & -0.01 & 0.00 & -0.01 & - & 0.92 & 0.94 & 0.94 & - \\
\$75,322--\$109,431 & -0.04 & 0.00 & -0.04 & - & 0.62 & 0.62 & 0.62 & - \\
\$109,432--\$162,799 & -0.02 & 0.00 & -0.02 & - & 0.80 & 0.80 & 0.80 & - \\
\$162,800 and higher & -0.22 & 0.00 & -0.21 & - & 0.01 & 0.01 & 0.01 & - \\
Non-landed immigrant & 0.01 & 0.00 & 0.01 & - & 0.86 & 0.87 & 0.87 & - \\
Immigration n.s. & 0.90 & 0.00 & 0.91 & - & 0.29 & 0.31 & 0.31 & - \\
NL & 0.14 & 0.00 & 0.14 & - & 0.20 & 0.20 & 0.20 & - \\
PEI & -0.03 & 0.00 & -0.03 & - & 0.78 & 0.78 & 0.78 & - \\
NS & -0.21 & 0.00 & -0.21 & - & 0.05 & 0.06 & 0.06 & - \\
NB & -0.33 & 0.00 & -0.32 & - & 0.01 & 0.00 & 0.00 & - \\
QC & -0.25 & 0.00 & -0.25 & - & 0.02 & 0.02 & 0.02 & - \\
ON & -0.16 & 0.00 & -0.16 & - & 0.06 & 0.06 & 0.06 & - \\
MB & -0.21 & 0.00 & -0.20 & - & 0.07 & 0.07 & 0.07 & - \\
SK & -0.01 & 0.00 & -0.01 & - & 0.90 & 0.90 & 0.90 & - \\
BC & 0.04 & 0.00 & 0.04 & - & 0.65 & 0.65 & 0.65 & - \\
\bottomrule
\end{tabular}
\end{center}
\footnotesize{Notes: $n = 17,031$. GLM, Lasso, DB and SI in the columns 2-5 denote the survey GLM, survey Lasso, debiased Lasso and SI one-step estimates of $\theta_{0(j)}$, respectively. The columns 6-9 report the p-values of the survey GLM, DB, $C(\alpha)$ and SI tests for $\theta_{0(j)}=0$, respectively. $``-"$ means ``not computed". n.s. and w.o.c.u. abbreviate ``not stated'' and ``without children under".}
\end{table}
\begin{table}[htbp]
\small
\caption{Point estimates of $\mathrm{AME}$ and test p-values for $H_0:\mathrm{AME}_{j}=0$ in Lasso Logistic Regression for Government Online Service Access.}
\begin{center}
\label{tab: AME}
\begin{tabular}{lLLLLLLL}
\toprule
\multicolumn{1}{l}{} & \multicolumn{3}{c}{Estimator} & \multicolumn{4}{c}{p-values}\\
\cmidrule(lr){2-4}\cmidrule(lr){5-8}
\multicolumn{1}{l}{Variable} & \multicolumn{1}{l}{GLM} & \multicolumn{1}{l}{DB} &
\multicolumn{1}{l}{SI} & \multicolumn{1}{l}{$t_{\mathrm{svy}}$} & \multicolumn{1}{l}{DB} & \multicolumn{1}{l}{$C(\alpha)$} &
\multicolumn{1}{l}{SI} \\
\midrule
Female & 0.02 & -0.01 & - & 0.43 & 0.63 & 0.61 & - \\
Aboriginal & 0.19 & 0.02 & - & 0.11 & 0.46 & 0.45 & - \\
Aboriginal n.s. & 0.01 & 0.16 & - & 0.65 & 0.16 & 0.16 & - \\
English & -0.22 & 0.06 & - & 0.01 & 0.46 & 0.44 & - \\
French & 0.08 & -0.03 & -0.13 & 0.29 & 0.70 & 0.78 & 1.00 \\
English and French & 0.08 & 0.10 & - & 0.00 & 0.23 & 0.22 & - \\
Language n.s. & 0.07 & -0.03 & - & 0.49 & 0.80 & 0.79 & - \\
Employed & 0.11 & 0.08 & 0.08 & 0.27 & 0.00 & 0.00 & 0.00 \\
Employment n.s. & -0.04 & 0.07 & - & 0.48 & 0.34 & 0.33 & - \\
High school or less & -0.06 & -0.11 & -0.11 & 0.00 & 0.00 & 0.00 & 1.00 \\
University degree & -0.01 & 0.08 & 0.08 & 0.61 & 0.00 & 0.00 & 0.00 \\
Education n.s. & -0.04 & -0.16 & - & 0.68 & 0.10 & 0.06 & - \\
Visible minority & -0.12 & 0.06 & 0.06 & 0.00 & 0.00 & 0.00 & 0.00 \\
Visible minority n.s. & 0.20 & 0.08 & - & 0.29 & 0.42 & 0.41 & - \\
Family household w.o.c.u 18 & 0.00 & -0.06 & -0.06 & 0.92 & 0.00 & 0.00 & 1.00 \\
Single & -0.01 & -0.15 & -0.14 & 0.62 & 0.00 & 0.00 & 1.00 \\
Other household type & 0.00 & -0.02 & - & 0.80 & 0.65 & 0.64 & - \\
Family n.s. & -0.05 & -0.04 & - & 0.01 & 0.50 & 0.48 & - \\
\$44,119 and less & -0.03 & 0.00 & - & 0.80 & 0.94 & 0.94 & - \\
\$75,322--\$109,431 & -0.05 & -0.01 & - & 0.07 & 0.63 & 0.62 & - \\
\$109,432--\$162,799 & -0.07 & 0.00 & - & 0.01 & 0.81 & 0.80 & - \\
\$162,800 and higher & 0.03 & -0.05 & - & 0.20 & 0.01 & 0.01 & - \\
Non-landed immigrant & 0.00 & 0.00 & - & 0.86 & 0.88 & 0.87 & - \\
Immigration n.s. & -0.05 & 0.21 & - & 0.05 & 0.31 & 0.31 & - \\
NL & -0.04 & 0.03 & - & 0.06 & 0.21 & 0.20 & - \\
PEI & -0.02 & -0.01 & - & 0.62 & 0.79 & 0.78 & - \\
NS & -0.01 & -0.05 & - & 0.78 & 0.07 & 0.06 & - \\
NB & -0.06 & -0.07 & - & 0.02 & 0.01 & 0.00 & - \\
QC & -0.16 & -0.05 & - & 0.00 & 0.03 & 0.02 & - \\
ON & 0.00 & -0.04 & - & 0.90 & 0.07 & 0.06 & - \\
MB & 0.08 & -0.04 & - & 0.00 & 0.08 & 0.07 & - \\
SK & 0.06 & 0.00 & - & 0.00 & 0.90 & 0.90 & - \\
BC & 0.09 & 0.01 & - & 0.33 & 0.66 & 0.65 & - \\
\bottomrule
\end{tabular}
\end{center}
\footnotesize{Notes: $n = 17,031$. GLM, DB and SI in the columns 2-4 denote the survey GLM, debiased Lasso and SI one-step estimates of $\mathrm{AME}_{j}$.
The columns 5-8 report the p-values of the survey GLM, DB, $C(\alpha)$ and SI tests for $\mathrm{AME}_{j}=0$, respectively.
The $C(\alpha)$ test p-values are identical to those reported in Table \ref{tab: coefs} (Lemma \ref{lem: Ca equiv.}).
The p-values of SI2 were identical to those of SI, thus not shown. $``-"$ means ``not computed".}
\end{table}
\section{Conclusion}\label{sec: conc}
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.
\newpage
\bibliographystyle{myagsm}
\bibliography{DigitalDivide}
\newpage