La experimentación es la base para probar hipótesis con rigor científico. En medicina, se utiliza para evaluar el efecto de nuevos tratamientos en los pacientes, mientras que en el mundo digital, gigantes tecnológicos como Amazonas, netflixy Úber Realiza miles de experimentos cada año para optimizar y mejorar sus plataformas.
Para experimentos a gran escala, se utiliza comúnmente la asignación aleatoria y se considera el “estándar de oro”. Con una cantidad sustancial de datos, la aleatoriedad tiende a producir grupos comparables, donde importantes factores previos al experimento están equilibrados y se cumple el supuesto de intercambiabilidad.
Sin embargo, cuando la muestra de un experimento es muy pequeña, la asignación aleatoria a menudo no logra crear grupos estadísticamente equivalentes. Entonces, ¿Cómo podemos dividir las unidades de manera eficiente entre los grupos de tratamiento y control?
Lo que aprenderás:
En esta publicación, explicaré un enfoque basado en optimización para construir grupos equivalentes para un experimento propuesto por Bertsimas et al. en Este artículo.
Con un ejemplo sencillo en Python aprenderemos a:
- Diseñar un experimento basado en optimización.
- Realizar inferencias utilizando técnicas bootstrap.
- Implemente el código en Python para sus propios experimentos.
Antes de profundizar en nuestro ejemplo y código Python, primero analicemos algunas ideas sobre los beneficios de utilizar un enfoque basado en optimización para diseñar experimentos.
La optimización hace que las estadísticas sean más precisas, lo que permite una inferencia más poderosa.
El enfoque basado en la optimización coincide con los grupos experimentales para minimizar las discrepancias masivas en medias y varianzas. Esto hace que las estadísticas sean mucho más precisas, concentrándolas estrechamente en torno a sus valores nominales y al mismo tiempo siendo estimaciones imparciales. Esta mayor precisión permite una inferencia más poderosa (usando el algoritmo bootstrap, que cubriremos más adelante).
Esto permite a los investigadores sacar conclusiones estadísticamente válidas con menos datos, lo que reduce los costos experimentales, un beneficio importante en disciplinas como la investigación oncológica, donde probar agentes quimioterapéuticos en modelos de cáncer de ratón es laborioso y costoso. Además, en comparación con otros métodos para tamaños de muestra pequeños, se ha demostrado que la optimización los supera, como veremos más adelante con experimentos simulados.
Ahora, ¡profundicemos en un ejemplo para ver cómo aplicar este enfoque usando Python!
Estudio de caso: un experimento farmacológico con 20 ratones
Consideremos un experimento con 20 ratones, en el que estamos interesados en estudiar el efecto de un fármaco sobre el crecimiento de tumores con diferentes tamaños iniciales de tumores. Supongamos que los tamaños iniciales de los tumores se distribuyen normalmente con una media de 200 mg y una desviación estándar de 300 mg (truncada para garantizar valores no negativos). Podemos generar la población de ratones con el siguiente código Python:
import numpy as np
import pandas as pd
import scipy.stats as statsdef create_experiment_data(n_mice, mu, sigma, seed):
lower, upper = 0, np.inf
initial_weights = stats.truncnorm(
(lower - mu) / sigma,
(upper - mu) / sigma,
loc=mu,
scale=sigma
).rvs(size=n_mice, random_state=seed)
return pd.DataFrame({
'mice': list(range(1, n_mice+1)),
'initial_weight': initial_weights,
})
tumor_data = create_experiment_data(n_mice=20, mu=200, sigma=300, seed=123)
print(tumor_data.head())
> mice initial_weight
0 1 424.736888
1 2 174.691035
2 3 141.016478
3 4 327.518749
4 5 442.239789
Ahora, necesitamos dividir a los 20 roedores en dos grupos de 10 cada uno: un grupo recibirá el tratamiento, mientras que el otro recibirá un placebo. Lo lograremos mediante la optimización.
También asumiremos que el crecimiento del tumor se observa durante un período de un día y sigue el modelo de Gompertz (detallado en Modelos matemáticos en la investigación del cáncer). Se supone que el tratamiento tiene un efecto determinista de reducir el tamaño del tumor en 250 mg.
Diseño experimental
Se pretende formular la creación de grupos equivalentes como un problema de optimización, donde el objetivo es minimizar la discrepancia tanto en la media como en la varianza. del peso inicial del tumor.
Para implementar esto, debemos seguir tres pasos:
Paso 1: Normalizar el peso inicial del tumor
Primero, se debe preprocesar toda la muestra y la métrica debe normalizarse para que tenga una media cero y una varianza unitaria:
mean = tumor_data['initial_weight'].mean()
std = tumor_data['initial_weight'].std()tumor_data['norm_initial_weight'] = (tumor_data['initial_weight'] - mean) / std
Paso 2: Crea los grupos usando la optimización.
A continuación, necesitamos implementar el modelo de optimización general que construye grupos m con k-sujetos cada uno, minimizando el discrepancia máxima entre dos grupos cualesquiera (una descripción completa de las variables del modelo se puede encontrar en el articulo) y pasándole la métrica normalizada:
El modelo matemático se puede implementar en Python usando la biblioteca ortools con el solucionador SCIP, de la siguiente manera:
from ortools.linear_solver import pywraplp
from typing import Unionclass SingleFeatureOptimizationModel:
"""
Implements the discrete optimization model proposed by Bertsimas et al. (2015) in "The Power of
Optimization Over Randomization in Designing Experiments Involving Small Samples".
See: https://doi.org/10.1287/opre.2015.1361.
"""
def __init__(self, data: pd.DataFrame, n_groups: int, units_per_group: int, metric: str, unit: str):
self.data = data.reset_index(drop=True)
self.parameters = {
'rho': 0.5,
'groups': range(n_groups),
'units': range(len(self.data)),
'unit_list': self.data[unit].tolist(),
'metric_list': self.data[metric].tolist(),
'units_per_group': units_per_group,
}
self._create_solver()
self._add_vars()
self._add_constraints()
self._set_objective()
def _create_solver(self):
self.model = pywraplp.Solver.CreateSolver("SCIP")
if self.model is None:
raise Exception("Failed to create SCIP solver")
def _add_vars(self):
self.d = self.model.NumVar(0, self.model.infinity(), "d")
self.x = {}
for i in self.parameters['units']:
for p in self.parameters['groups']:
self.x[i, p] = self.model.IntVar(0, 1, "")
def _set_objective(self):
self.model.Minimize(self.d)
def _add_constraints(self):
self._add_constraints_d_bounding()
self._add_constraint_group_size()
self._add_constraint_all_units_assigned()
def _add_constraints_d_bounding(self):
rho = self.parameters['rho']
for p in self.parameters['groups']:
for q in self.parameters['groups']:
if p < q:
self.model.Add(self.d >= self._mu(p) - self._mu(q) + rho * self._var(p) - rho * self._var(q))
self.model.Add(self.d >= self._mu(p) - self._mu(q) + rho * self._var(q) - rho * self._var(p))
self.model.Add(self.d >= self._mu(q) - self._mu(p) + rho * self._var(p) - rho * self._var(q))
self.model.Add(self.d >= self._mu(q) - self._mu(p) + rho * self._var(q) - rho * self._var(p))
def _add_constraint_group_size(self):
for p in self.parameters['groups']:
self.model.Add(
self.model.Sum([
self.x[i,p] for i in self.parameters['units']
]) == self.parameters['units_per_group']
)
def _add_constraint_all_units_assigned(self):
for i in self.parameters['units']:
self.model.Add(
self.model.Sum([
self.x[i,p] for p in self.parameters['groups']
]) == 1
)
def _add_contraint_symmetry(self):
for i in self.parameters['units']:
for p in self.parameters['units']:
if i < p:
self.model.Add(
self.x[i,p] == 0
)
def _mu(self, p):
mu = self.model.Sum([
(self.x[i,p] * self.parameters['metric_list'][i]) / self.parameters['units_per_group']
for i in self.parameters['units']
])
return mu
def _var(self, p):
var = self.model.Sum([
(self.x[i,p]*(self.parameters['metric_list'][i])**2) / self.parameters['units_per_group']
for i in self.parameters['units']
])
return var
def optimize(
self,
max_run_time: int = 60,
max_solution_gap: float = 0.05,
max_solutions: Union[int, None] = None,
num_threads: int = -1,
verbose: bool = False
):
"""
Runs the optimization model.
Args:
max_run_time: int
Maximum run time in minutes.
max_solution_gap: float
Maximum gap with the LP relaxation solution.
max_solutions: int
Maximum number of solutions until stop.
num_threads: int
Number of threads to use in solver.
verbose: bool
Whether to set the solver output.
Returns: str
The status of the solution.
"""
self.model.SetTimeLimit(max_run_time * 60 * 1000)
self.model.SetNumThreads(num_threads)
if verbose:
self.model.EnableOutput()
self.model.SetSolverSpecificParametersAsString(f"limits/gap = {max_solution_gap}")
self.model.SetSolverSpecificParametersAsString(f"limits/time = {max_run_time * 60}")
if max_solutions:
self.model.SetSolverSpecificParametersAsString(f"limits/solutions = {max_solutions}")
status = self.model.Solve()
if verbose:
if status == pywraplp.Solver.OPTIMAL:
print("Optimal Solution Found.")
elif status == pywraplp.Solver.FEASIBLE:
print("Feasible Solution Found.")
else:
print("Problem infeasible or unbounded.")
self._extract_solution()
return status
def _extract_solution(self):
tol = 0.01
self.assignment = {}
for i in self.parameters['units']:
for p in self.parameters['groups']:
if self.x[i,p].solution_value() > tol:
self.assignment.setdefault(p, []).append(self.parameters['unit_list'][i])
def get_groups_list(self):
return list(self.assignment.values())
model = SingleFeatureOptimizationModel(
data = tumor_data,
n_groups = 2,
units_per_group = 10,
unit = 'mice',
metric = 'norm_initial_weight',)
status = model.optimize()
optimized_groups = model.get_groups_list()
print(f"The optimized mice groups are: {optimized_groups}")
> The optimized mice groups are: [
[1, 4, 5, 6, 8, 12, 14, 16, 17, 18],
[2, 3, 7, 9, 10, 11, 13, 15, 19, 20]
]
Nota: El parámetro rho controla el equilibrio entre minimizar las discrepancias en el primer momento y el segundo momento y es elegido por el investigador. En nuestro ejemplo, hemos considerado rho es igual a 0,5.
Paso 3: Aleatorizar qué grupo recibe qué tratamiento
Por último, determinamos aleatoriamente qué grupo de ratones experimentales recibirá el fármaco y cuál recibirá el placebo:
import randomdef random_shuffle_groups(group_list, seed):
random.seed(seed)
random.shuffle(group_list)
return group_list
randomized_groups = random_shuffle_groups(optimized_groups, seed=123)
treatment_labels = ["Placebo", "Treatment"]
treatment_dict = {treatment_labels[i]: randomized_groups[i] for i in range(len(randomized_groups))}
print(f"The treatment assignment is: {treatment_dict}")
> The treatment assignment is: {
'Placebo': [2, 3, 7, 9, 10, 11, 13, 15, 19, 20],
'Treatment': [1, 4, 5, 6, 8, 12, 14, 16, 17, 18]
}
Veamos la calidad del resultado analizando la media y varianza de los pesos iniciales de los tumores en ambos grupos:
mice_assignment_dict = {inx: gp for gp, indices in treatment_dict.items() for inx in indices}
tumor_data['treatment'] = tumor_data['mice'].map(mice_assignment_dict)print(tumor_data.groupby('treatment').agg(
avg_initial_weight = ('initial_weight', 'mean'),
std_initial_weight = ('initial_weight', 'std'),
).round(2))
> avg_initial_weight std_initial_weight
treatment
Placebo 302.79 202.54
Treatment 303.61 162.12
Esta es la división del grupo con mínima discrepancia. en la media y varianza de los pesos de los tumores preexperimentales. Ahora realicemos el experimento y analicemos los resultados.
Simulando el crecimiento tumoral
Dada la asignación de tratamiento establecida, se simula el crecimiento del tumor durante un día utilizando el modelo de Gompertz, suponiendo un efecto de -250 mg para el grupo tratado:
import numpy as np
from scipy.integrate import odeint# Gompertz model parameters
a = 1
b = 5
# Critical weight
wc = 400
def gompex_growth(w, t):
"""
Gomp-ex differential equation model based on the initial weight.
"""
growth_rate = a + max(0, b * np.log(wc / w))
return w * growth_rate
def simulate_growth(w_initial, t_span):
"""
Simulate the tumor growth using the Gomp-ex model.
"""
return odeint(gompex_growth, w_initial, t_span).flatten()
def simulate_tumor_growth(data: pd.DataFrame, initial_weight: str, treatment_col: str, treatment_value: str, treatment_effect: float):
"""
Simulate the tumor growth experiment and return the dataset.
"""
t_span = np.linspace(0, 1, 2)
final_weights = np.array([simulate_growth(w, t_span)[-1] for w in data[initial_weight]])
experiment_data = data.copy()
mask_treatment = data[treatment_col] == treatment_value
experiment_data['final_weight'] = np.where(mask_treatment, final_weights + treatment_effect, final_weights)
return experiment_data.round(2)
experiment_data = simulate_tumor_growth(
data = tumor_data,
initial_weight = 'initial_weight',
treatment_col = 'treatment',
treatment_value = 'Treatment',
treatment_effect = -250
)
print(experiment_data.head())
> mice initial_weight norm_initial_weight treatment final_weight
0 1 424.74 0.68 Treatment 904.55
1 2 174.69 -0.72 Placebo 783.65
2 3 141.02 -0.91 Placebo 754.56
3 4 327.52 0.14 Treatment 696.60
4 5 442.24 0.78 Treatment 952.13
mask_tr = experiment_data.group == 'Treatment'
mean_tr = experiment_data[mask_tr]['final_weight'].mean()
mean_co = experiment_data[~mask_tr]['final_weight'].mean()
print(f"Mean difference between treatment and control: {round(mean_tr - mean_co)} mg")
> Mean difference between treatment and control: -260 mg
Ahora que tenemos los pesos finales de los tumores, observamos que el peso final promedio del tumor en el grupo de tratamiento es 260 mg menor que en el grupo de control. Sin embargo, para determinar si esta diferencia es estadísticamente significativa, debemos aplicar el siguiente mecanismo de arranque para calcular el valor p.
Inferencia bootstrap para diseño basado en optimización
“En un diseño basado en optimización, las estadísticas como la diferencia promedio entre grupos se vuelven más precisas pero ya no siguen las distribuciones habituales. Por lo tanto, se debe utilizar un método de inferencia de arranque para sacar conclusiones válidas”.
El método de inferencia boostrap propuesto por Bertsimas et al. (2015) Implica utilizar muestreo con reemplazo para construir la distribución de referencia de nuestro estimador. En cada iteración, la división del grupo se realiza mediante optimización y finalmente el valor p se deriva de la siguiente manera:
from tqdm import tqdm
from typing import Any, Listdef inference(data: pd.DataFrame, unit: str, outcome: str, treatment: str, treatment_value: Any = 1, n_bootstrap: int = 1000) -> pd.DataFrame:
"""
Estimates the p-value using bootstrap for two groups.
Parameters
-----------
data (pd.DataFrame): The experimental dataset with the observed outcome.
unit (str): The experimental unit column.
outcome (str): The outcome metric column.
treatment (str): The treatment column.
treatment_value (Any): The value referencing the treatment (other will be considered as control).
n_bootstrap (int): The number of random draws with replacement to use.
Returns
-----------
pd.DataFrame: The dataset with the results.
Raise
------------
ValueException: if there are more than two treatment values.
"""
responses = data[outcome].values
mask_tr = (data[treatment] == treatment_value).values
delta_obs = _compute_delta(responses[mask_tr], responses[~mask_tr])
deltas_B = _run_bootstrap(data, unit, outcome, n_bootstrap)
pvalue = _compute_pvalue(delta_obs, deltas_B)
output_data = pd.DataFrame({
'delta': [delta_obs],
'pvalue': [pvalue],
'n_bootstrap': [n_bootstrap],
'avg_delta_bootstrap': [np.round(np.mean(deltas_B), 2)],
'std_delta_bootstrap': [np.round(np.std(deltas_B), 2)]
})
return output_data
def _run_bootstrap(data: pd.DataFrame, unit: str, outcome: str, B: int = 1000) -> List[float]:
"""
Runs the bootstrap method and returns the bootstrapped deltas.
Parameters
-----------
data (pd.DataFrame): The dataframe from which sample with replacement.
outcome (str): The outcome metric observed in the experiment.
B (int): The number of random draws with replacement to perfrom.
Returns
-----------
List[float]: The list of bootstrap deltas.
"""
deltas_bootstrap = []
for i in tqdm(range(B), desc="Bootstrap Progress"):
sample_b = _random_draw_with_replacement(data, unit)
responses_b, mask_tr_b = _optimal_treatment_control_split(sample_b, unit, outcome, seed=i)
delta_b = _compute_delta(responses_b[mask_tr_b], responses_b[~mask_tr_b])
deltas_bootstrap.append(delta_b)
return deltas_bootstrap
def _compute_delta(response_tr, responses_co):
delta = np.mean(response_tr) - np.mean(responses_co)
return delta
def _compute_pvalue(obs_delta, bootstrap_deltas):
count_extreme = sum(1 for delta_b in bootstrap_deltas if abs(delta_b) >= abs(obs_delta))
p_value = (1 + count_extreme) / (1 + len(bootstrap_deltas))
return p_value
def _random_draw_with_replacement(data: pd.DataFrame, unit: str):
sample = data.sample(frac=1, replace=True)
sample[unit] = range(1, len(sample) + 1)
return sample
def _optimal_treatment_control_split(data: pd.DataFrame, unit: str, outcome: str, seed: int):
result = _sample(
data = data,
unit = unit,
normalized_feature = 'norm_initial_weight',
seed = seed
)
treatment_dict = {inx: gp for gp, indices in result.items() for inx in indices}
treatment = data[unit].map(treatment_dict)
mask_tr = (treatment == 'Treatment').values
responses = data[outcome].values
return responses, mask_tr
def _sample(data: pd.DataFrame, unit: str, normalized_feature: str, seed: int):
model = SingleFeatureOptimizationModel(
data,
n_groups = 2,
units_per_group = 10,
unit = unit,
metric = normalized_feature,
)
status = model.optimize()
optimized_groups = model.get_groups_list()
randomized_groups = random_shuffle_groups(optimized_groups, seed=seed)
treatment_labels = ["Placebo", "Treatment"]
return {treatment_labels[i]: randomized_groups[i] for i in range(len(randomized_groups))}
infer_result = inference(
data = experiment_data,
unit = 'mice',
outcome = 'final_weight',
treatment = 'group',
treatment_value = 'Treatment',
n_bootstrap = 1000
)print(infer_result)
> delta pvalue n_bootstrap avg_delta_bootstrap std_delta_bootstrap
0 -260.183 0.001998 1000 2.02 112.61
La diferencia observada de -260 mg entre los grupos es significativa al nivel de significancia del 5% (el valor de p inferior a 0,05). Por lo tanto, rechazamos la hipótesis nula de igualdad de medias y concluimos que el tratamiento tuvo un efecto estadísticamente significativo.
Podemos simular el experimento varias veces, generando poblaciones de ratones con diferentes pesos tumorales iniciales, extraídas de la misma distribución normal con una media de 200 mg y una desviación estándar de 300 mg.
Esto nos permite comparar el diseño basado en optimización con otros diseños experimentales. En el siguiente gráfico, comparo el enfoque de optimización con la asignación aleatoria simple y la asignación aleatoria estratificada (donde los estratos se crearon utilizando un algoritmo de k-medias basado en el peso inicial del tumor):
En su artículo, los autores también comparan el enfoque basado en la optimización con la nueva aleatorización y el emparejamiento por pares entre diferentes tamaños de efectos y tamaños de grupos. ¡Recomiendo encarecidamente leer el artículo completo si está interesado en explorar más los detalles!