\section{Introduction}

Repeated measurements data consists of observations from several experimental units, subjects, or groups under different conditions. 
This grouping or clustering of the individual responses into experimental units 
typically introduces dependencies: 
the different units are assumed to be independent, 
but there may be heterogeneity between units 
and correlation within units.

Mixed-effects models 
provide a powerful and flexible tool to analyze grouped data by incorporating fixed and random effects. Fixed effects are associated with the entire population, and random effects are associated with individual 
groups
and model the heterogeneity across them and the dependence structure within them~\citep{Pinheiro2000}. 
Linear mixed-effects models~\citep{Laird-Ware1982, Pinheiro2000, Verbeke-Molenberghs2009, Demidenko2004} 
 impose a linear relationship between all covariates and the response. 
Partially linear mixed-effects models~\citep{Zeger-Diggle1994}  extend the linear ones. 

We consider the partially linear mixed-effects model 
\begin{equation}\label{eq:initPLMM}
	\mathbf{Y}_i = \Xi\beta_0 + g(\mathbf{W}_i) + \mathbf{Z}_i\mathbf{b}_i+ \boldsymbol{\varepsilon}_i
\end{equation}
for groups $i\in\{1,\ldots,N\}$. There are $n_i$ observations per group $i$. 
The unobserved random variable $\mathbf{b}_i$, called random effect, introduces 
correlation
within its group $i$ because all $n_i$ observations within this group share $\mathbf{b}_i$. 
We make the assumption generally made that both the random effect $\mathbf{b}_i$ and the error term $\boldsymbol{\varepsilon}_i$ 
follow a Gaussian distribution~\citep{Pinheiro2000}. 
The matrices $\mathbf{Z}_i$ assigning the random effects to group-level observations are fixed. 
The linear covariables $\Xi$ and the nonparametric and potentially high-dimensional covariables $\mathbf{W}_i$ are observed and random, and they may be dependent. Furthermore, the nonparametric covariables may contain nonlinear transformations and interaction terms of the linear ones.
Please see Assumption~\ref{assumpt:Distribution} in Section~\ref{sect:ModelAndDML} for further details. 

Our aim is to estimate and make inference for the so-called fixed effect $\beta_0$ in~\eqref{eq:initPLMM} in the presence of a highly complex $g$ using general machine learning 
algorithms. 
The parametric component $\beta_0$ 
provides a simple summary of the covariate effects that are of main scientific interest. The nonparametric component $g$ enhances model flexibility  because time trends and further covariates with possibly nonlinear and interaction effects can be modeled nonparametrically. 
\\

Repeated measurements, or longitudinal, data is omnipresent
in empirical research.
For example, assume we want to study the effect of a treatment over time. Observing the same subjects repeatedly presents three main advantages over having cross-sectional data. 
First, subjects can serve as their own controls.  
Second, the between-subject variability is explicitly modeled and can be excluded from the experimental error. This yields more efficient estimators of the relevant model parameters. 
Third, data can be collected more reliably~\citep{Davies2002, Fitzmaurice2011}.
\\

Various approaches have been considered in the literature to estimate the nonparametric component 
$g$ in~\eqref{eq:initPLMM}: 
kernel methods~\citep{Hart-Wehrly1986, Zeger-Diggle1994, Taavoni2019b, Chen2017},
backfitting~\citep{Zeger-Diggle1994, Taavoni2019b}, 
spline methods~\citep{Rice-Silverman1991, Zhang2004, Qin-Zhu2007, Qin-Zhu2009, Li-Zhu2010, Kim2017, Aniley2019}, 
and local linear regression~\citep{Taavoni2019b, Liang2009}. 

Our aim is to make inference for $\beta_0$ in the presence of potentially highly complex effects of $\mathbf{W}_i$ on $\Xi$ and $\mathbf{Y}_i$. First, we adjust $\Xi$ and $\mathbf{Y}_i$ for $\mathbf{W}_i$ by regressing $\mathbf{W}_i$ out of them using  machine learning algorithms. These machine learning algorithms may yield biased results, especially if regularization methods are used, like for instance with the lasso~\citep{Tibshirani1996}. 
Second, we fit a linear mixed-effects model to these regression residuals to estimate $\beta_0$. Our estimator of $\beta_0$ converges at the parametric rate, follows a Gaussian distribution asymptotically, and is semiparametrically efficient.

We rely on the double machine learning 
framework of~\citet{Chernozhukov2018} to estimate $\beta_0$ using general machine learning algorithms. 
To the best of our knowledge, this is the first work to allow the nonparametric nuisance components of a 
partially linear mixed-effects model to be estimated with arbitrary machine learners like random forests~\citep{Breiman2001} or the lasso~\citep{Tibshirani1996, Buehlmann2011}. 
In contrast to the setting and proofs of~\citet{Chernozhukov2018}, we have dependent data and need to incorporate this accordingly. 

\citet{Chernozhukov2018} introduce double machine learning to estimate a low-dimensional parameter in the presence of nonparametric nuisance components using machine learning 
algorithms. 
This estimator converges at the parametric rate and is asymptotically Gaussian due to Neyman orthogonality and sample splitting with cross-fitting. 
We would like to remark that nonparametric components can be estimated without sample splitting and cross-fitting if the underlying function class satisfies some entropy conditions; 
see for instance~\citet{Geer-Mammen1997}. However, these regularity conditions limit the complexity of the function class, and 
machine learning algorithms usually do not satisfy them. Particularly, these conditions fail to hold if the dimension of the nonparametric variables increases with the sample size~\citep{Chernozhukov2018}. 


\subsection{Additional Literature}

Expositions and overviews of mixed-effects modeling techniques can be found in~\citet{Pinheiro1994, Davidian-Giltinan1995, Vonesh-Chinchilli1997, Pinheiro2000, Davidian-Giltinan2003}.
\\

