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.
73,401 characters · 27 sections · 36 citation commands
xtdml: Double Machine Learning Estimation to Static Panel Data Models with Fixed Effects in R
Double machine learning (DML) is a method for causal inference that leverages machine learning tools to flexibly model the functional form of the covariates while preserving valid statistical inference for the structural parameter (or treatment effect). Several works have contributed to the development of the DML method across a wide range of settings (i.e., treatment endogeneity, treatment heterogeneity, binary and continuous treatments) in cross‐sectional models chernozhukov2018,lewis2020,knaus2022double,huber2023,scaillet2025,chernozhukov2025, static and dynamic panel data models klosin2022,semenova2023,arganaraz2025,clarke2025,marquez2025, and difference-in-differences (DID) designs chang2020,haddad2024,zhang2025. The availability of \proglang{R} R packages for some of these methods has further contributed to the dissemination and practical adoption of DML in applied research.\footnote{A growing number of empirical studies have applied these DML methods to estimate treatment effects, including, for example, deryugina2019, langen2023, strittmatter2023, baiardi2024,baiardi2024plough, and li2025.}
This article introduces the \proglang{R} package \pkg{xtdml} xtdml, which implements the DML method for static panel data models with fixed effects as developed in clarke2025. The package fits partially linear panel regression (PLPR) models with low-dimensional fixed effects and treatment exogeneity, and permits flexible adjustment for high-dimensional confounding. To account for unobserved individual heterogeneity, the user can choose among the following panel data approaches:
The \pkg{xtdml} package builds on the object-oriented architecture of \pkg{DoubleML} DoubleML, which provides a comprehensive and user-friendly toolkit for the estimation of treatment effects in cross-sectional partially linear and interactive regression models within the DML framework.\footnote{The \pkg{xtdml} package can be viewed as the panel data counterpart of \pkg{DoubleML} for estimating partially linear regression models with DML, in the same way as \pkg{plm} plm is the panel data analogue of the \code{lm()} function in the \pkg{stats} stats package for linear regression. For more details on \pkg{DoubleML} see bach2024.} \pkg{xtdml} adopts the \pkg{R6} class system R6 and integrates with the \pkg{mlr3} mlr3 ecosystem to support state-of-the-art learners and tuning methods. The core model class \code{xtdml_plr} supports estimation within PLPR models with fixed effects by learning nuisance components through machine learning and constructing Neyman-orthogonal score functions based on adapted moment conditions. Common tasks implemented by \pkg{xtdml} include the estimation of the target parameter, the calculation of standard errors, t-tests, confidence intervals, and root mean squared errors (RMSEs). Similarly to \pkg{DoubleML}, the package offers users control over the estimation procedure, including the choice of learners for each nuisance function, hyperparameter tuning, the specification of the DML algorithm, and the selection of orthogonal score functions appropriate to the panel data setting. In contrast to the original \pkg{DoubleML}, \pkg{xtdml} introduces several features tailored to panel data including:
\pkg{xtdml} is part of a limited class of \proglang{R} packages that allow the estimation of longitudinal or panel data models with machine learning techniques. Other \proglang{R} packages include the following: \pkg{ablasso} ablasso combines the Arellano-Bond estimation method with the Least Absolute Shrinkage and Selection Operator (Lasso) for the estimation of structural parameters in dynamic linear panel models with lagged dependent variables and predetermined covariates as explanatory variables, and unobserved individual and time fixed effects chernozhukov2024. Among the class of causal estimators implemented in \pkg{causalweight} causalweight, \code{didcontDML} and \code{didcontDMLpanel} functions conduct the estimation of the average treatment effect on the treated using the DID methods for repeated cross-sections and panel data, respectively, while flexibly controlling for confounding via DML chang2020,haddad2024.\footnote{Similar to \pkg{causalweight}, \pkg{xtdml} can accommodate a classical DID design with two time periods and two groups (or a two-way fixed-effects setup) without staggered adoption by specifying \code{approach = "fd-exact"} in \code{xtdml_data_from_data_frame()} and using $Treat_i \times Post_t$ as the treatment variable. The advantage of \pkg{xtdml} over \pkg{causalweight} is the direct integration with the \pkg{mlr3} ecosystem, which provides access to a broad range of base learners.} The package \pkg{desla} desla supports the implementation of desparsified Lasso with both cross-sectional data van2014 and high-dimensional time series adamek2023,adamek2024. \pkg{REEMtrees} REEMtree fits a single regression tree with individual-specific effects sela2012, designed for predict models with longitudinal or clustered data.
The \pkg{hdm} package hdm implements double-selection procedures for high-dimensional cross-sectional settings, fitting Lasso regressions with non-Gaussian errors, with or without instrumental variables. Although originally designed for cross-sectional data, it can be adapted to panel data settings by manually transforming the data through conventional panel data transformations, such as within-group or first-differences. This is possible because the method relies on penalized and linear regressions, provided that the user specifies an extensive dictionary of covariates (such as high-order polynomials and interaction terms) to consistently estimate the parameter of interest. In contrast, the \pkg{DoubleML} toolkit for cross-sectional models cannot be used directly with panel data settings in the same way as \pkg{hdm} as it allows for nonlinear machine-learning learners, for which standard panel data transformations are not valid without further methodological and software adjustments, which are now implemented in \pkg{xtdml}.
The rest of the article is structured as follows. Section (ref) introduces the partially linear regression model with fixed effects and its transformed versions under three panel data approaches. Section (ref) presents the DML algorithm for panel data, and discusses its main elements by panel data approach: the learning procedure, the construction of Neyman orthogonal score functions, and block-k-fold sample splitting. Section (ref) discusses hyperparameter tuning and the transformation of the variables allowed in the package. Section (ref) provides a guide of the main functions and methods in the package with simulated data. Section (ref) illustrates the use of the package with real data. Last section concludes.
Consider the partially linear panel regression (PLPR) model with fixed effect $\alpha^*_i$
where $g_1$ is a non-linear nuisance function of $\vector X_{it}$ and $\mathbb{E}[U_{it}|D_{it},\vector X_{it},\alpha^*_i]=0$ and $\mathbb{E}[\alpha_i^*\mid D_{it},\vector X_{it}]\neq 0$. The target parameter $\theta_0$ is the average partial effect (APE) of continuous $D_{it}$ such that $d\theta_0=\mathbb{E}[Y_{it}(d)-Y_{it}(0)]$, or the average treatment effect (ATE) $\mathbb{E}[Y_{it}(1)-Y_{it}(0)]$ for binary treatments.
The partial-out PLPR (PO-PLPR) model is
where $l_1$ and $m_1$ are nuisance functions, $\alpha_i$ is a fixed effect, $\mathbb{E}[U_{it}\mid V_{it},\vector X_{it},\alpha_i]=0$, and $V_{it}$ the residual of a non-linear additive noise treatment model depending on fixed effect $\gamma_i$ and satisfying $\mathbb{E}[V_{it}|\vector X_{it},\gamma_i]=0$.
Models (ref)–(ref) cannot be directly estimated from the observed data due to the presence of unobserved individual heterogeneity, which is assumed to be correlated with the covariates and the treatment. Ignoring the panel structure and proceeding with a pooled regression without accounting for the fixed effects would make the double machine learning estimator inconsistent due to omitted variable bias. As standard practice in Econometrics, unobserved individual heterogeneity can either be eliminated through panel data transformations (through within-group or first-differences transformation) or explicitly modelled (correlated random effects).
In the following sections, we discuss the panel data approaches implemented in the \pkg{xtdml} package to address the presence of the unobserved individual heterogeneity. These are based on the panel data approaches adapted for DML by clarke2025: correlated random effects (CRE), first-differences exact (FD exact) transformation, and within-group approximation (WG approximation) transformation.
The first approach models unobserved individual heterogeneity through the Mundlak device mundlak1978, converting the fixed-effects model into a random-effects specification. Under additive separability
$\mathbb{E}[Y_{it}(0)\mid \vector X_{it}, \overline{X}_i]= \mathbb{E}[Y_{it}\mid\vector X_{it},\overline{X}_i]- \mathbb{E}[D_{it}\mid\vector X_{it},\overline{X}_i]\theta_0$.
The PO-PLPR model follows because $\mathbb{E}[Y_{it}\mid\vector X_{it},\overline{X}_i]=l_1(\vector X_{it})+\alpha_i$, where $l_1(\vector X_{it})=g_1(\vector X_{it})+m_1(\vector X_{it})\theta_0$ and $\alpha_i=\alpha^*_i+\gamma_i\theta_0$. Then, the PO-PLPR model for CRE with nuisance parameters $\widetilde{l}_1$ and $\widetilde{m}_1$ is
where $\mathbb{E}[U_{it} \mid V_{it},\vector X_{it},{\overline{\vector X}_i},a_i\big]=\mathbb{E}[V_{it}|\vector X_{it},{\overline{\vector X}_i},c_i]=0$, $\widetilde{l}_1(\vector X_{it},{\overline{\vector X}_i})=l_1(\vector X_{it}) + \omega_{\alpha}({\overline{\vector X}_i})$ and $\widetilde{m}_1(\vector X_{it},{\overline{\vector X}_i}) = m_1(\vector X_{it})+\omega_\gamma({\overline{\vector X}_i})$. The random effects, $a_i$ and $c_i$, capture autocorrelation between observations on the same individual. The nuisance parameters $\widetilde{l}_1$ and $\widetilde{m}_1$ can be learnt directly from $\mathbb{E}[Y_{it}\mid\vector X_{it}, {\overline{\vector X}_i}]$ and $\mathbb{E}[D_{it}\mid\vector X_{it}, {\overline{\vector X}_i}]$, respectively.
Similarly, the PLPR model for CRE is
where the nuisance parameter $\widetilde{g}_1$ can be learnt iteratively from $Y_{it} - \mathbb{E}[Y_{it}\mid\vector X_{it}, {\overline{\vector X}_i}]$.
The second type of approaches follows more conventional estimation techniques for panel data, where the unobserved individual heterogeneity is removed from the model by transforming the data. Let $Q$ be a panel data transformation operator, the within-group (WG) or time-demeaning transformation is $Q(W_{it})=W_{it}-\overline{W}_i$ for a generic random variable $W_{it}$, where $\overline{W}_{i}= T^{-1} \sum_{t=1}^T W_{it}$. The first-difference (FD) transformation is $Q(W_{it}) = W_{it}-W_{it-1}$ for $t=2,\dots,T$.
The transformed models for $Q(Y_{it})$ and $Q(V_{it})$ under PO-PLPR model (ref)-(ref) is
and under PLPR model (ref) is
which do not depend on fixed effects $\alpha_i$ and $\gamma_i$ because $Q(\alpha_i)=Q(\alpha^*_i)=Q(\gamma_i)=0$. The transformed nuisance functions $\{Q\big(l_1(\vector X_{it})\big),Q\big(m_1(\vector X_{it})\big)\}$ can be learnt directly from the data while $Q\big(g_1(\vector X_{it})\big)$ can be learnt iteratively from $Q(Y_{it}) - \mathbb{E}[Q(Y_{it})\mid Q(\vector X_{it})]$.
In the following section, we discuss how the learning procedure is implemented in \pkg{xtdml} from available data.
The objective is to make inferences on the target parameter $\theta_0$, given suitable predictions of the vector of nuisance parameters $\bm{\eta}_0=(l_0,m_0,g_0)$ obtained using algorithms from the \pkg{mlr3} ecosystem (e.g., \pkg{mlr3}, \pkg{mlr3learners}, \pkg{mlr3extralearners}). The estimation procedure implemented in \pkg{xtdml} accounts for (a) the presence of the unobserved individual heterogeneity correlated with the observed factors, and (b) possibly nonlinear functions of the covariates to obtain a consistent estimate of the structural (causal) parameter.
The DML procedure presented in this section relies on three key elements. First, the nuisance functions are predicted (or learnt) accounting for both (a) and (b). Second, Neyman orthogonal score functions are constructed based on moment conditions derived from models presented in Section (ref) to reduce the first-order bias introduced by regularized learners (low variance, but large bias). Third, the sample is divided into $K$ groups (or folds), where a part is used for predictions and the rest for estimating the target parameter to reduce overfitting bias (low bias, but large variance). The employed splitting scheme respects the longitudinal structure of the data by randomly assigning a subject with its entire time series to a unique fold $k$ (block-k-fold sample splitting).
Algorithm (ref) below summarizes the DML procedure implemented in \pkg{xtdml}, which adapts Algorithms 1 and 2 in chernozhukov2018 and bach2024 to longitudinal datasets. The argument \code{dml_procedure} in the \code{new()} method specifies the DML algorithm to implement (\code{dml_procedure = c('dml1', 'dml2')}, with \code{dml2} as the default).
The first stage of the DML procedure consists of learning the nuisance functions, $\bm{\eta}_0=(l_0,m_0,g_0)$, using machine learning algorithms while accounting for adjustment due to the longitudinal structure of the data. The \code{approach} argument in the \code{xtdml_data_from_data_frame} function allows the user to choose among \code{("fd-exact", "wg-approx", "cre", "pooled")}, where no default option is provided.\footnote{As previously mentioned, the pooled estimation approach implements the cross-sectional DML estimation with cluster-robust statistical inference, as in \pkg{DoubleML}. This option is recommended for repeated cross-sections.} The data is transformed according to the selected approach using the specified panel and time identifiers, and the prediction tasks for the nuisance parameters are adapted accordingly.
The prediction tasks (or learning stage) implemented in the \pkg{xtdml} package differ by panel data approach.
The CRE Approach (Mundlak Device). The procedure is implemented by selecting the argument \code{approx = "cre"} in the function \code{xtdml_data_from_data_frame}, which internally generates the individual means of the covariates by \code{panel_id}. The nuisance function
The WG Approximation Approach. The procedure is implemented by selecting the argument \code{approx = "wd-approx"}, which internally time-demeans all variables by \code{panel_id}. The nuisance function
With this type of panel approach, the approximation error may be large when the nuisance functions are highly nonlinear and non-smooth functions of the covariates.
The FD Exact Approach. The procedure is implemented by selecting the argument \code{approx = "fd-exact"} (default), which internally generates the first-order lags of the covariates by \code{panel_id}. The nuisance function
The second stage of the DML procedure requires the construction of Neyman orthogonal score functions with the predicted nuisance functions to consistently estimate the structural parameter.
The Neyman orthogonal score function for PLPR models has generic form
where $W_i=\{W_{it}: t=1,\dots,T\}$ with $W_{it} = \{Y_{it},D_{it},\vector X_{it}\}$; ${\vector V}^\perp_i=(V_{i1}^\perp, \dots, V^\perp_{it},\ldots, V^\perp_{iT})'$ contains the orthogonalized regressors chosen to ensure Neyman orthogonality; the row vector $\mathbf{r}_i=(r_{i1},\dots,r_{it},\dots,r_{iT})^\prime$ contains the residuals of the model (either (ref))-(ref)) or (ref))-(ref))); $\Sigma_0(\matrix X_i)=E[\mathbf{r}_i\mathbf{r}'_i\mid X_i]$ is the (potentially) heteroskedastic residual variance-covariance matrix; and $\matrix{X}_i=(\vector X_{i1},\dots,\vector X_{it},\dots,\vector X_{iT})$ is the appropriately chosen set of predictor variables.
The \code{new()} method allows the user to choose between two types of scores, the IV-type score and the partial-out (PO) score, which correspond to the untransformed Model (ref) and the untransformed Models (ref)–(ref), respectively.
For \code{score = "orth-PO"} (default), Equation (ref) takes the following form
For \code{score = "orth-IV"}), Equation (ref) becomes
The main computational distinction between the two score types is that \code{orth-PO} estimates $l_1$ directly, while \code{orth-IV} targets $g_1$ which is obtained iteratively via $l_1$ and $m_1$. As a result, \code{orth-IV} is more computationally intensive because it requires an additional machine-learning fit.
The DML estimator based on the Neyman-orthogonal score (10) solves the finite-sample analogue of the moment condition $\mathbb{E}[\psi^\perp(W_i; \theta_0, \bm{\eta}_0)]=0$ with respect to $\theta_0$, such that
and has closed-form solution
where
The DML estimator is $\sqrt{N}$-consistent for the population parameter $\theta_0$ with a normal limiting distribution as in chernozhukov2018 and approximate variance
where $J_0 =\mathbb{E}({\vector V}_i^{\perp'}{\vector v}_i)$. By default, the \pkg{xtdml} package computes the cluster-robust statistical inference at the \code{panel_id} level, or at the \code{cluster_cols} if specified otherwise.
The third key component of the DML procedure is sample splitting with cross-fitting to reduce the overfitting bias and restore efficiency.
The argument \code{draw_sample_splitting} in the \code{new()} method controls whether sample splitting is conducted. When \code{draw_sample_splitting = TRUE} (default), the sample $\mathcal{W}$ is divided in $k$ folds (\code{n_folds = 5} as default) of equal size $N_k\equiv |W_k| = N/K$.\footnote{For example, if the sample size is small, splitting it into many folds further reduces the information available in the estimation sample for estimation. In such cases, the complementary sample may differ substantially in its observed characteristics, effectively forcing the learner to extrapolate} To preserve the longitudinal structure of the data, each cross-sectional unit with its time series is randomly allocated in the same fold. \pkg{xtdml} uses the specified \code{panel_id} in the \code{xtdml_data_from_data_frame} function to randomly assign subjects into $k$ folds. This sampling strategy is referred as block-k-fold sampling.
For each fold $k$, $W_k$ serves as the main sample for estimating the target parameter $\theta_0$, while its complement $W_k^c$ is used to learn the nuisance functions $\bm{\eta}_0$. When the argument \code{apply_cross_fitting = TRUE} (default) in the \code{new()} method, cross-fitting is conducted: the roles of the folds are switched at every iteration so that each fold is used both as an estimation sample and as a complementary sample.
While Neyman orthogonality makes the DML estimator insensitive to minor prediction errors in the nuisance parameters, finding optimal configurations of hyperparameters remains an important task to bring the learned nuisance functions as close as possible to their population counterpart and reach state-of-the-art performance in treatment effect estimation machlanski2023,machlanski2024,bach2024hyper.
Hyperparameter tuning proceeds with various trials of different configurations of the hyperparameters to tune, and it is a computationally intensive task. Hyperparameter tuning in \pkg{xtdml} is conducted using the tuning methods provided in the \pkg{mlr3tuning} package. By default, the \code{tune()} method uses grid search bergstra2012 to compare different configurations of hyperparameter values. The optimizer (\code{tuner}) randomly searches among a specified number of different values to try per hyperparameter (resolution) and stops the optimization when the specified maximum number of evaluations is reached (terminator). Once the optimal configuration of hyperparameters is found, this is stored and then passed to the DML procedure after calling the \code{fit()} method.
The default tuning procedure in \pkg{xtdml} does not tune on folds (\code{tune_on_folds = FALSE}), and passes all data units into the tuning procedure. Note that the composition of the units in the $k$-th fold differs from the corresponding fold in the DML procedure. When tuning on folds (\code{tune_on_folds = TRUE}), only units in the training sample of fold $k$ are used for tuning and subsequently divided in additional $k$ folds by creating inner training and testing samples.
In practice, base learners require the data to be pre-processed in specific ways to enhance the learning performance. The \pkg{xtdml} package allows the user to choose among three types of data transformations to apply to the covariates $X$ when initializing the data environment, such as \code{transformX = ("no", "minmax", "poly")}. Selecting the option \code{"no"} (default) leaves the covariates untransformed; this option is recommended for tree-based learners. \code{"minmax"} applies the min-max normalization to the covariates, i.e., $\dot{X}=(X_{max}-X_{min})/X_{min}$; this option is recommended with neural networks. \code{"poly"} adds polynomials up to order three and interactions between all possible combinations of two and three variables, such as $\{X_{it,j}X_{it,k}, X_{it,j}^q : j,k\in \{1,\ldots,p\}, q=3\}$; this option is recommended for penalized regressors such as Lasso. Table (ref) shows how the data $\{Y_{it},D_{it},\vector X_{it}\}$ provided by the user is transformed internally by \pkg{xtdml} across panel data approaches.
This section outlines the main steps for estimating the structural (causal) parameter using the \pkg{xtdml} package, including:
The \pkg{xtdml} package is available for installation from CRAN:
After installation, the documentation of the package and available class methods can be accessed by typing:\footnote{It is not possible to run the examples by typing \code{example("xtdml")} because it is an abstract class intended to define the general structure and interface of the package. Examples can be run for the following classes and functions: \code{make_plpr_data} to generate a synthetic data frame for partially linear panel regression (PLPR) models; \code{xtdml_data_from_data_frame} to initialize the data environment; \code{xtdml_plr} to run the routine that estimates the PLPR model with DML.}
To demonstrate the use of the \pkg{xtdml} package, we use simulated data from the built-in function \code{make_plpr_data()}, based on the following data-generating process
where $U_{it} \sim \mathcal{N}(0,1)$, $V_{it} \sim \mathcal{N}(0,1)$, $\alpha_i = \rho A_i + \sqrt{1-\rho^2} B_i $ with $A_i\sim \mathcal{N}(3,3)$, $B_i\sim \mathcal{N}(0,1)$, and $ \gamma_i\sim \mathcal{N}(0,5)$. The covariates are distributed as $X_{it,p} \sim A_i + \mathcal{N}(0, 5)$, where $p$ is the number of covariates.
The nuisance functions are generated as
where $a_1=b_2=0.25$ and $a_2=b_1=0.5$, and $\mathds{1}(.)$ is the indicator operator.
The function \code{make_plpr_data} accepts several arguments, including: the number of cross-sectional observations (\code{n_obs}), number of time periods (\code{t_per}), number of covariates (\code{dim_x}) value of the treatment effect (\code{theta}), degree of correlation between the unobserved individual heterogeneity and the observable characteristics (\code{rho}). In this example, we generate a data set with 1,000 cross-sectional units observed over 10 time periods (for a total of 10,000 observations in the sample), and 20 control variables. The true treatment effect is 0.5; the parameter that governs the relationship between the unobserved individual heterogeneity and the covariates is 0.8.
A built-in example is available by typing \code{example("make_plpr_data")}, which generates a data frame under different parameter choices.
First, the user is required to initialize the data environment using the function \\ \code{xtdml_data_from_data_frame()}, which accepts the following arguments: data, set of raw covariates (\code{x_cols}), exogenous treatment (\code{d_cols}), output / dependent variable (\code{y_col}), panel identifier (\code{panel_id}), time identifier (\code{time_id}), and clustering variable (\code{cluster_cols}) if different from \code{panel_id}, panel data transformation approach (\code{approach}), type of transformation to apply to the set of raw covariates (\code{transformX}).\footnote{Note that panel-related arguments, such as \code{approach} and \code{transformX}, allow users to run the \pkg{xtdml} estimation procedures without performing any ex-ante data manipulation or panel transformations.}
In this example, we specify to prepare the data for the CRE estimation approach, and to leave the covariates untransformed (suitable for the tree-based method we specify in the next section).
In this example, the panel identifier and the clustering variable coincide, so specifying \code{cluster_vars = "id"} is unnecessary.
The output of the data environment can be print by typing:\footnote{A \code{plot()} method is not implemented for data objects because \code{xtdml_data_from_data_frame()} provides a “container” for storing and checking input data, later used in the DML procedure. Therefore, typing \code{obj_xtdml_data\$plot()} would return an error.}
Another built-in example is available by typing the code \code{example("xtdml_data_from_data_frame")} in the \proglang{R} console.
As an illustrative example, we use gradient boosting learner to predict the nuisance functions. We specify: the base learner (\code{xgboost}) for the two nuisance functions, the task to perform (\code{regr} as both the outcome and treatment variables are continuous), and the number of trees to grow (\code{nrounds = 100}).\footnote{The number of boosted trees is set to 100 to reduce the computational time, but it is preferable to be 1000 to enhance estimation accuracy.}
We now initialize the DML estimation environment.
At this stage, the model has not yet been fitted. Consequently, calling the \code{print()} or \code{summary()} methods will issue the warning message \code{"fit() not yet called"}, as no estimation results are available.
An additional built-in example for the classification and regression tree (\code{rpart}) can be run by typing the code \code{example("xtdml_plr")} in \proglang{R} console.
In this section, we proceed by setting up the tuning procedure for a subset of the hyperparameters. Tuning warning messages can be suppressed by typing:
The code below defines the hyperparameter search space and tuning settings, all of which can be customized by the user.
The model is then fitted using the selected (fixed and tuned) hyperparameters by calling the \code{fit()} method.
Calling the \code{print()} method displays the \code{xtdml_obj} object.
Alternatively, a concise summary of the estimated model can be displayed by calling the \code{summary()} method.\footnote{Note that the \code{plot()} method is not implemented in \pkg{xtdml}. Calling \code{xtdml_obj\$plot()} returns an warning message indicating that the plotting method is not supported for \pkg{xtdml} objects.}
Additional statistical inference can be conducted by calling the method \code{confint()} which calculates the 95% confidence intervals.
The lower and upper bounds of the confidence interval can be extracted separately typting:
The created object \code{xtdml_obj} stores the results after calling the \code{fit()} method. We can extract various estimates, such as the estimated coefficient of the treatment effect, cluster-robust standard error, p-value, root mean squared errors (RMSEs), selected hyperparameters, panel data metadata.
Regarding the quality of learners, we can extract the model RMSE and RMSE of the nuisance functions.
Panel-related metadata, such as the number of observations, number of unique subjects, and number of clusters used in the mode, can be extracted by calling:
We can also retrieve the selected (fixed and tuned) hyperparameters for each nuisance function, which are used in the prediction phase of DML procedure, using the dedicated method \code{xtdml_obj\$get_params()}.\footnote{The field \code{\$params} also provides access to selected hyperparameters of each base learner by the nuisance parameter, but it is not strictly equivalent to the \code{\$get_params()} method. That is, \code{\$params} stores the raw internal parameter lists for each learner, reflecting the object’s internal state, whereas \code{\$get_params()} is a public interface that returns validated and processed parameter values after tuning, performing internal checks and incorporating defaults.} A valid nuisance parameter must be specified in the method \code{xtdml_obj\$get_params()}, otherwise \proglang{R} issues an error message. For \code{score = "orth-PO"}, the valid arguments are \code{"ml_l"} and \code{"ml_m"}, which correspond to the nuisance functions $l_0(X) \equiv \mathbb{E}(Y | X)$ and $m_0(X) \equiv \mathbb{E}(D | X)$, respectively.
The hyperparameters used for predicting \code{"ml_l"} are:
The hyperparameters used for predicting \code{"ml_m"} are:
The stored results can be displayed in a customized table as follows.
Similarly, the analysis can be repeated using the IV-type orthogonal score. This requires setting \code{score = "orth-IV"} when calling \code{xtdml_plr\$new()}. For brevity, the corresponding code is omitted here but is provided in the replication package. The resulting estimation output is reported below.
To illustrate the use of the \pkg{xtdml} package, we use the \code{Produc} panel dataset, available in the \pkg{plm} package in \proglang{R}. The \code{Produc} dataframe contains information on 48 US states over 1970-1986 (17 years) in long format, including:
In this empirical example, we compare the effect of increased public investment (\code{pcap}) on state-level economic growth (\code{gsp}) using different base learners and panel data approaches within the DML framework for panel data. The nuisance functions, \code{ml_l} and \code{ml_m}, are learnt using three base learners from \pkg{mlr3} ecosystem: boosted trees (\code{xgboost}), cross-validated Lasso (\code{cv_glmet}), and neural networks (\code{nnet}). We select the partial-out score function (\code{score = 'orth-PO'}), 3 folds (\code{n_folds = 3}), with cross fitting (\code{apply_cross_fitting = TRUE}). This example allows us to compute clustered statistical inference at a level (the region) different from the panel identifier (the US states).
We first load data from \pkg{plm} package.
We then generate a binary treatment variable by defining treated units as those above the median value of \code{pcap}.
Generate the logarithm of the variables and define output variable, treatment variables, and controls to pass into the DML environment.
In the next sections, we report the code to estimate the structural parameter using the FD exact approach, and the estimation results obtained with the CRE and WG approximation approaches; detailed code is available in the replication package.
We use the \code{xgboost} learner available in \pkg{mlr3extralearners}, and fix the number of boosted trees to 100 to reduce the computational time.
We initialize the data environment selecting the options of estimating the model using the FD exact approach and raw variables without any transformation (\code{transformX = "no"}).
We initialize the DML estimation environment using the data object created above.
We define the ranges of the the hyperparameters to tune tune via grid search (\code{lambda} is the L2 regularization term on weights; and \code{max_depth} is the maximum depth of the tree), and the relative tuning options (\code{terminator} is the stop criterion of the tuning process after \code{n_evals} combinations of values to try, and \code{resolution} is the number of grid points generated per hyperparameter to tune).
We fit the model and print summary of results.
As a second learner, we use cross-validated Lasso (\code{cv_glmnet}), where the optimal value of the hyperparameter \code{lambda} is selected as the one that minimizes the cross-validated error.
We then initialize the data and estimation environments. We select the option to include a polynomials and their interactions of the raw covariates (\code{transformX = "poly"}) to satisfy Lasso weak sparsity assumption.
We start the DML estimation procedure and print the results.
The third learner we employ is a single-layer neural network with 100 maximum number of iterations (\code{maxit}), 1000 maximum allowable number of weights (\code{MaxNWts}), and no switch for tracing optimization (\code{trace = FALSE}).
We initialize the data environment and estimation environment. We apply a min-max transformation (\code{transformX = "minmax"}) to ensure that all raw covariates are on a comparable scale. This improves the stability and convergence of neural networks, which are sensitive to the relative magnitudes of the features.
We set the hyperparameters to tune: \code{size} being the number of units in the hidden layer, and \code{decay} being the parameter for weight decay. We then initialize the tuning task via grid search.
We now fit the DML model using tuned hyperparameters. We finally print the summary of the results.
We extract the DML objects relevant for statistical analysis and display them in a customized table.
The analysis is also repeated using alternative panel data approaches available in \pkg{xtdml}. For brevity, we omit the code needed to reproduce these additional results; the additional code is provided in the replication package.\footnote{The omitted code loops over the three panel approaches in the following order \code{approaches = c("cre", "wg-approx")}. Note that the classification task (\code{classif}) is possible only for the CRE approach, because the binary treatment becomes discrete after the WG and FD transformations. The analysis in this section can also be repeated choosing the Neyman-orthogonal IV-type score (\code{score = "orth-IV"}), or other base learners or ensembles from \pkg{mlr3}.}
The output for \code{approach = "cre"} is
The output for \code{approach = "wg-approx"} is
The overall running time for the code with the FD exact approach is 38.8 seconds, 36.79 seconds for CRE, and 35.44 seconds for WG approximation.
The DML results show notable variation in estimated effects across learners and panel data transformations. This reflects both the different nonlinear patterns that learners can accommodate and the modelling assumptions behind each panel data transformation.
The FD exact estimates are more stable across learners than CRE estimates, and noticeably smaller than those obtained with the WG approximation approach. For the CRE approximation approaches, the DML estimates obtained with gradient boosting have opposite signs than those produced by Lasso and neural network. This may reflect sensitivity of the learner to the small cross-sectional size ($N = 48$) and the limited number of boosting rounds ($B = 100$). In addition, the DML point estimates under FD exact and CRE are insignificant for Lasso and neural networks, unlike gradient boosting. Under the WG approximation, the effects are all statistically significant at 5% level.
Regarding the predictive performance, model RMSEs from the FD exact estimations are lower and more consistent across learners, indicating greater precision. The RMSEs of the nuisance functions $l$ and $m$ are also smaller for FD exact than for CRE and WG approximation, suggesting that learners more accurately capture the relationships involving $(X_{it}, X_{it-1})$ for FD exact than those based on $(X_{it}, \overline{X}_i)$ for CRE or $Q(X_{it})$ for WG approximation. Finally, neural networks consistently report the lowest RMSEs for $l$ and $m$, indicating its superior predictive accuracy than the other selected learners, particularly relative to gradient boosting, across all panel-data approaches.
Overall, the FD exact approach is generally more robust, requires no structural assumptions on individual effects, and produces more accurate predictions of the nuisance functions since the learners jointly predict the relationship between the covariates and their lags, which is a more demanding learning task than in CRE or WG settings.
This article introduces the R package \pkg{xtdml}, which currently implements partially linear panel regression models with high-dimensional confounding variables, low-dimensional fixed-effects and exogenous treatment within the double machine learning framework, as proposed by clarke2025. The article documents the functionalities of the \pkg{xtdml} package with a simulated example, and provides practical guidance for its use with real data when different panel data approaches and base learners are chosen.
The results in this paper were obtained using R version 4.5.1 (2025-06-13 ucrt) on a laptop equipped with a 12th Gen Intel Core i7-1260P processor with 12 cores and 16 threads (base clock speed 2.1 GHz), 16 GB of RAM (15.9 GB installed), and running a 64-bit Windows operating system on an x64-based processor. The following packages have been used to run the code: \pkg{xtdml} 0.1.11, \pkg{xgboost} 1.7.11.1, \pkg{nnet} 7.3-20, \pkg{glmnet} 4.1-10, \pkg{mlr3} 1.3.0, \pkg{mlr3misc} 0.19.0, \pkg{mlr3tuning} 1.5.0, \pkg{mlr3learners} 0.13.0, \pkg{mlr3measures} 1.2.0, \pkg{paradox} 1.0.1, \pkg{bbotk} 1.8.1, \pkg{R6} 2.6.1, \pkg{checkmate} 2.3.3, \pkg{clusterGeneration} 1.3.8, \pkg{MASS} 7.3-65, \pkg{MLmetrics} 1.1.3, \pkg{Rcpp} 1.1.0.
The author is grateful to Modesto Escobar Mercado for his assistance with the submission of the \pkg{xtdml} package to CRAN during his visit to the Institute for Social and Economic Research, and thanks Paul Clarke and Omar Hussein for their valuable comments. The author also acknowledges that this research was supported by the UK Economic and Social Research Council (ESRC) under award ES/S012486/1 (MiSoC).