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.
101,149 characters
A NEYMAN-ORTHOGONALIZATION APPROACH TO THE INCIDENTAL PARAMETER PROBLEM
\def\spacingset#1{\small\normalsize} \spacingset{1}
\if11
{
\title{\vspace{-1cm}{\bf \large A NEYMAN-ORTHOGONALIZATION APPROACH TO THE INCIDENTAL PARAMETER PROBLEM}\thanks{We are grateful to the Editor and referees, as well as Dmitry Arkhangelsky, Jin Hahn, Bo Honor\'e, Roger Moon, Whitney Newey, Andres Santos, Vira Semenova, and Vasilis Syrgkanis for comments and discussion. \newline Funded by the European Union (ERC-NETWORK-101044319) and by the French Government and the French National Research Agency under the Investissements d'Avenir program (ANR-17-EURE-0010). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. }}
\author{ St\'ephane Bonhomme\thanks{[email removed]} \\ {\small Department of Economics, University of Chicago} \and Koen Jochmans\thanks{[email removed]} \\ {\small Toulouse School of Economics, Universit\'e Toulouse Capitole} \and Martin Weidner\thanks{[email removed]} \\{\small Department of Economics and Nuffield College, University of Oxford}}
\date{\small February 2026}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\bf {\normalsize }}
\end{center}
\medskip
} \fi
\vspace{-1.18cm}
\begin{abstract}
\noindent A popular approach to perform inference on a target parameter in the presence of nuisance parameters is to construct estimating equations that are orthogonal to the nuisance parameters, in the sense that their expected first derivative is zero. Such first-order orthogonalization allows the estimator of the nuisance parameters to converge at a slower-than-parametric rate. It may, however, not suffice when the nuisance parameters are very imprecisely estimated. Leading examples are models for panel and network data that feature fixed effects. In this paper, we show how, in the conditional-likelihood setting, estimating equations can be constructed that are orthogonal to any chosen order $q$, in that their leading $q$ expected derivatives are zero. This yields estimators of target parameters that are unaffected by the presence of nuisance parameters to order $q$. In an empirical illustration, we apply our method to a fixed-effect model of team production.
\end{abstract}
\noindent
{\bf JEL Classification:} C13, C23, C55.
\medskip
\noindent
{\bf Keywords:} Neyman-orthogonality, incidental parameter, higher-order bias correction, networks.
\newpage
\onehalfspacing
\setcounter{equation}{0}
\section{Introduction}
Inference in the presence of nuisance parameters has received substantial attention. One fruitful way to proceed is to work with estimating equations that are orthogonal with respect to the nuisance parameters in the sense of \cite{Neyman1959}. Such equations underlie much of the results in semiparametric estimation (\citealp{Newey1994}) and are at the heart of recent advances on doubly-robust estimation and high-dimensional inference as discussed in, for example, \cite*{ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018}. A key finding is that Neyman-orthogonality permits the construction of asymptotically unbiased estimators that converge at the usual $n^{-\nicefrac{1}{2}}$-rate provided the nuisance parameter has a convergence rate that is faster than $n^{-\nicefrac{1}{4}}$, where $n$ is the sample size.
However, the faster-than-$n^{-\nicefrac{1}{4}}$ requirement often fails in problems where the dimension of the nuisance parameter is large relative to the sample size. Examples are panel data models with fixed effects, which are widely used in linear and nonlinear difference-in-differences settings. There, we observe $N$ units over $T$ periods of time and the model includes both common parameters and unit-specific nuisance parameters {(such as heterogeneous intercepts or slopes)}. The latter are estimated at the rate $T^{-\nicefrac{1}{2}}$. For an estimator of the former based on Neyman-orthogonalization to be successful we would therefore need that $T^{-\nicefrac{1}{2}} = o((N T)^{-\nicefrac{1}{4}})$, which translates into the requirement that $ N= o(T)$. This is usually not a realistic condition in microeconometric applications. In fact, under this requirement the standard fixed-effect estimator would permit asymptotically-valid inference. Consequently, (first-order) Neyman-orthogonalization does not solve the incidental parameter problem in panel data.\footnote{The problem is reminiscent of the poor performance of double machine-learning techniques in some settings, as recently documented by \cite{WuthrichZhu2021} and \cite{AngristFrandsen2022}.
A related problem where the conventional approach was formally shown to fail
is a nonlinear version of the judge-leniency design, see \citet{hahn2021problems}.
}
The issue can be even more severe in high-dimensional regressions on network data. In such settings, the convergence rate of the estimator of the nuisance parameter depends on the connectivity structure of the network (\citealp{JochmansWeidner2019}). Examples include the estimation of teacher value-added (\citealp{JacksonRockoffStaiger2014}), of the contributions of worker and firm heterogeneity to the variance of log wages and other covariance components (\citealp{AbowdKramarzMargolis1999}, \citealp{KlineSaggioSoelvsten2020}), as well as of complementarity patterns in team production (\citealp{AhmadpoodJones2019}, \citealp{Bonhomme2021}). Fixed effects in network-formation models are also poorly estimated, especially in the prevalent case where the network is sparse (see, e.g., \citealp{graham2020sparse}).
Motivated by these concerns, we focus on a higher-order generalization of Neyman-orthogonality that was proposed by \cite{MackeySyrgkanisZadik2018}. An estimating equation is Neyman-orthogonal to order $q$ when all $q$ leading derivatives with respect to the nuisance parameter have zero expectation. When $q=1$, this means that the expected Jacobian is zero, which recovers the conventional notion of Neyman-orthogonality (to order one). Working with estimating equations that are Neyman-orthogonal to order $q$, when combined with sample splitting, allows one to construct asymptotically-linear estimators when nuisance parameters are estimated at a rate {no slower} than $n^{-\nicefrac{1}{2(q+1)}}$. As an example, in the panel data problem, this reduces the bias from $O(T^{-1})$ down to $O(T^{-q})$, yielding valid inference under the requirement that $N=o(T^{2q-1})$. Combining orthogonalization with sample splitting (or cross-fitting) is important to achieve such an improvement, because orthogonalized estimating equations, by themselves, do not, in general, deliver estimators with improved sampling properties.
Working in the conditional-likelihood setting, we show how to construct estimating equations that are orthogonal to any chosen order. These estimating equations can be understood to be generalizations of the projected score of \cite{SmallMcLeish1989} and \cite{WatermanLindsay1996}.
They can also be seen as higher-order influence functions, as in \citet{RobinsLiTchetgenTchetgenvanderVaart2008} and \citet{vanderVaart2014}. Our approach applies to general low-dimensional target parameters that satisfy a moment restriction. This includes functions of the nuisance parameters such as average elasticities or other average effects. The conditional-likelihood framework allows us to orthogonalize a given estimating equation without introducing additional nuisance parameters. As is well known, this is not essential to achieve orthogonality to order one. However, avoiding such additional nuisance parameters turns out to be very helpful in enabling the construction of higher-order orthogonalized estimating equations.
Our approach relates to techniques to correct for bias in panel data (\citealp{HahnNewey2004}, \citealp{DhaeneJochmans2015b,DhaeneJochmans2015a}) and to the literature on small measurement error (\citealp{Chesher1991}, \citealp{evdokimov2023simple}). However, in contrast with these approaches, we do not restrict the nuisance parameters beyond the fact that they can be estimated at a certain rate.
We illustrate the usefulness of our approach in several examples and in an empirical application to the estimation of nonlinear regressions on network data; a problem for which, at present, no alternative solutions exist.
In this setting, we estimate a constant elasticity of substitution (CES) production function from the scientific output of research collaborations. As in \cite{AhmadpoodJones2019}, the production function depends on researcher-specific fixed effects.
Estimates of the parameters can be used to quantify the degree of complementarity among researchers within teams, and to compute the impact of counterfactual re-allocations in the spirit of earlier work by \citet{graham2014complementarity}.
This problem is difficult because in the data that we use (taken from \citealp{DuctorFafchampsGoyalvanderLeij2014} and concerning publications in economics on EconLit), the number of collaborations per researcher is quite low. A conventional estimator is thus likely to suffer from substantial bias. Our procedure uncovers the presence of complementarity among authors in the production of research articles. In a counterfactual exercise we also find that randomly pairing researchers would lead to a decrease in the average quality of articles. Our findings are corroborated in a simulation experiment targeted to our empirical application.
\setcounter{equation}{0}
\section{Problem statement and motivation\label{sec_mod}}
\subsection{Setup}
\label{subsec:Setup}
Let $Z_i=(Y_i,X_i)$ be random vectors, for $i=1,..,N$. We consider a setting where the conditional density function of $Y_i$ at $y$ given $X_i=x$, say $\ell(y\,|\, x;\theta_0,\eta_{i0})$, is known up to the parameters $\theta_0$ and $\eta_{i0}$. Throughout, we will treat $\eta_{10},\ldots,\eta_{N0}$ as nuisance parameters, and leave the marginal density of the conditioning variable, $\ell_{X_i}(x)$, unrestricted. We are interested in estimating a parameter $\mu_0$ that is defined through the moment condition
\begin{align}
\mathbb{E}\left(
\sum_{i=1}^N u(Z_i;\theta_0,\eta_{i0} ,\mu_0)\right) = 0,
\label{MainMoment}
\end{align}
where the expectations are over $Z_i$ under $\ell(y\,|\, x;\theta_0,\eta_{i0}) \,\ell_{X_i}(x)$. While the function $u$ could additionally depend on $i$, for instance in settings where the dimension of $\eta_{i0}$ differs across $i$, we omit this dependence for conciseness. We assume that, for all $i=1,\ldots,N$, $Z_i$ contains $n_i$ individual observations, and denote the total number of observations as $n=\sum_{i=1}^N n_i$. For example, in a balanced panel data setting with $N$ units and $T$ time periods, $Z_i$ is the time series of unit $i$'s observations, $n_i=T$ for all $i$, and $n=NT$.
Our setup accommodates different types of target parameters. As an example, we can set $\mu_0 = \theta_0$. In this case, using $u(z;\theta,\eta_i)$ as a shorthand for $u(z;\theta,\eta_i,\theta)$, one possibility is to use the score,
$$
u(z;\theta,\eta_i)
=
\frac{\partial \log \ell(y\,|\, x;\theta,\eta_i)}{\partial\theta}.
$$
More generally, the moment condition \eqref{MainMoment} defines the target parameter
$$\mu_0=\mu(\theta_0,\eta_{10},\ldots,\eta_{N0},\ell_{X_1},\ldots, \ell_{X_N}),$$ which can be a function of the
parameters $\theta_0$ and $\eta_{i0}$ describing the conditional distribution of $Y_i$ given $X_i$,
of the marginal distribution of $X_i$, and (implicitly) of the sample size. For example, we may be interested in an average effect of the form
$$
\mu_0 = \sum_{i=1}^N\int m(x;\theta_0,\eta_{i0})\ell_{X_i}(x) \, dx,
$$
where {$m$ is a known function.}
To illustrate the setup we will refer to two leading examples.
\paragraph{Example: Neyman-Scott model.}
Our first example is the well-known \cite{NeymanScott1948} model. Here,
\begin{equation}Y_{ij}=\eta_{i0} + \varepsilon_{ij},\quad \varepsilon_{ij}\sim \mathrm{iid}~{\cal{N}}\left(0,\sigma_0^2\right),\quad i=1,\ldots,N,\quad j=1,\ldots,T,\label{eq_neyman_scott}
\end{equation}
and the goal is to estimate $\theta_0 = \sigma_0^2$ in the presence of the nuisance parameters $\eta_{10},\ldots,\eta_{N0}$. Define, for all $i=1,\ldots,N$,
\begin{equation}\label{eq_NS}
u(Y_i;\sigma^2,\eta_i)=-\frac{T}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_{j=1}^T(Y_{ij}-\eta_i)^2,
\end{equation}
where $Y_i=(Y_{i1},\ldots,Y_{iT})^\top$ has dimension $n_i=T$, and the total number of observations is $n=NT$. It is well-known that the maximum-likelihood estimator of $\sigma_0^2$ is on average too small, suffering from bias $- \sigma_0^2 / T$. While in this panel data problem first-order orthogonality does not reduce the order of this bias, we demonstrate below that second-order orthogonalization fully removes it.
\paragraph{Example: CES production function.} Consider an environment where we observe workers producing output in $n$ teams of size 2. Moreover, let $k(j,1)$ and $k(j,2)$ denote the workers in team $j$, and write ${\cal{K}}=\{(k(j,1),k(j,2))\,:\, j=1,\ldots,n\}$ for the set of workers in all teams; note that a given worker may be part of multiple teams.
Consider a model for team production where team output is a CES aggregate of worker inputs (as in \citealp{AhmadpoodJones2019}),
\begin{equation}Y_{j}= \left(\frac{\eta_{k(j,1)0}^{\gamma_0}+\eta_{k(j,2)0}^{\gamma_0}}{2}\right)^{\frac{1}{\gamma_0}}\varepsilon_{j}^{\sigma_0},\quad \log \varepsilon_{j}\,|\, {\cal{K}}\sim \mathrm{iid}~{\cal{N}}\left(0,1\right),\quad j=1,\ldots,n.\label{eq_CES_size2}\end{equation}
In this model, one may be interested in estimating the substitution parameter $\gamma_0$ or the log error variance $\sigma_0^2$, average elasticities, or effects of counterfactual re-allocations of workers to teams, for example.{\footnote{The model relies on the assumption that the network ${\cal{K}}$ of co-workers is exogenous, i.e., independent of the shocks $\varepsilon_{j}$. Relaxing exogeneity through a parametric model of team formation is conceptually feasible within our likelihood approach, but we do not consider this extension here.}
To analyze this example we consider $N\leq n$ subsets of teams $j$, of size $n_i$ each. Let $Y_i$ denote the vector of team outcomes in subset $i$, let ${\cal{K}}_i$ denote the set of indicators for workers belonging to those teams, and let $\eta_i$ be the collection of all fixed effects of those workers. Finally, let $\theta=(\gamma,\sigma^2)^\top$. The scores with respect to $\gamma$ and $\sigma^2$ take the form $u(Y_i,{\cal{K}}_i;\theta,\eta_i)$. In contrast to our previous example, the theoretical literature on network models such as (\ref{eq_CES_size2}) is scarce, and to our knowledge no approach has as yet been developed for achieving bias reduction in such a setting.
In this example, a worker's fixed effect may appear in the nuisance parameter $\eta_i$ across multiple observations. While such ``overlapping fixed effects'' typically complicate the analysis and correction of incidental parameter bias, our orthogonalization and estimation methods straightforwardly accommodate this structure.
\subsection{The role of first-order orthogonality and its limitations\label{subsec_22}}
In the remainder of this section, we motivate our approach in a setting where one wishes to estimate $\mu_0=\theta_0$ based on a random sample $Z_1,\ldots, Z_n$, taking $u$ to be a univariate function, and $\eta_0$ {and $\theta_0$} to be scalar. Hence $N=1$, and $n_1=n$ is the total number of observations.
If $\mathbb{E}(\sum_{j=1}^nu(Z_j;\theta_0,\eta_0))=0$, a conventional estimator of {$\theta_0$, say $\widehat{\theta}$}, would be the solution to
$$
\sum_{j=1}^n u(Z_j;\theta,\widehat{\eta}) = 0,
$$
where $\widehat{\eta}$ is a consistent estimator of $\eta_0$ obtained in a preliminary step. However, it is well known that such a ``plug-in'' estimator is sensitive to the quality of the preliminary estimator $\widehat{\eta}$ used.
Assuming sufficient regularity, a standard argument based on a linearization around $\theta_0$ yields
$$
\left( \mathbb{E}\left(\frac{\partial u(Z_j;\theta_0,\eta_0)}{\partial \theta}\right) + o_P(1) \right)
\,
(\widehat{\theta}-\theta_0)
=
\frac{1}{n}
\sum_{j=1}^n u(Z_j;\theta_0,\widehat{\eta}),
$$
so that the sampling properties of $\widehat{\theta}-\theta_0$ are dictated by the sampling properties of the estimating equation.
We have
\begin{equation}\label{eq_expan_1}
\begin{split}
\frac{1}{n}
\sum_{j=1}^n u(Z_j;\theta_0,\widehat{\eta})
& =
\underset{(A)}{\underbrace{\frac{1}{n} \sum_{j=1}^n u(Z_j;\theta_0,\eta_0)}}
\\
&+
\underset{(B)}{\underbrace{ \left(
\frac{1}{n} \sum_{j=1}^n
\frac{\partial u(Z_j;\theta_0,\eta_0)}{\partial \eta}
-
\mathbb{E}\left( \frac{\partial u(Z_j;\theta_0,\eta_0)}{\partial \eta} \right)
\right)
\,
\left(\widehat{\eta}-\eta_0 \right)}}
\\
& +
\underset{(C)}{\underbrace{\mathbb{E}\left( \frac{\partial u(Z_j;\theta_0,\eta_0)}{\partial \eta} \right)
\,
\left(\widehat{\eta}-\eta_0 \right)}}
\\
&
+
O_P(\lvert \widehat{\eta}-\eta_0 \rvert^{2}).
\end{split}
\end{equation}
The (A) term in (\ref{eq_expan_1}) is a zero-mean sample average to which a standard central-limit theorem can be applied. Hence, it is generally $O_P(n^{-\nicefrac{1}{2}})$. The next two terms in the expansion capture the first-order effect of estimation noise in $\widehat{\eta}$. The (B) term can generally be ensured to be $o_P(n^{-\nicefrac{1}{2}})$. A generic approach to achieve this is to compute $\widehat{\eta}$ from data that are independent of $Z_1,\ldots, Z_n$, for example using sample splitting. In the case of (\ref{eq_expan_1}), (B) is the product of a sample average of zero-mean random variables---which is $O_P(n^{-\nicefrac{1}{2}})$---and an $o_P(1)$ term--- as $\widehat{\eta}$ is consistent for $\eta_0$---and, therefore, (B) is $o_P(n^{-\nicefrac{1}{2}})$.
The (C) term, however, features a non-random Jacobian that, in general, is non-zero. Hence, (C) is $O_P(\lvert\widehat{\eta}-\eta_0 \rvert)$, and will only be asymptotically negligible when $\widehat{\eta}$ is superconsistent for $\eta_0$, which is not usually the case.
Suppose now that $u$ is first-order orthogonal, in the sense that
\begin{equation} \label{eq:neyman1}
\mathbb{E}\left( \frac{\partial u(Z_j;\theta_0,\eta_0)}{\partial \eta} \right)
=
0.
\end{equation}
Then the (C) term vanishes from (\ref{eq_expan_1}) and we obtain
\begin{equation}
\begin{split}
\frac{1}{n}
\sum_{j=1}^n u(Z_j;\theta_0,\widehat{\eta})
= \frac{1}{n} \sum_{j=1}^n u(Z_j;\theta_0,\eta_0)
+
O_P(\lvert\widehat{\eta}-\eta_0 \rvert^2)
+ o_P(n^{-\nicefrac{1}{2}}).\label{eq_equ_1}
\end{split}
\end{equation}
The requirement that $\widehat{\eta}-\eta_0 = o_P(n^{-\nicefrac{1}{4}})$ then guarantees that the impact of the estimation error in $\widehat{\eta}$ on $\widehat{\theta}$ is asymptotically negligible. While a given function $u$ does not, in general, satisfy \eqref{eq:neyman1}, \cite{Neyman1959} proposed a general method to transform it into one that does. The resulting function is said to be Neyman-orthogonal.
Condition \eqref{eq:neyman1} has a long history in semiparametric estimation problems (\citealp{Bickel1982}, \citealp{Schick1986}, \citealp{Newey1994}). More recently, it has proved to be a fundamental ingredient in the literature on high-dimensional inference (see \citealp*{ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018}).
There are, however, instances where it is ineffective. To illustrate this it suffices to consider the simple panel data setting from the \cite{NeymanScott1948} problem.
{\paragraph{Example: Neyman-Scott model (continued).} In this problem it is easy to verify that
$$
\mathbb{E}\left( \frac{\partial u(Y_i;\sigma_0^2,\eta_{i0})}{\partial \eta_i} \right)=-\frac{1}{\sigma_0^4}\sum_{j=1}^T\mathbb{E}\left(Y_{ij}-\eta_{i0}\right)=0,
$$
and so the score is already first-order Neyman-orthogonal with respect to the fixed effects. Nevertheless, given preliminary estimators $\widehat{\eta}_1,\ldots, \widehat{\eta}_N$, and letting $\nu_i = \widehat{\eta}_i - \eta_{i0}$, the estimator
\begin{equation}
\widehat{\sigma}^2 =
\frac{1}{NT} \sum_{i=1}^N \sum_{j=1}^T (Y_{ij} - \widehat{\eta}_i)^2,\label{eq_sig2}
\end{equation}
has expectation
$
\sigma^2_0 - \nicefrac{2}{N} \sum_{i=1}^N \mathbb{E}(\overline{\varepsilon}_i \, \nu_i)
+
\nicefrac{1}{N}\sum_{i=1}^N \mathbb{E}(\nu_i^2),
$
for $\overline{\varepsilon}_i = \nicefrac{1}{T} \sum_{j=1}^T \varepsilon_{ij}$.
Thus, when using sample splitting, the bias is $\nicefrac{1}{N}\sum_{i=1}^N \mathbb{E}(\nu_i^2)$, the mean squared error of the preliminary estimator. With cross-fitting this is, at best, $O(T^{-1})$. Hence,
$
\sqrt{NT} (\widehat{\sigma}^2 - \sigma^2_0)
$
will not have a correctly-centered limit distribution unless $\nicefrac{N}{T}\rightarrow 0$. However, under this condition, the joint maximum-likelihood estimator of $\sigma^2_0$ and the fixed effects, too, is asymptotically unbiased. Hence, having a score that is Neyman-orthogonal, even when combined with sample splitting, does not suffice to resolve the incidental parameter problem in panel data problems.
}
\subsection{Higher-order orthogonality\label{subsec_23}}
To see how Neyman-orthogonality to a higher order can be helpful we now consider a further expansion of (\ref{eq_expan_1}). Again assuming sufficient regularity, we have, for any integer $q\geq 1$,
\begin{equation*} \label{eq_equ_2}
\begin{split}
\frac{1}{n}
\sum_{j=1}^n u(Z_j;\theta_0,\widehat{\eta})
& = \underset{(A)}{\underbrace{\frac{1}{n} \sum_{j=1}^n u(Z_j;\theta_0,\eta_0)}}
\\
& +
\underset{(B)}{\underbrace{\sum_{p=1}^q
\frac{1}{p!}
\left(
\frac{1}{n} \sum_{j=1}^n
\frac{\partial^p u(Z_j;\theta_0,\eta_0)}{\partial \eta^p}
-
\mathbb{E}\left( \frac{\partial^p u(Z_j;\theta_0,\eta_0)}{\partial \eta^p} \right)
\right)
\left(\widehat{\eta}-\eta_0 \right)^p}}
\\
& +
\underset{(C)}{\underbrace{\sum_{p=1}^q
\frac{1}{p!}
\mathbb{E}\left( \frac{\partial^p u(Z_j;\theta_0,\eta_0)}{\partial \eta^p} \right)
\,
\left(\widehat{\eta}-\eta_0 \right)^p}}
\\
&
+
O_P(\lvert \widehat{\eta}-\eta_0 \rvert^{q+1}).
\end{split}
\end{equation*}
Here, the (A) term is the same as before. Also, with sample splitting we can again ensure that the (B) term will be asymptotically negligible. On the other hand, if the function $u$ satisfies the higher-order orthogonality condition
\begin{equation} \label{eq:neymanq}
\mathbb{E}\left( \frac{\partial^p u(Z_j;\theta_0,\eta_0)}{\partial \eta^p} \right)
=
0
,\qquad
1\leq p \leq q,
\end{equation}
the (C) term is equal to zero, and so
\begin{equation}
\begin{split}
\frac{1}{n}
\sum_{j=1}^n u(Z_j;\theta_0,\widehat{\eta})
= \frac{1}{n} \sum_{j=1}^n u(Z_j;\theta_0,\eta_0)
+
O_P(\lvert\widehat{\eta}-\eta_0 \rvert^{q+1})
+ o_P(n^{-\nicefrac{1}{2}}).
\end{split} \label{eq_equ_2}
\end{equation}
Comparing (\ref{eq_equ_2}) to (\ref{eq_equ_1}) we see that the impact of estimation noise in $\widehat{\eta}$ on our estimator of $\theta_0$ has been reduced further. Moreover, for the impact of estimation error to be negligible, we now only require that $\lvert \widehat{\eta}-\eta_0 \rvert^{q+1} = o_P(n^{-\nicefrac{1}{2}})$. It then follows from standard results that, as $n\rightarrow\infty$,
$$
\sqrt{n}(\widehat{\theta}-\theta_0) \overset{d}{\rightarrow} {\cal N}(0,\Sigma_\theta)
$$
for some $\Sigma_\theta$, provided that
$$
\widehat{\eta}-\eta_0 = o_P(n^{-\nicefrac{1}{2(q+1)}}).
$$
The notion of $q$th-order Neyman-orthogonality as in (\ref{eq:neymanq}) was introduced by \citet{MackeySyrgkanisZadik2018}. In the context of our conditional-likelihood setup, we will give a general procedure to construct higher-order Neyman-orthogonal functions below.
{\paragraph{Example: Neyman-Scott model (continued)} In the model of \cite{NeymanScott1948},
$$\mathbb{E}\left( \frac{\partial^2 u(Y_i;\sigma_0^2,\eta_{i0})}{\partial \eta_i^2} \right)=\frac{T}{\sigma_0^4}\neq 0.$$ It thus follows that $u$ is not orthogonal to second order. Below we will show that a second-order Neyman-orthogonal score equation exists. This estimating equation no longer depends on a preliminary estimator of $\eta_{i0}$}, and its solution is the usual degrees-of-freedom corrected estimator
\begin{equation}
\widehat{\sigma}^2
=
\frac{1}{N(T-1)} \sum_{i=1}^N \sum_{j=1}^T (Y_{ij}-\overline{Y}_i)^2,\label{eq_sig2_hat}
\end{equation}
where $\overline{Y}_i = \nicefrac{1}{T} \sum_{j=1}^T Y_{ij}$. The estimator $\widehat{\sigma}^2$ is well-known to be fixed-$T$ consistent.
\setcounter{equation}{0}
\section{Estimation based on orthogonalized functions\label{sec_est}}
We now present our estimation approach in the general case where the target parameter $\mu_0$ may be equal to $\theta_0$ or may be a different parameter such as an average effect, and there are multiple, vector-valued nuisance parameters $\eta_{i0}$. We start by formally defining higher-order Neyman-orthogonality in this general setup and describe estimation based on higher-order Neyman-orthogonal moment functions. In the next section, we will then show how to construct such functions.
\subsection{Definition of higher-order orthogonality}
\label{subsec:DefHigherOrder}
Let $d_\eta$ be the dimension of $\eta$ and write $\eta=(\eta_1,\ldots,\eta_{d_\eta})$. For any non-negative integer $p$ and a vector of integers $m = (m_1,\ldots, m_p)$ satisfying $1\leq m_s \leq d_\eta$ for all $1\leq s\leq p$, define
\begin{equation}
D^m_{\eta} =
\frac{\partial^p }{\partial \eta_{m_1}\cdots \partial \eta_{m_p}}. \label{DefGenPartialDerivative}
\end{equation}
For a given $p$, there are $
d_p = \binom{d_\eta+p-1}{p}
$
unique such partial derivatives. Let
$
\nabla^{(p)}_\eta
$
be the vector operator of dimension $d_p$ that collects all these unique partial derivatives of order $p$. Finally, let $\nabla^q_\eta$ be the vector operator of dimension $\sum_{p=1}^q d_p$ obtained on stacking $\nabla^{(p)}_\eta$ for $p=1,\ldots,q$. Explicitly, we have
$$
\nabla^{q}_\eta
= \left(\begin{array}{c}
\nabla^{(1)}_\eta \\
\nabla^{(2)}_\eta \\
\vdots \\
\nabla^{(q)}_\eta
\end{array} \right)
= \left(\begin{array}{c}
\left[D^m_{\eta} \, :\, m \in \{1,\ldots,d_\eta\} \right]
\\
\left[D^m_{\eta} \, :\, m \in \{1,\ldots,d_\eta\}^2, \,
m_1 \leq m_2\right]
\\
\vdots \\
\left[D^m_{\eta} \, :\, m \in \{1,\ldots,d_\eta\}^q, \, m_1 \leq m_2 \leq \ldots \leq m_q\right]
\end{array} \right) .
$$
Neyman-orthogonality to order $q$ can now be defined as follows
(\citealp{MackeySyrgkanisZadik2018}).
\begin{definition}\label{def_ortho}
If the function $u$ satisfies
\begin{align}
\mathbb{E}\left[ \nabla^q_\eta \, u(Z;\theta_0,\eta_0,\mu_0) \right] = 0,
\label{NeymanOrthDef}
\end{align}
for some integer $q$, then we say that $u$ is Neyman-orthogonal to order $q$.
\end{definition}
\noindent
In this definition, {\it all} possible partial derivatives of $u(Z;\theta,\eta,\mu)$ with respect to $\eta$ up to order
$q$ have mean zero.
\subsection{Estimation}
\label{subsect:Estimation}
Let $\mu_0$ satisfy (\ref{MainMoment}) for a (possibly vector-valued) function $u$. We assume that $u$ is Neyman-orthogonal to order $q$ with respect to $\eta_i$, in the sense of Definition \ref{def_ortho}. Suppose that we have access to preliminary estimators $\widehat{\eta}_1,\ldots, \widehat{\eta}_N$ of the nuisance parameters that are independent of the data $Z_1,\ldots,Z_N$. If $\eta_{i0}$ is defined as the solution to a moment condition involving the same data, estimation based on sample-splitting, combined with cross-fitting (see, e.g., \citealp{NeweyRobins2017}), can be applied. When the observations are independent this is conventional. For situations where the data are dependent, modified sample-splitting strategies are available (see, e.g., \citealp{semenova2023inference}).
We estimate $\mu_0$ by the GMM estimator
\begin{equation}
\widehat{\mu}
=
\underset{\mu}{\mbox{argmin}}\, \left\|\sum_{i=1}^Nu(Z_i;\widehat{\theta},\widehat{\eta}_i,\mu)\right\|_{W},\label{muEstimation}
\end{equation}
where $W$ is a chosen symmetric positive-definite matrix, $\lVert u \rVert_W = \sqrt{u^\top W \, u}$, and $\widehat{\theta}$ is an estimator of $\theta_0$.
The estimator $\widehat{\theta}$ will depend on the problem at hand. If $\theta_0$ is defined through a moment condition of the form $\sum_{i=1}^N\mathbb{E}(\widetilde{u}(Z_i;\theta_0,\eta_{i0})) = 0$, for a function $\widetilde{u}$ that is Neyman-orthogonal to order $q$, then our framework can be applied and we can use
\begin{equation}\widehat{\theta}=\underset{\theta}{\mbox{argmin}}\, \left\|\sum_{i=1}^N\widetilde{u}(Z_i;\theta,\widehat{\eta}_i)\right\|_{\widetilde{W}},\label{thetaEstimation}\end{equation}
where $\widetilde W$ is again a chosen weight matrix. In this case, we may equally combine \eqref{muEstimation} and \eqref{thetaEstimation} into a single GMM estimation procedure.
In Section \ref{sec_asympt}, we provide conditions under which this approach yields estimators that are $n^{-\nicefrac{1}{2}}$-consistent and asymptotically normal, where $n$ is the total number of observations. We will impose two key conditions. The first one is that, although their number may increase with the sample size, the dimension of each $\eta_i$ remains bounded as $n$ tends to infinity. This imposes a suitable sense of sparsity in the relationship between the nuisance parameters and the outcomes. This condition is trivially satisfied in the panel data and network problems with fixed effects that we consider. The second key condition we impose is that the convergence rates of the preliminary estimates $\widehat{\eta}_i$ be faster than $n^{-\nicefrac{1}{2(q+1)}}$. This ensures that, after having orthogonalized to order $q$, any remainder terms are asymptotically negligible.
Our approach requires choosing an orthogonality order $q$. In practice, we recommend reporting estimates and standard errors for various values $q=1,2,...,q_{\rm max}$, each based on a $q$-orthogonal function $u_q^*$. For a given $q$ that satisfies the conditions of our theory (see Theorem \ref{th:Asymptotic} in Section \ref{sec_asympt}), the estimators based on $u_q^*$ and $u_{q+1}^*$ are both root-$n$ consistent and asymptotically normal. Given this, we propose to report as a diagnostic the difference in orthogonal moment functions $u^*_{q}-u^*_{q+1}$, suitably normalized. A large value of this statistic indicates that the order $q$ of orthogonality is likely too small. We describe this approach in Appendix \ref{App_q}, while leaving the formal construction of a selection method for $q$ and its impact on inference on $\mu_0$ to future work.
\setcounter{equation}{0}
\section{Achieving higher-order Neyman-orthogonality\label{sec_Construct}}
{
\subsection{Main result}
Let $u$ be a moment function. We now show how to construct an orthogonalized counterpart of $u$, which we call $u_q^*$, that is Neyman-orthogonal to order $q$, where $q\geq 1$ is any arbitrary order.
Recall the definition of the vector operators $\nabla^{(p)}_\eta$
and $\nabla^{q}_\eta$ in Section~\ref{subsec:DefHigherOrder}. It is convenient to
introduce the \cite{Bhattacharyya1946} basis $v_1,v_2,\ldots$, where
$$
v_p(z;\theta,\eta) = \frac{\nabla_{\eta}^{(p)} \ell(y\,|\, x;\theta,\eta)}{\ell(y\,|\, x;\theta,\eta)}.
$$
\cite{SmallMcLeish1994} discuss several properties of this basis. One important property for our purposes is that
\begin{equation} \label{eq:zeromean}
\mathbb{E}_{\theta,\eta}(v_p(Z;\theta;\eta) \,|\, X=x)
=
\int v_p(z;\theta,\eta) \, \ell(y\,|\, x;\theta,\eta) \, dy = 0
\end{equation}
for any $p$, so all elements of the Bhattacharyya basis have (conditional) mean equal to zero. In (\ref{eq:zeromean}), and throughout this section, $\mathbb{E}_{\theta,\eta}(\cdot \,|\, X=x) $ denotes the conditional expectation under $\ell(y\,|\, x;\theta,\eta)$.
The low-order basis functions are familiar from likelihood theory. For example,
\begin{equation*}
\begin{split}
v_1(z;\theta,\eta)
& =
\ \
\frac{\partial \log \ell(y\,|\, x;\theta,\eta)}{\partial \eta},
\\
v_2(z;\theta,\eta)
& =
\frac{\partial \log \ell(y\,|\, x;\theta,\eta)}{\partial \eta}\frac{\partial \log \ell(y\,|\, x;\theta,\eta)}{\partial \eta^{\top}}
+
\frac{\partial^2 \log \ell(y\,|\, x;\theta,\eta)}{\partial \eta\partial \eta^{\top}}.
\end{split}
\end{equation*}
The fact that these functions have mean zero follows from the unbiasedness of the score and from the information equality, respectively.
Stacking the leading $q$ basis functions, we obtain
$$
w_q(z;\theta,\eta) =\frac{\nabla_{\eta}^q \ell(y\,|\, x;\theta,\eta)}{ \ell(y\,|\, x;\theta,\eta)}
=
\left( \begin{array}{c}
v_1(z;\theta,\eta) \\
\vdots \\
v_q(z;\theta,\eta)
\end{array}
\right).
$$
The vectors $w_q$ are mean-zero ``generalized score functions''. The vector space spanning the Bhattacharyya basis at order $q$ is the tangent set of order $q$; see, e.g., \citet{vanderVaart2014}. While it is possible to achieve higher-order Neyman-orthogonality using other bases of functions, the Bhattacharyya basis delivers simple expressions through the use of Bartlett identities.
Next, let us define the matrices
$$
\varSigma_{w_qw_q}(x;\theta,\eta)
=
\mathbb{E}_{\theta,\eta}(w_q(Z;\theta,\eta)\, w_q(Z;\theta,\eta)^\top\,|\, X=x) ,
$$
and
$$
\varSigma_{w_qu}(x;\theta,\eta,\mu)
=
\mathbb{E}_{\theta,\eta}(w_q(Z;\theta,\eta)\, u(Z;\theta,\eta,\mu)^\top\,|\, X=x),
$$
which are, respectively, the (conditional) covariance matrix of the first $q$ members of the Bhattacharrya basis, and the covariance matrix
of the same $q$ basis functions with the vector function $u$.
Finally, let
$$
b_q(x;\theta,\eta,\mu) =
\nabla_{\eta}^q \,
\mathbb{E}_{\theta,\eta}(u(Z;\theta,\eta,\mu)^\top \vert X=x)
.
$$
Note that $b_q$ is zero when $u$ is the score for $\theta$, i.e.,
$\frac{\partial \log \ell(y\,|\, x;\theta,\eta)}{\partial\theta}$. In general, however, $b_q$ will be non-zero.
Here we
assume that $u(z;\theta,\eta,\mu)$ and $\ell(y\,|\, x;\theta,\eta)$
are sufficiently often differentiable in~$\eta$, and that the expectations in the definitions of
$\varSigma_{w_qw_q}$,
$\varSigma_{w_qu}$,
and $b_q$ are well-defined.
The proof of the following result is in Appendix \ref{App_proofs}.
\begin{theorem}\label{theo_neyman}
Suppose that $\varSigma_{w_qw_q}(x;\theta,\eta)$ is invertible and
let
$$
A(x;\theta,\eta,\mu)
=
\varSigma_{w_qw_q}(x;\theta,\eta)^{-1}
\left(
\varSigma_{w_q u}(x;\theta,\eta,\mu)
-
b_q(x;\theta,\eta,\mu)
\right).
$$
Then the function
$$
u_q^*(z;\theta,\eta,\mu)
=
u(z;\theta,\eta,\mu)
-
A(x;\theta,\eta,\mu)^\top \,
w_q(z;\theta,\eta)
$$
satisfies
$\mathbb{E}_{\theta,\eta} ( \nabla^q_\eta \, u_q^*(Z;\theta,\eta,\mu) \, \vert \, X=x ) = 0$.
This implies that $u_q^*$
is Neyman-orthogonal to order $q$, as defined above.
\end{theorem}
}
\noindent
Theorem \ref{theo_neyman} generalizes the projected-score construction of \citet{SmallMcLeish1989} and \citet{WatermanLindsay1996}. To see this, consider the case where $\mu_0=\theta_0$, and $u$ is the score function for $\theta$. Then $b_q=0$ and Theorem \ref{theo_neyman} yields
$$
u_q^*(z;\theta,\eta,\mu)
=
u(z;\theta,\eta,\mu)
-\left(\varSigma_{w_qw_q}(x;\theta,\eta)^{-1}
\varSigma_{w_q u}(x;\theta,\eta)\right)^\top \,
w_q(z;\theta,\eta),
$$
which is the projected score of order $q$. The projected score was originally developed as a tool to achieve E-ancillarity (\citealp{SmallMcLeish1988}) and to approximate the conditional score for $\theta$, when the latter exists (\citealp{WatermanLindsay1996}). It generalizes \cite{Neyman1959} in that $u_q^*$ is the (population) residual of a least-squares regression of $u$ on $w_q$; thus,
$$
\mathbb{E}_{\theta,\eta}
(
w_q(Z;\theta,\eta)
\,
u_q^*(Z;\theta,\eta,\mu)^\top
\vert X=x)
=0.
$$
While the fact that $u_q^*$ is Neyman-orthogonal is noted by \citet{WatermanLindsay1996} (although a link with Neyman's work is not made), it is not exploited. Moreover, unlike the conditional score, the projected score still depends on $\eta$, and it will generally not have improved properties over the score itself. As we highlight here, it is the combination of higher-order versions of Neyman-orthogonality with sample splitting that allows one to improve over working with the original score.
Theorem \ref{theo_neyman} covers more general estimating equations as well as more general parameters of interest, such as average elasticities or counterfactual quantities. Incorporating $b_q(x;\theta,\eta,\mu)$ into $u_q^*(z;\theta,\eta,\mu)$ is a key innovation that enables this generality. Observe that, in this case, we have that
$$
\mathbb{E}_{\theta,\eta}
(
w_q(Z;\theta,\eta)
\,
u_q^*(Z;\theta,\eta,\mu)^\top
\vert X=x)
=
b_q(x;\theta,\eta,\mu),
$$
thereby revealing $u_q^*$ to be an influence function of order $q$ per Equation (1.9) in \cite{vanderVaart2014}.
We note that Theorem \ref{theo_neyman} requires the matrix $\varSigma_{w_qw_q}(x;\theta,\eta)$ to be invertible. In the standard case of first-order Neyman-orthogonality this corresponds to non-singularity of the information matrix of the nuisance parameters. For higher-order Neyman-orthogonality this requirement imposes further restrictions. For example, in the team production model (\ref{eq_CES_size2}) with two-worker teams, the $2\times 2$ matrix $
\varSigma_{w_1w_1}({\cal{K}},\theta, \eta)$ is not invertible, as only the sum $\eta_{k(i,1)}^{\gamma}+\eta_{k(i,2)}^{\gamma}$ can be identified. In our application in Section \ref{sec_appli}, we will tackle this issue by combining data on teams of size 2 with single-author articles, and working with subsets $i$ of three teams each.{\footnote{As another example, consider the standard binary-choice panel data model
$$
\mathbb{P}_{\theta,\eta_i}(Y_{ij} = 1 \,|\, X_{i1},...,X_{iT})
=
\Phi(\eta_i + X_{ij}^\top \theta)
,\quad i=1,\ldots,N,\quad j=1,\ldots,T,
$$
for (conditionally-independent) binary outcomes $Y_{ij}$ and covariates $X_{ij}$. In this model, the rank of $\varSigma_{w_qw_q}(x;\theta,\eta)$ is bounded by $2^T$, so $\varSigma_{w_qw_q}(x;\theta,\eta)$ is singular for all $q>2^T$.}
\subsection{Intuition and discussion}
\noindent
To gain intuition into the construction in Theorem \ref{theo_neyman} it is useful to again consider the case where $u$ is a univariate function, the nuisance parameter is a scalar, and one wishes to estimate $\mu_0=\theta_0$.
\paragraph{First-order orthogonality.}
To relate our approach to the literature consider first $q=1$. Let
\begin{equation}
u_1^*(z;\theta,\eta) = u(z;\theta,\eta) - a_1(x;\theta,\eta) \, v_1(z;\theta,\eta),\label{eq_ustar}
\end{equation}
for some function $a_1$. Note that, by virtue of \eqref{eq:zeromean}, the term involving $v_1$ does not introduce any bias. We have
$$
\frac{\partial u_1^*(z;\theta,\eta)}{\partial\eta} =
\frac{\partial u(z;\theta,\eta)}{\partial\eta}
-
\frac{\partial a_1(x;\theta,\eta)}{\partial\eta}
\,
v_1(z;\theta,\eta)
-
a_1(x;\theta,\eta)
\,
\frac{\partial v_1(z;\theta,\eta)}{\partial\eta}.
$$
Take conditional expectations and exploit \eqref{eq:zeromean} to see that
$$
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial u_1^*(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
=
0
$$
if and only if
$$
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial u(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
-
a_1(x;\theta,\eta) \
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_1(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right) = 0.
$$
This is achieved by setting
\begin{equation}
a_1(x;\theta,\eta)
=
\left(
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_1(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
\right)^{-1}
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial u(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right).\label{eq_a1}
\end{equation}
Iterating expectations shows that the resulting function $u_1^*$ is Neyman-orthogonal to order $q=1$. By the information matrix equality we have
$$\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_1(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)=-\mathbb{E}_{\theta,\eta}
\left(
\left.
v_1(Z;\theta,\eta)^2
\right\rvert X=x
\right)=-\varSigma_{w_1w_1}(x;\theta,\eta),$$
and
$$\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial u(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)=-\mathbb{E}_{\theta,\eta}
\left(
\left.
v_1(Z;\theta,\eta)u(Z;\theta,\eta)
\right\rvert X=x
\right)=-\varSigma_{w_1u}(x;\theta,\eta),
$$
leading to the representation of the function $u_1^*$ as in the theorem.
The above derivation of \eqref{eq_a1} is well-known. Furthermore, it does not hinge on the likelihood structure. Indeed, recent work exploiting orthogonality, such as that surveyed in \cite*{ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018}, does so in the context of moment conditions. In our setup, as in \citeauthor{Neyman1959}'s (\citeyear{Neyman1959}) original work, the likelihood structure implies that $a_1$ is known up to the model parameters $\theta$ and $\eta$ (conditional on the regressors). Outside of this framework, in contrast, $a_1$ needs to be treated as an additional nuisance parameter. This is possible because, as $u_1^*$ is linear in $a_1$, it is automatically first-order Neyman-orthogonal to it by virtue of \eqref{eq:zeromean}. This logic, however, does not extend to higher order, as the implied system of equations becomes inconsistent, so that no solution exists, as we will see next.
\paragraph{Higher-order orthogonality.}
Let $q=2$, and again
consider a linear transformation of $u$, now involving the leading two Bhattacharyya basis functions. This gives
\begin{equation}u_{2}^*(z;\theta,\eta)
=
u(z;\theta,\eta)
-
\left(
\begin{array}{c}
a_{21}(x;\theta,\eta)
\\
a_{22}(x;\theta,\eta)
\end{array}
\right)^\top
\,
\left(
\begin{array}{c}
v_1(z;\theta,\eta)
\\
v_2(z;\theta,\eta)
\end{array}
\right).\label{eq_u2_star}
\end{equation}
Taking first-derivatives with respect to the nuisance parameter, and proceeding as in the first-order case, gives
$$
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial u(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
=
\left(
\begin{array}{c}
a_{21}(x;\theta,\eta)
\\
a_{22}(x;\theta,\eta)
\end{array}
\right)^\top
\,
\left(
\begin{array}{c}
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_1(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
\\
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_2(Z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
\end{array}
\right).
$$
Solving this equation for $a_{21}$ for given $a_{22}$ yields
\begin{equation}\label{eq_a21}
a_{21}(x;\theta,\eta)
=
a_1(x;\theta,\eta) - c_1(x;\theta,\eta) \, a_{22}(x;\theta,\eta),
\end{equation}
where $a_1$ is given by (\ref{eq_a1}) and
$$
c_1(x;\theta,\eta)
=
\left(
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_1(z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right)
\right)^{-1}
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial v_2(z;\theta,\eta)}{\partial\eta}
\right\rvert X=x
\right).
$$
The coefficient $c_1$ has the same form as $a_1$, except that it features $v_2$ instead of $u$. Moreover,
plugging (\ref{eq_a21}) back into (\ref{eq_u2_star}) yields
\begin{equation*}
u_{2}^*(z;\theta,\eta)
=
u_{1}^*(z;\theta,\eta) - a_{22}(x;\theta,\eta) \, v_2^*(z;\theta,\eta),
\end{equation*}
where
$
v_2^*(z;\theta,\eta)
=
v_2(z;\theta,\eta)
-
c_1(x;\theta,\eta) \,
v_1(z;\theta,\eta).
$
Note that $v_2^*$ is Neyman-orthogonal to order 1, that is,
$$
\mathbb{E}_{\theta,\eta}\left(\left. \frac{\partial v_2^*(Z;\theta,\eta)}{\partial\eta}\right\rvert X=x\right) = 0.
$$ It follows that $u_2^*$ is Neyman-orthogonal to order 1 for any $a_{22}$. We will now choose $a_{22} $ such that $u_2^*$ is Neyman-orthogonal to order 2.
Next, differentiating $u_2^*$ with respect to $\eta$ twice gives
\begin{equation*}
\begin{split}
\frac{\partial^2 u_2^*(z;\theta,\eta)}{\partial\eta^2}
& =
\frac{\partial^2 u_1^*(z;\theta,\eta)}{\partial\eta^2}
+
a_{22}(x;\theta,\eta) \,
\frac{\partial^2 v_2^*(z;\theta,\eta)}{\partial\eta^2}
\\
& +
\frac{\partial^2 a_{22}(x;\theta,\eta)}{\partial\eta^2}
\,
v_2^*(z;\theta,\eta)
+
2
\frac{\partial a_{22}(x;\theta,\eta)}{\partial\eta}
\,
\frac{\partial v_2^*(z;\theta,\eta)}{\partial\eta}.
\end{split}
\end{equation*}
Since $v_2^*$ has zero mean and is orthogonal to order 1, the terms involving the first and second derivative of $a_{22}$ drop out when taking expectations. It follows that $u_2^*$ in (\ref{eq_u2_star}) is Neyman-orthogonal to order 2 when one sets $a_{21}$ to its expression in (\ref{eq_a21}), and $a_{22}$ to
\begin{equation}
a_{22}(x;\theta,\eta)
=
\left(
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial^2 v_2^*(Z;\theta,\eta)}{\partial\eta^2}
\right\rvert X=x
\right)
\right)^{-1}
\mathbb{E}_{\theta,\eta}
\left(
\left.
\frac{\partial^2 u_1^*(Z;\theta,\eta)}{\partial\eta^2}
\right\rvert X=x
\right).\label{eq_a22}
\end{equation}
Note that this construction amounts to solving a system of linear equations. The fact that the solution in (\ref{eq_a21})--(\ref{eq_a22}) coincides with the expression in Theorem \ref{theo_neyman} may then again be verified by using Bartlett identities.
To appreciate the role of the likelihood structure in the above argument, suppose that $a_{21}$ and $a_{22}$ are not known up to the parameters $\theta$ and $\eta$. Then they are additional nuisance parameters and, thus, we require that all first- and second-order derivatives with respect to $(a_{21},a_{22})$, and $\eta$ be mean zero. The cross-derivatives between $(a_{21},a_{22})$ and $\eta$ are problematic, since having those to be mean zero would require
$$
\mathbb{E}_{\theta,\eta}\left(\left. \frac{\partial v_1(Z;\theta,\eta)}{\partial \eta} \right\vert X=x \right)=0,
\qquad
\mathbb{E}_{\theta,\eta}\left(\left. \frac{\partial v_2(Z;\theta,\eta)}{\partial \eta} \right\vert X=x \right) =0,
$$
which is not generally the case.
\setcounter{equation}{0}
\section{Examples\label{sec_ex}}
\subsection{Panel data models\label{subsec_panel}}
Consider an $N\times T$ panel data model with individual effects. Here, the likelihood factors across the cross-sectional observations and the likelihood contribution of unit $i$ takes the form
$$
\prod_{j=1}^T f(Y_{ij} \,|\, X_{ij}; \theta_{0}, \eta_{i0}).
$$
The maximum-likelihood estimator is well-known to suffer from a bias that is $O(T^{-1})$; see \cite{HahnNewey2004} and \cite{HahnKuersteiner2011} for derivations of this bias in static and dynamic models, respectively. Consider the estimation of $\theta_0$. The bias in the estimator comes from bias in the score stemming from estimation noise in the fixed effects. Taking $\eta_{i}$ to be scalar for notational simplicity, and letting $\widehat{\eta}_i$ be an estimator of $\eta_{i0}$, an expansion of the (normalized) score\footnote{In this discussion we work with the score divided by $T$, to facilitate the comparison with the panel data literature.}
$$
u(Z_i;\theta_0,\widehat{\eta}_i)
{=}
\frac{1}{T} \sum_{j=1}^T
\frac{\partial \log f(Y_{ij} \vert X_{ij};\theta_0,\widehat{\eta}_i)}{\partial\theta}
$$
yields
\begin{align*}
u(Z_i;\theta_0,\widehat{\eta}_i)=&
u(Z_i;\theta_0,\eta_{i0})
+
\frac{\partial u(Z_i;\theta_0,\eta_{i0})}{\partial \eta_{i}}
\,
(\widehat{\eta}_i - \eta_{i0})
+\frac{1}{2}
\frac{\partial^2 u(Z_i;\theta_0,\eta_{i0})}{\partial \eta_i^2}
\,
(\widehat{\eta}_i - \eta_{i0})^2\\&+o_P(\lvert\widehat{\eta}_i - \eta_{i0} \rvert^2).
\end{align*}
Taking expectations and re-arranging shows that
\begin{equation*}
\begin{split}
\mathbb{E}
(u(Z_i;\theta_0,\widehat{\eta}_i))
=&
\mathrm{cov}
\left(
\frac{\partial u(Z_i;\theta_0,\eta_{i0})}{\partial \eta_{i}},
\widehat{\eta}_i - \eta_{i0}
\right)
\\
& +
\mathbb{E}
\left(
\frac{\partial u(Z_i;\theta_0,\eta_{i0})}{\partial \eta_{i}}
\right)
\,
\mathbb{E}
(\widehat{\eta}_i - \eta_{i0})
\\
& +
\frac{1}{2}
\mathbb{E}
\left(
\frac{\partial^2 u(Z_i;\theta_0,\eta_{i0})}{\partial \eta_i^2}
\right)
\,
\mathbb{E} ((\widehat{\eta}_i - \eta_{i0})^2)
+
o(\mathbb{E}(\lvert \widehat{\eta}_i - \eta_{i0} \rvert^2)).
\end{split}
\end{equation*}
If we set $\widehat{\eta}_i = \widehat{\eta}_i(\theta_0) = \arg\max_{\eta} \prod_{j=1}^T \log f(Y_{ij} \,|\, X_{ij}; \theta_0, \eta)$, the maximum-likelihood estimator (MLE) given $\theta_0$, then each one of these terms is $O(T^{-1})$. If we use an estimator $\widehat{\eta}_i$ that is independent of the estimation sample, the first term disappears. However, the remaining terms, which capture the nonlinearity bias and variance in the estimator of $\eta_{i0}$, remain.
\cite{HahnNewey2004}, \cite{ArellanoHahn2007}, and \cite{DhaeneJochmans2015b,DhaeneJochmans2015a} present estimators of these terms based on the MLE that can be used to construct a bias-corrected estimator.
\cite{Lancaster2002} and \cite{woutersen2002robustness} integrate-out the fixed effects using a uniform prior after orthogonalizing to order 1 to obtain an estimator with bias $o(T^{-1})$; \cite{Arellano2003} presents an alternative derivation of the same result. First-order Neyman-orthogonality, by itself, does not suffice as it does not handle the third term in the expansion,
that is, it does not properly correct for the noise in the estimated fixed effects. \cite{LiLindsayWaterman2003}, building on \cite{WatermanLindsay1996}, show that their (second-order) projected score for $\theta$, when evaluated at $\widehat{\eta}_i(\theta)$, is a first-order unbiased estimating equation for $\theta$. Thus, here, a sample-splitting procedure is not needed to achieve bias reduction. This is a consequence of the (second- or higher-order) projected score being orthogonal to the influence function of $\widehat{\eta}_i(\theta)$. While interesting, it is not clear whether this property extends to higher-order projections or to other parameters of interest, such as average marginal effects.
More generally, with $\widehat{\eta}_i - \eta_{i0} = O_P(T^{-\nicefrac{1}{2}})$, the score admits a higher-order expansion of the form,
$$
\mathbb{E}
(u(Z_i;\theta_0,\widehat{\eta}_i))
=
\frac{B_1}{T}
+
\frac{B_2}{T^2} + \cdots + \frac{B_q}{T^q} + o(T^{-q})
$$
for constants $B_1,B_2,\ldots, B_q$.
The maximum-likelihood estimator has $B_1\neq 0$, in general, and so requires that $\nicefrac{N}{T}\rightarrow 0$ to be asymptotically unbiased. The approaches to bias correction mentioned above remove $B_1$ but not the remaining terms. Approaches that estimate and subsequently remove all $B_p$, $1\leq p \leq 1$, are given by \cite{DhaeneJochmans2015b,DhaeneJochmans2015a}. Likewise, an estimator based on Neyman-orthogonalization, combined with a sample-splitting estimator that uses preliminary estimators that satisfy $\widehat{\eta}_i - \eta_{i0} = O_P(T^{-\nicefrac{1}{2}})$, can be used to obtain the same result.
\paragraph{Example: Neyman-Scott model (continued).}
Recall that the (un-normalized) unit-specific score for $\sigma^2$ is given by (\ref{eq_NS}). The leading two elements of the Bhattacharyya basis for $\eta_i$ are
$$
v_1(Y_i;\sigma^2,\eta_i)
=
\sum_{j=1}^T \frac{Y_{ij}-\eta_i}{\sigma^2},
\qquad
v_2(Y_i;\sigma^2,\eta_i)
=
-
\frac{T}{\sigma^2}
+
\left( \sum_{j=1}^T \frac{Y_{ij}-\eta_i}{\sigma^2} \right)^2.
$$
We apply Theorem \ref{theo_neyman}. A small calculation yields $A(\sigma^2,\eta_i) = (0, \nicefrac{1}{2T})^\top$ and, after re-arranging,
$$
u_2^*(Y_i;\sigma^2,\eta_i)
=
\frac{1}{2\sigma^2}
\left(
\frac{\sum_{j=1}^T (Y_{ij}-\overline{Y}_i)^2}{\sigma^2}
-
(T-1)
\right),
$$
which does not depend on $\eta_i$. Summing over the cross-sectional units gives the second-order orthogonalized score equation for $\sigma^2$ as
$$
\sum_{i=1}^N
u_2^*(Y_i;\sigma^2,\eta_i)
=
\frac{1}{2\sigma^2}
\left(
\frac{\sum_{i=1}^N\sum_{j=1}^T (Y_{ij}-\overline{Y}_i){^2}}{\sigma^2}
-
N(T-1)
\right) = 0,
$$
which yields the degrees-of-freedom corrected estimator $\widehat{\sigma}^2$ in (\ref{eq_sig2_hat}).
Another parameter of interest in this problem is $\mu = \nicefrac{1}{N} \sum_{i=1}^N h(\eta_i)$, for $h$ a known function. This fits our framework, with
$$
u(Y_i; \sigma^2, \eta_i,\mu)
=
h(\eta_i) - \mu.
$$
The $q$-th orthogonalized counterpart to $u$ given by Theorem \ref{theo_neyman} is available in closed form, as
\begin{equation}u_{q}^*(Y_i; \sigma^2, \eta_i,\mu)=h(\eta_i) - \mu+\sum_{k=1}^q \sigma^{k}T^{-\frac{k}{2}}\frac{1}{k!}\nabla_{\eta}^{(k)}h(\eta_i)H_{k}\left(\frac{\sum_{j=1}^T(Y_{ij}-\eta_{i})}{ \sqrt{T}\sigma}\right),\label{eq_ustar_NS}\end{equation}
where $H_{k}$ denotes the $k$-th Hermite polynomial.
Interestingly, (\ref{eq_ustar_NS}) shows that, depending of the form of $h$, the variance of $u_{q}^*$ may converge or diverge as $q$ tends to infinity. Indeed, using a property of Hermite polynomials, the variance is
$$\mbox{Var}\left(u_q^*(Y_i; \sigma^2, \eta_i,\mu)\right)=\sum_{k=1}^q \sigma^{2k}T^{-k}\frac{1}{k!}\left[\nabla_{\eta}^{(k)}h(\eta_i)\right]^2.$$
It is instructive to consider the following three cases: (1) when $h(\eta_i)$ is a polynomial of order $K$ the series is stationary for $q\geq K$; (2) when $h(\eta_i)=\exp(\eta_i)$ the series converges; (3) in contrast, when $h(\eta_i)=\log(\eta_i)$ the series diverges. In the first two cases, using a large $q$ reduces bias without causing variance to diverge, while the third case presents a sharp trade-off between reduced bias and exploding variance as $q$ increases.
\subsection{Nonlinear network regression\label{sec_nonlin}}
Our next example is the nonlinear regression model with $d\geq 1$ outcomes,
\begin{equation}Y_i=m(X_i;\theta_0,\eta_{i0})+\sigma(X_i;\theta_0)\varepsilon_i,\quad \varepsilon_i\,|\, X\sim \mathrm{iid}~{\cal{N}}(0,I_d),\label{eq_network_nonlinreg}\end{equation}
where $m(x;\theta,\eta_i)$ is a $d\times 1$ vector, $\sigma(x;\theta)$ is an $d\times d$ diagonal matrix, and $m$ and $\sigma$ are known functions. We will show below that our CES production function example, in logarithms, fits into this framework.
For this model there are no analytical solutions for the orthogonalized estimators. We thus proceed numerically. To construct Neyman-orthogonal moment functions according to Theorem \ref{theo_neyman} we need to compute $\varSigma_{w_qw_q}(x;\theta,\eta)$, $\varSigma_{w_qu}(x;\theta,\eta,\mu)$, and $b_q(x;\theta,\eta,\mu)$, which involve higher-order derivatives of the conditional likelihood. To compute these derivatives, it is convenient to introduce the operator $\nabla_{m}^q$ that collects all derivatives with respect to
the $d$-vector $m$ up to order $q$. By the chain rule,
$$\nabla_{\eta_i}^q \ell(y\,|\, x;\theta,\eta_i)=M(x,\theta,\eta_i) \nabla_{m}^q \ell(y\,|\, x;\theta,\eta_i),$$
where the matrix $M$ has an analytical expression given by the multivariate Fa\`a di Bruno formula (\citealp{constantine1996multivariate}).
Given the matrix $M$ it is easy to compute $\varSigma_{w_qw_q}$, $\varSigma_{w_qu}$, and $b_q $. For example,
\begin{align*}
&\varSigma_{w_qw_q}(x;\theta,\eta_i)\\
&=M(x,\theta,\eta_i) \, \mathbb{E}_{\theta,\eta_i}\left(\frac{\nabla_{m}^q \ell(Y_i\,|\, X_i;\theta,\eta_i)}{\ell(Y_i\,|\, X_i;\theta,\eta_i)}\, \frac{\nabla_{m}^q \ell(Y_i\,|\, X_i;\theta,\eta_i)}{\ell(Y_i\,|\, X_i;\theta,\eta_i)}^\top\, \Bigg|\, X_i=x\right) M(x,\theta,\eta_i)^{\top},
\end{align*}
where the expectation on the right-hand can be readily computed by relying on formulas for moments of Hermite polynomials. We relegate further details to Appendix \ref{App_sec_implement}. In the next section we present simulations and an empirical application based on a version of (\ref{eq_network_nonlinreg}) designed to study team production.
\paragraph{Example: CES production function (continued).}
Consider the team production model
\begin{equation}Y_j=\beta_0(s_j)\left(\frac{1}{s_j}\sum_{r=1}^{s_j}\eta_{k(j,r)0}^{\gamma_0(s_j)}\right)^{\frac{1}{\gamma_0(s_j)}}\varepsilon_j^{\sigma_0(s_j)},\quad \log \varepsilon_j\,|\, {{\cal{K}}}\sim \mathrm{iid} ~{\cal{N}}\left(0,1\right),\label{eq_teams}\end{equation}
where $s_j$ is the size of team $j=1,...,n$, $(k(j,1),\ldots, k(j,s_j))$ are the $s_j$ workers in team $j$, and the set ${{\cal{K}}}=\{k(j,r)\,:\, r=1,\ldots,s_j,\, j=1,\ldots,n \}$ collects the workers in all teams. Model (\ref{eq_teams}) generalizes Model (\ref{eq_CES_size2}) by allowing for teams of varying sizes. Here we focus on teams of size 1 and 2, as in our application, and impose the normalization $\beta_0(1)=1$. For simplicity we will denote $\beta_0=\beta_0(2)$ and $\gamma_0=\gamma_0(2)$, which are the team size and substitution parameters, respectively, in teams of size 2.
We now explain how (\ref{eq_teams}) can be written as a special case of (\ref{eq_network_nonlinreg}), for a suitable choice of subsets of observations. To any team $j$ of size 2 involving workers $k$ and $k'$, we associate a team $j_1(j)$ of size 1 only involving worker $k$, and a team $j_2(j)$ of size 1 only involving worker $k'$. This construction results in $N$ subsets of three teams each. We then write the outcomes for these three teams, in logarithms, as
\begin{align}
\log Y_j&=\log \beta_0+\frac{1}{\gamma_0}\log\left(\frac{\eta_{k(j,1)0}^{\gamma_0}+\eta_{k(j,2)0}^{\gamma_0}}{2}\right)+\sigma_0(2)\log \varepsilon_j,\label{eq_subnet_1}\\
\log Y_{j_1(j)}&=\log\eta_{k(j,1)0}+\sigma_0(1)\log \varepsilon_{j_1(j)},\label{eq_subnet_2}\\
\log Y_{j_2(j)}&=\log\eta_{k(j,2)0}+\sigma_0(1)\log \varepsilon_{j_2(j)},\label{eq_subnet_3}
\end{align}
which takes the same form as (\ref{eq_network_nonlinreg}), for $d=3$, $\theta=\left(\beta_0,\gamma_0,\sigma_0^2(1),\sigma_0^2(2)\right)^\top$, $Y_i$ the vector of the three outcomes in (\ref{eq_subnet_1})--(\ref{eq_subnet_3}) for subset $i$, and $\eta_{i0}$ the $2\times 1$ vector of worker-specific effects in the corresponding teams.
\begin{remark}{(Implementation in other models)}
In models with discrete outcomes $Y$, one can express the matrices that feature in the expression for $A(x;\theta,\eta,\mu)$ in Theorem \ref{theo_neyman} in closed form, as sums over the support of $Y$. In models with continuous outcomes, one can proceed by simulation as follows, in the spirit of the ``reparameterization trick'' (\citealp{kingma2014auto}). Write $Y=g(X,U;\theta,\eta)$ where $U\,|\, X\sim F_U$ (for example, a standard multivariate Gaussian). Let $U^{(s)}$, $s=1,...,S$, be i.i.d. draws from $F_U$, and let $Y^{(s)}=g(x,U^{(s)};\theta,\eta)$ and $Z^{(s)}=(Y^{(s)},x)$. Assuming that $g$ is a smooth function of $\eta$ one can construct the simulation-based counterpart
\begin{align*}\widehat A(x;\theta,\eta,\mu)&=\left(\sum_{s=1}^Sw_q(Z^{(s)};\theta,\eta)w_q(Z^{(s)};\theta,\eta)^\top\right)^{-1}\\
&\times\left[\sum_{s=1}^Sw_q(Z^{(s)};\theta,\eta)u(Z^{(s)};\theta,\eta,\mu)^\top-\sum_{s=1}^S\nabla_{\eta}^qu\left(g(x,U^{(s)};\theta,\eta),x;\theta,\eta,\mu\right)^\top\right].\end{align*}
When focusing on $\theta_0$ instead of $\mu_0$, one can rely on the simpler expression
\begin{align*}\widehat A(x;\theta,\eta,\mu)=&\left(\sum_{s=1}^Sw_q(Z^{(s)};\theta,\eta)w_q(Z^{(s)};\theta,\eta)^\top\right)^{-1}\sum_{s=1}^Sw_q(Z^{(s)};\theta,\eta)u(Z^{(s)};\theta,\eta)^\top,\end{align*}
and $u_q^*$ can be obtained by regressing $u(Z^{(s)};\theta,\eta)$ on $w_q(Z^{(s)};\theta,\eta)$. We leave the study of the impact of a finite number $S$ of draws on inference about $\mu_0$ and $\theta_0$ to future work.
\end{remark}
\setcounter{equation}{0}
\section{Application to team production\label{sec_appli}}
\subsection{Model, data, and implementation}
We wish to estimate the parameters of the team production model in \eqref{eq_subnet_1}--\eqref{eq_subnet_3}. We will be especially interested in estimating the substitution parameter $\gamma$, which drives the nature of complementarities in teams of size 2, and the team size parameter $\beta$, which reflects the premium (or penalty) associated with working together relative to working alone. In addition to estimating production-function parameters, we will also report estimates of a counterfactual random re-allocation of workers to teams. Under random assignment, average output in teams of size 2 can be written as
\begin{equation}\mathbb{E}^{\rm rand}(Y_j)=\frac{2}{n_2(n_2-1)}\sum_{k_1<k_2}{\beta_0}\left[\frac{1}{2}\left({\eta_{k_10}}^{{\gamma_0}}+{\eta_{k_20}}^{{\gamma_0}}\right)\right]^{\frac{1}{{\gamma_0}}}\exp\left(\frac{1}{2}{\sigma_0^2(2)}\right),\label{eq_rand_alloc}
\end{equation}
where $n_2$ denotes the number of teams of size 2. As this quantity is an average over the worker fixed effects, it can be orthogonalized with respect to them using our approach.
\citet{AhmadpoodJones2019} consider model (\ref{eq_teams}) without the error term $\varepsilon_j$. Here our goal is to address the statistical challenge caused by the presence of a large number of possibly imprecisely estimated fixed effects. An alternative would be to specify a distribution for author heterogeneity conditional on the team network (i.e., for all the $\eta_{i0}$'s conditional on ${{\cal{K}}}$), as in \citet{Bonhomme2021}. An advantage of such a procedure is that, under correct specification, estimates are consistent even in poorly connected networks. This random-effect approach requires, however, to model how authors sort and collaborate in teams. Our approach avoids the need to do so. On the other hand, a fixed-effect approach requires that the author effects can be consistently estimated. In less well-connected networks, the convergence rate will be slower. Orthogonalization to a higher-order allows us to reduce the impact of estimation noise.
We look at the production of academic work in economics. We use data from \citet{DuctorFafchampsGoyalvanderLeij2014}, drawn from the EconLit database. These data contain a large collection of articles, indicated by their ID, together with author identifiers and a measure of {journal quality} proposed by \citet{kodrzycki2006new}. This measure is a ranking between 0 and 100, which we net of multiplicative time effects and will use as our outcome variable. In order to mitigate the variation in author fixed effects over time while ensuring sufficiently many collaborations, we restrict the sample to articles published between 1990 and 1999, written either alone or with a single co-author.{\footnote{We have also estimated the model on the entire sample, which ranges from 1970 to 1999. While the estimates of the substitution parameter differ somewhat in this larger sample, they are also less stable due to the fact that connectivity is lower. The estimates of the other parameters are similar to the ones on the 1990-1999 subsample.} We only include authors who produced {at least two sole-authored articles} during the sampling period.
Our sample contains 91,626 articles, 10\% of which are co-authored, and 16,408 authors. Average journal quality differs greatly across authors, with the 10th percentile of the quality measure being 0.4, the median being 0.9, and the 90th percentile being 8.5. The between-author variance in journal quality is 42\% of the overall variance. The distribution of journal quality, in turn, is skewed to the right, with a median of 0.6, a 90th percentile of 12, and a 99th percentile of 52. The number of publications per author varies substantially, with a 10th percentile of 2, a median of 4, and a 90th percentile of 13.
To implement our approach, we construct subsets of three papers, one co-authored ($j$) and two sole-authored ($j_1(j),j_2(j)$), as described in (\ref{eq_subnet_1})--(\ref{eq_subnet_3}). The score for $\theta$ based on subset $i$ then involves the three teams $j$, $j_1(j)$, and $j_2(j)$. Proceeding in this way is helpful as it limits the dimension of the parameter $\eta_i$ to two. This is not only in line with the assumptions we make in deriving asymptotics, but also helpful in terms of computation. Moreover, it reduces the number of derivatives that need to be computed. The number of derivatives nevertheless remains substantial, as we need to compute $9$ derivatives at order 2, $19$ at order 3, and $55$ at order 5, for example. Yet, using the computational remarks from Section \ref{sec_nonlin}, this can be implemented quite fast.
Finally, we exploit the network structure of the data to perform our sample splitting. For every worker, we construct a preliminary estimator of her fixed effect (in logs) as the average quality of her single-authored papers, except for one that we select at random and use later in estimation. This strategy is feasible due to our sample restriction. For each subset $i$ of three teams, we then stack the two worker fixed effects together to form our preliminary estimate $ \widehat{\eta}_i$. We next estimate the parameters $\beta_0,\gamma_0,\sigma_0^2(1),\sigma_0^2(2)$ on the sample from which all these single-authored articles have been removed. In the present case, $\widetilde{u}$ in (\ref{thetaEstimation}) has four components that correspond to the score with respect to all the parameters, and the weight matrix $\widetilde{W}$ is irrelevant since the problem is just-identified. In order to limit the variability due to the choice of split, we average parameter estimates across 100 random splits, through cross-fitting. The bias in the parameter estimates takes a complex form due to the team network environment. In Appendix \ref{AppC} we assess the ability of our orthogonalization approach to alleviate this bias in a Monte Carlo simulation.
\subsection{Empirical estimates}
Table \ref{tab_appli} shows the estimates of $\beta_0$, $\gamma_0$, $\sigma_0^2(2)$, and $\sigma_0^2(1)$ for various estimators. These are the plug-in estimator based on the preliminary estimates $\widehat{\eta}_i$ and six estimators based on Neyman-orthogonalized moments, for $1\leq q \leq 6$. In addition to point estimates, we report standard errors based on the parametric bootstrap.\footnote{Bootstrap replications are based on Neyman-orthogonalized estimates of $\beta_0$, $\gamma_0$, $\sigma_0^2(2)$, and $\sigma_0^2(1)$ to order $q=6$, together with the sample-split estimates $\widehat{\eta}_i$ of author effects. Within each bootstrap replication, we cross-fit the estimates 10 times. Results are based on 200 bootstrap replications.}
\begin{table}
\caption{Estimation results\label{tab_appli}}
\begin{center}
\begin{tabular}{|c||cccc|c|}\hline\hline
& Substitution ${\gamma}$ & Team size ${\beta}$ & Variance ${\sigma^2(2)}$ & Variance ${\sigma^2(1)}$ & Diagnostic \\ \hline\hline
Plug-in & $\underset{\small (0.0466)}{0.1267}$ & $\underset{\small (0.0217)}{1.2893}$ & $\underset{\small (0.0249)}{1.6232}$ & $\underset{\small (0.0204)}{1.6395}$ & - \\
$q=1$ & $\underset{\small (0.2123)}{-1.8905}$ & $\underset{\small (0.0279)}{1.3260}$ & $\underset{\small (0.0264)}{1.6793}$ & $\underset{\small (0.0265)}{1.7330}$ & $\underset{[0.0000]}{1285.7}$\\
$q=2$ & $\underset{\small (0.2367)}{0.7678}$ & $\underset{\small (0.0359)}{1.2956}$ & $\underset{\small (0.0260)}{1.4339}$ & $\underset{\small (0.0260)}{1.4683}$ &$\underset{[0.0000]}{63.846}$\\
$q=3$ & $\underset{\small (0.2034)}{0.4702}$ & $\underset{\small (0.0369)}{1.2845}$ & $\underset{\small (0.0254)}{1.4399}$ & $\underset{\small (0.0236)}{1.4407}$ &$\underset{[0.0006]}{19.490}$\\
$q=4$ & $\underset{\small (0.1763)}{0.4176}$ & $\underset{\small (0.0362)}{1.2838}$ & $\underset{\small (0.0254)}{1.4348}$ & $\underset{\small (0.0232)}{1.4214}$ & $\underset{[0.3126]}{4.7622}$ \\
$q=5$ & $\underset{\small (0.1730)}{0.4143}$ & $\underset{\small (0.0361)}{1.2836}$ & $\underset{\small (0.0254)}{1.4327}$ & $\underset{\small (0.0231)}{1.4180}$ & $\underset{[0.7890]}{1.7096}$ \\
$q=6$ & $\underset{\small (0.1770)}{0.4145}$ & $\underset{\small (0.0365)}{1.2836}$ & $\underset{\small (0.0254)}{1.4315}$ & $\underset{\small (0.0230)}{1.4165}$ & - \\ \hline\hline
\end{tabular}
\end{center}
\raggedright
\par\textit{{\footnotesize Notes: Point estimates based on $q$-ordered orthogonalized estimators, cross-fitted estimates (100 splits). Parametric bootstrap standard errors in parentheses (200 replications). In the last column we report the diagnostic for the orthogonality order $q$ (see Appendix \ref{App_q}, bootstrapped with 200 replications), with p-values based on the $\chi^2(4)$ distribution in brackets.}}
\end{table}
Starting with the substitution parameter $\gamma$, the uncorrected estimate is $0.13$, which is close to the Cobb-Douglas case. The value of the first-order Neyman-orthogonalized estimate is quite different. However, since the preliminary estimates of the author fixed effects are based on very few observations, we do not expect this estimator to adequately correct for bias. This is confirmed by the fact that all other Neyman-orthogonalized estimates, for $q\in\{2,\ldots,6\}$, range between $0.41$ and $0.77$, which is higher than the plug-in estimate, and very different from the first-order orthogonalized estimate. Relative to the plug-in, the orthogonalized estimates with $q\geq 2$ all indicate somewhat less complementarity between authors in team production. Notice the stability of estimates for larger values of $q$. A substitution parameter $\gamma=0.4$ corresponds to the case of imperfect complements; see Figure \ref{fig_prodf} in Appendix \ref{App_fig} for a graphical illustration.
Turning to the other parameters, the estimates of the team size parameter $\beta$ are virtually unaffected by the orthogonalization. This suggests the bias is limited for this parameter. Its value is close to $1.3$, implying that producing a paper with a co-author increases the paper's quality to some extent. Next, the log-error variance $\sigma^2(2)$ in teams of two coauthors is larger when using plug-in estimates ($1.6$) than when using orthogonalization with $q\geq 2$ ($1.4$), suggesting that the plug-in and first-order corrected estimates are biased upward. Lastly, the variance $\sigma^2(1)$ in teams of a single author is also larger under the plug-in estimator.
An interesting feature of Table \ref{tab_appli} is that point-estimates and standard errors appear to converge as $q$ increases. This suggests the absence of a sharp bias-variance trade-off as a function of $q$. As we have seen in the case of the Neyman-Scott example in Subsection \ref{subsec_panel}, this phenomenon is specific to the model and target parameter of interest, and the variance may converge or diverge as $q$ grows depending on the context.
One potential explanation for convergence in the present setting is that Model (\ref{eq_teams}) implies some non-trivial restrictions on the parameters $\gamma,\beta,\sigma^2(1),\sigma^2(2)$ that do not depend on the author-specific effects $\eta_i$, as we show in Appendix \ref{AppD}. Moment functions that are independent of $\eta$ are common in panel data settings, obtained e.g. by differencing, quasi-differencing, or functional differencing. In some models, such as models with discrete outcomes, such functions may not exist. In Appendix \ref{AppD} we exploit two types of restrictions as robustness checks. Our findings suggest that, while those restrictions seem broadly consistent with the higher-order orthogonal estimates reported in Table \ref{tab_appli}, using them directly for estimation may lead to very imprecise estimates. In contrast, Table \ref{tab_appli} suggests that higher-order Neyman-orthogonality may be a successful approach to approximate such functions in settings where those exist.
In the last column of Table \ref{tab_appli} we report our diagnostic statistic for the orthogonality order $q$, described in Appendix \ref{App_q}.\footnote{The statistic in Appendix \ref{App_q} depends on a variance matrix estimate $\widehat V$, which we compute based on the parametric bootstrap (200 replications).} A value larger than the 95-th quantile of the $\chi^2(4)$ distribution (9.488) should be interpreted as suggesting that $q$ is too low. The values of the statistic reported in the table, together with the associated p-values, suggest that values of $q\leq 3$ are too low, and values $q\geq 4$ are sufficiently large for biases to be asymptotically negligible.
\begin{table}[tbp]
\caption{Empirical estimates: average output\label{tab_appli2}}
\begin{center}
\begin{tabular}{|c||cc|}\hline\hline
& Average output & Counterfactual\\\hline\hline
Plug-in & $\underset{\small (0.2492)}{8.4434}$ & $\underset{\small (0.1949)}{7.2048}$ \\
$q=1$ & $\underset{\small (0.4837)}{6.8366}$ & $\underset{\small (0.3936)}{5.4218}$ \\
$q=2$ & $\underset{\small (0.5096)}{9.0903}$ & $\underset{\small (0.7626)}{8.7391}$ \\
$q=3$ & $\underset{\small (0.4452)}{6.1290}$ & $\underset{\small (0.4354)}{5.4619}$ \\
$q=4$ & $\underset{\small (0.4726)}{7.4412}$ & $\underset{\small (0.6227)}{6.7197}$ \\
$q=5$ & $\underset{\small (0.3889)}{6.9758}$ & $\underset{\small (0.3896)}{6.2202}$ \\
$q=6$ & $\underset{\small (0.3194)}{7.1168}$ & $\underset{\small (0.4225)}{6.3867}$ \\ \hline\hline
\end{tabular}
\end{center}
\raggedright
\par\textit{{\footnotesize Notes: Average output (the value in the data is $6.9995$, standard error $0.2309$), and counterfactual average output in a random allocation. Point estimates based on orthogonalized estimators to order $q$, cross-fitted estimates (100 splits). Parametric bootstrap standard errors in parentheses (200 replications).}}
\end{table}
Lastly, we report estimates of average journal quality in a counterfactual scenario where authors are randomly assigned across teams of two co-authors, see (\ref{eq_rand_alloc}). The first column in Table \ref{tab_appli2} shows estimates of the average output in the empirical allocation. This quantity can be estimated without bias as the sample mean of the journal quality variable, which is equal to $7.0$. We see that the plug-in estimate is $8.4$, larger than the empirical value. In comparison, Neyman-orthogonalized estimates for $q\geq 3$ range between $6.1$ and $7.4$, and estimates for $q=5$ and $q=6$ are closest to the empirical value. The second column in Table \ref{tab_appli2} shows estimates of average article quality under random assignment of authors to teams, using the plug-in method and Neyman-orthogonalized estimates to order $q\geq 1$.\footnote{To speed up computation, we approximate (\ref{eq_rand_alloc}) using a random subset of 1000 authors, for each random sample split (and each bootstrap replication).} The estimates vary with the order of orthogonalization. When taking $q\geq 4$, estimates range between $6.2$ and $6.7$. In addition, comparing the two columns of Table \ref{tab_appli2} shows that, irrespective of the order of orthogonalization, the estimates of average output are lower in the counterfactual scenario where workers are randomly allocated across teams.
The main takeaway from Table \ref{tab_appli2} is that randomly allocating authors among teams would tend to lower average paper quality. This is due to two economic forces. The first one is complementarity in production, as reflected by estimates of $\gamma$ lower than $1$. The second force is positive sorting. Indeed, the preliminary estimates of worker fixed effects are positively correlated within teams in the data. In the presence of complementarity, decreasing assortative matching leads to lower output, which is what we find in Table \ref{tab_appli2}.
\setcounter{equation}{0}
\section{Asymptotic properties\label{sec_asympt}}
In this section, we show that, under higher-order orthogonality,
the estimators $\widehat{\theta}$ and $\widehat{\mu}$
introduced in Section~\ref{subsect:Estimation} are $\sqrt{n}$-consistent and asymptotically normal
under appropriate assumptions, even if the convergence rate of $\widehat{\eta}_i$ is slower than $\sqrt{n}$. We focus on deriving the asymptotic distribution of $\widehat \mu$, assuming that we have already worked out the corresponding asymptotic result of $\widehat \theta$. However, the corresponding theory for $\widehat \theta$ is actually a special case of our results for $\widehat \mu$, where $\theta$ is dropped from the arguments, $\mu$ is replaced by $\theta$, and $u$ is replaced by $\widetilde u$. Thus, our focus on $\widehat \mu$ is without loss of generality.
\subsection{Notation}
For the presentation of the asymptotic theory, it is useful to be explicit about which parameters depend on the sample size and which ones do not. Recall that $n$ is the total number of observations in $(Z_1,...,Z_N)$, where each $Z_i$ comprises $n_i$ observations. In the asymptotic sequence, we let $N$ and $n_i$ depend on $n$, although we do not explicitly indicate this dependence. For example, in a panel data model, our assumptions allow both $N$ and $T$ to grow as the number $NT$ of observations tends to infinity.
To indicate the dependence on the sample size, we will write $\eta_n$ and $\mu_n$ instead of $\eta$ and $\mu$ in this section. While the dimension of $\mu$ is not changing with $n$, the true parameter $\mu_{0,n}$ is implicitly defined as the solution
of $\sum_{i=1}^N\mathbb{E}_{\theta_0,\eta_{0,n}}\left( u(Z_i; \theta_0, \eta_{0,n,i},\mu)\right) = 0$, which may depend on $n$. By contrast, the parameter $\theta$ and its true value $\theta_0$ are independent of $n$.
Remember also that $n=\sum_{i=1}^N n_i$, and note that if the observations within each unit $i$
are independent, then we have $\ell(y_i\,|\, x_i;\theta,\eta_{n,i}) = \prod_{j=1}^{n_i}
\ell(y_{ij}\,|\, x_{ij};\theta,\eta_{n,i})$. Hence, in the case of the score for $\theta$,
\begin{align*}
u(Z_i; \theta, \eta_{n,i}) =
\sum_{j=1}^{n_i} \frac{\partial \log\ell(y_{ij}\,|\, x_{ij};\theta,\eta_{n,i})}{\partial\theta}.
\end{align*}
More generally, whenever $n_i \to \infty$ we expect that $u$ scales linearly
with $n_i$, {explaining the scaling of $u(Z_i; \theta, \eta_{n,i})$ and of various other terms in Assumption~\ref{ass:asymptotic1} below.}
\subsection{A useful lemma}
With this notation in hand, we now state our first assumption.
\begin{assumption}
\label{ass:asymptotic1}
\phantom{a}
\begin{enumerate}[(i)]
\item We have
$\left[\frac 1 n \sum_{i=1}^N
\frac{\partial u^\top (Z_i; \widehat \theta,\widehat \eta_{n,i},\widehat\mu_n )} {\partial \mu}
\right] W \left[\frac 1 {\sqrt{n}} \sum_{i=1}^N u(Z_i; \widehat \theta,\widehat \eta_{n,i},\widehat\mu_n ) \right] = o_P(1)$, for some non-random symmetric
positive definite weight matrix $W$.
\item As $n\rightarrow \infty$, $(\widehat \theta,\widehat \eta_n , \widehat \mu_n)$ is contained
in a convex neighborhood ${\cal B}_n$ of $(\theta_0,\eta_{0,n},\mu_{0,n})$. Let ${\cal B}_{n,i}$ be the convex neighborhood of $(\theta_0,\eta_{0,n,i},\mu_{0,n})$ obtained by intersecting ${\cal B}_n$ with the parameter parameter subspace for observation $i$.
\item
$\max_i \dim(\eta_{n,i})=O(1)$.
\item For every $i$, the function $u(Z_i,\theta,\eta_{n,i},\mu)$ is $(q+1)$ times continuously differentiable
in the parameters $(\theta,\eta_{n,i},\mu)$, and we assume that for all its components all the partial derivatives of $u(Z_i;\theta,\eta_{n,i},\mu)$ up to order $(q+1)$ are bounded in
absolute value
by $n_iC_{n,i}(Z_i) \geq 0$, uniformly in the neighborhood ${\cal B}_{n,i}$, such that
$\frac{1}{n} \sum_{i=1}^N n_i\mathbb{E} \left[ C_{n,i}(Z_i)^2 \right] = O(1)$.
\item
$\widehat \mu_n - \mu_{0,n} = o_P(1)$ and $\frac 1 n \sum_{i=1}^N n_i \mathbb{E}\left( \left\|
\widehat \eta_{n,i}-\eta_{0,n,i} \right\|^{2(q+1)}\right) = o(n^{-1})$.
\item $\widehat \theta = \theta_0+\frac 1 n \sum_{i=1}^N \psi_{n,i} + o_P(n^{-1/2})$, where $\mathbb{E} (\psi_{n,i}) = 0$ and $\frac{1}{n} \sum_{i=1}^N \mathbb{E} \left(\left\| \psi_{n,i} \right\|^2 \right)=O(1)$.
\item The probability limits
$$G_\mu = \operatorname*{plim}_{n \rightarrow \infty} \frac 1 {n} \sum_{i=1}^N \frac{\partial u(Z_i; \theta_0, \eta_{0,n,i} , \mu_{0,n})}{\partial \mu^{\top}},\quad G_\theta = \operatorname*{plim}_{n \rightarrow \infty} \frac 1 {n} \sum_{i=1}^N \frac{\partial u(Z_i; \theta_0, \eta_{0,n,i} , \mu_{0,n})}{\partial \theta^{\top}}$$
exist, and ${\rm rank}(G_\mu)={\rm dim}(\mu)$.
\end{enumerate}
\end{assumption}
Part $(i)$ in Assumption \ref{ass:asymptotic1} is satisfied if $\widehat\mu_n$ is computed using GMM, see (\ref{muEstimation}). In Part $(ii)$, the neighborhood ${\cal B}_n$ depends on the sample size $n$, partly because the number of nuisance parameters of $\eta_{n,i}$ generally depends on $n$. Part $(iii)$ assumes that the maximal dimension of $ \eta_{n,i}$ is bounded as $n \rightarrow \infty$. Part $(iv)$ requires the derivatives of the moment functions (properly rescaled) to be suitably bounded. The first half of Part $(v)$ is a high-level
consistency assumption for $\widehat \mu_n$, which
can be justified by guaranteeing that the
objective function in \eqref{muEstimation} converges
uniformly to a population counterpart that has a
unique minimum at $\mu_0$. The second half of Part $(v)$ is the rate requirement on the preliminary estimates $\widehat{\eta}_{n,i}$, imposing a rate faster than $n^{-\nicefrac{1}{2(q+1)}}$. Part $(vi)$ requires $\widehat{\theta}$ to be asymptotically linear, in particular requiring $ \widehat \theta - \theta_0 = O_P(n^{-1/2})$. In the case where $\mu_{0,n}=\theta_0$ this condition is not needed. Lastly, Part $(vii)$ assumes existence of Jacobian matrices and a rank condition.
In the statement of the following lemma, $D^m_{\eta_{n,i}}$ denote the derivative operator with respect to $\eta_{n,i}$.
\begin{lemma}
\label{lemma:Expansion}
Under Assumption~\ref{ass:asymptotic1} we have
\begin{align*}
&\sqrt{n} \left( \widehat \mu_n - \mu_{0,n} \right)
\\
& \; \; = - \left( G^{\top}_{\mu} \, W \, G_\mu \right)^{-1} \, G^{\top}_{\mu} \,W \left\{ \frac 1 {\sqrt{n}} \sum_{i=1}^N
\Big[ u(Z_i; \theta_0, \eta_{0,n,i} , \mu_{0,n}) + G_{\theta} \, \psi_{n,i} \Big]
+ R_n
\right\} + o_P(1),
\end{align*}
where
\begin{align*}
R_{n} = \frac 1 {\sqrt{n}} \sum_{i=1}^N \sum_{m \in {\mathbb K}_{q,n,i}} \frac{1}{m!} \left[ D^m_{\eta_{n,i}} u(Z_i; \theta_0, \eta_{0,n,i} , \mu_{0,n}) \right] \left(\widehat \eta_{n,i} - \eta_{0,n,i}\right)^{m} ,
\end{align*}
and ${\mathbb K}_{q,n,i} = \left\{ m \in \mathbb{Z}^{{\rm dim}(\eta_{n,i})} \, : \, 1\leq
\sum_{r=1}^{{\rm dim}(\eta_{n,i})} m_r
\leq q \right\}$.
\end{lemma}
\subsection{Main result}
We are now in position to establish the main result of this section, which concerns root-$n$ consistency and asymptotic normality of estimators based on orthogonal equations. For this, we first state our second assumption.
\begin{assumption}
\label{ass:asymptotic2}
\phantom{a}
\begin{enumerate}[(i)]
\item The moment function
$u(Z_i; \theta, \eta_{n,i} , \mu)$ is Neyman-orthogonal to order $q$, and furthermore
$\sum_{i=1}^N\mathbb{E} \left( u(Z_i; \theta_0, \eta_{0,n,i} , \mu_{0,n})\right)=0$.
\item
$\widehat \eta_{n,i}$ are independent of $(Z_1,\ldots,Z_N)$ for all $i$.
\item
The $Z_1,\ldots,Z_N$ are independent across $i$.
\item
$\xi_{n,i}= u(Z_i; \theta_0, \eta_{0,n,i} , \mu_{0,n}) + G_{\theta} \, \psi_{n,i}$ satisfies Lindeberg's condition,\footnote{
That is, for any $\epsilon > 0$, $\frac{1}{s_n^2} \sum_{i=1}^N \mathbb{E}\left[ \xi_{n,i}^2 \cdot \mathbbm{1}(|\xi_{n,i}| > \epsilon s_n)\right] \to 0$ as $n \to \infty$, where $s_n^2 = \sum_{i=1}^N \mathrm{Var}(\xi_{n,i})$ and $\mathbbm{1}$ is the indicator function.}
and the following probability limit exists:
$$
V_\xi = \operatorname*{plim}_{n \rightarrow \infty} \frac 1 {n} \sum_{i=1}^N {\rm Var}\left( \xi_{n,i} \right) .
$$
\end{enumerate}
\end{assumption}
Part $(i)$ in Assumption \ref{ass:asymptotic2} requires $u$ to be Neyman-orthogonal in the sense of Definition \ref{def_ortho}. Part $(ii)$ requires the preliminary estimates to be independent from the estimation sample. With independent observations, this can be achieved by sample splitting. Part $(iii)$ imposes independence between the $Z_i$'s. We impose this assumption to simplify the presentation. It is straightforward to modify the variance expression in Theorem \ref{th:Asymptotic} below to account for particular forms of dependence (e.g., clustered) by using an appropriate expression for the matrix $V_{\xi}$ introduced in Part $(iv)$.
The following theorem provides an asymptotic characterization of $\widehat{\mu}_{n}$.
\begin{theorem}
\label{th:Asymptotic}
Let Assumptions~\ref{ass:asymptotic1} and \ref{ass:asymptotic2} hold with the same value of $q \in \{1,2,3,\ldots\}$. Then we have
\begin{align*}
\sqrt{n} \left( \widehat \mu_n - \mu_{0,n} \right) &\overset{d}{\rightarrow} {\cal N}\big(0, \, \left( G^{\top}_{\mu} \, W \, G_\mu \right)^{-1}
G^{\top}_{\mu} \,W \, V_\xi \, W G_\mu \left( G^{\top}_{\mu} \, W \, G_\mu \right)^{-1} \big) .
\end{align*}
\end{theorem}
Note that, although we leave the dependence on $q$ implicit in Theorem \ref{th:Asymptotic}, the asymptotic variance does depend on the order of orthogonality $q$ that the moment function satisfies. Note also that $ \widehat \mu_n$ in the theorem is based on a single set of preliminary estimates $\widehat{\eta}_i$. The variability caused by the use of a single sample split can be mitigated through the use of cross-fitting, as we do in the application.
\section{Final remarks}
In this paper we show how to construct higher-order Neyman-orthogonal moment functions in conditional-likelihood models. We use these functions, together with sample splitting, to reduce bias in estimation. Our application suggests that our higher-order corrections can be effective in network settings with fixed effects. An area of application is to double/debiased machine learning with fixed effects, where the nuisance parameters contains some components, such as low-dimensional functions, for which first-order orthogonality may suffice. However, for such applications it is important to extend the approach to non-likelihood models. As $q$ increases, orthogonalization imposes growing demands on the likelihood structure -- to achieve first-order orthogonality it is sufficient for the score to have mean zero, while to achieve second-order orthogonality our approach requires the information identity to hold, for example. {This reflects a trade-off between the robustness to nuisance parameters that our method achieves and robustness to model misspecification.} We are working on a strategy to construct orthogonal functions in semi-parametric models defined by moment conditions.
\clearpage
\begin{thebibliography}{}
\bibitem[\protect\citeauthoryear{Abowd, Kramarz, and Margolis}{Abowd, Kramarz
and Margolis}{1999}]{AbowdKramarzMargolis1999}
Abowd, J.~M., F.~Kramarz, and D.~N. Margolis (1999).
\newblock High wage workers and high wage firms.
\newblock {\em Econometrica\/}~{\em 67}, 251--333.
\bibitem[\protect\citeauthoryear{Ahmadpoor and Jones}{Ahmadpoor and
Jones}{2019}]{AhmadpoodJones2019}
Ahmadpoor, M. and B.~F. Jones (2019).
\newblock Decoding team and individual impact in science and invention.
\newblock {\em Proceedings of the National Academy of Sciences\/}~{\em 116},
13885--13890.
\bibitem[\protect\citeauthoryear{Andrews, Gill, Schank, and Upward}{Andrews,
Gill, Schank and Upward}{2008}]{AndrewsMartynGillSchankUpward}
Andrews, M.~J., L.~Gill, T.~Schank, and R.~Upward (2008).
\newblock High wage workers and low wage firms: negative assortative matching
or limited mobility bias?
\newblock {\em Journal of the Royal Statistical Society: Series A\/}~{\em 171},
673--697.
\bibitem[\protect\citeauthoryear{Angrist and Frandsen}{Angrist and
Frandsen}{2022}]{AngristFrandsen2022}
Angrist, J.~D. and B.~Frandsen (2022).
\newblock Machine labor.
\newblock {\em Journal of Labor Economics\/}~{\em 40}, S97--S140.
\bibitem[\protect\citeauthoryear{Arellano}{Arellano}{2003}]{Arellano2003}
Arellano, M. (2003).
\newblock Discrete choices with panel data.
\newblock {\em Investigaciones Economicas\/}~{\em XXVII}, 423--458.
\bibitem[\protect\citeauthoryear{Arellano and Bonhomme}{Arellano and
Bonhomme}{2012}]{arellano2012identifying}
Arellano, M. and S.~Bonhomme (2012).
\newblock Identifying distributional characteristics in random coefficients
panel data models.
\newblock {\em The Review of Economic Studies\/}~{\em 79}, 987--1020.
\bibitem[\protect\citeauthoryear{Arellano and Hahn}{Arellano and
Hahn}{2007}]{ArellanoHahn2007}
Arellano, M. and J.~Hahn (2007).
\newblock Understanding bias in nonlinear panel models: Some recent
developments.
\newblock In R.~Blundell, W.~K. Newey, and T.~Persson (Eds.), {\em Advances In
Economics and Econometrics}, Volume III. Econometric Society: Cambridge
University Press.
\bibitem[\protect\citeauthoryear{Bhattacharyya}{Bhattacharyya}{1946}]{Bhattacharyya1946}
Bhattacharyya, A. (1946).
\newblock On some analogues of the amount of information and their use in
statistical estimation.
\newblock {\em Sankhy{\=a}\/}~{\em 8}, 1--14.
\bibitem[\protect\citeauthoryear{Bickel}{Bickel}{1982}]{Bickel1982}
Bickel, P. (1982).
\newblock On adaptive estimation.
\newblock {\em Annals of Statistics\/}~{\em 10}, 647--671.
\bibitem[\protect\citeauthoryear{Bonhomme}{Bonhomme}{2021}]{Bonhomme2021}
Bonhomme, S. (2021).
\newblock Teams: Heterogeneity, sorting, and complementarity.
\newblock Mimeo.
\bibitem[\protect\citeauthoryear{Chernozhukov, Chetverikov, Demirer, Duflo,
Hansen, Newey, and Robins}{Chernozhukov
et~al.}{2018}]{ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018}
Chernozhukov, V., D.~Chetverikov, M.~Demirer, E.~Duflo, C.~Hansen, W.~Newey,
and J.~Robins (2018).
\newblock Double/debiased machine learning for treatment and structural
parameters.
\newblock {\em Econometrics Journal\/}~{\em 21}, C1--C68.
\bibitem[\protect\citeauthoryear{Chesher}{Chesher}{1991}]{Chesher1991}
Chesher, A. (1991).
\newblock The effect of measurement error.
\newblock {\em Biometrika\/}~{\em 78}, 451--462.
\bibitem[\protect\citeauthoryear{Constantine and Savits}{Constantine and
Savits}{1996}]{constantine1996multivariate}
Constantine, G. and T.~Savits (1996).
\newblock A multivariate {F}aa di {B}runo formula with applications.
\newblock {\em Transactions of the American Mathematical Society\/}~{\em 348},
503--520.
\bibitem[\protect\citeauthoryear{Dhaene and Jochmans}{Dhaene and
Jochmans}{2015a}]{DhaeneJochmans2015b}
Dhaene, G. and K.~Jochmans (2015a).
\newblock Profile-score adjustments for incidental-parameter problems.
\newblock Mimeo.
\bibitem[\protect\citeauthoryear{Dhaene and Jochmans}{Dhaene and
Jochmans}{2015b}]{DhaeneJochmans2015a}
Dhaene, G. and K.~Jochmans (2015b).
\newblock Split-panel jackknife estimation of fixed-effect models.
\newblock {\em The Review of Economic Studies\/}~{\em 82}, 991--1030.
\bibitem[\protect\citeauthoryear{Dhaene and Jochmans}{Dhaene and
Jochmans}{2016}]{DhaeneJochmans2016}
Dhaene, G. and K.~Jochmans (2016).
\newblock Likelihood inference in an autoregression with fixed effects.
\newblock {\em Econometric Theory\/}~{\em 32}, 1178--1215.
\bibitem[\protect\citeauthoryear{Ductor, Fafchamps, Goyal, and Van~der
Leij}{Ductor, Fafchamps, Goyal and Van~der
Leij}{2014}]{DuctorFafchampsGoyalvanderLeij2014}
Ductor, L., M.~Fafchamps, S.~Goyal, and M.~J. Van~der Leij (2014).
\newblock Social networks and research output.
\newblock {\em Review of Economics and Statistics\/}~{\em 96}, 936--948.
\bibitem[\protect\citeauthoryear{Evdokimov and Zeleneev}{Evdokimov and
Zeleneev}{2023}]{evdokimov2023simple}
Evdokimov, K.~S. and A.~Zeleneev (2023).
\newblock Simple estimation of semiparametric models with measurement errors.
\newblock {\em arXiv preprint arXiv:2306.14311\/}.
\bibitem[\protect\citeauthoryear{Ghazal and Neudecker}{Ghazal and
Neudecker}{2000}]{ghazal2000second}
Ghazal, G.~A. and H.~Neudecker (2000).
\newblock On second-order and fourth-order moments of jointly distributed
random matrices: a survey.
\newblock {\em Linear Algebra and its Applications\/}~{\em 321}, 61--93.
\bibitem[\protect\citeauthoryear{Graham}{Graham}{2020}]{graham2020sparse}
Graham, B.~S. (2020).
\newblock Sparse network asymptotics for logistic regression.
\newblock {\em arXiv preprint arXiv:2010.04703\/}.
\bibitem[\protect\citeauthoryear{Graham, Imbens, and Ridder}{Graham, Imbens and
Ridder}{2014}]{graham2014complementarity}
Graham, B.~S., G.~W. Imbens, and G.~Ridder (2014).
\newblock Complementarity and aggregate implications of assortative matching: A
nonparametric analysis.
\newblock {\em Quantitative Economics\/}~{\em 5}, 29--66.
\bibitem[\protect\citeauthoryear{Hahn and Hausman}{Hahn and
Hausman}{2021}]{hahn2021problems}
Hahn, J. and J.~Hausman (2021).
\newblock Problems with the control variable approach in achieving unbiased
estimates in nonlinear models in the presence of many instruments.
\newblock {\em Journal of Quantitative Economics\/}~{\em 19}, 39--58.
\bibitem[\protect\citeauthoryear{Hahn and Kuersteiner}{Hahn and
Kuersteiner}{2011}]{HahnKuersteiner2011}
Hahn, J. and G.~Kuersteiner (2011).
\newblock Bias reduction for dynamic nonlinear panel models with fixed effects.
\newblock {\em Econometric Theory\/}~{\em 27}, 1152--1191.
\bibitem[\protect\citeauthoryear{Hahn and Newey}{Hahn and
Newey}{2004}]{HahnNewey2004}
Hahn, J. and W.~K. Newey (2004).
\newblock Jackknife and analytical bias reduction for nonlinear panel models.
\newblock {\em Econometrica\/}~{\em 72}, 1295--1319.
\bibitem[\protect\citeauthoryear{Jackson, Rockof\mbox{}f, and Staiger}{Jackson,
Rockof\mbox{}f and Staiger}{2014}]{JacksonRockoffStaiger2014}
Jackson, C.~K., J.~E. Rockof\mbox{}f, and D.~O. Staiger (2014).
\newblock Teacher ef\mbox{}fects and teacher related policies.
\newblock {\em Annual Review of Economics\/}~{\em 6}, 801--825.
\bibitem[\protect\citeauthoryear{Jochmans and Weidner}{Jochmans and
Weidner}{2019}]{JochmansWeidner2019}
Jochmans, K. and M.~Weidner (2019).
\newblock Fixed-effect regressions on network data.
\newblock {\em Econometrica\/}~{\em 87}, 1543--1560.
\bibitem[\protect\citeauthoryear{Kingma and Welling}{Kingma and
Welling}{2014}]{kingma2014auto}
Kingma, D.~P. and M.~Welling (2014).
\newblock Auto-encoding variational bayes.
\newblock {\em stat\/}~{\em 1050}, 1.
\bibitem[\protect\citeauthoryear{Kline, Saggio, and S{\o}lvsten}{Kline, Saggio
and S{\o}lvsten}{2020}]{KlineSaggioSoelvsten2020}
Kline, P., R.~Saggio, and M.~S{\o}lvsten (2020).
\newblock Leave-out estimation of variance components.
\newblock {\em Econometrica\/}~{\em 88}, 1859--1898.
\bibitem[\protect\citeauthoryear{Kodrzycki and Yu}{Kodrzycki and
Yu}{2006}]{kodrzycki2006new}
Kodrzycki, Y.~K. and P.~Yu (2006).
\newblock New approaches to ranking economics journals.
\newblock {\em The BE Journal of Economic Analysis \& Policy\/}~{\em 5}.
\bibitem[\protect\citeauthoryear{Lancaster}{Lancaster}{2002}]{Lancaster2002}
Lancaster, T. (2002).
\newblock Orthogonal parameters and panel data.
\newblock {\em Review of Economic Studies\/}~{\em 69}, 647--666.
\bibitem[\protect\citeauthoryear{Li, Lindsay, and Waterman}{Li, Lindsay and
Waterman}{2003}]{LiLindsayWaterman2003}
Li, H., B.~Lindsay, and R.~Waterman (2003).
\newblock Efficiency of projected score methods in rectangular array
asymptotics.
\newblock {\em Journal of the Royal Statistical Society, Series B\/}~{\em 65},
191--208.
\bibitem[\protect\citeauthoryear{Mackey, Syrgkanis, and Zadik}{Mackey,
Syrgkanis and Zadik}{2018}]{MackeySyrgkanisZadik2018}
Mackey, L., V.~Syrgkanis, and I.~Zadik (2018).
\newblock Orthogonal machine learning: Power and limitations.
\newblock In {\em International Conference on Machine Learning}, pp.\
3375--3383. PMLR.
\bibitem[\protect\citeauthoryear{Magnus and Neudecker}{Magnus and
Neudecker}{1979}]{magnus1979commutation}
Magnus, J.~R. and H.~Neudecker (1979).
\newblock The commutation matrix: some properties and applications.
\newblock {\em Annals of Statistics\/}~{\em 7}, 381--394.
\bibitem[\protect\citeauthoryear{Magnus and Neudecker}{Magnus and
Neudecker}{1980}]{magnus1980elimination}
Magnus, J.~R. and H.~Neudecker (1980).
\newblock The elimination matrix: some lemmas and applications.
\newblock {\em SIAM Journal on Algebraic Discrete Methods\/}~{\em 1\/}(4),
422--449.
\bibitem[\protect\citeauthoryear{McLeish and Small}{McLeish and
Small}{1994}]{SmallMcLeish1994}
McLeish, D.~L. and C.~G. Small (1994).
\newblock {\em Hilbert Space Methods in Probability and Statistical Inference}.
\newblock Wiley NY.
\bibitem[\protect\citeauthoryear{Newey}{Newey}{1994}]{Newey1994}
Newey, W.~K. (1994).
\newblock The asymptotic variance of semiparametric estimators.
\newblock {\em Econometrica\/}~{\em 62}, 1349--1382.
\bibitem[\protect\citeauthoryear{Newey and Robins}{Newey and
Robins}{2017}]{NeweyRobins2017}
Newey, W.~K. and J.~M. Robins (2017).
\newblock Cross-fitting and fast remainder rates for semiparametric estimation.
\newblock Mimeo.
\bibitem[\protect\citeauthoryear{Neyman}{Neyman}{1959}]{Neyman1959}
Neyman, J. (1959).
\newblock Optimal asymptotic tests of composite hypotheses.
\newblock In U.~Grenander (Ed.), {\em Probability and Statistics}, pp.\
416--444. Wiley NY.
\bibitem[\protect\citeauthoryear{Neyman and Scott}{Neyman and
Scott}{1948}]{NeymanScott1948}
Neyman, J. and E.~L. Scott (1948).
\newblock Consistent estimates based on partially consistent observations.
\newblock {\em Econometrica\/}~{\em 16}, 1--32.
\bibitem[\protect\citeauthoryear{Robins, Li, Tchetgen~Tchetgen, and van~der
Vaart}{Robins, Li, Tchetgen~Tchetgen and van~der
Vaart}{2008}]{RobinsLiTchetgenTchetgenvanderVaart2008}
Robins, J., L.~Li, E.~Tchetgen~Tchetgen, and A.~van~der Vaart (2008).
\newblock Higher order influence functions and minimax estimation of nonlinear
functionals.
\newblock In {\em Probability and Statistics: Essays in Honor of David
A.~Freedman}, pp.\ 335--421. Institute of Mathematical Statistics.
\bibitem[\protect\citeauthoryear{Schick}{Schick}{1986}]{Schick1986}
Schick, A. (1986).
\newblock On asymptotically efficient estimation in semiparametric models.
\newblock {\em Annals of Statistics\/}~{\em 14}, 1139--1151.
\bibitem[\protect\citeauthoryear{Semenova, Goldman, Chernozhukov, and
Taddy}{Semenova, Goldman, Chernozhukov and
Taddy}{2023}]{semenova2023inference}
Semenova, V., M.~Goldman, V.~Chernozhukov, and M.~Taddy (2023).
\newblock Inference on heterogeneous treatment effects in high-dimensional
dynamic panels under weak dependence.
\newblock {\em Quantitative Economics\/}~{\em 14}, 471--510.
\bibitem[\protect\citeauthoryear{Small and McLeish}{Small and
McLeish}{1988}]{SmallMcLeish1988}
Small, C.~G. and D.~L. McLeish (1988).
\newblock Generalizations of ancillarity, completeness, and sufficiency in an
inference function space.
\newblock {\em Annals of Statistics\/}~{\em 16}, 534--551.
\bibitem[\protect\citeauthoryear{Small and McLeish}{Small and
McLeish}{1989}]{SmallMcLeish1989}
Small, C.~G. and D.~L. McLeish (1989).
\newblock Projection as a method for increasing sensitivity and eliminating
nuisance parameters.
\newblock {\em Biometrika\/}~{\em 76}, 693--703.
\bibitem[\protect\citeauthoryear{van~der Vaart}{van~der
Vaart}{2014}]{vanderVaart2014}
van~der Vaart, A. (2014).
\newblock Higher order tangent spaces and influence functions.
\newblock {\em Statistical Science\/}~{\em 29}, 679--686.
\bibitem[\protect\citeauthoryear{Waterman and Lindsay}{Waterman and
Lindsay}{1996}]{WatermanLindsay1996}
Waterman, R.~P. and B.~G. Lindsay (1996).
\newblock Projected score methods for approximating conditional scores.
\newblock {\em Biometrika\/}~{\em 83}, 1--13.
\bibitem[\protect\citeauthoryear{Woutersen}{Woutersen}{2002}]{woutersen2002robustness}
Woutersen, T. (2002).
\newblock Robustness against incidental parameters.
\newblock Technical report, Research Report.
\bibitem[\protect\citeauthoryear{W{\"u}thrich and Zhu}{W{\"u}thrich and
Zhu}{2021}]{WuthrichZhu2021}
W{\"u}thrich, K. and Y.~Zhu (2021).
\newblock Omitted variable bias of {L}asso-based inference methods: {A} finite
sample analysis.
\newblock Forthcoming in {\it Review of Economics and Statistics}.
\end{thebibliography}
\clearpage