\citet{Zhang1998} consider partially linear
mixed-effects models and estimate the nonparametric component with natural cubic splines. They treat the smoothing parameter as an extra variance component that is jointly estimated with the other variance components of the model. 
\citet{Masci2019} consider partially linear
mixed-effects models for unsupervised classification with discrete random effects. 
\citet{Schelldorfer2011} consider  high-dimensional linear mixed-effects models where the number of fixed effects coefficients may be much larger than the overall sample size. 
\citet{Taavoni2019a} employ a regularization approach in generalized 
partially linear mixed-effects models using regression splines to approximate the nonparametric component. 
\citet{gamm4} use penalized regression splines where the penalized components are treated as random effects.
\\

The unobserved random variables in the partially linear mixed-effects model~\eqref{eq:initPLMM} 
are assumed to follow a Gaussian distribution. 
\citet{Taavoni2021} introduce multivariate $t$ 
partially linear
mixed-effects models for longitudinal data. They consider $t$-distributed random effects to account for outliers in the data. 
\citet[Chapter 4]{Fahrmeir2011} relax the assumption of Gaussian random effects in generalized linear mixed models. They consider nonparametric Dirichlet processes and Dirichlet process mixture priors for the random effects.
\citet[Chapter 3]{Ohinata2012} consider partially linear mixed-effects models and make no distributional assumptions for the random terms, 
and the nonparametric component is estimated with kernel methods.
\citet{Lu2016} consider a partially linear mixed-effects model that is nonparametric in time and that features asymmetrically distributed errors  and missing data.

Furthermore, methods have been developed to analyze repeated measurements data that are robust to outliers. 
\citet{Guoyou2008} consider robust estimating equations and estimate the nonparametric component with a regression spline.
\citet{Tang2015} consider median-based regression methods in a partially linear model with longitudinal data to account for highly skewed responses. 
\citet{Lin2018} present an estimation technique in partially linear models for longitudinal data that is doubly robust in the sense that  it simultaneously accounts for missing responses and mismeasured covariates. 
\\

It is prespecified in 
the  partially linear mixed-effects model~\eqref{eq:initPLMM} which covariates are modeled with random effects. 
Simultaneous variable selection for fixed effects variables and random effects has been developed by~\citet{Bondell2010, Ibrahim2011}. They use penalized likelihood approaches.  
\citet{Li-Zhu2010} use a nonparametric test to test the existence of random effects in partially linear mixed-effects models. 
\citet{Zhang2020} propose a variable selection procedure for the linear covariates of a generalized partially linear model with longitudinal data. 
\\

\textit{Outline of the Paper.}
Section~\ref{sect:ModelAndDML} presents our double machine learning estimator of the linear coefficient in a partially linear mixed-effects model. Section~\ref{sect:numerical-experiments} presents our numerical results. 
\\

\textit{Notation.}
We denote by $\indset{N}$ the set $\{1,2,\ldots,N\}$. We add the probability law as a subscript to the probability operator $\Prob$ and the expectation operator $\E$ whenever we want to emphasize the corresponding dependence.
We denote the $L^p(P)$ norm by $\normP{\cdot}{p}$ and the Euclidean or  operator
norm  by $\norm{\cdot}$, depending on the context. 
We implicitly assume that given expectations and conditional expectations exist. We denote by $\stackrel{d}{\rightarrow}$ convergence in distribution.
The symbol $\independent$ denotes independence of random variables.
We denote by $\mathds{1}_{n}$ the $n\times n$ identity matrix and omit the subscript $n$ if we do not want to emphasize the dimension. 
We denote the $d$-variate Gaussian distribution by $\mathcal{N}_d$.


\section{Model Formulation and the Double Machine Learning Estimator}\label{sect:ModelAndDML}

We consider repeated measurements data that is grouped according to experimental units or subjects. 
This grouping structure introduces dependency in the data. The individual experimental units or groups are assumed to be independent, but there may be some 
between-group heterogeneity and within-group correlation.
We consider the partially linear mixed-effects model 
\begin{equation}\label{eq:PLMM}
	\mathbf{Y}_i = \Xi\beta_0 + g(\mathbf{W}_i) + \mathbf{Z}_i\mathbf{b}_i+ \boldsymbol{\varepsilon}_i, \quad i\in\indset{N}
\end{equation}
for groups $i$ as in~\eqref{eq:initPLMM} to model the 
between-group heterogeneity and within-group correlation with random effects.
We have $n_i$ observations per group that are concatenated row-wise into $\mathbf{Y}_i\in\R^{n_i}$, $\Xi\in\R^{n_i\times d}$, and $\mathbf{W}_i\in\R^{n_i\times v}$. 
The nonparametric  variables may be high-dimensional, but $d$ is fixed.
Both $\Xi$ and $\mathbf{W}_i$ are random. 
The  $\Xi$ and  $\mathbf{W}_i$ belonging to the same group $i$ may be dependent.
For groups $i\neq j$, we assume $\Xi\independent\mathbf{X}_j$, $\mathbf{W}_i\independent\mathbf{W}_j$, and $\Xi\independent\mathbf{W}_j$.
We assume that $\mathbf{Z}_i\in\R^{n_i\times q}$ is fixed. 
The random variable $\mathbf{b}_i\in\R^q$ denotes a group-specific vector of random regression coefficients that is assumed to follow a Gaussian distribution. 
The dimension $q$ of the random effects model is fixed.
Also the error terms are assumed to follow a Gaussian distribution as is commonly used 
in a mixed-effects models framework~\citep{Pinheiro2000}. 
All groups $i$ share the common linear coefficient $\beta_0$ and the potentially complex function $g\colon\R^{v}\rightarrow\R$. The function $g$ is applied row-wise to $\mathbf{W}_i$, denoted by $g(\mathbf{W}_i)$.

