diff --git a/doubleml/double_ml_data.py b/doubleml/double_ml_data.py index 04ecb1e32..fdee739dd 100644 --- a/doubleml/double_ml_data.py +++ b/doubleml/double_ml_data.py @@ -110,7 +110,7 @@ class DoubleMLData(DoubleMLBaseData): Default is ``None``. s_col : None or str - The selection variable (only relevant/used for SSM Estimatiors). + The score or selection variable (only relevant/used for RDD or SSM Estimatiors). Default is ``None``. use_other_treat_as_covariate : bool @@ -182,7 +182,7 @@ def _data_summary_str(self): if self.t_col is not None: data_summary += f'Time variable: {self.t_col}\n' if self.s_col is not None: - data_summary += f'Selection variable: {self.s_col}\n' + data_summary += f'Score/Selection variable: {self.s_col}\n' data_summary += f'No. Observations: {self.n_obs}\n' return data_summary @@ -212,7 +212,7 @@ def from_arrays(cls, x, y, d, z=None, t=None, s=None, use_other_treat_as_covaria Default is ``None``. s : :class:`numpy.ndarray` - Array of the selection variable (only relevant/used for SSM models). + Array of the score or selection variable (only relevant/used for RDD and SSM models). Default is ``None``. use_other_treat_as_covariate : bool @@ -351,7 +351,7 @@ def t(self): @property def s(self): """ - Array of selection variable. + Array of score or selection variable. """ if self.s_col is not None: return self._s.values @@ -538,7 +538,7 @@ def t_col(self, value): @property def s_col(self): """ - The selection variable. + The score or selection variable. """ return self._s_col @@ -547,10 +547,10 @@ def s_col(self, value): reset_value = hasattr(self, '_s_col') if value is not None: if not isinstance(value, str): - raise TypeError('The selection variable s_col must be of str type (or None). ' + raise TypeError('The score or selection variable s_col must be of str type (or None). ' f'{str(value)} of type {str(type(value))} was passed.') if value not in self.all_variables: - raise ValueError('Invalid selection variable s_col. ' + raise ValueError('Invalid score or selection variable s_col. ' f'{value} is no data column.') self._s_col = value else: @@ -725,24 +725,24 @@ def _check_disjoint_sets_t_s(self): if self.s_col is not None: s_col_set = {self.s_col} if not s_col_set.isdisjoint(x_cols_set): - raise ValueError(f'{str(self.s_col)} cannot be set as selection variable ``s_col`` and covariate in ' + raise ValueError(f'{str(self.s_col)} cannot be set as score or selection variable ``s_col`` and covariate in ' '``x_cols``.') if not s_col_set.isdisjoint(d_cols_set): - raise ValueError(f'{str(self.s_col)} cannot be set as selection variable ``s_col`` and treatment variable in ' - '``d_cols``.') + raise ValueError(f'{str(self.s_col)} cannot be set as score or selection variable ``s_col`` and treatment ' + 'variable in ``d_cols``.') if not s_col_set.isdisjoint(y_col_set): - raise ValueError(f'{str(self.s_col)} cannot be set as selection variable ``s_col`` and outcome variable ' - '``y_col``.') + raise ValueError(f'{str(self.s_col)} cannot be set as score or selection variable ``s_col`` and outcome ' + 'variable ``y_col``.') if self.z_cols is not None: z_cols_set = set(self.z_cols) if not s_col_set.isdisjoint(z_cols_set): - raise ValueError(f'{str(self.s_col)} cannot be set as selection variable ``s_col`` and instrumental ' - 'variable in ``z_cols``.') + raise ValueError(f'{str(self.s_col)} cannot be set as score or selection variable ``s_col`` and ' + 'instrumental variable in ``z_cols``.') if self.t_col is not None: t_col_set = {self.t_col} if not s_col_set.isdisjoint(t_col_set): - raise ValueError(f'{str(self.s_col)} cannot be set as selection variable ``s_col`` and time variable ' - '``t_col``.') + raise ValueError(f'{str(self.s_col)} cannot be set as score or selection variable ``s_col`` and time ' + 'variable ``t_col``.') class DoubleMLClusterData(DoubleMLData): @@ -780,7 +780,7 @@ class DoubleMLClusterData(DoubleMLData): Default is ``None``. s_col : None or str - The selection variable (only relevant/used for SSM Estimatiors). + The score or selection variable (only relevant/used for RDD and SSM Estimatiors). Default is ``None``. use_other_treat_as_covariate : bool @@ -854,7 +854,7 @@ def _data_summary_str(self): if self.t_col is not None: data_summary += f'Time variable: {self.t_col}\n' if self.s_col is not None: - data_summary += f'Selection variable: {self.s_col}\n' + data_summary += f'Score/Selection variable: {self.s_col}\n' data_summary += f'No. Observations: {self.n_obs}\n' return data_summary @@ -888,7 +888,7 @@ def from_arrays(cls, x, y, d, cluster_vars, z=None, t=None, s=None, use_other_tr Default is ``None``. s : :class:`numpy.ndarray` - Array of the selection variable (only relevant/used for SSM models). + Array of the score or selection variable (only relevant/used for RDD or SSM models). Default is ``None``. use_other_treat_as_covariate : bool @@ -1039,7 +1039,7 @@ def _check_disjoint_sets_cluster_cols(self): 'cluster variable in ``cluster_cols``.') if self.s_col is not None: if not s_col_set.isdisjoint(cluster_cols_set): - raise ValueError(f'{str(self.s_col)} cannot be set as selection variable ``s_col`` and ' + raise ValueError(f'{str(self.s_col)} cannot be set as score or selection variable ``s_col`` and ' 'cluster variable in ``cluster_cols``.') def _set_cluster_vars(self): diff --git a/doubleml/rdd/__init__.py b/doubleml/rdd/__init__.py new file mode 100644 index 000000000..99190a943 --- /dev/null +++ b/doubleml/rdd/__init__.py @@ -0,0 +1,9 @@ +""" +The :mod:`doubleml.rdd` module implements double machine learning estimates for regression discontinuity designs. +""" + +from .rdd import RDFlex + +__all__ = [ + "RDFlex", +] diff --git a/doubleml/rdd/datasets/__init__.py b/doubleml/rdd/datasets/__init__.py new file mode 100644 index 000000000..c54c7827a --- /dev/null +++ b/doubleml/rdd/datasets/__init__.py @@ -0,0 +1,9 @@ +""" +The :mod:`doubleml.rdd.datasets` module implements data generating processes for regression discontinuity designs. +""" + +from .simple_dgp import make_simple_rdd_data + +__all__ = [ + "make_simple_rdd_data", +] diff --git a/doubleml/rdd/datasets/simple_dgp.py b/doubleml/rdd/datasets/simple_dgp.py new file mode 100644 index 000000000..0964ec585 --- /dev/null +++ b/doubleml/rdd/datasets/simple_dgp.py @@ -0,0 +1,103 @@ +import numpy as np +from numpy.polynomial.polynomial import Polynomial + + +def make_simple_rdd_data(n_obs=5000, p=4, fuzzy=True, binary_outcome=False, **kwargs): + """ + Generates synthetic data for a regression discontinuity design (RDD) analysis. + + .. math:: + Y_0 &= g_0 + g_{cov} + \\epsilon_0 \\ + Y_1 &= g_1 + g_{cov} + \\epsilon_1 \\ + g_0 &= 0.1 \\cdot \\text{score}^2 \\ + g_1 &= \tau + 0.1 \\cdot \\text{score}^2 - 0.5 \\cdot \\text{score}^2 \\ + g_{cov} &= \\sum_{i=1}^{\text{dim\\_x}} \text{Polynomial}(X_i) \\ + \\epsilon_0, \\epsilon_1 &\\sim \\mathcal{N}(0, 0.2^2) + + Parameters + ---------- + n_obs : int + Number of observations to generate. Default is 5000. + + p : int + Degree of the polynomial for covariates. Default is 4. + + fuzzy : bool + If True, generates data for a fuzzy RDD. Default is True. + + binary_outcome : bool + If True, generates binary outcomes. Default is False. + + **kwargs : Additional keyword arguments. + cutoff : float + The cutoff value for the score. Default is 0.0. + dim_x : int + The number of independent covariates. Default is 3. + a : float + Factor to control interaction of score and covariates to the outcome equation. Default is 0.0. + tau : float + Parameter to control the true effect in the generated data at the given cutoff. Default is 1.0. + + Returns + ------- + dict: A dictionary containing the generated data with keys: + 'score' (np.ndarray): The running variable. + 'X' (np.ndarray): The independent covariates. + 'Y0' (np.ndarray): The potential outcomes without treatment. + 'Y1' (np.ndarray): The potential outcomes with treatment. + 'intended_treatment' (np.ndarray): The intended treatment assignment. + """ + + cutoff = kwargs.get('cutoff', 0.0) + dim_x = kwargs.get('dim_x', 3) + a = kwargs.get('a', 0.0) + tau = kwargs.get('tau', 1.0) + + score = np.random.normal(size=n_obs) + # independent covariates + X = np.random.uniform(size=(n_obs, dim_x), low=-1, high=1) + + # Create polynomials of covariates + if p == 0: + covs = np.zeros((n_obs, 1)) + else: + covs = np.column_stack([Polynomial(np.arange(p + 1))(X[:, i]) for i in range(X.shape[1])]) + g_cov = np.sum(covs, axis=1) + + g0 = 0.1 * score**2 + g1 = tau + 0.1 * score**2 - 0.5 * score**2 + a * np.sum(X, axis=1) * score + + eps_scale = 0.2 + # potential outcomes with independent errors + if not binary_outcome: + Y0 = g0 + g_cov + np.random.normal(size=n_obs, scale=eps_scale) + Y1 = g1 + g_cov + np.random.normal(size=n_obs, scale=eps_scale) + else: + p_Y0 = 1 / (1 + np.exp(-1.0 * (g0 + g_cov))) + p_Y1 = 1 / (1 + np.exp(-1.0 * (g1 + g_cov))) + Y0 = np.random.binomial(n=1, p=p_Y0, size=n_obs) + Y1 = np.random.binomial(n=1, p=p_Y1, size=n_obs) + + intended_treatment = (score >= cutoff).astype(int) + if fuzzy: + prob = 0.3 + 0.4 * intended_treatment + 0.01 * score**2 - 0.02 * score**2 * intended_treatment + 0.2 * g_cov + prob = np.clip(prob, 0.0, 1.0) + D = np.random.binomial(n=1, p=prob, size=n_obs) + else: + D = intended_treatment + + D = D.astype(int) + Y = Y0 * (1 - D) + Y1 * D + + oracle_values = { + 'Y0': Y0, + 'Y1': Y1, + } + res_dict = { + 'score': score, + 'Y': Y, + 'D': D, + 'X': X, + 'oracle_values': oracle_values + } + return res_dict diff --git a/doubleml/rdd/rdd.py b/doubleml/rdd/rdd.py new file mode 100644 index 000000000..1bbe68306 --- /dev/null +++ b/doubleml/rdd/rdd.py @@ -0,0 +1,600 @@ +import warnings +import numpy as np +import pandas as pd +from collections.abc import Callable + +from scipy.stats import norm +from rdrobust import rdrobust, rdbwselect + +from sklearn.base import clone +from sklearn.utils.multiclass import type_of_target + +from doubleml import DoubleMLData +from doubleml.double_ml import DoubleML +from doubleml.utils.resampling import DoubleMLResampling +from doubleml.utils._checks import _check_resampling_specification, _check_supports_sample_weights + + +class RDFlex(): + """Flexible adjustment with double machine learning for regression discontinuity designs + + Parameters + ---------- + obj_dml_data : :class:`DoubleMLData` object + The :class:`DoubleMLData` object providing the data and specifying the variables for the causal model. + + ml_g : estimator implementing ``fit()`` and ``predict()`` + A machine learner implementing ``fit()`` and ``predict()`` methods and support ``sample_weights`` (e.g. + :py:class:`sklearn.ensemble.RandomForestRegressor`) for the nuisance functions + :math:`g_0^{\\pm}(X) = E[Y|\\text{score}=\\text{cutoff}^{\\pm}, X]`. The adjustment function is then + defined as :math:`\\eta_0(X) = (g_0^{+}(X) + g_0^{-}(X))/2`. + + ml_m : classifier implementing ``fit()`` and ``predict_proba()`` or None + A machine learner implementing ``fit()`` and ``predict_proba()`` methods and support ``sample_weights``(e.g. + :py:class:`sklearn.ensemble.RandomForestClassifier`) for the nuisance functions + :math:`m_0^{\\pm}(X) = E[D|\\text{score}=\\text{cutoff}^{\\pm}, X]`. The adjustment function is then + defined as :math:`\\eta_0(X) = (m_0^{+}(X) + m_0^{-}(X))/2`. + Or None, in case of a non-fuzzy design. + Default is ``None``. + + fuzzy : bool + Indicates whether to fit a fuzzy or a sharp design. + That is if the intended treatment defined by the cutoff can diverge from the actual treatment given + with ``obj_dml_data.d``. + Default is ``False``. + + n_folds : int + Number of folds. + Default is ``5``. + + n_rep : int + Number of repetitons for the sample splitting. + Default is ``1``. + + cutoff : float or int + A float or intspecifying the cutoff in the score. + Default is ``0``. + + h_fs : float or None + Initial bandwidth in the first stage estimation. If ``None``, then the optimal bandwidth without + covariates will be used. + Default is ``None``. + + fs_specification : str + Specification of the first stage regression. The options are ``cutoff``, ``cutoff and score`` and + ``interacted cutoff and score``. + Default is ``cutoff``. + + fs_kernel : str + Kernel for the first stage estimation. ``uniform``, ``triangular`` and ``epanechnikov``are supported. + Default is ``triangular``. + + **kwargs : kwargs + Key-worded arguments that are not used within RDFlex but directly handed to rdrobust. + + Examples + -------- + + Notes + ----- + + """ + + def __init__(self, + obj_dml_data, + ml_g, + ml_m=None, + fuzzy=False, + cutoff=0, + n_folds=5, + n_rep=1, + h_fs=None, + fs_specification="cutoff", + fs_kernel="triangular", + **kwargs): + + self._check_data(obj_dml_data, cutoff) + self._dml_data = obj_dml_data + + self._score = self._dml_data.s - cutoff + self._cutoff = cutoff + self._intendend_treatment = (self._score >= 0).astype(bool) + self._fuzzy = fuzzy + + if not fuzzy and any(self._dml_data.d != self._intendend_treatment): + warnings.warn('A sharp RD design is being estimated, but the data indicate that the design is fuzzy.') + + self._check_and_set_learner(ml_g, ml_m) + + _check_resampling_specification(n_folds, n_rep) + self._n_folds = n_folds + self._n_rep = n_rep + + if h_fs is None: + fuzzy = self._dml_data.d if self._fuzzy else None + self._h_fs = rdbwselect(y=obj_dml_data.y, + x=self._score, + fuzzy=fuzzy).bws.values.flatten().max() + else: + if not isinstance(h_fs, (float)): + raise TypeError("Initial bandwidth 'h_fs' has to be a float. " + f'Object of type {str(type(h_fs))} passed.') + self._h_fs = h_fs + + self._fs_specification = self._check_fs_specification(fs_specification) + self._fs_kernel_function, self._fs_kernel_name = self._check_and_set_kernel(fs_kernel) + self._w = self._calc_weights(kernel=self._fs_kernel_function, h=self.h_fs) + + self._check_effect_sign() + + # TODO: Add further input checks + self.kwargs = kwargs + + self._smpls = DoubleMLResampling(n_folds=self.n_folds, n_rep=self.n_rep, n_obs=obj_dml_data.n_obs, + stratify=obj_dml_data.d).split_samples() + + self._M_Y, self._M_D, self._h, self._rdd_obj, \ + self._all_coef, self._all_se, self._all_ci = self._initialize_arrays() + + # Initialize all properties to None + self._coef = None + self._se = None + self._ci = None + self._N_h = None + self._final_h = None + self._all_cis = None + self._i_rep = None + + def __str__(self): + if np.any(~np.isnan(self._M_Y[:, 0])): + method_names = ["Conventional", "Robust"] + lines = [ + "Method Coef. S.E. t-stat P>|t| 95% CI", + "-------------------------------------------------------------------------", + ] + + for i, name in enumerate(method_names): + if name == "Conventional": + line = ( + f"{name:<18}" + f"{self.coef[i]:<10.3f}" + f"{self.se[i]:<10.3f}" + f"{self.t_stat[i]:<9.3f}" + f"{self.pval[i]:<11.3e}" + f"[{self.ci[i, 0]:.3f}, {self.ci[i, 1]:.3f}]" + ) + else: + assert name == "Robust" + # Access robust values from index 2 as specified + line = ( + f"{name:<17}" + " - - " + f"{self.t_stat[2]:<9.3f}" + f"{self.pval[2]:<11.3e}" + f"[{self.ci[2, 0]:.3f}, {self.ci[2, 1]:.3f}]" + ) + + lines.append(line) + result = "\n".join(lines) + + additional_info = ( + "\nDesign Type: " + ("Fuzzy" if self.fuzzy else "Sharp") + + f"\nCutoff: {self.cutoff}" + + f"\nFirst Stage Kernel: {self.fs_kernel}" + + f"\nFinal Bandwidth: {self.h}" + ) + + return result + additional_info + + else: + return "DoubleML RDFlex Object. Run `.fit()` for estimation." + + @property + def fuzzy(self): + """ + Indicates whether the design is fuzzy or not. + """ + return self._fuzzy + + @property + def n_folds(self): + """ + Number of folds. + """ + return self._n_folds + + @property + def n_rep(self): + """ + Number of repetitions for the sample splitting. + """ + return self._n_rep + + @property + def h_fs(self): + """ + Initial bandwidth in the first stage estimation. + """ + return self._h_fs + + @property + def h(self): + """ + Array of final bandwidths in the last stage estimation (shape (``n_rep``,)). + """ + return self._h + + @property + def fs_kernel(self): + """ + Kernel for the first stage estimation. + """ + return self._fs_kernel_name + + @property + def w(self): + """ + Weights for the first stage estimation. + """ + return self._w + + @property + def cutoff(self): + """ + Cutoff at which the treatment effect is estimated. + """ + return self._cutoff + + @property + def coef(self): + """ + Estimates for the causal parameter after calling :meth:`fit`. + """ + return self._coef + + @property + def se(self): + """ + Standard errors for the causal parameter(s) after calling :meth:`fit`. + """ + return self._se + + @property + def t_stat(self): + """ + t-statistics for the causal parameter(s) after calling :meth:`fit`. + """ + t_stat = self.coef / self.se + return t_stat + + @property + def pval(self): + """ + p-values for the causal parameter(s) after calling :meth:`fit`. + """ + pval = 2 * norm.cdf(-np.abs(self.t_stat)) + return pval + + @property + def ci(self): + """ + Confidence intervals for the causal parameter(s) after calling :meth:`fit`. + """ + return self._ci + + @property + def all_coef(self): + """ + Estimates of the causal parameter(s) for the ``n_rep`` different sample splits after calling :meth:`fit`. + """ + return self._all_coef + + @property + def all_se(self): + """ + Standard errors of the causal parameter(s) for the ``n_rep`` different sample splits after calling :meth:`fit`. + """ + return self._all_se + + def fit(self, n_iterations=2): + """ + Estimate RDFlex model. + + Parameters + ---------- + + n_iterations : int + Number of iterations for the iterative bandwidth fitting. + Default is ``2``. + + Returns + ------- + self : object + """ + + self._check_iterations(n_iterations) + + # set variables for readablitity + Y = self._dml_data.y + D = self._dml_data.d + for i_rep in range(self.n_rep): + self._i_rep = i_rep + + # set initial weights smpls + weights = self.w + + for iteration in range(n_iterations): + eta_Y = self._fit_nuisance_model(outcome=Y, estimator_name="ml_g", + weights=weights, smpls=self._smpls[i_rep]) + self._M_Y[:, i_rep] = Y - eta_Y + + if self.fuzzy: + eta_D = self._fit_nuisance_model(outcome=D, estimator_name="ml_m", + weights=weights, smpls=self._smpls[i_rep]) + self._M_D[:, i_rep] = D - eta_D + + # update weights via iterative bandwidth fitting + if iteration < (n_iterations - 1): + h, b, weights = self._update_weights() + else: + if n_iterations == 1: + h = None + b = None + + rdd_res = self._fit_rdd(h=h, b=b) + self._set_coefs(rdd_res, h) + + self.aggregate_over_splits() + + return self + + def confint(self, level=0.95): + """ + Confidence intervals for RDFlex models. + + Parameters + ---------- + level : float + The confidence level. + Default is ``0.95``. + + Returns + ------- + df_ci : pd.DataFrame + A data frame with the confidence interval(s). + """ + if not isinstance(level, float): + raise TypeError('The confidence level must be of float type. ' + f'{str(level)} of type {str(type(level))} was passed.') + if (level <= 0) | (level >= 1): + raise ValueError('The confidence level must be in (0,1). ' + f'{str(level)} was passed.') + + # compute critical values + alpha = 1 - level + percentages = np.array([alpha / 2, 1. - alpha / 2]) + + critical_values = np.repeat(norm.ppf(percentages[1]), self._n_rep) + + # compute all cis over repetitions (shape: n_coef x 2 x n_rep) + self._all_cis = np.stack( + (self.all_coef - self.all_se * critical_values, + self.all_coef + self.all_se * critical_values), + axis=1) + ci = np.median(self._all_cis, axis=2) + df_ci = pd.DataFrame(ci, columns=['{:.1f} %'.format(i * 100) for i in percentages], + index=['Conventional', 'Bias-Corrected', 'Robust']) + + return df_ci + + def _fit_nuisance_model(self, outcome, estimator_name, weights, smpls): + + # Include transformation of score and cutoff if necessary + if self._fs_specification == "cutoff": + Z = self._intendend_treatment # instrument for treatment + Z_left = np.zeros_like(Z) + Z_right = np.ones_like(Z) + elif self._fs_specification == "cutoff and score": + Z = np.column_stack((self._intendend_treatment, self._score)) + Z_left = np.zeros_like(Z) + Z_right = np.column_stack((np.ones_like(self._intendend_treatment), np.zeros_like(self._score))) + else: + assert self._fs_specification == "interacted cutoff and score" + Z = np.column_stack((self._intendend_treatment, self._intendend_treatment * self._score, self._score)) + Z_left = np.zeros_like(Z) + Z_right = np.column_stack((np.ones_like(self._intendend_treatment), np.zeros_like(self._score), + np.zeros_like(self._score))) + + X = self._dml_data.x + ZX = np.column_stack((Z, X)) + ZX_left = np.column_stack((Z_left, X)) + ZX_right = np.column_stack((Z_right, X)) + + mu_left, mu_right = np.full_like(outcome, fill_value=np.nan), np.full_like(outcome, fill_value=np.nan) + + for train_index, test_index in smpls: + estimator = clone(self._learner[estimator_name]) + estimator.fit(ZX[train_index], outcome[train_index], sample_weight=weights[train_index]) + + if self._predict_method[estimator_name] == "predict": + mu_left[test_index] = estimator.predict(ZX_left[test_index]) + mu_right[test_index] = estimator.predict(ZX_right[test_index]) + else: + assert self._predict_method[estimator_name] == "predict_proba" + mu_left[test_index] = estimator.predict_proba(ZX_left[test_index])[:, 1] + mu_right[test_index] = estimator.predict_proba(ZX_right[test_index])[:, 1] + + return (mu_left + mu_right)/2 + + def _update_weights(self): + rdd_res = self._fit_rdd() + # TODO: "h", "b" features "left" and "right" + h = rdd_res.bws.loc["h"].max() + b = rdd_res.bws.loc["b"].max() + weights = self._calc_weights(kernel=self._fs_kernel_function, h=h) + + return h, b, weights + + def _fit_rdd(self, h=None, b=None): + if self.fuzzy: + rdd_res = rdrobust(y=self._M_Y[:, self._i_rep], x=self._score, + fuzzy=self._M_D[:, self._i_rep], h=h, b=b, **self.kwargs) + else: + rdd_res = rdrobust(y=self._M_Y[:, self._i_rep], x=self._score, + h=h, b=b, **self.kwargs) + return rdd_res + + def _set_coefs(self, rdd_res, h): + self._h[self._i_rep] = h + self._all_coef[:, self._i_rep] = rdd_res.coef.values.flatten() + self._all_se[:, self._i_rep] = rdd_res.se.values.flatten() + self._all_ci[:, :, self._i_rep] = rdd_res.ci.values + self._rdd_obj[self._i_rep] = rdd_res + + def _calc_weights(self, kernel, h): + weights = kernel(self._score, h) + return weights + + def _initialize_arrays(self): + M_Y = np.full(shape=(self._dml_data.n_obs, self.n_rep), fill_value=np.nan) + M_D = np.full(shape=(self._dml_data.n_obs, self.n_rep), fill_value=np.nan) + h = np.full(shape=self.n_rep, fill_value=np.nan) + rdd_obj = [None] * self.n_rep + all_coef = np.full(shape=(3, self.n_rep), fill_value=np.nan) + all_se = np.full(shape=(3, self.n_rep), fill_value=np.nan) + all_ci = np.full(shape=(3, 2, self.n_rep), fill_value=np.nan) + + return M_Y, M_D, h, rdd_obj, all_coef, all_se, all_ci + + def _check_data(self, obj_dml_data, cutoff): + if not isinstance(obj_dml_data, DoubleMLData): + raise TypeError('The data must be of DoubleMLData type. ' + f'{str(obj_dml_data)} of type {str(type(obj_dml_data))} was passed.') + + # score checks + if obj_dml_data.s_col is None: + raise ValueError('Incompatible data. ' + + 'Score variable has not been set. ') + is_continuous = (type_of_target(obj_dml_data.s) == 'continuous') + if not is_continuous: + raise ValueError('Incompatible data. ' + + 'Score variable has to be continuous. ') + + if not isinstance(cutoff, (int, float)): + raise TypeError('Cutoff value has to be a float or int. ' + f'Object of type {str(type(cutoff))} passed.') + if not (obj_dml_data.s.min() <= cutoff <= obj_dml_data.s.max()): + raise ValueError('Cutoff value is not within the range of the score variable. ') + + # treatment checks + one_treat = (obj_dml_data.n_treat == 1) + binary_treat = (type_of_target(obj_dml_data.d) == 'binary') + zero_one_treat = np.all((np.power(obj_dml_data.d, 2) - obj_dml_data.d) == 0) + if not (one_treat & binary_treat & zero_one_treat): + raise ValueError('Incompatible data. ' + 'To fit an RDFlex model with DML ' + 'exactly one binary variable with values 0 and 1 ' + 'needs to be specified as treatment variable.') + + # instrument checks + if obj_dml_data.z_cols is not None: + raise ValueError('Incompatible data. ' + + ' and '.join(obj_dml_data.z_cols) + + ' have been set as instrumental variable(s). ') + + def _check_and_set_learner(self, ml_g, ml_m): + # check ml_g + ml_g_is_classifier = DoubleML._check_learner(ml_g, 'ml_g', regressor=True, classifier=True) + _check_supports_sample_weights(ml_g, 'ml_g') + self._learner = {'ml_g': ml_g} + if ml_g_is_classifier: + if self._dml_data.binary_outcome: + self._predict_method = {'ml_g': 'predict_proba'} + else: + raise ValueError(f'The ml_g learner {str(ml_g)} was identified as classifier ' + 'but the outcome variable is not binary with values 0 and 1.') + else: + self._predict_method = {'ml_g': 'predict'} + + # check ml_m + if self._fuzzy: + if ml_m is not None: + _ = DoubleML._check_learner(ml_m, 'ml_m', regressor=False, classifier=True) + _check_supports_sample_weights(ml_m, 'ml_m') + + self._learner['ml_m'] = ml_m + self._predict_method['ml_m'] = 'predict_proba' + else: + raise ValueError('Fuzzy design requires a classifier ml_m for treatment assignment.') + + else: + if ml_m is not None: + warnings.warn(('A learner ml_m has been provided for for a sharp design but will be ignored. ' + 'A learner ml_m is not required for estimation.')) + + def _check_and_set_kernel(self, fs_kernel): + if not isinstance(fs_kernel, (str, Callable)): + raise TypeError('fs_kernel must be either a string or a callable. ' + f'{str(fs_kernel)} of type {str(type(fs_kernel))} was passed.') + + kernel_functions = { + "uniform": lambda x, h: np.array(np.abs(x) <= h, dtype=float), + "triangular": lambda x, h: np.array(np.maximum(0, (h - np.abs(x)) / h), dtype=float), + "epanechnikov": lambda x, h: np.array(np.where(np.abs(x) < h, .75 * (1 - np.square(x / h)), 0), dtype=float) + } + + if isinstance(fs_kernel, str): + fs_kernel = fs_kernel.casefold() + if fs_kernel not in kernel_functions: + raise ValueError(f"Invalid kernel '{fs_kernel}'. Valid kernels are {list(kernel_functions.keys())}.") + + kernel_function = kernel_functions[fs_kernel] + kernel_name = fs_kernel + + else: + assert callable(fs_kernel) + kernel_function = fs_kernel + kernel_name = 'custom_kernel' + + return kernel_function, kernel_name + + def _check_fs_specification(self, fs_specification): + if not isinstance(fs_specification, str): + raise TypeError("fs_specification must be a string. " + f'{str(fs_specification)} of type {str(type(fs_specification))} was passed.') + expected_specifications = ["cutoff", "cutoff and score", "interacted cutoff and score"] + if fs_specification not in expected_specifications: + raise ValueError(f"Invalid fs_specification '{fs_specification}'. " + f"Valid specifications are {expected_specifications}.") + return fs_specification + + def _check_iterations(self, n_iterations): + """Validate the number of iterations.""" + if not isinstance(n_iterations, int): + raise TypeError('The number of iterations for the iterative bandwidth fitting must be of int type. ' + f'{str(n_iterations)} of type {str(type(n_iterations))} was passed.') + if n_iterations < 1: + raise ValueError('The number of iterations for the iterative bandwidth fitting has to be positive. ' + f'{str(n_iterations)} was passed.') + + def _check_effect_sign(self, tolerance=1e-6): + d_left, d_right = self._dml_data.d[self._score < 0], self._dml_data.d[self._score > 0] + w_left, w_right = self._w[self._score < 0], self._w[self._score > 0] + treatment_prob_difference = np.average(d_left, weights=w_left) - np.average(d_right, weights=w_right) + if treatment_prob_difference > tolerance: + warnings.warn("Treatment probability within bandwidth left from cutoff higher than right from cutoff.\n" + "Treatment assignment might be based on the wrong side of the cutoff.") + + def aggregate_over_splits(self): + var_scaling_factors = np.array([np.sum(res.N_h) for res in self._rdd_obj]) + + self._coef = np.median(self.all_coef, 1) + coefs_deviations = np.square(self.all_coef - self._coef.reshape(-1, 1)) + scaled_coefs_deviations = np.divide(coefs_deviations, var_scaling_factors.reshape(1, -1)) + + rescaled_variances = np.median(np.square(self.all_se) + scaled_coefs_deviations, 1) + self._se = np.sqrt(rescaled_variances) + + self._ci = np.median(self._all_ci, axis=2) + self._N_h = np.median([res.N_h for res in self._rdd_obj], axis=0) + self._final_h = np.median([res.bws for res in self._rdd_obj], axis=0)[0, 0] diff --git a/doubleml/rdd/tests/__init__.py b/doubleml/rdd/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/doubleml/rdd/tests/conftest.py b/doubleml/rdd/tests/conftest.py new file mode 100644 index 000000000..1297572f4 --- /dev/null +++ b/doubleml/rdd/tests/conftest.py @@ -0,0 +1,129 @@ +import pytest + +import numpy as np +import pandas as pd + +from doubleml.rdd.datasets import make_simple_rdd_data +from doubleml import DoubleMLData +from doubleml.rdd import RDFlex + +from rdrobust import rdrobust + +from sklearn.dummy import DummyRegressor, DummyClassifier + + +DATA_SIZE = 500 + +ml_g_dummy = DummyRegressor(strategy='constant', constant=0) +ml_m_dummy = DummyClassifier(strategy='constant', constant=0) + + +@pytest.fixture(scope='module') +def predict_dummy(): + """ + - make predictions using rd-flex with constant model + - make predictions using rdrobust as a reference + """ + def _predict_dummy(data: DoubleMLData, cutoff, alpha, n_rep, p, fs_specification, ml_g=ml_g_dummy): + dml_rdflex = RDFlex( + data, + ml_g=ml_g, + ml_m=ml_m_dummy, + cutoff=cutoff, + n_rep=n_rep, + p=p, + fs_specification=fs_specification + ) + dml_rdflex.fit(n_iterations=1) + ci_manual = dml_rdflex.confint(level=1-alpha) + + rdrobust_model = rdrobust( + y=data.y, + x=data.s, + c=cutoff, + level=100*(1-alpha), + p=p + ) + + reference = { + 'model': rdrobust_model, + 'coef': rdrobust_model.coef.values.flatten(), + 'se': rdrobust_model.se.values.flatten(), + 'ci': rdrobust_model.ci.values + } + + actual = { + 'model': dml_rdflex, + 'coef': dml_rdflex.coef, + 'se': dml_rdflex.se, + 'ci': ci_manual, + } + return reference, actual + + return _predict_dummy + + +def defier_mask(fuzzy, data, actual_cutoff): + if fuzzy == 'left': + # right defiers (not treated even if score suggested it + return (data['D'] == 0) & (data['score'] >= actual_cutoff) + elif fuzzy == 'right': + # left defiers (treated even if score not suggested it + return (data['D'] == 1) & (data['score'] < actual_cutoff) + elif fuzzy in ['both', 'none']: + return None + raise ValueError(f'Invalid type of fuzzyness {fuzzy}') + + +def generate_data( + n_obs: int, + fuzzy: str, + cutoff: float, + binary_outcome: bool = False +): + data = make_simple_rdd_data( + n_obs=n_obs, + fuzzy=fuzzy in ['both', 'left', 'right'], + cutoff=cutoff, + binary_outcome=binary_outcome, + ) + + mask = defier_mask(fuzzy, data, cutoff) + if mask is not None: + data = {k: v[~mask] for k, v in data.items() if k != 'oracle_values'} + + columns = ['y', 'd', 'score'] + ['x' + str(i) for i in range(data['X'].shape[1])] + df = pd.DataFrame( + np.column_stack((data['Y'], data['D'], data['score'], data['X'])), + columns=columns + ) + return DoubleMLData(df, y_col='y', d_cols='d', s_col='score') + + +@pytest.fixture(scope='module') +def rdd_sharp_data(): + def _rdd_sharp_data(cutoff, binary_outcome=False): + return generate_data(n_obs=DATA_SIZE, fuzzy='none', cutoff=cutoff, binary_outcome=binary_outcome) + return _rdd_sharp_data + + +@pytest.fixture(scope='module') +def rdd_fuzzy_data(): + def _rdd_fuzzy_data(cutoff, binary_outcome=False): + return generate_data(n_obs=DATA_SIZE, fuzzy='both', cutoff=cutoff, binary_outcome=binary_outcome) + return _rdd_fuzzy_data + + +@pytest.fixture(scope='module') +def rdd_fuzzy_left_data(): + def _rdd_fuzzy_left_data(cutoff, binary_outcome=False): + return generate_data(n_obs=DATA_SIZE, fuzzy='left', cutoff=cutoff, binary_outcome=binary_outcome) + return _rdd_fuzzy_left_data + + +@pytest.fixture(scope='module') +def rdd_fuzzy_right_data(): + def _rdd_fuzzy_right_data(cutoff, binary_outcome=False): + data = generate_data(n_obs=DATA_SIZE, fuzzy='left', cutoff=cutoff, binary_outcome=binary_outcome) + return data + return _rdd_fuzzy_right_data diff --git a/doubleml/rdd/tests/test_rdd_classifier.py b/doubleml/rdd/tests/test_rdd_classifier.py new file mode 100644 index 000000000..166ab192c --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_classifier.py @@ -0,0 +1,29 @@ +import pytest +import pandas as pd +import numpy as np + +from doubleml.rdd.datasets import make_simple_rdd_data +import doubleml as dml +from doubleml.rdd import RDFlex + +from sklearn.linear_model import LogisticRegression + +np.random.seed(3141) + +n_obs = 300 +data = make_simple_rdd_data(n_obs=n_obs) +data["Y_bin"] = (data["Y"] < np.median(data["Y"])).astype("int") + + +df = pd.DataFrame( + np.column_stack((data['Y_bin'], data['D'], data['score'], data['X'])), + columns=['y', 'd', 'score'] + ['x' + str(i) for i in range(data['X'].shape[1])] +) +dml_data = dml.DoubleMLData(df, y_col='y', d_cols='d', s_col='score') + +dml_rdflex = RDFlex(dml_data, ml_g=LogisticRegression(), ml_m=LogisticRegression(), fuzzy=True) + + +@pytest.mark.ci +def test_rdd_classifier(): + dml_rdflex.fit() diff --git a/doubleml/rdd/tests/test_rdd_classifier_fuzzy.py b/doubleml/rdd/tests/test_rdd_classifier_fuzzy.py new file mode 100644 index 000000000..18ea2ce12 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_classifier_fuzzy.py @@ -0,0 +1,98 @@ +""" +Dummy test using fixed learner for fuzzy data +""" +import pytest +import numpy as np +from sklearn.dummy import DummyClassifier + +ml_g_dummy = DummyClassifier(strategy='constant', constant=0) + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_fuzzy_data, cutoff): + return rdd_fuzzy_data(cutoff, binary_outcome=True) + + +@pytest.fixture(scope='module') +def data_zero(rdd_fuzzy_data): + return rdd_fuzzy_data(0.0, binary_outcome=True) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_classifier_fuzzy_left.py b/doubleml/rdd/tests/test_rdd_classifier_fuzzy_left.py new file mode 100644 index 000000000..50235399d --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_classifier_fuzzy_left.py @@ -0,0 +1,98 @@ +""" +Dummy test using fixed learner for left sided fuzzy data +""" +import pytest +import numpy as np +from sklearn.dummy import DummyClassifier + +ml_g_dummy = DummyClassifier(strategy='constant', constant=0) + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_fuzzy_left_data, cutoff): + return rdd_fuzzy_left_data(cutoff, binary_outcome=True) + + +@pytest.fixture(scope='module') +def data_zero(rdd_fuzzy_left_data): + return rdd_fuzzy_left_data(0.0, binary_outcome=True) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_classifier_fuzzy_right.py b/doubleml/rdd/tests/test_rdd_classifier_fuzzy_right.py new file mode 100644 index 000000000..21d1b5c39 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_classifier_fuzzy_right.py @@ -0,0 +1,98 @@ +""" +Dummy test using fixed learner for right sided fuzzy data +""" +import pytest +import numpy as np +from sklearn.dummy import DummyClassifier + +ml_g_dummy = DummyClassifier(strategy='constant', constant=0) + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_fuzzy_right_data, cutoff): + return rdd_fuzzy_right_data(cutoff, binary_outcome=True) + + +@pytest.fixture(scope='module') +def data_zero(rdd_fuzzy_right_data): + return rdd_fuzzy_right_data(0.0, binary_outcome=True) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_classifier_sharp.py b/doubleml/rdd/tests/test_rdd_classifier_sharp.py new file mode 100644 index 000000000..2bd8faa40 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_classifier_sharp.py @@ -0,0 +1,98 @@ +""" +Dummy test using fixed learner for sharp data +""" +import pytest +import numpy as np +from sklearn.dummy import DummyClassifier + +ml_g_dummy = DummyClassifier(strategy='constant', constant=0) + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_sharp_data, cutoff): + return rdd_sharp_data(cutoff, binary_outcome=True) + + +@pytest.fixture(scope='module') +def data_zero(rdd_sharp_data): + return rdd_sharp_data(0.0, binary_outcome=True) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification, ml_g=ml_g_dummy + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_default_values.py b/doubleml/rdd/tests/test_rdd_default_values.py new file mode 100644 index 000000000..b6e7eafc5 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_default_values.py @@ -0,0 +1,35 @@ +import pytest +import numpy as np +import pandas as pd + +import doubleml as dml +from doubleml.rdd import RDFlex +from doubleml.rdd.datasets import make_simple_rdd_data + +from sklearn.linear_model import Lasso, LogisticRegression + +np.random.seed(3141) + +n_obs = 300 +data = make_simple_rdd_data(n_obs=n_obs, fuzzy=False) +df = pd.DataFrame( + np.column_stack((data['Y'], data['D'], data['score'], data['X'])), + columns=['y', 'd', 'score'] + ['x' + str(i) for i in range(data['X'].shape[1])] +) +dml_data = dml.DoubleMLData(df, y_col='y', d_cols='d', s_col='score') + +dml_rdflex = RDFlex(dml_data, ml_g=Lasso(), ml_m=LogisticRegression()) + + +def _assert_resampling_default_settings(dml_obj): + assert dml_obj.n_folds == 5 + assert dml_obj.n_rep == 1 + assert dml_obj.cutoff == 0 + assert dml_obj.h_fs is not None + assert dml_obj.fs_kernel == "triangular" + assert dml_obj.fuzzy is False + + +@pytest.mark.ci +def test_rdd_defaults(): + _assert_resampling_default_settings(dml_rdflex) diff --git a/doubleml/rdd/tests/test_rdd_exceptions.py b/doubleml/rdd/tests/test_rdd_exceptions.py new file mode 100644 index 000000000..3defef5da --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_exceptions.py @@ -0,0 +1,224 @@ +import pytest +import pandas as pd +import numpy as np +import copy + +from doubleml import DoubleMLData +from doubleml.rdd.datasets import make_simple_rdd_data +from doubleml.rdd import RDFlex + +from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin +from sklearn.linear_model import Lasso, LogisticRegression + +n = 500 +data = make_simple_rdd_data(n_obs=n, fuzzy=False) +df = pd.DataFrame( + np.column_stack((data['Y'], data['D'], data['score'], data['X'])), + columns=['y', 'd', 'score'] + ['x' + str(i) for i in range(data['X'].shape[1])] +) + +dml_data = DoubleMLData(df, y_col='y', d_cols='d', s_col='score') + +ml_g = Lasso() +ml_m = LogisticRegression() + + +# dummy learners for testing +class DummyRegressorNoSampleWeight(BaseEstimator, RegressorMixin): + """ + A dummy regressor that predicts the mean of the target values, + and does not support sample weights. + """ + def fit(self, X, y): + self.mean_ = np.mean(y) + return self + + def predict(self, X): + return np.full(shape=(X.shape[0],), fill_value=self.mean_) + + +class DummyClassifierNoSampleWeight(BaseEstimator, ClassifierMixin): + """ + A dummy classifier that predicts the most frequent class, + and does not support sample weights. + """ + def fit(self, X, y): + self.classes_, self.counts_ = np.unique(y, return_counts=True) + self.most_frequent_ = self.classes_[np.argmax(self.counts_)] + return self + + def predict(self, X): + return np.full(shape=(X.shape[0],), fill_value=self.most_frequent_) + + def predict_proba(self, X): + return np.column_stack( + (np.full(shape=(X.shape[0],), fill_value=1), + np.full(shape=(X.shape[0],), fill_value=0)) + ) + + +@pytest.mark.ci +def test_rdd_exception_data(): + # DoubleMLData + msg = r"The data must be of DoubleMLData type. \[\] of type was passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex([], ml_g) + + # score column + msg = 'Incompatible data. Score variable has not been set. ' + with pytest.raises(ValueError, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data._s_col = None + _ = RDFlex(tmp_dml_data, ml_g) + msg = 'Incompatible data. Score variable has to be continuous. ' + with pytest.raises(ValueError, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data._s = tmp_dml_data._d + _ = RDFlex(tmp_dml_data, ml_g) + + # existing instruments + msg = r'Incompatible data. x0 have been set as instrumental variable\(s\). ' + with pytest.raises(ValueError, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data._z_cols = ['x0'] + _ = RDFlex(tmp_dml_data, ml_g) + + # treatment exceptions + msg = ('Incompatible data. ' + 'To fit an RDFlex model with DML ' + 'exactly one binary variable with values 0 and 1 ' + 'needs to be specified as treatment variable.') + # multiple treatment variables + with pytest.raises(ValueError, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data._d_cols = ['d', 'x0'] + _ = RDFlex(tmp_dml_data, ml_g) + # non-binary treatment + with pytest.raises(ValueError, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data.x_cols = ['x1'] # reset x to only x1 to enable setting d to x0 + tmp_dml_data.d_cols = ['x0'] + _ = RDFlex(tmp_dml_data, ml_g) + + +@pytest.mark.ci +def test_rdd_exception_cutoff(): + msg = "Cutoff value has to be a float or int. Object of type passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, cutoff=[200]) + + msg = 'Cutoff value is not within the range of the score variable. ' + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, cutoff=200) + + +@pytest.mark.ci +def test_rdd_warning_fuzzy(): + msg = 'A sharp RD design is being estimated, but the data indicate that the design is fuzzy.' + with pytest.warns(UserWarning, match=msg): + _ = RDFlex(dml_data, ml_g, cutoff=0.1) + + +@pytest.mark.ci +def test_rdd_warning_treatment_assignment(): + msg = ("Treatment probability within bandwidth left from cutoff higher than right from cutoff.\n" + "Treatment assignment might be based on the wrong side of the cutoff.") + with pytest.warns(UserWarning, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data._s = -1.0*tmp_dml_data._s + _ = RDFlex(tmp_dml_data, ml_g, ml_m, fuzzy=True) + + +@pytest.mark.ci +def test_rdd_exception_learner(): + + # ml_g + msg = (r'The ml_g learner LogisticRegression\(\) was identified as classifier but the outcome variable is not' + ' binary with values 0 and 1.') + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g=LogisticRegression()) + msg = (r"The ml_g learner DummyRegressorNoSampleWeight\(\) does not support sample weights. Please choose a learner" + " that supports sample weights.") + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g=DummyRegressorNoSampleWeight(), ml_m=ml_m) + + # ml_m + msg = r'Invalid learner provided for ml_m: Lasso\(\) has no method .predict_proba\(\).' + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m=Lasso(), fuzzy=True) + msg = 'Fuzzy design requires a classifier ml_m for treatment assignment.' + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, fuzzy=True) + msg = (r"The ml_m learner DummyClassifierNoSampleWeight\(\) does not support sample weights. Please choose a learner" + " that supports sample weights.") + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m=DummyClassifierNoSampleWeight(), fuzzy=True) + + msg = ('A learner ml_m has been provided for for a sharp design but will be ignored. ' + 'A learner ml_m is not required for estimation.') + with pytest.warns(UserWarning, match=msg): + tmp_dml_data = copy.deepcopy(dml_data) + tmp_dml_data._data['sharp_d'] = (tmp_dml_data.s >= 0) + tmp_dml_data.d_cols = 'sharp_d' + _ = RDFlex(tmp_dml_data, ml_g, ml_m, fuzzy=False) + + +@pytest.mark.ci +def test_rdd_exception_resampling(): + # n_folds + msg = r"The number of folds must be of int type. \[1\] of type was passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, n_folds=[1]) + msg = 'The number of folds greater or equal to 2. 1 was passed.' + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, n_folds=1) + + # n_rep + msg = r"The number of repetitions for the sample splitting must be of int type. \[0\] of type was passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, n_rep=[0]) + msg = 'The number of repetitions for the sample splitting has to be positive. 0 was passed.' + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, n_rep=0) + + +@pytest.mark.ci +def test_rdd_exception_kernel(): + msg = "fs_kernel must be either a string or a callable. 2 of type was passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, fs_kernel=2) + msg = r"Invalid kernel 'rbf'. Valid kernels are \['uniform', 'triangular', 'epanechnikov'\]." + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, fs_kernel='rbf') + + +@pytest.mark.ci +def test_rdd_exception_h_fs(): + msg = "Initial bandwidth 'h_fs' has to be a float. Object of type passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, h_fs=1) + + +@pytest.mark.ci +def test_rdd_exception_fs_specification(): + msg = "fs_specification must be a string. 1 of type was passed." + with pytest.raises(TypeError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, fs_specification=1) + + msg = ("Invalid fs_specification 'local_constant'. " + r"Valid specifications are \['cutoff', 'cutoff and score', 'interacted cutoff and score'\].") + with pytest.raises(ValueError, match=msg): + _ = RDFlex(dml_data, ml_g, ml_m, fs_specification='local_constant') + + +@pytest.mark.ci +def test_rdd_exception_fit(): + rdd_model = RDFlex(dml_data, ml_g, ml_m) + msg = (r"The number of iterations for the iterative bandwidth fitting must be of int type. \[0\] of type " + "was passed.") + with pytest.raises(TypeError, match=msg): + rdd_model.fit(n_iterations=[0]) + + msg = 'The number of iterations for the iterative bandwidth fitting has to be positive. 0 was passed.' + with pytest.raises(ValueError, match=msg): + rdd_model.fit(n_iterations=0) diff --git a/doubleml/rdd/tests/test_rdd_fuzzy.py b/doubleml/rdd/tests/test_rdd_fuzzy.py new file mode 100644 index 000000000..6203abcf0 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_fuzzy.py @@ -0,0 +1,95 @@ +""" +Dummy test using fixed learner for fuzzy data +""" +import pytest +import numpy as np + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_fuzzy_data, cutoff): + return rdd_fuzzy_data(cutoff) + + +@pytest.fixture(scope='module') +def data_zero(rdd_fuzzy_data): + return rdd_fuzzy_data(0.0) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_fuzzy_left.py b/doubleml/rdd/tests/test_rdd_fuzzy_left.py new file mode 100644 index 000000000..2c9f887a0 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_fuzzy_left.py @@ -0,0 +1,95 @@ +""" +Dummy test using fixed learner for left sided fuzzy data +""" +import pytest +import numpy as np + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_fuzzy_left_data, cutoff): + return rdd_fuzzy_left_data(cutoff) + + +@pytest.fixture(scope='module') +def data_zero(rdd_fuzzy_left_data): + return rdd_fuzzy_left_data(0.0) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_fuzzy_right.py b/doubleml/rdd/tests/test_rdd_fuzzy_right.py new file mode 100644 index 000000000..f89d9d2a3 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_fuzzy_right.py @@ -0,0 +1,95 @@ +""" +Dummy test using fixed learner for right sided fuzzy data +""" +import pytest +import numpy as np + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_fuzzy_right_data, cutoff): + return rdd_fuzzy_right_data(cutoff) + + +@pytest.fixture(scope='module') +def data_zero(rdd_fuzzy_right_data): + return rdd_fuzzy_right_data(0.0) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/rdd/tests/test_rdd_return_types.py b/doubleml/rdd/tests/test_rdd_return_types.py new file mode 100644 index 000000000..418bc0c43 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_return_types.py @@ -0,0 +1,56 @@ +import pytest +import numpy as np +import pandas as pd + +import doubleml as dml +from doubleml.rdd import RDFlex +from doubleml.rdd.datasets import make_simple_rdd_data + +from sklearn.linear_model import Lasso, LogisticRegression + +np.random.seed(3141) + +n_obs = 300 +data = make_simple_rdd_data(n_obs=n_obs) +df = pd.DataFrame( + np.column_stack((data['Y'], data['D'], data['score'], data['X'])), + columns=['y', 'd', 'score'] + ['x' + str(i) for i in range(data['X'].shape[1])] +) +dml_data = dml.DoubleMLData(df, y_col='y', d_cols='d', s_col='score') + +dml_rdflex = RDFlex(dml_data, ml_g=Lasso(), ml_m=LogisticRegression()) + + +def _assert_return_types(dml_obj): + assert isinstance(dml_obj.n_folds, int) + assert isinstance(dml_obj.n_rep, int) + assert (isinstance(dml_obj.cutoff, float) | isinstance(dml_obj.cutoff, int)) + assert isinstance(dml_obj.fuzzy, bool) + assert isinstance(dml_obj.fs_kernel, str) + assert isinstance(dml_obj.w, np.ndarray) + assert isinstance(dml_obj.h, np.ndarray) + assert dml_obj.w.shape == (n_obs,) + assert dml_obj.w.dtype == float + assert isinstance(dml_obj.__str__(), str) + + +def _assert_return_types_after_fit(dml_obj): + assert isinstance(dml_obj.fit(), RDFlex) + assert isinstance(dml_obj.__str__(), str) + assert isinstance(dml_obj.n_folds, int) + assert isinstance(dml_obj.n_rep, int) + assert (isinstance(dml_obj.cutoff, float) | isinstance(dml_obj.cutoff, int)) + assert isinstance(dml_obj.fuzzy, bool) + assert isinstance(dml_obj.fs_kernel, str) + assert isinstance(dml_obj.w, np.ndarray) + assert isinstance(dml_obj.h, np.ndarray) + assert dml_obj.w.shape == (n_obs,) + assert dml_obj.w.dtype == float + assert isinstance(dml_obj.confint(), pd.DataFrame) + # TODO: Add Coefficient tests + + +@pytest.mark.ci +def test_rdd_returntypes(): + _assert_return_types(dml_rdflex) + _assert_return_types_after_fit(dml_rdflex) diff --git a/doubleml/rdd/tests/test_rdd_sharp.py b/doubleml/rdd/tests/test_rdd_sharp.py new file mode 100644 index 000000000..686276e99 --- /dev/null +++ b/doubleml/rdd/tests/test_rdd_sharp.py @@ -0,0 +1,95 @@ +""" +Dummy test using fixed learner for sharp data +""" +import pytest +import numpy as np + + +@pytest.fixture(scope='module', + params=[-0.2, 0.0, 0.4]) +def cutoff(request): + return request.param + + +@pytest.fixture(scope='module', + params=[0.05, 0.1]) +def alpha(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 4]) +def n_rep(request): + return request.param + + +@pytest.fixture(scope='module', + params=[1, 2]) +def p(request): + return request.param + + +@pytest.fixture(scope='module', + params=["cutoff", "cutoff and score", "interacted cutoff and score"]) +def fs_specification(request): + return request.param + + +@pytest.fixture(scope='module') +def data(rdd_sharp_data, cutoff): + return rdd_sharp_data(cutoff) + + +@pytest.fixture(scope='module') +def data_zero(rdd_sharp_data): + return rdd_sharp_data(0.0) + + +@pytest.fixture(scope='module') +def predict_placebo(predict_dummy, data_zero, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data_zero, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.fixture(scope='module') +def predict_nonplacebo(predict_dummy, data, cutoff, alpha, p, n_rep, fs_specification): + return predict_dummy( + data, cutoff=cutoff, alpha=alpha, n_rep=n_rep, p=p, fs_specification=fs_specification + ) + + +@pytest.mark.ci +def test_rdd_placebo_coef(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_coef(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['coef'], reference['coef'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_se(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_se(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['se'], reference['se'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_placebo_ci(predict_placebo): + reference, actual = predict_placebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) + + +@pytest.mark.ci +def test_rdd_nonplacebo_ci(predict_nonplacebo): + reference, actual = predict_nonplacebo + assert np.allclose(actual['ci'], reference['ci'], rtol=1e-9, atol=1e-4) diff --git a/doubleml/tests/test_dml_data.py b/doubleml/tests/test_dml_data.py index dd8e43f97..f9575d56a 100644 --- a/doubleml/tests/test_dml_data.py +++ b/doubleml/tests/test_dml_data.py @@ -450,11 +450,11 @@ def test_s_col_setter(): dml_data.s_col = 's_new' assert np.array_equal(dml_data.s, s_comp) - msg = r'Invalid selection variable s_col. a13 is no data column.' + msg = r'Invalid score or selection variable s_col. a13 is no data column.' with pytest.raises(ValueError, match=msg): dml_data.s_col = 'a13' - msg = (r'The selection variable s_col must be of str type \(or None\). ' + msg = (r'The score or selection variable s_col must be of str type \(or None\). ' "5 of type was passed.") with pytest.raises(TypeError, match=msg): dml_data.s_col = 5 @@ -607,19 +607,19 @@ def test_disjoint_sets(): with pytest.raises(ValueError, match=msg): _ = DoubleMLData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1', 'xx2'], z_cols='zz', t_col='zz') - msg = 'xx2 cannot be set as selection variable ``s_col`` and covariate in ``x_cols``.' + msg = 'xx2 cannot be set as score or selection variable ``s_col`` and covariate in ``x_cols``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1', 'xx2'], s_col='xx2') - msg = 'dd1 cannot be set as selection variable ``s_col`` and treatment variable in ``d_cols``.' + msg = 'dd1 cannot be set as score or selection variable ``s_col`` and treatment variable in ``d_cols``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1', 'xx2'], s_col='dd1') - msg = 'yy cannot be set as selection variable ``s_col`` and outcome variable ``y_col``.' + msg = 'yy cannot be set as score or selection variable ``s_col`` and outcome variable ``y_col``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1', 'xx2'], s_col='yy') - msg = 'zz cannot be set as selection variable ``s_col`` and instrumental variable in ``z_cols``.' + msg = 'zz cannot be set as score or selection variable ``s_col`` and instrumental variable in ``z_cols``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1', 'xx2'], z_cols='zz', s_col='zz') - msg = 'tt cannot be set as selection variable ``s_col`` and time variable ``t_col``.' + msg = 'tt cannot be set as score or selection variable ``s_col`` and time variable ``t_col``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1', 'xx2'], t_col='tt', s_col='tt') @@ -642,7 +642,7 @@ def test_disjoint_sets(): msg = 'xx2 cannot be set as time variable ``t_col`` and cluster variable in ``cluster_cols``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLClusterData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1'], t_col='xx2', cluster_cols='xx2') - msg = 'xx2 cannot be set as selection variable ``s_col`` and cluster variable in ``cluster_cols``.' + msg = 'xx2 cannot be set as score or selection variable ``s_col`` and cluster variable in ``cluster_cols``.' with pytest.raises(ValueError, match=msg): _ = DoubleMLClusterData(df, y_col='yy', d_cols=['dd1'], x_cols=['xx1'], s_col='xx2', cluster_cols='xx2') diff --git a/doubleml/utils/__init__.py b/doubleml/utils/__init__.py index 0439d692b..5281435ec 100644 --- a/doubleml/utils/__init__.py +++ b/doubleml/utils/__init__.py @@ -8,6 +8,7 @@ from .blp import DoubleMLBLP from .policytree import DoubleMLPolicyTree from .gain_statistics import gain_statistics +from .global_learner import GlobalClassifier, GlobalRegressor __all__ = [ "DMLDummyRegressor", @@ -16,5 +17,7 @@ "DoubleMLClusterResampling", "DoubleMLBLP", "DoubleMLPolicyTree", - "gain_statistics" + "gain_statistics", + "GlobalClassifier", + "GlobalRegressor" ] diff --git a/doubleml/utils/_checks.py b/doubleml/utils/_checks.py index d7d2881ed..8a546cb15 100644 --- a/doubleml/utils/_checks.py +++ b/doubleml/utils/_checks.py @@ -1,5 +1,6 @@ import numpy as np import warnings +import inspect from sklearn.utils.multiclass import type_of_target @@ -361,12 +362,30 @@ def _check_framework_compatibility(dml_framework_1, dml_framework_2, check_treat if dml_framework_1._is_cluster_data != dml_framework_2._is_cluster_data: raise ValueError('The cluster structure in DoubleMLFrameworks must be the same. ' f'Got {str(dml_framework_1._is_cluster_data)} and {str(dml_framework_2._is_cluster_data)}.') + return def _check_set(x): return {x} if x is not None else {} +def _check_resampling_specification(n_folds, n_rep): + if not isinstance(n_folds, int): + raise TypeError('The number of folds must be of int type. ' + f'{str(n_folds)} of type {str(type(n_folds))} was passed.') + if n_folds < 2: + raise ValueError('The number of folds greater or equal to 2. ' + f'{str(n_folds)} was passed.') + + if not isinstance(n_rep, int): + raise TypeError('The number of repetitions for the sample splitting must be of int type. ' + f'{str(n_rep)} of type {str(type(n_rep))} was passed.') + if n_rep < 1: + raise ValueError('The number of repetitions for the sample splitting has to be positive. ' + f'{str(n_rep)} was passed.') + return + + def _check_cluster_partitions(smpls, values): test_indices = np.concatenate([test_index for test_index in smpls]) if len(test_indices) != len(values): @@ -474,3 +493,10 @@ def _check_sample_splitting(all_smpls, all_smpls_cluster, dml_data, is_cluster_d smpls_cluster = None return smpls, smpls_cluster, n_rep, n_folds + + +def _check_supports_sample_weights(learner, learner_name): + if not ('sample_weight' in inspect.signature(learner.fit).parameters): + raise ValueError(f"The {learner_name} learner {str(learner)} does not support sample weights. " + "Please choose a learner that supports sample weights.") + return diff --git a/doubleml/utils/global_learner.py b/doubleml/utils/global_learner.py new file mode 100644 index 000000000..1a38062a3 --- /dev/null +++ b/doubleml/utils/global_learner.py @@ -0,0 +1,112 @@ +from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin, is_regressor, is_classifier, clone + +from sklearn.utils.multiclass import unique_labels + + +class GlobalRegressor(BaseEstimator, RegressorMixin): + """ + A global regressor that ignores the attribute `sample_weight` when being fit to ensure a global fit. + + Parameters + ---------- + base_estimator: regressor implementing ``fit()`` and ``predict()`` + Regressor that is used when ``fit()`` ``predict()`` and ``predict_proba()`` are being called. + """ + def __init__(self, base_estimator): + + if not is_regressor(base_estimator): + raise ValueError(f'base_estimator must be a regressor. Got {base_estimator.__class__.__name__} instead.') + + self.base_estimator = base_estimator + + def fit(self, X, y, sample_weight=None): + """ + Fits the regressor provided in ``base_estimator``. Ignores ``sample_weight``. + + Parameters + ---------- + X: array-like of shape (n_samples, n_features) + Training data. + + y: array-like of shape (n_samples,) or (n_samples, n_targets) + Target values. + + sample_weight: array-like of shape (n_samples,). + Individual weights for each sample. Ignored. + """ + self._fitted_learner = clone(self.base_estimator) + self._fitted_learner.fit(X, y) + + return self + + def predict(self, X): + """ + Predicts using the regressor provided in ``base_estimator``. + + Parameters + ---------- + X: array-like of shape (n_samples, n_features) + Samples. + """ + return self._fitted_learner.predict(X) + + +class GlobalClassifier(BaseEstimator, ClassifierMixin): + """ + A global classifier that ignores the attribute ``sample_weight`` when being fit to ensure a global fit. + + Parameters + ---------- + base_estimator: classifier implementing ``fit()`` and ``predict_proba()`` + Classifier that is used when ``fit()``, ``predict()`` and ``predict_proba()`` are being called. + """ + def __init__(self, base_estimator): + + if not is_classifier(base_estimator): + raise ValueError(f'base_estimator must be a classifier. Got {base_estimator.__class__.__name__} instead.') + + self.base_estimator = base_estimator + + def fit(self, X, y, sample_weight=None): + """ + Fits the classifier provided in ``base_estimator``. Ignores ``sample_weight``. + + Parameters + ---------- + X: array-like of shape (n_samples, n_features) + Training data. + + y: array-like of shape (n_samples,) or (n_samples, n_targets) + Target classes. + + sample_weight: array-like of shape (n_samples,). + Individual weights for each sample. Ignored. + """ + self.classes_ = unique_labels(y) + self._fitted_learner = clone(self.base_estimator) + self._fitted_learner.fit(X, y) + + return self + + def predict(self, X): + """ + Predicts using the classifier provided in ``base_estimator``. + + Parameters + ---------- + X: array-like of shape (n_samples, n_features) + Samples. + """ + return self._fitted_learner.predict(X) + + def predict_proba(self, X): + """ + Probability estimates using the classifier provided in ``base_estimator``. + The returned estimates for all classes are ordered by the label of classes. + + Parameters + ---------- + X: array-like of shape (n_samples, n_features) + Samples to be scored. + """ + return self._fitted_learner.predict_proba(X) diff --git a/doubleml/utils/tests/test_exceptions_global_learners.py b/doubleml/utils/tests/test_exceptions_global_learners.py new file mode 100644 index 000000000..ccd393222 --- /dev/null +++ b/doubleml/utils/tests/test_exceptions_global_learners.py @@ -0,0 +1,18 @@ +import pytest +from sklearn.linear_model import LogisticRegression, LinearRegression + +from doubleml.utils import GlobalRegressor, GlobalClassifier + + +@pytest.mark.ci +def test_global_regressor_input(): + msg = "base_estimator must be a regressor. Got LogisticRegression instead." + with pytest.raises(ValueError, match=msg): + _ = GlobalRegressor(base_estimator=LogisticRegression(random_state=42)) + + +@pytest.mark.ci +def test_global_classifier_input(): + msg = "base_estimator must be a classifier. Got LinearRegression instead." + with pytest.raises(ValueError, match=msg): + _ = GlobalClassifier(base_estimator=LinearRegression()) diff --git a/doubleml/utils/tests/test_global_learners.py b/doubleml/utils/tests/test_global_learners.py new file mode 100644 index 000000000..9cae0941d --- /dev/null +++ b/doubleml/utils/tests/test_global_learners.py @@ -0,0 +1,121 @@ +import pytest +import numpy as np +from doubleml.utils import GlobalRegressor, GlobalClassifier +from sklearn.base import clone +from sklearn.linear_model import LinearRegression, LogisticRegression +from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier + + +@pytest.fixture(scope='module', + params=[LinearRegression(), + RandomForestRegressor(n_estimators=10, max_depth=2, random_state=42)]) +def regressor(request): + return request.param + + +@pytest.fixture(scope='module', + params=[LogisticRegression(random_state=42), + RandomForestClassifier(n_estimators=10, max_depth=2, random_state=42)]) +def classifier(request): + return request.param + + +@pytest.fixture(scope="module") +def gl_fixture(regressor, classifier): + + global_reg = GlobalRegressor(base_estimator=regressor) + weighted_reg = clone(regressor) + unweighted_reg = clone(regressor) + + global_clas = GlobalClassifier(base_estimator=classifier) + weighted_clas = clone(classifier) + unweighted_clas = clone(classifier) + + X = np.random.normal(0, 1, size=(100, 10)) + y_con = np.random.normal(0, 1, size=(100)) + y_cat = np.random.binomial(1, 0.5, size=(100)) + sample_weight = np.random.random(size=(100)) + + # fit models + global_reg.fit(y=y_con, X=X, sample_weight=sample_weight) + weighted_reg.fit(y=y_con, X=X, sample_weight=sample_weight) + unweighted_reg.fit(y=y_con, X=X) + + global_clas.fit(y=y_cat, X=X, sample_weight=sample_weight) + weighted_clas.fit(y=y_cat, X=X, sample_weight=sample_weight) + unweighted_clas.fit(y=y_cat, X=X) + + global_reg_pred = global_reg.predict(X=X) + weighted_reg_pred = weighted_reg.predict(X=X) + unweighted_reg_pred = unweighted_reg.predict(X=X) + + global_clas_pred = global_clas.predict(X=X) + weighted_clas_pred = weighted_clas.predict(X=X) + unweighted_clas_pred = unweighted_clas.predict(X=X) + + global_clas_pred_proba = global_clas.predict(X=X) + weighted_clas_pred_proba = weighted_clas.predict(X=X) + unweighted_clas_pred_proba = unweighted_clas.predict(X=X) + + result_dict = { + "GlobalRegressor": global_reg, + "WeightedRegressor": weighted_reg, + "UnweightedRegressor": unweighted_reg, + "GlobalClassifier": global_clas, + "WeightedClassifier": weighted_clas, + "UnweightedClassifier": unweighted_clas, + "X": X, + "y_con": y_con, + "y_cat": y_cat, + "sample_weight": sample_weight, + "global_reg_pred": global_reg_pred, + "weighted_reg_pred": weighted_reg_pred, + "unweighted_reg_pred": unweighted_reg_pred, + "global_clas_pred": global_clas_pred, + "weighted_clas_pred": weighted_clas_pred, + "unweighted_clas_pred": unweighted_clas_pred, + "global_clas_pred_proba": global_clas_pred_proba, + "weighted_clas_pred_proba": weighted_clas_pred_proba, + "unweighted_clas_pred_proba": unweighted_clas_pred_proba + } + + return result_dict + + +@pytest.mark.ci +def test_predict(gl_fixture): + assert np.allclose(gl_fixture["global_reg_pred"], gl_fixture["unweighted_reg_pred"]) + assert np.allclose(gl_fixture["global_clas_pred"], gl_fixture["unweighted_clas_pred"]) + + assert not np.allclose(gl_fixture["global_reg_pred"], gl_fixture["weighted_reg_pred"]) + assert not np.allclose(gl_fixture["global_clas_pred"], gl_fixture["weighted_clas_pred"]) + + +@pytest.mark.ci +def test_predict_proba(gl_fixture): + assert np.allclose(gl_fixture["global_clas_pred_proba"], gl_fixture["unweighted_clas_pred_proba"]) + assert not np.allclose(gl_fixture["global_clas_pred_proba"], gl_fixture["weighted_clas_pred_proba"]) + + +@pytest.mark.ci +def test_clone(gl_fixture): + clone_reg = None + clone_clas = None + + try: + clone_reg = clone(gl_fixture["GlobalRegressor"]) + clone_clas = clone(gl_fixture["GlobalClassifier"]) + except Exception as e: + pytest.fail(f"clone() raised an exception:\n{str(e)}\n") + + # Test if they still work cloned + clone_reg.fit(y=gl_fixture["y_con"], X=gl_fixture["X"], sample_weight=gl_fixture["sample_weight"]) + clone_clas.fit(y=gl_fixture["y_cat"], X=gl_fixture["X"], sample_weight=gl_fixture["sample_weight"]) + + pred_global_reg = gl_fixture["GlobalRegressor"].predict(X=gl_fixture["X"]) + pred_global_clas = gl_fixture["GlobalClassifier"].predict_proba(X=gl_fixture["X"]) + pred_clone_reg = clone_reg.predict(X=gl_fixture["X"]) + pred_clone_clas = clone_clas.predict_proba(X=gl_fixture["X"]) + + np.allclose(pred_global_reg, pred_clone_reg) + np.allclose(pred_global_clas, pred_clone_clas) diff --git a/requirements.txt b/requirements.txt index 8aa7f9a95..99f30f016 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,5 @@ scipy scikit-learn statsmodels plotly -matplotlib \ No newline at end of file +matplotlib +rdrobust \ No newline at end of file