Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
101,182 characters · 22 sections · 67 citation commands
A NEYMAN-ORTHOGONALIZATION APPROACH TO THE INCIDENTAL PARAMETER PROBLEM
\def\spacingset#1 \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\bf JEL Classification:} C13, C23, C55.
{\bf Keywords:} Neyman-orthogonality, incidental parameter, higher-order bias correction, networks.
\onehalfspacing
\setcounter{equation}{0}
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 Neyman1959. Such equations underlie much of the results in semiparametric estimation (Newey1994) and are at the heart of recent advances on doubly-robust estimation and high-dimensional inference as discussed in, for example, 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 WuthrichZhu2021 and AngristFrandsen2022. A related problem where the conventional approach was formally shown to fail is a nonlinear version of the judge-leniency design, see 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 (JochmansWeidner2019). Examples include the estimation of teacher value-added (JacksonRockoffStaiger2014), of the contributions of worker and firm heterogeneity to the variance of log wages and other covariance components (AbowdKramarzMargolis1999, KlineSaggioSoelvsten2020), as well as of complementarity patterns in team production (AhmadpoodJones2019, Bonhomme2021). Fixed effects in network-formation models are also poorly estimated, especially in the prevalent case where the network is sparse (see, e.g., graham2020sparse).
Motivated by these concerns, we focus on a higher-order generalization of Neyman-orthogonality that was proposed by 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 SmallMcLeish1989 and WatermanLindsay1996. They can also be seen as higher-order influence functions, as in RobinsLiTchetgenTchetgenvanderVaart2008 and 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 (HahnNewey2004, DhaeneJochmans2015b,DhaeneJochmans2015a) and to the literature on small measurement error (Chesher1991, 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 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 graham2014complementarity.
This problem is difficult because in the data that we use (taken from 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}
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
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 (ref) 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 NeymanScott1948 model. Here,
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$,
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 AhmadpoodJones2019),
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)) 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.
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
The (A) term in ((ref)) 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)), (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
Then the (C) term vanishes from ((ref)) and we obtain
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 (ref), Neyman1959 proposed a general method to transform it into one that does. The resulting function is said to be Neyman-orthogonal.
Condition (ref) has a long history in semiparametric estimation problems (Bickel1982, Schick1986, 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 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
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. }
To see how Neyman-orthogonality to a higher order can be helpful we now consider a further expansion of ((ref)). Again assuming sufficient regularity, we have, for any integer $q\geq 1$,
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
the (C) term is equal to zero, and so
Comparing ((ref)) to ((ref)) 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)) was introduced by 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 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
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}
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.
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
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(
\right) = \left(
\right) . $$
Neyman-orthogonality to order $q$ can now be defined as follows (MackeySyrgkanisZadik2018).
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.
Let $\mu_0$ satisfy ((ref)) 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). 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., 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., semenova2023inference).
We estimate $\mu_0$ by the GMM estimator
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
where $\widetilde W$ is again a chosen weight matrix. In this case, we may equally combine (ref) and (ref) into a single GMM estimation procedure.
In Section (ref), 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) in Section (ref)), 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), 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}
{
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). It is convenient to introduce the 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)}. $$ SmallMcLeish1994 discuss several properties of this basis. One important property for our purposes is that
for any $p$, so all elements of the Bhattacharyya basis have (conditional) mean equal to zero. In ((ref)), 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,
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(
\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., 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).
}
Theorem (ref) generalizes the projected-score construction of SmallMcLeish1989 and 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) 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 (SmallMcLeish1988) and to approximate the conditional score for $\theta$, when the latter exists (WatermanLindsay1996). It generalizes 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 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) 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 vanderVaart2014.
We note that Theorem (ref) 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)) 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), 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$.}
To gain intuition into the construction in Theorem (ref) 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
for some function $a_1$. Note that, by virtue of (ref), 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 (ref) 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
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 (ref) is well-known. Furthermore, it does not hinge on the likelihood structure. Indeed, recent work exploiting orthogonality, such as that surveyed in ChernozhukovChetverikovDemirerDufloHansenNeweyRobins2018, does so in the context of moment conditions. In our setup, as in Neyman1959's (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 (ref). 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
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(
\right)^\top \, \left(
\right). $$ Solving this equation for $a_{21}$ for given $a_{22}$ yields
where $a_1$ is given by ((ref)) 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)) back into ((ref)) yields
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
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)) is Neyman-orthogonal to order 2 when one sets $a_{21}$ to its expression in ((ref)), and $a_{22}$ to
Note that this construction amounts to solving a system of linear equations. The fact that the solution in ((ref))--((ref)) coincides with the expression in Theorem (ref) 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}
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 HahnNewey2004 and 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
Taking expectations and re-arranging shows that
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. HahnNewey2004, ArellanoHahn2007, and DhaeneJochmans2015b,DhaeneJochmans2015a present estimators of these terms based on the MLE that can be used to construct a bias-corrected estimator.
Lancaster2002 and woutersen2002robustness integrate-out the fixed effects using a uniform prior after orthogonalizing to order 1 to obtain an estimator with bias $o(T^{-1})$; 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. LiLindsayWaterman2003, building on 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 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)). 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). 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)).
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) is available in closed form, as
where $H_{k}$ denotes the $k$-th Hermite polynomial.
Interestingly, ((ref)) 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.
Our next example is the nonlinear regression model with $d\geq 1$ outcomes,
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) 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 (constantine1996multivariate). Given the matrix $M$ it is easy to compute $\varSigma_{w_qw_q}$, $\varSigma_{w_qu}$, and $b_q $. For example,
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). In the next section we present simulations and an empirical application based on a version of ((ref)) designed to study team production.
\paragraph{Example: CES production function (continued).}
Consider the team production model
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)) generalizes Model ((ref)) 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)) can be written as a special case of ((ref)), 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
which takes the same form as ((ref)), 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))--((ref)) for subset $i$, and $\eta_{i0}$ the $2\times 1$ vector of worker-specific effects in the corresponding teams.
\setcounter{equation}{0}
We wish to estimate the parameters of the team production model in (ref)--(ref). 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
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.
AhmadpoodJones2019 consider model ((ref)) 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 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 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 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))--((ref)). 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), 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)) 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) we assess the ability of our orthogonalization approach to alleviate this bias in a Monte Carlo simulation.
Table (ref) 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.}
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) in Appendix (ref) 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) 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), 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)) 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). 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) 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), using them directly for estimation may lead to very imprecise estimates. In contrast, Table (ref) 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) we report our diagnostic statistic for the orthogonality order $q$, described in Appendix (ref).\footnote{The statistic in Appendix (ref) 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.
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)). The first column in Table (ref) 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) 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)) 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) 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) 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).
\setcounter{equation}{0}
In this section, we show that, under higher-order orthogonality, the estimators $\widehat{\theta}$ and $\widehat{\mu}$ introduced in Section (ref) 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.
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$,
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) below.}
With this notation in hand, we now state our first assumption.
Part $(i)$ in Assumption (ref) is satisfied if $\widehat\mu_n$ is computed using GMM, see ((ref)). 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 (ref) 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}$.
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.
Part $(i)$ in Assumption (ref) requires $u$ to be Neyman-orthogonal in the sense of Definition (ref). 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) 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}$.
Note that, although we leave the dependence on $q$ implicit in Theorem (ref), 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.
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.