We denote the total number of observations by $\NN_T:=\sum_{i=1}^{N}n_i$. 
We assume that the numbers $n_i$ of within-group observations are uniformly upper bounded by $n_{\mathrm{max}}<\infty$. 
Asymptotically, the number of groups, $N$,  goes to infinity. 

Our distributional and independency assumptions are summarized as follows:

\begin{assumptions}\label{assumpt:Distribution} Consider the partially linear mixed-effects model~\eqref{eq:PLMM}. We assume that there is some $\sigma_0>0$ and some symmetric positive definite matrix 
$\Gamma_{0}\in\R^{q\times q}$ 
such that the following conditions hold.
\begin{enumerate}[label={\theassumptions.\arabic*}]
	\item \label{assumpt:D1} 
		The random effects $\mathbf{b}_1,\ldots,\mathbf{b}_{\NN}$ are independent and identically distributed\ $\mathcal{N}_q(\boldsymbol{0}, \Gamma_{0})$.
	\item \label{assumpt:D2}
		The error terms $\boldsymbol{\varepsilon}_1,\ldots,\boldsymbol{\varepsilon}_{\NN}$ are independent
		and follow a Gaussian distribution, $\boldsymbol{\varepsilon}_i\sim\mathcal{N}_{n_i}(\boldsymbol{0}, \sigma_0^2\mathds{1}_{n_i})$ for $i\in\indset{N}$, with the common variance component $\sigma_0^2$. 
	\item \label{assumpt:D3}
		The variables $\mathbf{b}_1,\ldots,\mathbf{b}_{\NN}, \boldsymbol{\varepsilon}_1,\ldots,\boldsymbol{\varepsilon}_{\NN}$ are independent. 
	\item \label{assumpt:D4-2}
		For all $i, j \in\indset{N}$, $i\neq j$, we have $(\mathbf{b}_i,\boldsymbol{\varepsilon}_i)\independent(\mathbf{W}_i, \Xi)$ and $(\mathbf{b}_i,\boldsymbol{\varepsilon}_i)\independent(\mathbf{W}_j, \mathbf{X}_j)$.
	\item \label{assumpt:D6}
	         For all $i, j \in\indset{N}$, $i\neq j$, we have $\Xi\independent\mathbf{X}_j$, $\mathbf{W}_i\independent\mathbf{W}_j$, and $\Xi\independent\mathbf{W}_j$.
\end{enumerate}
\end{assumptions}

We would like to remark that the distribution of the error terms $\boldsymbol{\varepsilon}_i$ in Assumption~\ref{assumpt:D2} can be generalized to $\boldsymbol{\varepsilon}_i\sim\mathcal{N}_{n_i}(\boldsymbol{0},\sigma_0^2\Lambda_i(\boldsymbol{\lambda}))$, where $\Lambda_i(\boldsymbol{\lambda})\in\R^{n_i\timesn_i}$ is a symmetric positive definite matrix parametrized by some finite-dimensional parameter vector $\boldsymbol{\lambda}$ that all groups have in common. 
For the sake of notational simplicity, we restrict ourselves to Assumption~\ref{assumpt:D2}. 

Moreover, we may consider stochastic random effects matrices $\mathbf{Z}_i$. Alternatively, the nonparametric variables $\mathbf{W}_i$ may be part of the random effects matrix. In this case, 
we consider 
the random effects matrix $\widetilde{\mathbf{Z}}_i = \zeta(\mathbf{Z}_i,\mathbf{W}_i)$ for some known function $\zeta$ in~\eqref{eq:PLMM} instead of $\mathbf{Z}_i$.  
Please see Section~\ref{sec:nonfixedZi} in the appendix for further details. 
For simplicity, we restrict ourselves to fixed random effects matrices $\mathbf{Z}_i$ that are disjoint from $\mathbf{W}_i$. 
\\

The unknown parameters in our model are $\beta_0$, $\Gamma_{0}$, and $\sigma_0$. 
Our aim is to estimate 
$\beta_0$ and 
make inference for it.
Although the variance parameters $\Gamma_{0}$ and $\sigma_0$ need to be estimated consistently to construct an estimator of $\beta_0$, it is not our goal to perform inference for them.

\subsection{The Double Machine Learning Fixed-Effects Estimator}

Subsequently, we describe our estimator of $\beta_0$ in~\eqref{eq:PLMM}. 
To motivate our procedure, we first consider the population version with the residual terms
\begin{displaymath}
	\mathbf{R}_{\Xi} := \Xi - \E[\Xi | \mathbf{W}_i] \quad\textrm{and}\quad
	\mathbf{R}_{\Yi} := \mathbf{Y}_i - \E[\mathbf{Y}_i | \mathbf{W}_i] \quad\textrm{for}\quad
	i\in\indset{N}
\end{displaymath}
that adjust $\Xi$ and $\mathbf{Y}_i$ for $\mathbf{W}_i$. On this adjusted level, we have the linear mixed-effects model
\begin{equation}\label{eq:LMM}
	\mathbf{R}_{\Yi} = \mathbf{R}_{\Xi} \beta_0 + \mathbf{Z}_i\mathbf{b}_i + \boldsymbol{\varepsilon}_i, \quad i\in\indset{N}
\end{equation}
due to~\eqref{eq:PLMM} and Assumption~\ref{assumpt:D4-2}. 
In particular, the adjusted and grouped responses in this model are independent in the sense that we have $\mathbf{R}_{\Yi} \independent \mathbf{R}_{\Yj}$ for $i\neq j$. 
The strategy now is to first estimate the residuals with machine learning algorithms and then use 
linear mixed model techniques to infer $\beta_0$. This is done with sample splitting and cross-fitting, and the details are described next.
\\

Let us define $\Sigma_{0} := \sigma_0^{-2}\Gamma_{0}$ 
and $\mathbf{V}_{0, i}:=(\mathbf{Z}_i\Sigma_{0}\mathbf{Z}_i^T + \mathds{1}_{n_i})$ so that we have
\begin{equation}\label{eq:resNormal}
	(\mathbf{R}_{\Yi} | \mathbf{W}_i,\Xi) \sim \mathcal{N}_{n_i}\big(\mathbf{R}_{\Xi}\beta_0, \sigma_0^2\mathbf{V}_{0, i}\big). 
\end{equation}
We assume that there exist functions $m_X^0\colon\R^{v}\rightarrow\R^d$ and $m_Y^0\colon\R^{v}\rightarrow\R$ that we can apply row-wise to $\mathbf{W}_i$ to have
$\E[\Xi | \mathbf{W}_i] = m_X^0(\mathbf{W}_i)$ and $\E[\mathbf{Y}_i | \mathbf{W}_i] = m_Y^0(\mathbf{W}_i)$. 
In particular, $m_X^0$ and $m_Y^0$ do not depend on the grouping index $i$. 
Let $\eta^0:= (m_X^0,m_Y^0)$ denote the true unknown nuisance parameter. 
Let us denote by $\theta_0 := (\beta_0, \sigma_0^2,\Sigma_{0})$ the complete true unknown parameter vector and by $\theta := (\beta,\sigma^2,\Sigma)$ and $\mathbf{V}_{i}:=\mathbf{Z}_i\Sigma\mathbf{Z}_i^T+\mathds{1}_{n_i}$ respective general parameters. 
The log-likelihood of group $i$ is given by 
\begin{equation}\label{eq:log-likelihood}
\begin{array}{rcl}
	\libig{\theta, \eta^0} &=& -\frac{n_i}{2} \log(2\pi) - \frac{n_i}{2}\log(\sigma^2) - \frac{1}{2}\log\big( \det(\mathbf{V}_{i})\big) \\
	&&\quad- \frac{1}{2\sigma^2} (\mathbf{R}_{\Yi}-\mathbf{R}_{\Xi}\beta)^T\mathbf{V}_{i}^{-1}(\mathbf{R}_{\Yi}-\mathbf{R}_{\Xi}\beta)  - \log \big(p(\mathbf{W}_i,\Xi)\big),
	\end{array}
\end{equation}
where $p(\mathbf{W}_i,\Xi)$ denotes the joint density of $\mathbf{W}_i$ and $\Xi$. We assume that $p(\mathbf{W}_i,\Xi)$ does not depend on $\theta$. 
The true nuisance parameter $\eta^0$ in the log-likelihood~\eqref{eq:log-likelihood} is unknown and estimated with machine learning algorithms (see below).
Denote by $\eta := (m_X,m_Y)$ some general nuisance parameter. The terms that adjust $\Xi$ and $\mathbf{Y}_i$ for $\mathbf{W}_i$ with this general nuisance parameter are given by $\Xi-m_X(\mathbf{W}_i)$ and $\mathbf{Y}_i-m_Y(\mathbf{W}_i)$. Up to additive constants that do not depend on $\theta$ and $\eta$, we thus consider maximum likelihood estimation with the likelihood
\begin{displaymath}
	\begin{array}{rl}
	&\ell_i(\theta,\eta) 
	= -\frac{n_i}{2}\log(\sigma^2) - \frac{1}{2}\log\big(\det(\mathbf{V}_{i}) \big)\\
	&\quad- \frac{1}{2\sigma^2} \Big(\mathbf{Y}_i-m_Y(\mathbf{W}_i)-\big(\Xi-m_X(\mathbf{W}_i)\big)\beta\Big)^T\mathbf{V}_{i}^{-1}\Big(\mathbf{Y}_i-m_Y(\mathbf{W}_i)-\big(\Xi-m_X(\mathbf{W}_i)\big)\beta\Big),
	\end{array}
\end{displaymath}
which is a function of both the finite-dimensional parameter $\theta$ and the infinite-dimensional
nuisance parameter $\eta$. 
\\

Our estimator of $\beta_0$ is constructed as follows using double machine learning. 
First, we estimate $\eta^0$ with machine learning 
algorithms and plug these estimators into the estimating equations for $\theta_0$, equation~\eqref{eq:setToZero} below, 
to obtain an estimator for $\beta_0$. This procedure is done with sample splitting and cross-fitting as explained next.

Consider repeated measurements from $N$ experimental units, subjects, or groups as in~\eqref{eq:PLMM}. Denote by $\mathbf{S}_i :=(\mathbf{W}_i,\Xi,\mathbf{Z}_i,\mathbf{Y}_i)$ the observations of group $i$.
First, we split the group indices $\indset{N}$ into $K\ge 2$ disjoint sets $I_1, \ldots, I_{K}$ of approximately equal size; 
please see Section~\ref{sect:AssumptionsDefinitions} in the appendix for further details. 

For each $k\in\indset{K}$, we estimate the conditional expectations $m_X^0(W)$ and $m_Y^0(W)$ with data from $I_{\kk}^c$. We call the resulting estimators  $\hat m_X^{I_{\kk}^c}$ and $\hat m_Y^{I_{\kk}^c}$, respectively. 
Then, the adjustments $\widehat{\mathbf{R}}_{\Xi}^{I_{\kk}}:=\Xi-\hat m_X^{I_{\kk}^c}(\mathbf{W}_i)$, and $\widehat{\mathbf{R}}_{\mathbf{Y}_i}^{I_{\kk}}:=\mathbf{Y}_i-\hat m_Y^{I_{\kk}^c}(\mathbf{W}_i)$ for $i\inI_{\kk}$ are evaluated on $I_{\kk}$, the complement of $I_{\kk}^c$. 
Let $\hat\eta^{I_{\kk}^c} := (\hat m_X^{I_{\kk}^c}, \hat m_Y^{I_{\kk}^c})$ denote the estimated nuisance parameter.
Consider the score function
	$\psi(\mathbf{S}_i; \theta,\eta) := \nabla_{\theta}\ell_i(\theta,\eta)$, 
where $\nabla_{\theta}$ denotes the gradient with respect to $\theta$ interpreted as a vector.
On each set $I_{\kk}$, we consider an estimator $\hat\theta_{\kk} = (\hat\beta_{\kk},\hat\sigma_{\kk}^2, \hat{\Sigma}_{\kk})$ of $\theta_0$ that, approximately, in the sense of Assumption~\ref{assumpt:Theta2} in the appendix, solves
\begin{equation}\label{eq:setToZero}
	\frac{1}{\nn_{T, k}}\sum_{i\inI_{\kk}} \psibig{\mathbf{S}_i;\hat\theta_{\kk},\hat\eta^{I_{\kk}^c}} 
	= \frac{1}{\nn_{T, k}}\sum_{i\inI_{\kk}} \nabla_{\theta}\ell_i(\theta,\eta)
	\stackrel{!}{=}  \boldsymbol{0},
\end{equation}
where $\nn_{T, k} := \sum_{i\inI_{\kk}}n_i$ denotes the total number of observations from experimental units that belong to the set $I_{\kk}$.
These $K$ estimators $\hat\theta_{\kk}$ for $k\in\indset{K}$ are assembled to form the final cross-fitting estimator
\begin{equation}\label{eq:betahat}
	\hat\beta := \frac{1}{K} \sum_{k=1}^{K} \hat\beta_{\kk}
\end{equation}
of $\beta_0$. 
We remark that one can simply use linear mixed model computation and software to compute $\hat\beta_{\kk}$ based on the estimated residuals $\widehat{\mathbf{R}}^{I_{\kk}}$. 
The estimator $\hat\beta$ fundamentally depends on the particular sample split. To alleviate this effect, the overall procedure may be repeated $\mathcal{S}$ times~\citep{Chernozhukov2018}. The $\mathcal{S}$ point estimators are aggregated by the median, and an additional term accounting for the random splits is added to the variance estimator of $\hat\beta$; please see 
Algorithm~\ref{algo:Summary} that presents the complete procedure.

\begin{algorithm}[h!]
	\SetKwInOut{Input}{Input}
   	\SetKwInOut{Output}{Output}


 \Input{$N$ independent and identically distributed\  observations $\{\mathbf{S}_i=(\mathbf{W}_i,\Xi, \mathbf{Z}_i, \mathbf{Y}_i)\}_{i\in\indset{N}}$ from the model~\eqref{eq:PLMM} satisfying Assumption~\ref{assumpt:Distribution}, a natural number $K$, a natural number $\mathcal{S}$.}
 \Output{An estimator of $\beta_0$ in~\eqref{eq:PLMM} together with its estimated asymptotic variance.}
 
 \For{$s\in\indset{\mathcal{S}}$}
 {
 Split the grouped observation index set $\indset{N}$ into $K$ sets $I_1, \ldots, I_{K}$ of approximately equal size.
 
 \For{$k\inK$}
 {
 Compute the conditional expectation estimators $\hat m_X^{I_{\kk}^c}$ and $\hat m_Y^{I_{\kk}^c}$ with some  
 machine learning algorithm and data from $I_{\kk}^c$. 
 
 Evaluate the adjustments $\widehat{\mathbf{R}}_{\Xi}^{I_{\kk}}=\Xi-\hat m_X^{I_{\kk}^c}(\mathbf{W}_i)$ and $\widehat{\mathbf{R}}_{\mathbf{Y}_i}^{I_{\kk}}=\mathbf{Y}_i-\hat m_Y^{I_{\kk}^c}(\mathbf{W}_i)$ for $i\inI_{\kk}$. 
 
Compute $\hat\theta_{\kk, s} = (\hat\beta_{\kk, s}, \hat\sigma_{\kk, s}^2, \hat\Sigma_{\kk, s})$ using, for instance, linear mixed model techniques. 
 }

Compute $\hat\beta_{s} = \frac{1}{K}\sum_{k=1}^{K}\hat\beta_{\kk, s}$ as an approximate solution to~\eqref{eq:setToZero}.

Compute an estimate $\hat T_{0,s}$ of the asymptotic variance-covariance matrix $T_0$ in Theorem~\ref{thm:asymptoticGaussian}. 
 }
 
 Compute $\hat\beta = \mathrm{median}_{s\in\indset{\mathcal{S}}}(\hat\beta_{s})$.
 
 Estimate $T_0$ by $\hat T_0 = \mathrm{median}_{s\in\indset{\mathcal{S}}}(\hat T_{0,s} + (\hat\beta-\hat\beta_{s})(\hat\beta-\hat\beta_{s})^T)$.\label{algo:varianceCorrection}

 \caption{Double machine learning in a partially linear mixed-effects model with repeated measurements.}\label{algo:Summary}
\end{algorithm}


\subsection{Theoretical Properties of the Fixed-Effects Estimator}

The estimator $\hat\beta$ as in~\eqref{eq:betahat} converges at the parametric rate, $N^{-1/2}$, and is asymptotically Gaussian distributed and semiparametrically efficient. 

\begin{theorem}\label{thm:asymptoticGaussian}
Consider grouped observations $\{\mathbf{S}_i=(\mathbf{W}_i,\Xi,\mathbf{Y}_i)\}_{i\in\indset{N}}$ from the partially linear mixed-effects model~\eqref{eq:PLMM} that satisfy Assumption~\ref{assumpt:Distribution} such that $p(\mathbf{W}_i,\Xi)$ does not depend on $\theta$. 
Let $\NN_T:=\sum_{i=1}^{N}n_i$ denote the total number of unit-level observations. 
Furthermore, suppose the assumptions in Section~\ref{sect:AssumptionsDefinitions} in the appendix hold,
and consider the symmetric positive-definite matrix $T_0$ given in Assumption~\ref{assumpt:regularity7} in the appendix. 
Then, $\hat\beta$ as in~\eqref{eq:betahat} concentrates in a $1/\surd{\NN_{T}}$ neighborhood of $\beta_0$ and is centered Gaussian, namely
\begin{equation}\label{eq:thmEquation}
	\surd{\NN_{T}}T_0^{\frac{1}{2}}(\hat\beta-\beta_0)
	\stackrel{d}{\rightarrow}\mathcal{N}_d(\boldsymbol{0}, \mathds{1}_{d}) \quad (N\rightarrow\infty), 
\end{equation}
and semiparametrically efficient. 
The convergence in~\eqref{eq:thmEquation} is in fact uniformly over the law $P$ of $\{\mathbf{S}_i=(\mathbf{W}_i,\Xi,\mathbf{Y}_i)\}_{i\in\indset{N}}$. 
\end{theorem}

Please see Section~\ref{sect:asymptoticDistribution} in the appendix for a proof of Theorem~\ref{thm:asymptoticGaussian}. 
Our proof builds on~\citet{Chernozhukov2018}, but we have to take into account the correlation within units that is introduced by the random effects.

The inverse asymptotic variance-covariance matrix $T_0$ can be consistently estimated;
see Lemma~\ref{lem:multiplierMatrix} in the appendix. 
The estimator $\hat\beta$ is semiparametrically efficient because the score function comes from the log-likelihood of our data and because $\eta^0$ solves a concentrating-out equation for fixed $\theta$; see~\citet{Chernozhukov2018, Newey1994b}. 

The assumptions in Section~\ref{sect:AssumptionsDefinitions} of the appendix specify regularity conditions and  required 
convergence rates of the machine learning estimators. The 
machine learning errors need to satisfy the product relationship 
\begin{displaymath}
	\normP{m_X^0(W)-\hat m_X^{I_{\kk}^c}(W)}{2}\big(\normP{m_Y^0(W)-\hat m_Y^{I_{\kk}^c}(W)}{2}+ \normP{m_X^0(W)-\hat m_X^{I_{\kk}^c}(W)}{2} \big)\llN^{-\frac{1}{2}}.
\end{displaymath}
This bound requires that only the products of the  
machine learning estimation errors 
$\normP{m_X^0(W)-\hat m_X^{I_{\kk}^c}(W)}{2}$ and $\normP{m_Y^0(W)-\hat m_Y^{I_{\kk}^c}(W)}{2}$
but not the individual ones need to vanish at a rate smaller than
$N^{-1/2}$. In particular, the individual estimation errors  may vanish at the  rate 
smaller than $N^{-1/4}$. This is achieved by many machine learning  
methods (cf. 
\citet{Chernozhukov2018}):
$\ell_1$-penalized and related methods in a variety of sparse models
\citep{Bickel2009, Buehlmann2011, Belloni2011, Belloni-Chernozhukov2011, Belloni2012, Belloni-Chernozhukov2013}, forward selection in sparse models
\citep{Kozbur2020}, $L_2$-boosting in sparse linear models
\citep{Luo2016}, a class of regression trees and random forests 
\citep{Wager2016}, and neural networks \citep{Chen1999}. 

We note that so-called Neyman orthogonality makes score functions insensitive to inserting potentially biased machine learning 
estimators of the nuisance parameters. A score function is Neyman orthogonal if its Gateaux derivative vanishes at the true $\theta_0$ and the true $\eta^0$. In particular, Neyman orthogonality is a first-order property. The product relationship of the machine learning 
estimating errors described above is used to bound second-order terms. We refer to Section~\ref{sect:asymptoticDistribution} in the appendix for more details. 


\section{Numerical Experiments}\label{sect:numerical-experiments}

We apply our method to an empirical and a pseudorandom dataset and in a simulation study. 
Our implementation is available in the \textsf{R}-package \texttt{dmlalg}~\citep{dmlalg}. 


\subsection{Empirical Analysis: CD4 Cell Count Data}

Subsequently, we apply our method to longitudinal CD4 cell counts data collected from human immunodeficiency virus (HIV) seroconverters. This data has previously been analyzed by~\citet{Zeger-Diggle1994} and is available in the \textsf{R}-package \texttt{jmcm}~\citep{jmcm} as \texttt{aids}. 
It contains $2376$ observations of CD4 cell counts measured on $369$ subjects. The data was collected during a period ranging from $3$ years before to $6$ years after seroconversion. The number of observations per subject ranges from $1$ to $12$, but for most subjects,  $4$ to $10$ observations are available. 
Please see~\citet{Zeger-Diggle1994} for more details on this dataset. 

Apart from time, five other covariates are measured: the age at seroconversion in years (age), the smoking status measured by the number of cigarette packs consumed per day (smoking), a binary variable indicating drug use (drugs), the number of sex partners (sex), and the depression status measured on the Center for Epidemiologic Studies Depression (CESD) scale (cesd), where higher CESD values indicate the presence of more depression symptoms. 

We incorporate a random intercept per person. 
Furthermore, we consider a square-root transformation of the CD4 cell counts to reduce the skewness of this variable as proposed by~\citet{Zeger-Diggle1994}. 
The CD4 counts are our response. The covariates that are of scientific interest are considered as $X$'s, 
and the remaining covariates are considered as $W$'s in the partially linear mixed-effects model~\eqref{eq:PLMM}.
The effect of time is modeled nonparametrically, but there are several options to model the other covariates. 
Other models than partially linear mixed-effects model
have also been considered in the literature to analyze this dataset. 
For instance, \citet{Fan2000} consider a functional linear model where the linear coefficients are a function of the time. 
\\

We consider two partially linear mixed-effects models for this dataset.
First, we incorporate all covariates except time linearly. 
Most approaches in the literature considering a partially linear mixed-effects model for this data that model time nonparametrically
report that sex and cesd are significant and that either smoking or drugs is significant as well; see for instance~\citet{Zeger-Diggle1994, Taavoni2019b, Wang2011}. 
\citet{Guoyou2008} develop a robust estimation method 
for longitudinal data and estimate 
nonlinear effects from time with regression splines. 
With the CD4 dataset, 
They find that smoking and cesd are significant. 

We apply our method with $K=2$ sample splits, $\mathcal{S}=100$ repetitions of splitting the data, and learn the conditional expectations with random forests that consist of $500$ trees whose minimal node size is $5$. 

Like~\citet{Guoyou2008},
we conclude that smoking and cesd are significant; please see the first row of Table~\ref{tab:empirical} for a more precise account of our findings.
Therefore,
we can expect
that our method implicitly performs robust estimation. Apart from sex, our point estimators are larger or of about the same size in absolute value as what~\citet{Guoyou2008} obtain. This suggests that our method incorporates potentially less bias. However, apart from age, the standard deviations are slightly larger with our method. 
This can be expected because random forests are more complex than the regression splines~\citet{Guoyou2008} employ. 
\\

\begingroup
\begin{table}
\footnotesize
\centering
\begin{tabular}{| L{2.5cm} | C{1.9cm}  | C{1.9cm} | C{1.9cm} | C{1.9cm} | C{2.2cm} |}
\hline
&  \textbf{age} & \textbf{smoking} & \textbf{drugs} & \textbf{sex} & \textbf{cesd}\\
\hline
\hline
$W = (\mathrm{time})$ &  $0.004$ ($0.027$) & $0.752$ ($0.123$) & $0.704$ ($0.360$) & $0.001$ ($0.043$)    &   $-0.042$ ($0.015$)\\
\hline
$W = (\mathrm{time}, \mathrm{age}, \mathrm{sex})$ & - & $0.620$ ($0.126$) & $0.602$ ($0.335$) &- & $-0.047$ ($0.015$) \\
\hline
\hline
\citet{Zeger-Diggle1994} 
& $0.037$ ($0.18$) & $0.27$ ($0.15$) & $0.37$ ($0.31$) & $0.10$ ($0.038$) &  $-0.058$ ($0.015$) \\
\hline
\citet{Taavoni2019b} 
& $1.5\cdot 10^{-17}$ ($3.5\cdot 10^{-17}$)& $0.152$ ($0.208$) & $0.130$ ($0.071$) &  $0.0184$ ($0.0039$) & $-0.0141$ ($0.0061$)  \\
\hline
\citet{Wang2011} 
& $0.010$ ($0.033$) & $0.549$ ($0.144$) & $0.584$ ($0.331$) & $0.080$ ($0.038$) & $-0.045$ ($0.013$)\\
\hline
\citet{Guoyou2008} 
& $0.006$ ($0.038$) &  $0.538$ ($0.136$) & $0.637$ ($0.350$)& $0.066$ ($0.040$) & $-0.042$ ($0.015$)\\
\hline
\end{tabular}
\caption{\label{tab:empirical}Estimates of the linear coefficient and its standard deviation in parentheses with our method 
for nonparametrically adjusting for time (first row) and for time, age, and sex (second row).
The remaining rows display the results from~\citet[Section 5]{Zeger-Diggle1994}, \citet[Table 1, ``Kernel'']{Taavoni2019b}, \citet[Table 2, ``Semiparametric efficient scenario I'']{Wang2011}, and \citet[Table 5, ``Robust'']{Guoyou2008}, respectively. 
}
\end{table}
\endgroup

We consider a second estimation approach where we model the variables time, age, and sex nonparametrically and allow them to interact. 
It is conceivable that these variables are not (causally) influenced by smoking, drugs, and cesd and that they are therefore exogenous. 
The variables smoking, drugs, and cesd are modeled linearly,  
and they are considered as treatment variables. Some direct causal effect interpretations are possible if one is willing to assume, for instance, that the nonparametric adjustment variables are causal parents of the linear variables or the response. 
However, we do not pursue this line of thought further.

We estimate the conditional expectations given the three nonparametric variables time, age, and sex again with random forests that consist of $500$ trees whose minimal node size is $5$ and use $K=2$ and $\mathcal{S}=100$ in Algorithm~\ref{algo:Summary}. 
We again find that smoking and cesd are significant; please see the second row of Table~\ref{tab:empirical}. 
This cannot be expected a priori because this second model incorporates more complex adjustments, which can lead to less significant variables. 


\subsection{Pseudorandom Simulation Study: CD4 Cell Count Data}\label{sect:pseudorandom}

Second, we consider the CD4 cell count data from the previous subsection and perform a pseudorandom simulation study. The variables smoking, drugs, and cesd are modeled linearly and the variables time, age, and sex nonparametrically. 
We condition on these six variables in our simulation. That is, they are the same in all repetitions. 
The function $g$ in~\eqref{eq:PLMM} is chosen as a regression tree that we built beforehand. 
We let $\beta_0=(0.62, 0.6, -0.05)^T$, where the first component corresponds to smoking, the second one to drugs, and the last one to cesd, consider a  standard deviation of the random intercept per subject of $4.36$, and a standard deviation of the error term of $4.35$. These are the point estimates of the respective quantities obtained in the previous subsection. 

Our fitting procedure uses random forests consisting of $500$ trees whose minimal node size is $5$ to estimate the conditional expectations, and we use $K=2$ and $\mathcal{S}=10$ in Algorithm~\ref{algo:Summary}. We perform $5000$ simulation runs.
We compare the performance of our method with that of the spline-based function \texttt{gamm4} from the package~\texttt{gamm4}~\citep{gamm4} for the statistical software \textsf{R}~\citep{R}.
This method represents the nonlinear part of the model by smooth additive functions and estimates them by penalized regression splines. The penalized components are treated as random effects and the unpenalized components 
as fixed. 

The results are displayed in Figure~\ref{fig:pseudorandom}. 
With our method,  \texttt{mmdml}, the two-sided confidence intervals for $\beta_0$ are of about the same length  but achieve a coverage that is closer to the nominal $95\%$ level than with \texttt{gamm4}. The \texttt{gamm4} method largely undercovers the packs component of $\beta_0$, which can be explained by the incorporated bias.

\begin{figure}[]
	\centering
	\caption[]{\label{fig:pseudorandom} 
	Coverage and  length of two-sided confidence intervals at significance level 5\% and  bias for our method, \texttt{mmdml}, and 
	\texttt{gamm4}. 
	In the coverage plot, solid dots represent point estimators, and circles represent $95\%$ confidence bands with respect to the $5000$ simulation runs.
	The confidence interval length and bias are displayed with boxplots without outliers.}
	\includegraphics[width=0.75\textwidth]{pics/pseudoRand_paper_pics__25_January_2022__17_35_42__pseudorandom__Xforest_Yforest_S10_reps5000_FINAL.pdf}
\end{figure}


\subsection{Simulation Study}\label{sect:simul-study}

We consider a partially linear mixed-effects model with $q=3$ 
random effects and where $\beta_0$ is one-dimensional.
Every subject has their own random intercept term and a nested random effect with two levels.
Thus, the random effects structure is more complex than in the previous two subsections because these models only used a random intercept. 
We compare three data generating mechanisms: 
One where the function $g$ is nonsmooth and the number of observations per group is balanced, 
one where the function $g$ is smooth and the number of observations per group is balanced, 
and one where the function $g$ is nonsmooth and the number of observations per group is unbalanced;
please see Section~\ref{sect:dataSimulation} in the appendix for more details.

We estimate the nonparametric nuisance components, that is, the conditional expectations, with
random forests consisting of $500$ trees whose minimal node size is $5$. 
Furthermore, we use $K=2$ and $\mathcal{S}=10$ in Algorithm~\ref{algo:Summary}. 

We perform $1000$ simulation runs and consider different numbers of groups $N$. As in the previous subsection, we compare the performance of our method with  
\texttt{gamm4}.

The results are displayed in Figure~\ref{fig:simulation}. Our method, \texttt{mmdml}, highly outperforms \texttt{gamm4} in terms of coverage for nonsmooth $g$ because the coverage of \texttt{gamm4} equals $0$ due to its substantial bias. Our method overcovers slightly due to the correction factor that results from the $\mathcal{S}$ repetitions. However, this correction factor is highly recommended in practice.  
With smooth $g$, \texttt{gamm4} is closer to the nominal coverage and has shorter confidence intervals than our method. Because the underlying model is smooth and additive, a spline-based estimator is better suited. In all scenarios,  our method outputs longer confidence intervals than \texttt{gamm4} because we use random forests; consistent with theory, the difference in absolute value decreases though when $N$ increases.

\begin{figure}[]
	\centering
	\caption[]{\label{fig:simulation} 
	Coverage and median length of two-sided confidence intervals for $\beta_0$ at significance level 5\% (true $\beta_0 = 0.5$) and median bias
	for three data generating scenarios for our method, \texttt{mmdml}, and 
	\texttt{gamm4}. 
	The shaded regions in the coverage plot represent $95\%$ confidence bands with respect to the $1000$ simulation runs. 
	The dots in the coverage and bias plot are jittered, but neither are their interconnecting lines nor their confidence bands.}
	\includegraphics[width=\textwidth]{pics/coverage_paper_pics__26_January_2022__99_99_99__simulation__Xforest_Yforest_S10_reps1000_FINAL_merged_1.pdf}
\end{figure}


\section{Conclusion}\label{sect:conclusion}

Our aim was to develop inference for the linear coefficient $\beta_0$ of a partially linear mixed-effects model 
that includes a linear term and potentially complex nonparametric terms.
Such models can be used to describe heterogenous and correlated data that feature 
some grouping structure, which may result from taking repeated measurements. 
Traditionally, spline or kernel approaches are used to cope with the nonparametric part of such a model. We presented a scheme that uses the double machine learning  
framework of~\citet{Chernozhukov2018} to estimate any 
nonparametric components with arbitrary machine learning algorithms. This allowed us to consider complex nonparametric components with interaction structures and high-dimensional variables. 

Our proposed method is as follows. First, the nonparametric variables are regressed out from the response and the linear variables. This step adjusts the response and the linear variables  for the nonparametric variables and may be performed with any machine learning algorithm.
The adjusted variables satisfy a linear mixed-effects model, where the linear coefficient $\beta_0$ can be estimated with standard linear 
mixed-effects
techniques. We showed that the estimator of $\beta_0$ 
asymptotically follows a Gaussian distribution, converges at the parametric rate, and is semiparametrically efficient.
This asymptotic result 
allows us to perform inference for $\beta_0$. 

Empirical experiments demonstrated the performance of our proposed method.  
We conducted an empirical and pseudorandom data analysis and a  simulation study. 
The simulation study and the pseudorandom experiment confirmed the effectiveness of our method in terms of coverage, length of confidence intervals, and estimation bias
compared to a penalized regression spline approach relying on additive models. 
In the empirical experiment, we analyzed longitudinal CD4 cell counts data collected from 
HIV-infected individuals. 
In the literature, most methods only incorporate the time component nonparametrically to analyze this dataset.
Because we estimate nonparametric components with machine learning algorithms, we can allow several variables to enter the model nonlinearly, and we can allow these variables to interact. 
A comparison of our results with the literature suggests that our method  
may perform robust estimation. 

Implementations of our method are available in the \textsf{R}-package \texttt{dmlalg}~\citep{dmlalg}. 

    
\section*{Acknowledgements}

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 786461).

\phantomsection
\addcontentsline{toc}{section}{References}
