import numpy as np
import pandas as pd
from scipy.stats import t as t_dist


def simulate_exp_growth(
    group_name: str,
    duration: int,
    initial_tumor_volume: float,
    group_growth_rate_day_pct: float,
    growth_rate_cv: float,
    residual_cv: float,
    n_animals: int,
    max_tumor_volume: float,
    df=None
) -> pd.DataFrame:
    """
    Simulates tumor growth data for a group of animals over a given duration.

    Args:
        group_name (str): Name of the group (e.g., "control", "therapy_a").
        duration (int): Duration of the study in days.
        initial_tumor_volume (float): Initial tumor volume (mm^3).
        group_growth_rate_day_pct (float): Average daily growth rate of the group (%).
        growth_rate_cv (float): Coefficient of variation for the growth rate within the group.
        residual_cv (float): Coefficient of variation for the measurement error.
        n_animals (int): Number of animals in the group.
        max_tumor_volume (float): Maximum allowed tumor volume before removal (mm^3).
        df (int, optional): Degrees of freedom for the t-distribution. If None,
        defaults to n_animals - 1.

    Returns:
        pandas.DataFrame: Simulated data with columns: day, animal, group, tv (Tumor Volume).
    """

    # Initialize an empty list to store DataFrames for each animal in the group
    group_data = []

    # Loop over each animal in the group
    for i in range(n_animals):
        # Create an array of time points (days) when measurements are taken, starting from
        # day 0 up to the end of the study (`duration`). Measurements are taken every 3 days.
        days = np.arange(0, duration + 1, 3)

        # Simulate individual growth rates
        #
        # individual_growth_rate = np.random.normal(
        #     group_growth_rate_day_pct / 100,
        #     group_growth_rate_day_pct * growth_rate_cv / 100,
        #     len(days) - 1
        # )
        #
        # Since we're simulating tumor growth data, and the sample size is usually small
        # (below 15 animals per group), using the t distribution could be a reasonable choice.
        # It will introduce slightly more variability into the simulated data, reflecting the
        # uncertainty associated with small samples and the stochastic nature of tumor growth.

        # Shift and scale the t-distribution for individual growth rate
        df = n_animals - 1 if df is None else df  # Default df if not provided
        individual_growth_rate = t_dist.rvs(
            df=df,
            loc=group_growth_rate_day_pct / 100,
            scale=group_growth_rate_day_pct * growth_rate_cv / 100,
            size=len(days) - 1
        )

        # Note: if the growth rates are strictly positive and tend to be right-skewed (some
        # tumors grow much faster than others), the gamma distribution can be a good option.

        # Initialize an empty array to store the tumor volumes over time
        # ln_tumor_volume = np.empty(len(days))
        # we're now working with log-normal tumor volumes directly (instead of log-transformed values)
        tumor_volume = np.empty(len(days))

        # Simulate the initial log-transformed tumor volume
        #
        # If the initial tumor sizes are tightly controlled and have low variability,
        # a normal distribution might be suitable. Log-normal or truncated normal
        # distributions could also be considered.
        tumor_volume[0] = np.random.normal(
            loc=initial_tumor_volume,
            scale=initial_tumor_volume * residual_cv
        )

        # Iterate over time points (days) to simulate tumor growth
        for t in range(1, len(days)):
            # ln_tumor_volume[t] = np.random.normal(
            #     ln_tumor_volume[t - 1] + individual_growth_rate[t - 1] * days[t],
            #     np.sqrt(np.log(residual_cv**2 + 1))
            # )
            #
            # the log-normal distribution is a popular choice for modeling tumor volumes
            # over time, as it captures the right-skewed nature of tumor growth data.
            #
            # Calculate mean and standard deviation of the log-normal distribution
            # Apply exponential growth to the previous tumor volume (not log-transformed)
            mean = tumor_volume[t-1] * \
                np.exp(individual_growth_rate[t-1] * days[t])
            # Calculate the mean and standard deviation of the underlying normal
            # distribution for log-normal
            mean_log = np.log(mean)
            # With std_log = sqrt(ln(CV² + 1))
            std_log = np.sqrt(np.log(residual_cv**2 + 1))
            # Simulate tumor volume at time t directly from the log-normal distribution
            # The `np.random.lognormal` function takes the mean and standard deviation
            # of the underlying normal distribution as input. Remember, if a variable X
            # follows a log-normal distribution, then the logarithm of that variable,
            # ln(X), follows a normal distribution.
            tumor_volume[t] = np.random.lognormal(
                mean=mean_log,
                sigma=std_log,
            )
            # Considere Gompertz Growth Model that describes sigmoidal growth patterns,
            # which are often observed in tumors.

        # # Transform log-tumor volumes back to the original scale (tumor volumes in mm^3).
        # tumor_volume = np.exp(ln_tumor_volume)

        # Find indices where tumor volume exceeds the maximum allowed.
        welfare_missing = np.where(tumor_volume > max_tumor_volume)[0]
        if len(welfare_missing) > 1:
            # Exclude from its first occurrence until the next timepoint for survival .
            tumor_volume[welfare_missing[1:]] = np.nan

        # Create a DataFrame for the current animal's data
        animal_data = pd.DataFrame(
            {
                "day":    days,
                "animal": f"{group_name}_{i + 1}",
                "group":  group_name,
                "tv":     tumor_volume,
            }
        )

        # Add the animal's data to the list of group data
        group_data.append(animal_data)

    # Concatenate all animal dataframes into one for the entire group
    # ignore_index=True resets the index of the resulting DataFrame.
    return pd.concat(group_data, ignore_index=True)

import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm

# Initiating values for tumor_volume and residual_cv
tumor_volume = 600
residual_cv = 0.25

# Calculate parameters for the underlying normal distribution of the log-normal
mean_log = np.log(tumor_volume)
std_log = np.sqrt(np.log(residual_cv**2 + 1))

# Calculate the mean and SD of the equivalent normal distribution
mean_normal = np.exp(mean_log + std_log**2 / 2)
std_normal = np.sqrt((np.exp(std_log**2) - 1) *
                     np.exp(2*mean_log + std_log**2))
# When we use std = initial_tumor_volume * residual_cv, we're calculating the standard deviation
# directly from the initial tumor volume and the coefficient of variation. This is correct if
# residual_cv is the CV of the normal distribution. However, in the original simulation setup,
# residual_cv represents the coefficient of variation of the log-normal distribution,
# not the normal distribution.

# Create x values for plotting
x = np.linspace(0, 2000, 1000)

# Calculate PDF values for normal and log-normal distributions
normal_pdf = norm.pdf(x, loc=mean_normal, scale=std_normal)
lognormal_pdf = lognorm.pdf(x, s=std_log, scale=np.exp(mean_log))

# Plot the PDFs
plt.plot(x, normal_pdf, label='Normal distribution')
plt.plot(x, lognormal_pdf, label='Log-normal distribution')
plt.xlabel('Tumor volume')
# plt.xscale("log")
plt.ylabel('Probability density')
plt.title('Normal vs. log-normal distribution of tumor volumes')
plt.legend()
plt.grid(axis='y', alpha=0.75)

# Set a seed for reproducibility
np.random.seed(444)

# Simulate the vehicle (control) group data
vehicle_group = simulate_exp_growth(
    group_name="control",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=5,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Simulate data for the single-agent therapy A group
therapy_a_group = simulate_exp_growth(
    group_name="therapy_a",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=4,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Simulate data for the single-agent therapy B group
therapy_b_group = simulate_exp_growth(
    group_name="therapy_b",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=3,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Simulate data for the combination therapy A+B group
therapy_ab_group = simulate_exp_growth(
    group_name="therapy_combo_ab",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=1.3,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Combine all group data into a single DataFrame
TGI_dataset_test = pd.concat(
    [
        vehicle_group,
        therapy_a_group,
        therapy_b_group,
        therapy_ab_group
    ],
    ignore_index=True
)

# Let's have a look at the dataset structure and values
print(TGI_dataset_test.head())

   day     animal    group          tv
0    0  control_1  control  103.815784
1    3  control_1  control  128.940403
2    6  control_1  control  176.041978
3    9  control_1  control  209.166393
4   12  control_1  control  442.556814

# For comparison with the R tutorial we can import the original TGI.dataset
TGI_dataset = pd.read_csv("TGI_dataset.csv")
TGI_dataset.columns = [col.lower() for col in TGI_dataset.columns]
for col in TGI_dataset.select_dtypes(include='object').columns:
    TGI_dataset[col] = TGI_dataset[col].astype(str).str.lower()
TGI_dataset.head(6)

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.simplefilter(action='ignore', category=DeprecationWarning)

def summary_stats_tgi(dataset: pd.DataFrame) -> pd.DataFrame:
    """Calculates geometric mean, standard error, and 95% confidence intervals 
        for tumor volumes per group per day.

    Args:
        dataset (pd.DataFrame): TGI dataset with columns 'day', 'group', and 'tv' (tumor volume).

    Returns:
        pd.DataFrame: Summary statistics with columns 'day', 'group', 'est' (geometric mean), 
        'se' (standard error), 'lb95_est' (lower 95% CI bound), and 'ub95_est' (upper 95% CI bound).
    """

    # Split the dataset by group and day
    grouped_data = dataset.groupby(['group', 'day'])

    # Apply the summary statistics calculation to each group
    summary_list = []
    for (group, day), group_df in grouped_data:
        # First, we take the natural logarithm of the tumor volumes
        # Drop missing values (NaN)
        log_tumor_volumes = np.log(group_df['tv'].dropna())

        # Analyze the log-transformed values using standard normal distribution techniques
        # To get back to the original scale (tumor volumes), we exponentiate the mean and confidence
        # interval limits. Since the exponential is the inverse of the logarithm, this back-transformation
        # gives us the geometric mean and the corresponding confidence interval on the original scale.
        geometric_mean = np.exp(np.mean(log_tumor_volumes))
        # In R, the default behavior of the var() function is to use a delta degrees of freedom (ddof) of 1.
        # This means it calculates the sample variance, which is considered an unbiased estimator of the
        # population variance. NumPy's np.var function, by default, uses ddof = 0 (population variance),
        # but not pandas's var.
        se = np.sqrt(log_tumor_volumes.var() / len(log_tumor_volumes))
        lb95 = geometric_mean * np.exp(norm.ppf(0.025) * se)
        ub95 = geometric_mean * np.exp(norm.ppf(0.975) * se)

        summary_list.append({
            "day": day,
            "group": group,
            "est": geometric_mean,
            "se": se,
            "lb95_est": lb95,
            "ub95_est": ub95
        })

    # Combine results into a DataFrame
    return pd.DataFrame(summary_list)

TGI_stats = summary_stats_tgi(dataset=TGI_dataset)
print(TGI_stats)

    day             group          est        se     lb95_est     ub95_est
0     0           control   101.536731  0.034077    94.976633   108.549939
1     3           control   132.261775  0.061715   117.193124   149.267948
2     6           control   184.300154  0.064576   162.389861   209.166672
3     9           control   281.463182  0.102086   230.422992   343.809105
4    12           control   477.111360  0.125217   373.280079   609.824264
5    15           control  1169.973866  0.153067   866.732241  1579.309945
6    18           control  2750.176428  0.195191  1875.917983  4031.876902
7     0         therapy_a   103.594835  0.067338    90.786160   118.210638
8     3         therapy_a   105.617074  0.100664    86.705978   128.652794
9     6         therapy_a   136.011443  0.113316   108.923209   169.836279
10    9         therapy_a   195.629951  0.127840   152.271268   251.334860
11   12         therapy_a   312.601850  0.145064   235.240926   415.403554
12   15         therapy_a   551.161945  0.132813   424.843493   715.038583
13   18         therapy_a  1137.049884  0.128526   883.849561  1462.785632
14    0         therapy_b    98.736314  0.061430    87.536227   111.369431
15    3         therapy_b   108.043137  0.105557    87.851067   132.876241
16    6         therapy_b   147.498501  0.108302   119.289214   182.378666
17    9         therapy_b   199.556878  0.133412   153.640644   259.195397
18   12         therapy_b   270.856421  0.164438   196.231703   373.860083
19   15         therapy_b   443.544059  0.191394   304.804928   645.433567
20   18         therapy_b   677.167458  0.182431   473.599020   968.236306
21    0  therapy_combo_ab    96.364433  0.062295    85.288652   108.878542
22    3  therapy_combo_ab    97.289897  0.081205    82.974796   114.074690
23    6  therapy_combo_ab    99.155402  0.132551    76.469620   128.571237
24    9  therapy_combo_ab   111.844686  0.154870    82.563880   151.509761
25   12  therapy_combo_ab   136.879392  0.176595    96.832279   193.488869
26   15  therapy_combo_ab   178.214612  0.170720   127.533986   249.035170
27   18  therapy_combo_ab   225.100813  0.195905   153.328348   330.469719

TGI_dataset.groupby(['group', 'day']).describe()

from scipy.stats import gmean

TGI_dataset.pivot_table(
    values='tv',
    index='day',
    columns='group',
    aggfunc=gmean)

def plot_tgi_stats(
    summ_stats: pd.DataFrame,
    ln_scale: bool = False,
    colormap: str = "viridis"
) -> None:
    """Plots tumor growth curves with 95% confidence intervals.

    Args:
        summ_stats (pd.DataFrame): Dataframe with summary statistics (output of `summary_stats_tgi`).
        ln_scale (bool, optional): If True, plot on natural log scale. Defaults to False.
        colormap (str, optional): Colormap to use for group colors. Defaults to "viridis".
    """

    # Map groups to colors and symbols
    groups = summ_stats["group"].unique()
    n_groups = len(groups)

    # Map groups to colors from the colormap, and markers
    colors = plt.get_cmap(colormap)(np.linspace(0.25, 1, n_groups))

    # Different markers for each group
    markers = ['o', 's', '^', 'D', '*', 'v', 'x', 'p']

    # Set up plot
    fig, ax = plt.subplots()

    # Plot each group
    for i, group in enumerate(groups):
        group_data = summ_stats[summ_stats["group"] == group]
        ax.errorbar(
            group_data["day"],
            group_data["est"],
            yerr=[
                group_data["est"] - group_data["lb95_est"],
                group_data["ub95_est"] - group_data["est"]],
            fmt=markers[i],
            color=colors[i],
            label=group,
            capsize=5)

        ax.plot(
            group_data["day"],
            group_data["est"],
            color=colors[i],
            linewidth=2)

    # Formatting
    ax.set_xlabel("Days", fontsize=12)
    # Ensure x-axis ticks are integers
    ax.set_xticks(group_data['day'].unique())
    ax.set_xticklabels(group_data['day'].unique(), fontsize=10)
    if ln_scale:
        ax.set_yscale("log")
    ax.set_ylabel("Tumor volume (mm³)", fontsize=12)
    ax.set_title("Tumor growth curves", fontsize=14)
    ax.grid(which='major', axis='y', linestyle='--')
    ax.legend(loc='upper left', fontsize=10)

plot_tgi_stats(TGI_stats, ln_scale=False, colormap='viridis')

import seaborn as sns

# Create line plot of tumor volume over time for each group, using geometric mean
ax = sns.lineplot(
    # Input DataFrame containing tumor growth data, not the summary statistics
    data=TGI_dataset,
    x='day',                       # Column name for time points (x-axis)
    y='tv',                        # Column name for tumor volumes (y-axis)
    # Column name for grouping variable (different lines/colors for each group)
    hue='group',
    # Column name for assigning different line styles to groups (automatic)
    style='group',
    # Disable automatic dashes for different line styles (we'll use markers)
    dashes=False,
    markers=["o", "s", "^", "D"],  # Manually assign markers for each group
    estimator=gmean,               # Use geometric mean as the estimator for central tendency
    # Show 95% confidence intervals as error bars
    errorbar=('ci', 95),
    palette='viridis',             # Color palette to use for different groups
    # Show confidence intervals as filled bands around the lines
    err_style='band')

# Optionally plot on a logarithmic scale
# plt.yscale('log')
# Ensure x-axis ticks are integers
ax.set_xticks(TGI_dataset['day'].unique())
ax.set_xticklabels(TGI_dataset['day'].unique())
ax.set_xlabel("Days", fontsize=12)
ax.set_ylabel("Tumor volume (mm³)", fontsize=12)
ax.set_title("Tumor growth curves", fontsize=14)
ax.grid(which='major', axis='y', linestyle='--')
ax.legend(loc='upper left', fontsize=10);

def plot_individual_curves(
    dataset: pd.DataFrame,
    colormap: str = "viridis",
    max_cols: int = 2
) -> None:
    """Plots individual tumor growth curves for each animal, with one subplot per group.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'tv' (tumor volume), 'group', and 'animal'.
        colormap (str, optional): Colormap to use for group colors. Defaults to "viridis".
        max_cols (int, optional): Maximum number of columns in the plot grid. Defaults to 2.
    """

    groups = dataset["group"].unique()
    n_groups = len(groups)

    # Map groups to colors from the colormap
    colors = plt.get_cmap(colormap)(np.linspace(0.25, 1, n_groups))

    # Calculate number of rows and columns for subplots
    num_cols = min(n_groups, max_cols)
    num_rows = (n_groups + num_cols - 1) // num_cols

    # Find the maximum tumor volume across all groups
    max_tumor_volume = dataset['tv'].max()

    # Create the figure and axes
    fig, axes = plt.subplots(
        nrows=num_rows,
        ncols=num_cols,
        # Adjust figsize based on num_cols
        figsize=(3 * num_cols, 3 * num_rows),
        sharey=True,
    )

    # Flatten the axes array for easier iteration
    axes = axes.ravel()

    # Plot individual curves for each group
    for i, group in enumerate(groups):
        ax = axes[i]
        group_data = dataset[dataset["group"] == group]

        for animal in group_data["animal"].unique():
            animal_data = group_data[group_data["animal"] == animal]
            ax.plot(
                animal_data["day"],
                animal_data["tv"],
                color=colors[i],
                label=animal,
                lw=2)

        # Formatting
        ax.set_title(group)
        ax.set_xlabel("Days")
        ax.set_ylabel("Tumor volume (mm3)")
        ax.set_ylim(0, max_tumor_volume * 1.1)  # Set y-axis limits

    # Hide any unused subplots if the number of groups is not even
    for j in range(n_groups, num_cols * num_rows):
        axes[j].set_visible(False)

    # Add overall title and adjust layout
    plt.suptitle("Individual tumor growth curves")
    plt.tight_layout()
    plt.show()

# Plot the current TGI_dataset
plot_individual_curves(TGI_dataset_test)

def calculate_tgi_from_data(
    dataset: pd.DataFrame,
    type: str = 'TGI',
    time_points: list = None,
    baseline: int = 0,
    treatment_name: str = None,
    control_name: str = "control"
) -> pd.DataFrame:
    """Calculates either TGI or TGI_delta from a tumor growth dataset.

    Args:
        dataset (pd.DataFrame): Dataframe with columns 'day', 'group', and 'tv' (Tumor Volume).
        type (str, optional): Types of TGI to calculate ('TGI' or 'TGI_delta'). Defaults to "TGI".
        time_points (list, optional): List of days at which to calculate TGI. If None, calculates for all days.
        baseline (int, optional): Baseline day for TGI_delta calculation. Required if 'TGI_delta' is in `type`.
        treatment_name (str, optional): Name of the treatment group. Defaults to None (calculates 
        for all groups vs. control).
        control_name (str, optional): Name of the vehicle (control) group. Defaults to "control".

    Returns:
        pd.DataFrame: DataFrame containing TGI or TGI_delta results.

    Raises:
        ValueError: If no treatment_name provided and no other group than control found in the dataset
    """

    # Ensure treatment_name is specified
    if not treatment_name:  # not None returns True
        non_control_groups = dataset[dataset['group']
                                     != control_name]['group'].unique()
        if len(non_control_groups) == 0:
            raise ValueError("No treatment groups found in the dataset.")
        # otherwise we take the first non-control group from the list as the treatment group
        treatment_name = non_control_groups[0]

    # Add baseline for calculating TGI_delta
    time_points = [baseline] + \
        time_points if type == "TGI_delta" else time_points

    # Filter control and treatment data
    control_data = dataset[(dataset["group"] == control_name) & (
        dataset["day"].isin(time_points))]
    treatment_data = dataset[(dataset["group"] == treatment_name) & (
        dataset["day"].isin(time_points))]

    # Calculate means
    control_means = control_data.groupby("day")["tv"].mean().reset_index()
    treatment_means = treatment_data.groupby("day")["tv"].mean().reset_index()

    # Calculate TGI or TGI_delta
    if type == "TGI":
        result_data = control_means.merge(
            treatment_means, on="day", suffixes=("_control", "_treatment"))
        result_data["tgi_pct"] = (
            1 - (result_data["tv_treatment"] / result_data["tv_control"])) * 100
        result_data["contrast"] = f"{treatment_name} vs. {control_name}"
        return result_data[["contrast", "day", "tgi_pct"]]

    elif type == "TGI_delta":
        control_baseline = control_data[control_data["day"] == baseline]
        treatment_baseline = treatment_data[treatment_data["day"] == baseline]

        # Calculate mean baseline tumor volumes for both groups
        control_baseline_mean = control_baseline["tv"].mean()
        treatment_baseline_mean = treatment_baseline["tv"].mean()

        # Adjust tumor volumes relative to their respective baseline means
        control_means["tv_control_adjusted"] = control_means["tv"] - \
            control_baseline_mean
        treatment_means["tv_treatment_adjusted"] = treatment_means["tv"] - \
            treatment_baseline_mean

        # Merge and calculate TGI_delta
        result_data = control_means.merge(treatment_means, on="day")
        result_data["tgi_delta_pct"] = (
            1 - (result_data["tv_treatment_adjusted"] / result_data["tv_control_adjusted"])) * 100
        result_data["contrast"] = f"{treatment_name} vs. {control_name}"

        # mask to exclude baseline from the report
        ix = result_data['day'] == baseline
        return result_data.loc[~ix, ["contrast", "day", "tgi_delta_pct"]]

    else:
        raise ValueError(
            "Invalid type. Please specify either 'TGI' or 'TGI_delta'.")

# Example 1 - Single time point TGI. Imagine we want to assess the efficacy of
# therapy_b compared to the control group at day 18 of the study.
calculate_tgi_from_data(
    dataset=TGI_dataset,
    type='TGI',
    time_points=[18, 21, 24, 27],  # it doesn't hurt if days don't exist
    treatment_name="therapy_b",
    control_name='control')

# Example 2 - TGI-delta over multiple time points. If we're interested in
# how the effect of therapy_B changes over time, we can calculate TGI-delta
# at multiple time points relative to a baseline (e.g., day 0).
calculate_tgi_from_data(
    dataset=TGI_dataset,
    type='TGI_delta',
    time_points=[6, 12, 18],
    # baseline=0, # default
    treatment_name="therapy_b",
    # control_name='control'  # default name of the control group
)

# we apply arithmetic 'mean' and scipy gmean to tumor volumes
# grouped by treatment group and day
TGI_dataset.groupby(['group', 'day'])['tv'].agg(['mean', gmean])
#TGI_dataset.groupby(['group', 'day']).agg({'tv': ['mean', gmean]})

# TGI at multiple time-points
# use .xs to access values from each group in the hierarchical index
(
    TGI_dataset
    .groupby(['group', 'day'])
    .tv  # equivalent to ['tv']
    .mean()
    .transform(lambda x: 1 - x / x.xs('control', level=0))
    # pull out the outer-most index and push it up into columns
    .unstack(level=0)
)

# Adjusted geometric means of the tumor volume
(
    TGI_dataset
    .groupby(['group', 'day'])
    .agg(adjusted_geometric_mean=('tv', gmean))  # named aggregation
    .transform(lambda x: x - x.xs(0, level=1))
    .unstack(level='group')
)

# TGI-delta (arithmetic mean) over all groups and time points
(
    TGI_dataset
    .groupby(['group', 'day'])
    ['tv']
    .agg(TGI_delta_arithmetic_mean='mean')
    .transform(lambda x: x - x.xs(0, level=1))
    .transform(lambda x: 1 - x/x.xs('control', level=0))
    .unstack(level='group')
)

Data
|
|-- Categorical (compare frequencies)
|   |-- Chi-squared tests (trend is not important)
|   `-- Cochran-Armitage (trend is mportant)
|
|-- Continuous (typical TGI)
|   |
|   |-- Single time-point (usually last day of study)
|   |   |-- T-tests or ANOVAs (with/without adjustment for covariates)
|   |   `-- Bliss or HSA (analysis of synergy)
|   |
|   `-- Multiple time-points (several assessments throughout study)
|       |-- Growth rates or AUCs (summarizes into a few parameters)
|       |-- T-tests, ANOVAs, or ANCOVAs (with/without covariates, assessment at each time point)
|   |   `-- Bliss or HSA (analysis of synergy)
|
`-- Time-to-event (modeling events)
    |-- Log-rank tests (non-parametric)
    `-- Cox Proportional Hazards models (semi-parametric with/without adjustment for covariates)

from typing import List
import statsmodels.formula.api as smf

def tgi_single_time(
    dataset: pd.DataFrame,
    time_point: int,
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    baseline_adjusted: bool = False,
    baseline: int = 0,
    p_tails: int = 1,
    correction_method: str = 'none',
) -> pd.DataFrame:
    """Estimates Tumor Growth Inhibition (TGI) at a single time point using a linear model.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', and 'tv' (tumor volume).
        time_point (int): Day at which to calculate TGI.
        treatments_numerator (List[str]): List of treatment group names (numerators in TGI calculation).
        treatment_denominator (str): Name of the comparison group (denominator in TGI calculation).
        Defaults to 'control'.
        baseline_adjusted (bool, optional): If True, adjust for baseline tumor volume. Defaults to False.
        baseline (int, optional): Baseline day for adjustment. Defaults to 0.
        p_tails (int, optional): Number of tails for the hypothesis test (1 for one-sided, 2 for two-sided).
        Defaults to 1.
        correction_method (str, optional): Method for multiple comparisons correction ('none', 'bonferroni', 
        'sidak'). Defaults to 'none'.

    Returns:
        pd.DataFrame: DataFrame containing TGI estimate, confidence intervals, and P-value.
    """

    # Filter data for the relevant time points and groups
    tv_data = dataset[
        (dataset['day'].isin([baseline, time_point]))
        &
        (dataset['group'].isin(treatments_numerator + [treatment_denominator]))
    ].copy()

    # Log-transform tumor volume (TV) and add lnTV column
    tv_data['log_tv'] = np.log(tv_data['tv'])

    # Baseline adjustment (if applicable)
    if baseline_adjusted:
        tv_data['log_tv_baseline'] = (
            tv_data
            .groupby('animal')  # we pick for each animal
            .log_tv             # the log_tv value
            .transform('first') # of the first day, i.e., day 0
        )

    # We filter the time point we want to analyze
    tv_data = tv_data[tv_data['day'] == time_point]

    # Build and fit the linear model
    if baseline_adjusted:
        model = smf.ols(
            formula=f"log_tv ~ log_tv_baseline + C(group, Treatment(reference='{
                treatment_denominator}'))",
            data=tv_data)
    else:
        model = smf.ols(
            formula=f"log_tv ~ C(group, Treatment(reference='{
                treatment_denominator}'))",
            data=tv_data)
    est = model.fit()

    degrees_of_freedom = est.df_resid

    # Extract model estimates and calculate TGI
    # print(est.summary2())
    results = est.summary2().tables[1]  # returns a true DataFrame

    # Create a list to store the final results for all treatments
    model_estimates_list = []
    for treatment in treatments_numerator:
        # Extract TGI, CI, and P-value for the current treatment
        ix = f"C(group, Treatment(reference='{
            treatment_denominator}'))[T.{treatment}]"
        tgi = (1 - np.exp(results.loc[ix, "Coef."]))

        # Calculate confidence intervals and p-value
        std_err = results.loc[ix, "Std.Err."]

        if p_tails == 1:
            lb95_tgi = np.nan
            ub95_tgi = (
                1 - np.exp(results.loc[ix, "Coef."] + t_dist.ppf(0.95, degrees_of_freedom) * std_err))
            pvalue = t_dist.cdf(results.loc[ix, "t"], degrees_of_freedom)
        else:
            lb95_tgi = (
                1 - np.exp(results.loc[ix, "Coef."] + t_dist.ppf(0.025, degrees_of_freedom) * std_err))
            ub95_tgi = (
                1 - np.exp(results.loc[ix, "Coef."] + t_dist.ppf(0.975, degrees_of_freedom) * std_err))
            pvalue = 2 * (1 - t_dist.cdf(abs(results.loc[ix, "t"]), degrees_of_freedom))

        # Adjust p-values if multiple comparison correction is specified
        if correction_method == 'bonferroni':
            pvalue *= len(treatments_numerator)
        elif correction_method == 'sidak':
            pvalue = 1 - (1 - pvalue) ** len(treatments_numerator)

        # Create a dataframe of the results for the current treatment
        model_estimates = pd.DataFrame({
            "contrast": [f"{treatment} vs. {treatment_denominator}"],
            "day": [time_point],
            "tgi": [tgi],
            "lb95_tgi": [lb95_tgi],
            "ub95_tgi": [ub95_tgi],
            "pvalue": [pvalue],
            "df": [degrees_of_freedom],
            "baseline_adj": [baseline_adjusted],
            "correction": [correction_method],
        })

        # Append results for this treatment to the list
        model_estimates_list.append(model_estimates)

    # Combine results into a single DataFrame
    return pd.concat(model_estimates_list, ignore_index=True)

tgi_single_time(
    TGI_dataset,
    time_point=18,
    treatments_numerator=['therapy_a', 'therapy_b'],
    treatment_denominator='control',  # default
    baseline_adjusted=True,
    baseline=0, # default
    p_tails=1,  # default
)

tgi_single_time(
    TGI_dataset,
    time_point=18,
    treatments_numerator=['therapy_a', 'therapy_b'])

tgi_single_time(
    TGI_dataset,
    time_point=12,
    treatments_numerator=['therapy_combo_ab'],
    treatment_denominator='therapy_b',
    p_tails=2)

tgi_single_time(
    TGI_dataset,
    time_point=12,
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    correction_method='bonferroni')

tgi_single_time(
    TGI_dataset,
    time_point=12,
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    correction_method=None)

# Example dataset for this demonstration
tv_data = TGI_dataset[
    (TGI_dataset['group'].isin(['therapy_b', 'control']))
].copy()

# Log transformation
tv_data['log_tv'] = np.log(tv_data['tv'])

# Putting the baseline tumor volume in column
tv_data['log_tv_baseline'] = (
    tv_data
    .groupby('animal')  # we pick for each animal
    .log_tv             # the log_tv value
    .transform('first') # of the first day, i.e., day 0
)

# Keeping only the data point we want to analyze
tv_data = tv_data[tv_data['day'] == 18]

tv_data.head()

# First model that reproduces the default behavior of `tgi_single_time` (adjusted)
smf.ols(
    formula="log_tv ~ log_tv_baseline + C(group, Treatment(reference='control'))",
    data=tv_data).fit().summary2().tables[1]

# Second model with intercept removed
smf.ols(
    formula="log_tv ~ log_tv_baseline + group - 1",
    data=tv_data).fit().summary2().tables[1]

# Precalculated individual TGI-delta
# Example dataset for this demonstration (fresh copy)
tv_data = TGI_dataset[
    (TGI_dataset['group'].isin(['therapy_b', 'control']))
].copy()

# Putting the baseline tumor volume in column
tv_data['tv_baseline'] = (
    tv_data
    .groupby('animal')  # we pick for each animal
    .tv                 # the raw tv value
    .transform('first') # of the first day, i.e., day 0
)

tv_data['tv_adjusted'] = tv_data['tv'] - tv_data['tv_baseline']

# Keeping only the data point we want to analyze
tv_data = tv_data[tv_data['day'] == 18]

# Pseudocount to avoid division by zero
espilon = 1e-8
# Log transformation
tv_data['log_tv_adjusted'] = np.log(tv_data['tv_adjusted']+espilon)


smf.ols(formula="log_tv_adjusted ~ group",
        data=tv_data).fit().summary2().tables[1]

# Keeping only the data point we want to analyze
tv_data_all = TGI_dataset[TGI_dataset['day'] == 18].copy()

# Log transformation
tv_data_all['log_tv'] = np.log(tv_data_all['tv'])

# Create a weight matrix
# Add a small value to avoid zero variances
epsilon = 1e-8  # A small constant to prevent division by zero

# We might want to use the population variance (ddof = 0) if we consider
# each group to represent its own population
weights = 1.0 / (tv_data_all.groupby('group')
                 ['log_tv'].transform('var', ddof=0) + epsilon)

# Fit weighted least squares model
smf.wls(
    formula="log_tv ~ group",
    data=tv_data_all,
    weights=weights).fit().summary2().tables[1]

from scipy.stats import dunnett

def tgi_multiple_time(
    dataset: pd.DataFrame,
    time_points: List[int],
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    baseline_adjusted: bool = False,
    baseline: int = 0,
    p_tails: int = 1,
    correction_method: str = 'none',
) -> pd.DataFrame:
    """
    Estimates Tumor Growth Inhibition (TGI) at multiple time points using linear models.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', and 'tv' (tumor volume).
        time_points (List[int]): List of days at which to calculate TGI. Must have at least 2 time points.
        treatments_numerator (List[str]): List of treatment group names.
        treatment_denominator (str, optional): Name of the comparison group. Defaults to 'control'.
        baseline_adjusted (bool, optional): If True, adjust for baseline tumor volume. Defaults to False.
        baseline (int, optional): Baseline day for adjustment. Defaults to 0.
        p_tails (int, optional): Number of tails for hypothesis test (1 for one-sided, 2 for two-sided). Defaults to 1.
        correction_method (str, optional): Method for multiple comparisons correction ('none', 'bonferroni', 'sidak', 'dunnett'). 
        Defaults to 'none'.

    Returns:
        pd.DataFrame: DataFrame containing TGI estimates, confidence intervals, and P-values for each treatment at each time point.
    """

    if len(time_points) < 2:
        raise ValueError(
            "Please provide at least 2 time points for TGI calculation across multiple time points.")

    # Filter data for relevant groups and time points (including baseline if adjusting)
    relevant_days = time_points + \
        [baseline] if baseline_adjusted else time_points
    tv_data = dataset[
        (dataset['day'].isin(relevant_days))
        &
        (dataset['group'].isin(treatments_numerator + [treatment_denominator]))
    ].copy()

    # Log-transform tumor volume
    tv_data['log_tv'] = np.log(tv_data['tv'])

    # Baseline adjustment
    if baseline_adjusted:
        tv_data['log_tv_baseline'] = tv_data.groupby(
            'animal')['log_tv'].transform('first')

    # Convert Day to categorical, not simply integers, with correct order
    tv_data['day'] = pd.Categorical(
        tv_data['day'], categories=sorted(time_points), ordered=True)

    # Results container
    model_estimates_list = []

    # Iterate over time points
    for time_point in time_points:
        # Filter data for the current time point
        data_time_point = tv_data[tv_data['day'] == time_point]

        # Fit the linear model
        if baseline_adjusted:
            model = smf.ols(
                formula=f"log_tv ~ log_tv_baseline + C(group, Treatment(reference='{treatment_denominator}'))",
                data=data_time_point)
        else:
            model = smf.ols(
                formula=f"log_tv ~ C(group, Treatment(reference='{treatment_denominator}'))",
                data=data_time_point)
        est = model.fit()

        # Extract model estimates and calculate TGI for each treatment group
        results = est.summary2().tables[1]
        # print(results)

        for treatment in treatments_numerator:
            # Extract TGI, CI, and P-value for the current treatment
            ix = f"C(group, Treatment(reference='{
                treatment_denominator}'))[T.{treatment}]"
            tgi = (1 - np.exp(results.loc[ix, "Coef."]))

            # Calculate confidence intervals and P-value
            std_err = results.loc[ix, "Std.Err."]
            degrees_of_freedom = model.df_resid

            if p_tails == 1:
                lb95_tgi = np.nan
                ub95_tgi = (
                    1 - np.exp(results.loc[ix, "Coef."] + t_dist.ppf(0.95, degrees_of_freedom) * std_err))
                pvalue = t_dist.cdf(results.loc[ix, "t"], degrees_of_freedom)

            else:
                lb95_tgi = (
                    1 - np.exp(results.loc[ix, "Coef."] + t_dist.ppf(0.025, degrees_of_freedom) * std_err))
                ub95_tgi = (
                    1 - np.exp(results.loc[ix, "Coef."] + t_dist.ppf(0.975, degrees_of_freedom) * std_err))
                pvalue = 2 * \
                    (1 -
                     t_dist.cdf(abs(results.loc[ix, "t"]), degrees_of_freedom))

            # Adjust p-values if multiple comparison correction is specified
            if correction_method == 'bonferroni':
                pvalue *= len(treatments_numerator)
            elif correction_method == 'sidak':
                pvalue = 1 - (1 - pvalue) ** len(treatments_numerator)

            # Create a dataframe of the results for the current treatment
            model_estimates = pd.DataFrame({
                "contrast": [f"{treatment} vs. {treatment_denominator}"],
                "day": [time_point],
                "tgi": [tgi],
                "lb95_tgi": [lb95_tgi],
                "ub95_tgi": [ub95_tgi],
                "pvalue": [pvalue],
                "df": [degrees_of_freedom],
                "baseline_adj": [baseline_adjusted],
                "correction": [correction_method],
            })

            # Append results for this treatment to the list
            model_estimates_list.append(model_estimates)

    # Combine results into a single DataFrame
    all_results = pd.concat(model_estimates_list, ignore_index=True)

    # scipy.stats >= 1.11 includes Dunnett's test for multiple comparisons
    # Apply Dunnett's test
    # TO DO: DOUBLE CHECK THE RESULTS OF THE DUNNETT'S CORRECTION
    # TO DO: IMPLEMENT DUNNETT'S CORRECTION IN `tgi_single_time`?
    if correction_method == 'dunnett':
        control_data = tv_data[tv_data['group'] ==
                               treatment_denominator]['log_tv'].values
        treatment_data = [tv_data[tv_data['group'] == trt]
                          ['log_tv'].values for trt in treatments_numerator]

        dunnett_result = dunnett(*treatment_data, control=control_data)

        # Update P-values in the results DataFrame
        for i, trt in enumerate(treatments_numerator):
            all_results.loc[all_results["contrast"] == f"{trt} vs. {
                treatment_denominator}", "pvalue"] = dunnett_result.pvalue[i]

    # Combine results into a single DataFrame
    return all_results

tgi_multiple_time(
    dataset=TGI_dataset,
    time_points=[0, 3, 6, 9, 12, 15, 18],
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    treatment_denominator='control',)

tgi_multiple_time(
    dataset=TGI_dataset,
    time_points=[3, 6, 9, 12, 15, 18],
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    treatment_denominator='control',
    baseline_adjusted=True,
    baseline=0,
    correction_method='dunnett',)

tgi_multiple_time(
    dataset=TGI_dataset,
    time_points=[3, 6, 9, 12, 15, 18],
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    treatment_denominator='control',
    baseline_adjusted=True,
    baseline=0,
    correction_method='bonferroni',)

tgi_multiple_time(
    dataset=TGI_dataset,
    time_points=[3, 6, 9, 12, 15, 18],
    treatments_numerator=['therapy_a'],
    treatment_denominator='therapy_b',
    baseline_adjusted=True,
    baseline=0,
    p_tails=2,)

plot_tgi_stats(TGI_stats, ln_scale=True, colormap='viridis')

import warnings

from statsmodels.tools.sm_exceptions import ConvergenceWarning
warnings.simplefilter('ignore', ConvergenceWarning)
# These are warnings that we show as a way of alerting you that we may be in a non-standard situation. Most likely, one of
# the variance parameters is converging to zero. If we get a model result and the parameter estimates and standard errors are
# finite, we can just use them and ignore all the warnings. Note that inference for variance parameters is not meaningful when
# the point estimate is at or near zero.

def tumor_growth_curves_exponential(
    dataset: pd.DataFrame,
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    time_points: List[int] = None,
    p_tails: int = 1
) -> pd.DataFrame:
    """
    Estimates and compares tumor growth rates (TGR) between treatment groups using a linear mixed-effects model.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', and 'tv' (tumor volume).
        treatments_numerator (List[str]): List of treatment group names.
        treatment_denominator (str, optional): Name of the comparison group. Defaults to 'control'.
        time_points (List[int], optional): List of days at which to calculate TGR. If None, uses all unique days.
        p_tails (int, optional): 1 for one-sided test (tumor reduction vs. control), 2 for two-sided test. Defaults to 1.

    Returns:
        pd.DataFrame: DataFrame containing TGR estimates, confidence intervals, and P-values for each contrast.
    """

    # Filter the dataset for relevant groups and time points
    included_groups = treatments_numerator + [treatment_denominator]

    if not time_points:
        # uses all unique days in the dataset
        time_points = dataset['day'].unique()
    tv_data = dataset[
        (dataset['group'].isin(included_groups))
        &
        (dataset['day'].isin(time_points))
    ].copy()

    # Log-transform tumor volume
    tv_data['log_tv'] = np.log(tv_data['tv'])

    # Convert 'day' to numeric and 'group' to categorical variable (good practice)
    tv_data['day'] = pd.to_numeric(tv_data['day'])
    tv_data['group'] = pd.Categorical(
        tv_data['group'], categories=[treatment_denominator] + treatments_numerator)

    # Fit mixed effects model (similar to lme in R)
    model = smf.mixedlm(
        # equivalent to R formula 'lnTV ~ Group * Day, random = ~ 1 + Day | Animal'
        # log_tv modeled as a function of group, day, and their interaction
        f"log_tv ~ C(group, Treatment(reference='{
            treatment_denominator}'))*day",
        data=tv_data,
        # random effects
        # 'animal' is the grouping variable for the random effects
        groups=tv_data['animal'],
        re_formula='~day',  # including a random slope for day for each animal
        # the mixedlm function automatically includes a random intercept term
    )

    est = model.fit(
        reml=True,  # default
        method=["lbfgs"]
    )

    # Extract model estimates for the interaction terms (group:day)
    # print(est.summary())
    results = est.summary().tables[1]
    # Need to trick a little bit the output
    for col in results.columns:
        results[col] = pd.to_numeric(results[col], errors='coerce')
    # print(results)

    # Results container
    model_estimates_list = []

    for treatment in treatments_numerator:
        # Extract TGI, CI, and P-value for the interaction terms (group:day)
        ix = f"C(group, Treatment(reference='{
            treatment_denominator}'))[T.{treatment}]:day"
        tgr = results.loc[ix, "Coef."]

        # Calculate confidence intervals and P-value
        std_err = results.loc[ix, "Std.Err."]

        if p_tails == 1:
            lb95_tgr = np.nan
            ub95_tgr = results.loc[ix, "Coef."] + norm.ppf(0.95) * std_err
            # Although the t-distribution is used in R
            pvalue = norm.cdf(results.loc[ix, "z"])
        # we just pickup the results from statsmodels as it (normal distribution)
        else:
            lb95_tgr = results.loc[ix, "[0.025"]
            ub95_tgr = results.loc[ix, "0.975]"]
            pvalue = results.loc[ix, "P>|z|"]

        # Create a dataframe of the results for the current treatment
        model_estimates = pd.DataFrame({
            "contrast": [f"{treatment} vs. {treatment_denominator}"],
            "time_window": [f"{tv_data['day'].min()}-{tv_data['day'].max()}"],
            "tgr": [tgr],
            "lb95_tgr": [lb95_tgr],
            "ub95_tgr": [ub95_tgr],
            "pvalue": [pvalue],
        })

        # Append results for this treatment to the list
        model_estimates_list.append(model_estimates)

    # Combine results into a single DataFrame
    return pd.concat(model_estimates_list, ignore_index=True)

tumor_growth_curves_exponential(
    dataset=TGI_dataset,
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    treatment_denominator='control',      # default
    time_points=[0, 3, 6, 9, 12, 15, 18], # default all time-points analyzed
    p_tails=1,                            # default
)

tumor_growth_curves_exponential(
    dataset=TGI_dataset,
    treatments_numerator=['therapy_a'],
    treatment_denominator='therapy_b',
    time_points=[9, 12, 15, 18],
    p_tails=2)

tumor_growth_curves_exponential(
    dataset=TGI_dataset,
    treatments_numerator=['therapy_a'],
    treatment_denominator='therapy_b',
    time_points=[0, 3, 6, 9, 12, 15, 18],
    p_tails=2)

def simulate_quad_exp_growth(
    group_name: str,
    duration: int,
    initial_tumor_volume: float,
    group_growth_rate_day_pct: float,
    growth_rate_cv: float,
    group_decay_rate_day_pct: float,
    decay_rate_cv: float,
    residual_cv: float,
    n_animals: int,
    max_tumor_volume: float,
    df=None,
) -> pd.DataFrame:
    """
    Simulates tumor growth data with quadratic exponential growth for a group of animals.

    Args:
        group_name (str): Name of the group (e.g., "control", "therapy_a").
        duration (int): Duration of the study in days.
        initial_tumor_volume (float): Initial tumor volume (mm^3).
        group_growth_rate_day_perc (float): Average daily growth rate of the group (%).
        growth_rate_cv (float): Coefficient of variation for the growth rate within the group.
        group_decay_rate_day_perc (float): Average daily decay rate of the group (%).
        decay_rate_cv (float): Coefficient of variation for the decay rate within the group.
        residual_cv (float): Coefficient of variation for the measurement error.
        n_animals (int): Number of animals in the group.
        max_tumor_volume (float): Maximum allowed tumor volume before removal (mm^3).
        df (int, optional): Degrees of freedom for the t-distribution. If None,
        defaults to n_animals - 1.

    Returns:
        pd.DataFrame: Simulated data with columns: day, animal, group, tv (Tumor Volume).
    """

    # Initialize an empty list to store DataFrames for each animal in the group
    group_data = []

    # Loop over each animal in the group
    for i in range(n_animals):
        # Create an array of time points (days) when measurements are taken, starting from
        # day 0 up to the end of the study (`duration`). Measurements are taken every 3 days.
        days = np.arange(0, duration + 1, 3)

        # Sample individual growth and decay rates
        #
        # growth_rate_individual = np.random.normal(
        #     group_growth_rate_day_pct / 100,
        #     group_growth_rate_day_pct * growth_rate_cv / 100,
        #     len(days) - 1
        # )
        #
        # Since we're simulating tumor growth data, and the sample size is usually small
        # (below 15 animals per group), using the t distribution could be a reasonable choice.
        # It will introduce slightly more variability into the simulated data, reflecting the
        # uncertainty associated with small samples and the stochastic nature of tumor growth.

        # Shift and scale the t-distribution for individual growth rate
        df = n_animals - 1 if df is None else df  # Default df if not provided
        growth_rate_individual = t_dist.rvs(
            df=df,
            loc=group_growth_rate_day_pct / 100,
            scale=group_growth_rate_day_pct * growth_rate_cv / 100,
            size=len(days) - 1
        )

        # Same for the average daily decay rate
        # decay_rate_individual = np.random.normal(
        #     -group_decay_rate_day_pct / 100,  # Note the negative sign for decay
        #     group_decay_rate_day_pct * decay_rate_cv / 100,
        #     len(days) - 1
        # )
        decay_rate_individual = t_dist.rvs(
            df=df,
            loc=-group_decay_rate_day_pct / 100,  # Note the negative sign for decay
            scale=group_decay_rate_day_pct * decay_rate_cv / 100,
            size=len(days) - 1
        )

        # Note: if the growth rates are strictly positive and tend to be right-skewed (some
        # tumors grow much faster than others), the gamma distribution can be a good option.

        # Initialize an empty array to store the tumor volumes over time
        ln_tumor_volume = np.empty(len(days))

        # ln.tumor.volume[1] <- rnorm(n = 1, mean = log(initial.tumor.volume), sd = sqrt(log((residual.cv^2) + 1)))
        ln_tumor_volume[0] = np.random.normal(
            loc=np.log(initial_tumor_volume),
            scale=np.sqrt(np.log(residual_cv**2 + 1))
        )

        # Iterate over time points (days) to simulate tumor growth
        for t in range(1, len(days)):
            ln_tumor_volume[t] = np.random.normal(
                ln_tumor_volume[t - 1] + growth_rate_individual[t - 1] *
                days[t] + decay_rate_individual[t - 1] * days[t]**2,
                np.sqrt(np.log(residual_cv**2 + 1))
            )
            # Considere Gompertz Growth Model that describes sigmoidal growth patterns,
            # which are often observed in tumors.

        # # Transform log-tumor volumes back to the original scale (tumor volumes in mm^3).
        tumor_volume = np.exp(ln_tumor_volume)

        # Find indices where tumor volume exceeds the maximum allowed.
        welfare_missing = np.where(tumor_volume > max_tumor_volume)[0]
        if len(welfare_missing) > 1:
            # Exclude from its first occurrence until the next timepoint for survival .
            tumor_volume[welfare_missing[1:]] = np.nan

        # Create a DataFrame for the current animal's data
        animal_data = pd.DataFrame(
            {
                "day":    days,
                "animal": f"{group_name}_{i + 1}",
                "group":  group_name,
                "tv":     tumor_volume,
            }
        )

        # Add the animal's data to the list of group data
        group_data.append(animal_data)

    # Concatenate all animal dataframes into one for the entire group
    # ignore_index=True resets the index of the resulting DataFrame.
    return pd.concat(group_data, ignore_index=True)

# Set a seed for reproducibility
np.random.seed(1234)

# Simulate the vehicle (control) group data
vehicle_group_quad = simulate_quad_exp_growth(
    group_name="control",
    duration=24,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=5,
    growth_rate_cv=0.1,
    group_decay_rate_day_pct=.2,
    decay_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20,
    max_tumor_volume=4000)

# Simulate data for the single-agent therapy A group
therapy_a_group_quad = simulate_quad_exp_growth(
    group_name="therapy_a",
    duration=24,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=4.5,
    growth_rate_cv=0.1,
    group_decay_rate_day_pct=0.125,
    decay_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20,
    max_tumor_volume=4000)

# Simulate data for the single-agent therapy B group
therapy_b_group_quad = simulate_quad_exp_growth(
    group_name="therapy_b",
    duration=24,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=3.5,
    growth_rate_cv=0.1,
    group_decay_rate_day_pct=0.1,
    decay_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20,
    max_tumor_volume=4000)

# Simulate data for the combination therapy A+B group
therapy_ab_group_quad = simulate_quad_exp_growth(
    group_name="therapy_combo_ab",
    duration=24,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=1,
    growth_rate_cv=0.1,
    group_decay_rate_day_pct=0,
    decay_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20,
    max_tumor_volume=4000)

# Combine all group data into a single DataFrame
TGI_dataset_quad = pd.concat(
    [
        vehicle_group_quad,
        therapy_a_group_quad,
        therapy_b_group_quad,
        therapy_ab_group_quad
    ],
    ignore_index=True
)

# Let's have a look at the summary statistics dataset structure and values
print(summary_stats_tgi(TGI_dataset_quad).head(18))

    day      group          est        se    lb95_est     ub95_est
0     0    control    91.691684  0.055588   82.226772   102.246077
1     3    control   111.463101  0.062560   98.600641   126.003470
2     6    control   131.565475  0.096492  108.894927   158.955744
3     9    control   169.299692  0.121101  133.528760   214.653276
4    12    control   249.098542  0.130870  192.741077   321.934922
5    15    control   320.900927  0.146368  240.869553   427.523544
6    18    control   414.639052  0.170672  296.752212   579.357244
7    21    control   478.331367  0.208414  317.926447   719.666133
8    24    control   482.961511  0.221178  313.072613   745.040643
9     0  therapy_a   104.284857  0.050122   94.527387   115.049528
10    3  therapy_a   118.469700  0.080468  101.184204   138.708111
11    6  therapy_a   162.688813  0.085912  137.476671   192.524665
12    9  therapy_a   217.878030  0.108315  176.203841   269.408632
13   12  therapy_a   310.574403  0.114910  247.944348   389.024637
14   15  therapy_a   467.623324  0.122894  367.526349   594.982030
15   18  therapy_a   690.926564  0.131255  534.204138   893.627516
16   21  therapy_a  1032.256951  0.158114  757.182426  1407.262473
17   24  therapy_a  1372.125118  0.176029  971.755894  1937.448850

# For comparison with the R tutorial we can import the original TGI.dataset
TGI_dataset_quad = pd.read_csv("TGI_dataset_quadratic.csv")
TGI_dataset_quad.columns = [col.lower() for col in TGI_dataset_quad.columns]
for col in TGI_dataset_quad.select_dtypes(include='object').columns:
    TGI_dataset_quad[col] = TGI_dataset_quad[col].astype(str).str.lower()
TGI_dataset.head(6)

plot_tgi_stats(summary_stats_tgi(TGI_dataset_quad))

#############################################################################
# USE THIS FUNCTION WITH CAUTION ! IT DOESN'T GIVE THE SAME RESULTS AS THE  #
# ORIGINAL R SNIPPET, THOUGH THE R SNIPPET FAILS TO SUMMARIZE PROPERLY THE  #
# TGR AND TDR VALUES, AS WELL AS THE P VALUES AND THE EXPECTED CONCLUSIONS. #
#############################################################################

def tumor_growth_curves_quadratic_exponential(
    dataset: pd.DataFrame,
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    time_points: List[int] = None,
    p_tails: int = 1
) -> pd.DataFrame:
    """
    Estimates and compares tumor growth rates (TGR) and tumor decay rates (TDR) between treatment groups using a linear mixed-effects model
    with a quadratic term for time.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', and 'tv' (tumor volume).
        treatments_numerator (List[str]): List of treatment group names.
        treatment_denominator (str, optional): Name of the comparison group. Defaults to 'control'.
        time_points (List[int], optional): List of days at which to calculate TGR. If None, uses all unique days.
        p_tails (int, optional): 1 for one-sided test (tumor reduction vs. control), 2 for two-sided test. Defaults to 1.

    Returns:
        pd.DataFrame: DataFrame containing TGR and TDR estimates, with their confidence intervals, and P-values for each contrast.
    """

    # Filter the dataset for relevant groups and time points
    included_groups = treatments_numerator + [treatment_denominator]

    if not time_points:
        # uses all unique days in the dataset
        time_points = dataset['day'].unique()
    tv_data = dataset[
        (dataset['group'].isin(included_groups))
        &
        (dataset['day'].isin(time_points))
    ].copy()

    # Log-transform tumor volume
    tv_data['log_tv'] = np.log(tv_data['tv'])

    # Convert 'day' to numeric and 'group' to categorical variable (good practice)
    tv_data['day'] = pd.to_numeric(tv_data['day'])
    tv_data['group'] = pd.Categorical(
        tv_data['group'], categories=[treatment_denominator] + treatments_numerator)

    # Generate orthogonal polynomial features for Day (degree 2)
    poly_features = np.polynomial.polynomial.polyvander(tv_data['day'], deg=2)
    tv_data['day_linear'] = poly_features[:, 1]  # Linear term
    tv_data['day_quadratic'] = poly_features[:, 2]  # Quadratic term

    # Fit a linear mixed-effects model with a quadratic time component and random effects for both intercept and slope
    model = smf.mixedlm(
        # equivalent to R formula 'lnTV ~ Group * poly(Day, 2), random = ~ 1 + poly(Day, 2) | Animal'
        f"log_tv ~ C(group, Treatment(reference='{
            treatment_denominator}')) * (day_linear + day_quadratic)",
        data=tv_data,
        groups=tv_data["animal"],
        re_formula="day_linear + day_quadratic",
        missing='drop')

    est = model.fit(
        reml=True,  # default
        method=["lbfgs"]
    )

    # Extract model estimates for the interaction terms (group:day)
    # print(est.summary())
    results = est.summary().tables[1]
    # Need to trick a little bit the output
    for col in results.columns:
        results[col] = pd.to_numeric(results[col], errors='coerce')

    # Results container
    model_estimates_list = []

    for treatment in treatments_numerator:
        # Extract TGI, TGR, and their CI, and P-value related to the interaction terms with
        # the quadratic day component ('poly(day, 2)')
        ix_tgr = f"C(group, Treatment(reference='{treatment_denominator}'))[T.{
            treatment}]:day_linear"
        ix_tdr = f"C(group, Treatment(reference='{treatment_denominator}'))[T.{
            treatment}]:day_quadratic"

        tgr = results.loc[ix_tgr, "Coef."]
        tdr = results.loc[ix_tdr, "Coef."]

        # Calculate confidence intervals and P-values for TGR and TDR
        std_err_tgr = results.loc[ix_tgr, "Std.Err."]
        std_err_tdr = results.loc[ix_tdr, "Std.Err."]

        if p_tails == 1:
            lb95_tgr = np.nan
            ub95_tgr = results.loc[ix_tgr, "Coef."] + \
                norm.ppf(0.95) * std_err_tgr
            # Although the t-distribution is used in R
            pvalue_tgr = norm.cdf(results.loc[ix_tgr, "z"])
            lb95_tdr = np.nan
            ub95_tdr = results.loc[ix_tdr, "Coef."] + \
                norm.ppf(0.95) * std_err_tdr
            # Although the t-distribution is used in R
            pvalue_tdr = norm.cdf(results.loc[ix_tdr, "z"])

        # we just pickup the results from statsmodels as it (normal distribution)
        else:
            lb95_tgr = results.loc[ix_tgr, "[0.025"]
            ub95_tgr = results.loc[ix_tgr, "0.975]"]
            pvalue_tgr = results.loc[ix_tgr, "P>|z|"]
            lb95_tdr = results.loc[ix_tdr, "[0.025"]
            ub95_tdr = results.loc[ix_tgr, "0.975]"]
            pvalue_tdr = results.loc[ix_tdr, "P>|z|"]

        # Create a dataframe of the results for the current treatment
        model_estimates = pd.DataFrame({
            "contrast": [f"{treatment} vs. {treatment_denominator}"],
            "time_window": [f"{tv_data['day'].min()}-{tv_data['day'].max()}"],
            "tgr": [tgr],
            "lb95_tgr": [lb95_tgr],
            "ub95_tgr": [ub95_tgr],
            "pvalue_tgr": [pvalue_tgr],
            "tdr": [tdr],
            "lb95_tdr": [lb95_tdr],
            "ub95_tdr": [ub95_tdr],
            "pvalue_tdr": [pvalue_tdr],
        })

        # Append results for this treatment to the list
        model_estimates_list.append(model_estimates)

    # Combine results into a single DataFrame
    return pd.concat(model_estimates_list, ignore_index=True)

tumor_growth_curves_quadratic_exponential(
    TGI_dataset_quad,
    treatments_numerator=['therapy_a'],)

tumor_growth_curves_quadratic_exponential(
    TGI_dataset_quad,
    treatments_numerator=['therapy_a'],
    p_tails=2)

def simulate_gompertz_growth(
        group_name: str,
        duration: int,
        initial_tumor_volume: float,
        maximum_tumor_volume: float,
        group_growth_rate_day_pct: float,
        growth_rate_cv: float,
        residual_cv: float,
        n_animals: int,
        df=None
) -> pd.DataFrame:
    """
    Simulates tumor growth data for a group of animals over a given duration, using the Gompertzian model 
    with log-normal error.

    Args:
        group_name (str): Name of the group (e.g., "control", "therapy_a").
        duration (int): Duration of the study in days.
        initial_tumor_volume (float): Initial tumor volume (mm^3).
        maximum_tumor_volume (float): Maximum theoretical tumor volume (mm^3).
        group_growth_rate_day_pct (float): Average initial daily growth rate of the group (%).
        growth_rate_cv (float): Coefficient of variation for the growth rate within the group.
        residual_cv (float): Coefficient of variation for the measurement error.
        n_animals (int): Number of animals in the group.
        df (int, optional): Degrees of freedom for the t-distribution (growth rate variation). If None, 
        defaults to n_animals - 1.

    Returns:
        pd.DataFrame: Simulated data with columns: day, animal, group, tv (Tumor Volume).
    """

    group_data = []

    for i in range(n_animals):
        days = np.arange(0, duration + 1, 3)

        # Sample individual growth rates from a t-distribution
        df = n_animals - 1 if df is None else df

        # growth.rate.individual <- rnorm(n = (length(days) - 1), mean = group.growth.rate.day.perc / 100,
        # sd = group.growth.rate.day.perc * growth.rate.cv / 100)
        growth_rate_individual = t_dist.rvs(
            df=df,
            loc=group_growth_rate_day_pct / 100,
            scale=group_growth_rate_day_pct * growth_rate_cv / 100,
            size=len(days) - 1
        )

        # Calculate carrying capacity on the log scale (log of maximum tumor volume)
        carrying_capacity = np.log(maximum_tumor_volume)

        ln_tumor_volume = np.empty(len(days))
        # Day 0 simulation using normal distribution
        # ln.tumor.volume[1] <- rnorm(n = 1, mean = log(initial.tumor.volume), sd = sqrt(log((residual.cv^2) + 1)))
        ln_tumor_volume[0] = np.random.normal(
            loc=np.log(initial_tumor_volume),
            scale=np.sqrt(np.log(residual_cv**2 + 1))
        )

        # Simulate tumor volumes using the Gompertzian growth equation
        for t in range(1, len(days)):
            # growth.term <- (ln.tumor.volume[1] - carrying.capacity) * exp(-growth.rate.individual[t - 1] * days[t])
            # ln.tumor.volume[t] <- rnorm(n = 1, mean = carrying.capacity + growth.term, sd = sqrt(log((residual.cv^2) + 1)))
            # Calculate mean and standard deviation of the underlying normal distribution for log-normal
            growth_term = (ln_tumor_volume[0] - carrying_capacity) * \
                np.exp(-growth_rate_individual[t-1] * days[t])

            ln_tumor_volume[t] = np.random.normal(
                loc=carrying_capacity + growth_term,
                scale=np.sqrt(np.log(residual_cv**2 + 1)),
            )

        # Store the data for the current animal in a DataFrame
        animal_data = pd.DataFrame({
            "day": days,
            "animal": f"{group_name}_{i + 1}",
            "group": group_name,
            # Retrieve the tumor volume in linear scale
            "tv": np.exp(ln_tumor_volume),
        })
        group_data.append(animal_data)
    # Combine all animal dataframes into one for the entire group.
    return pd.concat(group_data, ignore_index=True)

# Set random seed for reproducibility
np.random.seed(1)

# Simulate the control group (Vehicle) data
control = simulate_gompertz_growth(
    group_name="control",
    duration=60,
    initial_tumor_volume=100,
    maximum_tumor_volume=6000,
    group_growth_rate_day_pct=6,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20)

# Simulate data for Therapy_A group
therapy_a = simulate_gompertz_growth(
    group_name="therapy_a",
    duration=60,
    initial_tumor_volume=100,
    maximum_tumor_volume=6000,
    group_growth_rate_day_pct=4,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20)

# Simulate data for Therapy_B group
therapy_b = simulate_gompertz_growth(
    group_name="therapy_b",
    duration=60,
    initial_tumor_volume=100,
    maximum_tumor_volume=6000,
    group_growth_rate_day_pct=3,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20)

# Simulate data for Therapy_Combo_AB group
therapy_combo_ab = simulate_gompertz_growth(
    group_name="therapy_combo_ab",
    duration=60,
    initial_tumor_volume=100,
    maximum_tumor_volume=6000,
    group_growth_rate_day_pct=2,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=20)

# Combine all group data into a single DataFrame
TGI_dataset_gompertz = pd.concat(
    [control, therapy_a, therapy_b, therapy_combo_ab], ignore_index=True)

TGI_dataset_gompertz.head(6)

plot_tgi_stats(summary_stats_tgi(TGI_dataset_gompertz))

# For comparison with the R tutorial we can import the original TGI.dataset
TGI_dataset_gompertz = pd.read_csv("TGI_dataset_gompertz.csv")
TGI_dataset_gompertz.columns = [
    col.lower() for col in TGI_dataset_gompertz.columns]
for col in TGI_dataset_gompertz.select_dtypes(include='object').columns:
    TGI_dataset_gompertz[col] = TGI_dataset_gompertz[col].astype(
        str).str.lower()
TGI_dataset_gompertz.head(6)

def tumor_growth_curves_gompertzian(
    dataset: pd.DataFrame,
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    p_tails: int = 1
) -> pd.DataFrame:
    """
    Estimates and compares Gompertzian parameters (lambda) between treatment groups using a linear mixed-effects model.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', 'animal', and 'tv' (tumor volume).
        treatments_numerator (List[str]): List of treatment group names (numerator in the comparison).
        treatment_denominator (str, optional): Name of the comparison group (denominator). Defaults to 'control'.
        p_tails (int, optional): 1 for one-sided test (tumor reduction), 2 for two-sided test. Defaults to 1.

    Returns:
        pd.DataFrame: DataFrame containing Gompertzian parameter estimates, confidence intervals, and p-values for each contrast.
    """

    # Filter the dataset for relevant groups and time points
    included_groups = treatments_numerator + [treatment_denominator]
    tv_data = dataset[dataset['group'].isin(included_groups)].copy()

    # Create Gompertz transformed variable
    # Pseudocount to avoid division by zero
    espilon = 1e-8
    tv_data["y_cross"] = np.log(
        np.log(tv_data['tv'].max() / (tv_data['tv'] + epsilon) + epsilon))
    # Handle undefined values (-inf) by replacing them with a very small negative number
    # tv_data["y_cross"] = tv_data["y_cross"].replace(-np.inf, tv_data["y_cross"][np.isfinite(tv_data["y_cross"])].min())
    tv_data["y_cross"] = tv_data["y_cross"].replace(-np.inf, espilon)

    # Ensure 'day' is numeric and convert 'group' to categorical variable (good practice)
    tv_data['day'] = pd.to_numeric(tv_data['day'])
    tv_data['group'] = pd.Categorical(
        tv_data['group'], categories=[treatment_denominator] + treatments_numerator)

    # Fit a linear mixed-effects model
    model = smf.mixedlm(
        # equivalent to R formula 'y.cross ~ Group * Day, random = ~ 1 + Day | Animal'
        f"y_cross ~ C(group, Treatment(reference='{
            treatment_denominator}')) * day",
        data=tv_data,
        groups=tv_data["animal"],
        re_formula="day",
        missing='drop')

    est = model.fit(
        reml=True,  # default
        # method=["lbfgs"]
    )

    # Extract model estimates for the interaction terms (group:day)
    # print(est.summary())
    results = est.summary().tables[1]
    # Need to trick a little bit the output
    for col in results.columns:
        results[col] = pd.to_numeric(results[col], errors='coerce')

    # Results container
    model_estimates_list = []

    for treatment in treatments_numerator:
        # Extract model estimates
        ix = f"C(group, Treatment(reference='{
            treatment_denominator}'))[T.{treatment}]:day"

        Gompertz_par = results.loc[ix, "Coef."]

        # Calculate confidence intervals and P-values for TGR and TDR
        std_err = results.loc[ix, "Std.Err."]

        if p_tails == 1:
            lb95 = np.nan
            ub95 = Gompertz_par + norm.ppf(0.95) * std_err
            # Although the t-distribution is used in R
            pvalue = 1 - norm.cdf(results.loc[ix, "z"])

        # we just pickup the results from statsmodels as it (normal distribution)
        else:
            lb95 = results.loc[ix, "[0.025"]
            ub95 = results.loc[ix, "0.975]"]
            pvalue = results.loc[ix, "P>|z|"]

        # Create a dataframe of the results for the current treatment
        model_estimates = pd.DataFrame({
            "contrast": [f"{treatment} vs. {treatment_denominator}"],
            "Gompertz_par": [Gompertz_par],
            "lb95": [lb95],
            "ub95": [ub95],
            "pvalue": [pvalue],
        })

        # Append results for this treatment to the list
        model_estimates_list.append(model_estimates)

    # Combine results into a single DataFrame
    return pd.concat(model_estimates_list, ignore_index=True)

tumor_growth_curves_gompertzian(
    TGI_dataset_gompertz,
    treatments_numerator=['therapy_a', 'therapy_b'])

tumor_growth_curves_gompertzian(
    TGI_dataset_gompertz,
    treatments_numerator=['therapy_a'],
    treatment_denominator='therapy_b',
    p_tails=2)

from scipy.integrate import trapezoid as trapz
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.trapezoid.html

def tumor_auc(
    dataset: pd.DataFrame,
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    p_tails: int = 1
) -> pd.DataFrame:
    """
    Calculates and compares area under the curve (AUC) of tumor growth curves between treatment groups.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', 'animal', and 'tv' (tumor volume).
        treatments_numerator (List[str]): List of treatment group names.
        treatment_denominator (str, optional): Name of the comparison group. Defaults to 'control'.
        p_tails (int, optional): 1 for one-sided test (tumor reduction), 2 for two-sided test. Defaults to 1.

    Returns:
        pd.DataFrame: DataFrame with columns 'contrast', 'AUC_diff', 'lb95_AUC_diff', 'ub95_AUC_diff', and 'pvalue'.
    """

    # Filter relevant data
    tv_data = dataset[dataset['group'].isin(
        treatments_numerator + [treatment_denominator])].copy()

    # Log-transform tumor volume (tv) and add log_tv column
    tv_data['log_tv'] = np.log(tv_data['tv'])

    # Extract baseline (first day) data
    tv_data['log_tv_baseline'] = tv_data.groupby(
        'animal')['log_tv'].transform('first')
    # remove baseline from further calculations for standardization, to avoid bias,
    # and focusing on treatment's effect on tumor growth dynamics
    tv_data = tv_data.query('day != day.min()').copy()

    # Calculate height for AUC (subtracting baseline and setting negative values to 0)
    tv_data['height'] = tv_data['log_tv'] - tv_data['log_tv_baseline']
    tv_data['height'] = tv_data['height'].apply(
        lambda x: max(0, x))  # Ensure non-negative heights

    # Calculate AUC for each animal using trapezoidal rule
    def calculate_auc(x):
        x = x.sort_values('day')
        return trapz(x['height'], x['day'])

    auc_data = (
        tv_data
        .groupby('animal')
        .apply(
            calculate_auc,
            include_groups=False
        )
        .reset_index(name='AUC')
    )

    # Merge AUC data with original data to get group information
    auc_data = auc_data.merge(
        tv_data[['animal', 'group']], on='animal', how='left').drop_duplicates()

    # Convert group to factor
    auc_data['group'] = pd.Categorical(
        auc_data['group'],
        categories=[treatment_denominator] + treatments_numerator)

    # Fit a linear model to compare AUCs between groups
    model = smf.ols('AUC ~ group', data=auc_data)
    est = model.fit()
    degrees_of_freedom = est.df_resid

    # Extract model estimates for the interaction terms (group:day)
    # print(est.summary2())
    results = est.summary2().tables[1]

    # Results container
    model_estimates_list = []

    for treatment in treatments_numerator:
        # Extract model estimates
        ix = f"group[T.{treatment}]"

        Gompertz_par = results.loc[ix, "Coef."]

        # Calculate confidence intervals and P-values for TGR and TDR
        std_err = results.loc[ix, "Std.Err."]

        if p_tails == 1:
            lb95 = np.nan
            ub95 = Gompertz_par +t_dist.ppf(0.95, df=degrees_of_freedom) * std_err
            pvalue = t_dist.cdf(results.loc[ix, "t"], df=degrees_of_freedom)

        # we just pickup the results from statsmodels as it (normal distribution)
        else:
            lb95 = results.loc[ix, "[0.025"]
            ub95 = results.loc[ix, "0.975]"]
            pvalue = results.loc[ix, "P>|t|"]

        # Create a dataframe of the results for the current treatment
        model_estimates = pd.DataFrame({
            "contrast": [f"{treatment} vs. {treatment_denominator}"],
            "AUC_diff": [Gompertz_par],
            "lb95_AUC_diff": [lb95],
            "ub95_AUC_diff": [ub95],
            "pvalue": [pvalue],
        })

        # Append results for this treatment to the list
        model_estimates_list.append(model_estimates)

    # Combine results into a single DataFrame
    return pd.concat(model_estimates_list, ignore_index=True)

tumor_auc(
    dataset=TGI_dataset,
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    treatment_denominator='control',
    p_tails=1)

tumor_auc(
    dataset=TGI_dataset,
    treatments_numerator=['therapy_a'],
    treatment_denominator='therapy_b',
    p_tails=2)

from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

def tte_analysis(
    dataset: pd.DataFrame,
    event_threshold: float,
    treatments_numerator: List[str],
    treatment_denominator: str = 'control',
    method: str = 'cox',
    plotKM: bool = False,
    colormap: str = "viridis"
) -> pd.DataFrame:
    """Performs time-to-event analysis using either CoxPH regression or log-rank test.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', 'animal', and 'tv' (tumor volume).
        event_threshold (float): Tumor volume threshold defining the event of interest.
        treatments_numerator (List[str]): List of treatment group names.
        treatment_denominator (str, optional): Name of the comparison group. Defaults to 'control'.
        method (str, optional): Method for TTE analysis ('cox' or 'logrank'). Defaults to 'cox'.
        plotKM (bool, optional): Whether to plot Kaplan-Meier curves. Defaults to False.
        colormap (str, optional): Colormap to use for group colors. Defaults to "viridis".
    Returns:
        pd.DataFrame: DataFrame containing results of the analysis.
    """

    # Filter for relevant groups
    included_groups = treatments_numerator + [treatment_denominator]
    tv_data = dataset[dataset['group'].isin(included_groups)].copy()

    # Ensure 'day' is numeric and convert 'group' to categorical variable (good practice)
    tv_data['day'] = pd.to_numeric(tv_data['day'])
    tv_data['group'] = pd.Categorical(
        tv_data['group'], categories=[treatment_denominator] + treatments_numerator)

    # The goal is to produce a table compatible with lifelines, such as we can pick up
    # the time of the first occurence of tumor volume exceeding the threshold ('event'),
    # or the latest study day if this never happens ('censoring'), e.g.,
    #           0  3  6  9 12 15 18
    # animal_a  -  -  -  -  -  x  x
    #                          ^
    # animal_b  -  -  -  -  -  -  -
    #                             ^

    # Create a Series of observed 'events' based on tumor volume threshold
    event = (
        tv_data
        .groupby('animal')
        .tv
        .apply(lambda x: (x > event_threshold).any())
    ).rename('event')

    # Create a Series of duration to event/censoring for each animal
    duration = (
        tv_data
        .set_index('day')  # for getting the index via idxmax
        .groupby('animal')  # for each animal
        .tv                # working with tumor volumes
        .apply(            # if event detected, take the first date
            lambda x: (x > event_threshold).idxmax() - tv_data['day'].min()
            # in the case any event is in the Series
            if (x > event_threshold).any()
            # otherwise take the last date
            else tv_data['day'].max() - tv_data['day'].min()
        )
    ).rename('duration')

    # Create lifelines dataFrame using both Series and group information
    tte_data = pd.merge(left=duration, right=event, on='animal', how='right').merge(
        tv_data.groupby("animal")["group"].first(), on='animal')

    # Results container
    model_estimates_list = []

    match method:
        # Cox Proportional Hazards Model
        case 'cox':
            # Fit the Cox model
            cph = CoxPHFitter()
            # Fit cph to data
            cph.fit(
                df=tte_data,
                duration_col="duration",
                event_col="event",
                formula='group',
            )

            # We pick up the columns of interest from the results table
            results = cph.summary[[
                'coef', 'exp(coef)', 'se(coef)', 'z', 'p', 'exp(coef) lower 95%', 'exp(coef) upper 95%']]
            results = results.rename(
                columns={
                    'covariate': 'contrast',
                    'exp(coef)': 'hr',
                    'exp(coef) lower 95%': 'lb95_hr',
                    'exp(coef) upper 95%': 'up95_hr',
                    'p': 'pvalue'
                })
            results.insert(0, "contrast", results.index.str.replace(
                r'^group\[T\.([^\]]+)\]$', rf'\1 vs. {treatment_denominator}', regex=True))
            # r'^group\[T\.|\]$', '', regex=True) + f" vs. {treatment_denominator}")

            model_estimates_list.append(results)
            # Differences in underlying numerical libraries and floating-point arithmetic between
            # R and Python can lead to minor variations in the calculations and results.

        # Log-Rank Test
        case "logrank":
            # For each pair of treatment tested
            for group in treatments_numerator:
                ix = tte_data['group'] == group

                # We pick up the columns of interest from the results table
                results = logrank_test(
                    durations_A=tte_data[ix]['duration'],
                    durations_B=tte_data[tte_data['group']
                                         == treatment_denominator]['duration'],
                    event_observed_A=tte_data[ix]['event'],
                    event_observed_B=tte_data[tte_data['group']
                                              == treatment_denominator]['event'],
                    alpha=.95).summary[['test_statistic', 'p']]

                results = results.rename(
                    columns={
                        'test_statistic': 'chisq',
                        'p': 'pvalue'
                    })
                results.insert(0, "contrast", f"{group} vs. {
                               treatment_denominator}")

                model_estimates_list.append(results)

        case _:
            raise ValueError("Invalid method. Choose 'cox' or 'logrank'.")

    # Plot Kaplan-Meier Curves (if requested)
    if plotKM:
        # Set up figure
        fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 5))
        # Map groups to colors from the colormap, and markers
        num_groups = len(groups := tte_data['group'].unique())
        colors = plt.get_cmap(colormap)(np.linspace(0.25, 1, num_groups))

        # Create KaplanMeierFitter object
        kmf = KaplanMeierFitter()

        for i, group in enumerate(groups):
            kmf = KaplanMeierFitter()
            ix = tte_data['group'] == group
            # Fit the data to the model
            kmf.fit(
                durations=tte_data[ix]['duration'],
                event_observed=tte_data[ix]['event'],
                label=group)

            kmf.plot(
                show_censors=True,
                ci_show=False,
                color=colors[i],
                lw=3,
                ax=ax)

        # Formatting
        ax.set_xlabel("Days", fontsize=12)
        ax.set_ylabel(r"$1 - \text{P(event)}$", fontsize=12)
        ax.set_yticks([0, .25, .5, .75, 1])
        ax.set_title("Kaplan-Meier analysis", fontsize=14)
        ax.grid(which='major', axis='y', linestyle='--')
        ax.legend(loc='lower left', fontsize=10)

    # Return (combined) results
    return pd.concat(model_estimates_list, ignore_index=True)

tte_analysis(
    dataset=TGI_dataset_gompertz,
    event_threshold=2000,
    treatments_numerator=['therapy_a', 'therapy_b'],
    treatment_denominator='control',
    method='cox',
    plotKM=True,
    colormap='viridis',)

tte_analysis(
    dataset=TGI_dataset_gompertz,
    event_threshold=2000,
    treatments_numerator=['therapy_a', 'therapy_b'],
    treatment_denominator='control',
    method='logrank',
    plotKM=False,)

from scipy.stats import fisher_exact
from statsmodels.stats.contingency_tables import Table

def response_analysis(
        dataset: pd.DataFrame,
        treatments_numerator: str,
        treatment_denominator: str,
        day_fin: int,
        day_ini: int = 0,
        pd_cutoff: float = 0.5,
        pr_cutoff: float = 0.3,
        sd_include: str = 'non-responder',
        p_tails: int = 2
    ) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Analyzes categorical response data from a TGI study using RECIST-like criteria.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', 'animal', and 'tv' (tumor volume).
        treatments_numerator (List[str]): List of treatment group names.
        treatment_denominator (str): Name of the control (vehicle) group.
        day_fin (int): Day (final) for assessing response.
        day_ini (int): Baseline/initial day for tumor volume measurement.
        pd_cutoff (float, optional): Threshold for Progressive Disease (%). Defaults to 0.5.
        pr_cutoff (float, optional): Threshold for Partial Response (%). Defaults to 0.3.
        sd_include (str, optional): How to handle stable disease ('non-responder', 'responder', or 'exclude'). Defaults to 'non-responder'.
        p_tails (int, optional): 1 for one-sided, 2 for two-sided test. Defaults to 2.

    Returns:
        Tuple[pd.DataFrame, pd.DataFrame]: A tuple containing three DataFrames:
            - The first DataFrame include a 4-level (PD, SD, PR, CR) + binary (responder, non-responder) response representation of the data.
            - The second DataFrame contains p-values from statistical tests (Cochran-Armitage trend and Fisher's exact).

        Additional tests, such as the chi-squared test for independence, can be performed using the first DataFrame returned.
    """

    # Filter for relevant groups
    included_groups = treatments_numerator + [treatment_denominator]
    tv_data = dataset[dataset['group'].isin(included_groups)].copy()

    # Ensure 'day' is numeric and convert 'group' to categorical variable
    tv_data['day'] = pd.to_numeric(tv_data['day'])
    tv_data['group'] = pd.Categorical(
        tv_data['group'], categories=[treatment_denominator] + treatments_numerator)
    
    # Filter data for the specified initial and final days
    filtered_data = tv_data[
        (tv_data['group'].isin(included_groups))
        & (tv_data['day'].isin([day_ini, day_fin]))
    ].copy()

    #  Calculate the relative change in tumor volume from day_ini to day_fin
    response_data = filtered_data.pivot(
        index=['animal', 'group'],
        columns='day',
        values='tv')

    response_data['relative_change'] = response_data.loc[:, day_fin] / response_data.loc[:, day_ini]
    
    # Categorize responses using RECIST-like criteria
    response_data['response'] = pd.cut(
        response_data['relative_change'],
        bins=[-np.inf, 1-pr_cutoff, 1+pd_cutoff, np.inf],
        labels=['pr', 'sd', 'pd'],
        right=False  # Include lower bound, exclude upper bound
    )

    # If the final tumor volume is <= 10, then it's categorized as CR
    response_data['response'] = response_data['response'].cat.add_categories('cr')
    ix = response_data[day_fin] <= 10
    response_data.loc[ix, 'response'] = 'cr'

    response_data['response'] = response_data['response'].cat.reorder_categories(
        ['pd', 'sd', 'pr', 'cr'])

    # Create a dataFrame with both response Series and group information
    response = response_data.reset_index(level='group').loc[:, ['group', 'response']]

    # Create separate columns for 'responder' and 'non_responder' based on SD_include
    match sd_include:
        case 'non-responder':
            response['responder'] = (
                response['response'].isin(['pr', 'cr']))
            response['non_responder'] = (
                response['response'].isin(['pd', 'sd']))
        case 'responder':
            response['responder'] = (
                response['response'].isin(['pr', 'sd', 'cr']))
            response['non_responder'] = (
                response['response'] == 'pd')
        case 'exclude':
            response['responder'] = (
                response['response'].isin(['pr', 'cr']))
            response['non_responder'] = (
                response['response'] == 'pd')
        case _:
            raise ValueError(
                "Invalid sd_include value. Choose from 'non-responder', 'responder', or 'exclude'.")

    # We drop rows with neither responder or non_responder
    response = response.drop(response[~response['responder'] & ~response['non_responder']].index)
    
    # Remove the column name generated from the pivot
    response.columns.name = ''

    # Print the 4-level and binary contingency tables
    print("Contingency table (quaternary/4-level response):")
    print(
        rt4:=pd.crosstab(
            response['group'],
            response['response'],
            dropna=False  # Include all categories even with zero counts
        )
    )

    print("\nContingency table (binary response):")
    print(
        rt2:=pd.crosstab(
            response['group'],
            response['responder'],
            dropna=False)
    )

    # Initialize lists to store p-values (initialize with NaN to handle cases where the test is not performed)
    resp_pv = np.full(len(treatments_numerator), np.nan)
    p_fisher = np.full(len(treatments_numerator), np.nan)

    # Iterate over treatment groups (start from 0 since enumerate starts from 0)
    for i, treatment in enumerate(treatments_numerator):
        # Check if there are enough non-zero cells in the contingency table
        control_row = rt4.loc[treatment_denominator]
        treatment_row = rt4.loc[treatment]

        # Check for at least two non-zero column sums
        if np.sum((control_row + treatment_row) != 0) >= 2:
            # Perform Cochran-Armitage Test
            result = Table(
                rt4.loc[[treatment_denominator, treatment]],
                shift_zeros=False
            ).test_ordinal_association(
                col_scores=np.array([1, 2, 3, 4])
            )
            p_trend = result.pvalue

            # Adjust p-value for one-sided tests
            if (p_tails == 1) & (p_trend != 1):
                p_trend = p_trend / 2 if result.statistic > 0 else 1 - p_trend / 2

            resp_pv[i] = p_trend

        # Check if both the control and treatment groups have at least one observation in any of the response categories
        if (control_row.sum() != 0) and (treatment_row.sum() != 0):
            # Perform Fisher's exact test on the 2x2 table
            _, p_fisher_val = fisher_exact(rt2.loc[[treatment_denominator, treatment]])

            # Adjust p-value for one-sided tests
            if (p_tails == 1) & (p_fisher_val != 1):
                odds_ratio, _ = fisher_exact(rt2.loc[[treatment_denominator, treatment]], alternative='greater')
                p_fisher_val = p_fisher_val / 2 if odds_ratio > 1 else 1 - p_fisher_val / 2

            p_fisher[i] = p_fisher_val

    # Create DataFrame for results
    res = pd.DataFrame({
        'contrast': [f"{treatments_numerator[i]} vs. {treatment_denominator}" for i in range(len(treatments_numerator))],
        "time_window": f"{day_ini}-{day_fin}",
        'p.trend': resp_pv,
        'p.fisher': p_fisher
    })

    return response, res

response_table, res = response_analysis(
    TGI_dataset,
    treatments_numerator=['therapy_a', 'therapy_b', 'therapy_combo_ab'],
    treatment_denominator='control',
    day_fin=6,
    sd_include='responder')

print(res)

Contingency table (quaternary/4-level response):
response          pd  sd  pr  cr
group                           
control           10   5   0   0
therapy_a          4  11   0   0
therapy_b          9   6   0   0
therapy_combo_ab   2  10   3   0

Contingency table (binary response):
responder         False  True 
group                         
control              10      5
therapy_a             4     11
therapy_b             9      6
therapy_combo_ab      2     13
                       contrast time_window   p.trend  p.fisher
0         therapy_a vs. control         0-6  0.030860  0.065595
1         therapy_b vs. control         0-6  0.709521  1.000000
2  therapy_combo_ab vs. control         0-6  0.002044  0.007775

import pingouin as pg

pg.chi2_independence(
    response_table[response_table['group'].isin(['control', 'therapy_a'])],
    x='group',
    y='responder',
    # correction=True,
)

c:\Users\Sébastien\Desktop\Bonato_tutorial_python\.env\Lib\site-packages\pingouin\contingency.py:151: UserWarning: Low count on observed frequencies.
  warnings.warn(f"Low count on {name} frequencies.")

(responder  False  True 
 group                  
 control      7.0    8.0
 therapy_a    7.0    8.0,
 responder  False  True 
 group                  
 control      9.5    5.5
 therapy_a    4.5   10.5,
                  test    lambda      chi2  dof      pval    cramer     power
 0             pearson  1.000000  3.348214  1.0  0.067278  0.334077  0.448299
 1        cressie-read  0.666667  3.362745  1.0  0.066687  0.334801  0.449867
 2      log-likelihood  0.000000  3.414737  1.0  0.064617  0.337379  0.455457
 3       freeman-tukey -0.500000  3.474878  1.0  0.062307  0.340337  0.461881
 4  mod-log-likelihood -1.000000  3.554471  1.0  0.059385  0.344213  0.470313
 5              neyman -2.000000  3.778385  1.0  0.051919  0.354889  0.493602)

rt = response_table[response_table['group'].isin(['control', 'therapy_combo_ab'])]

print(
    Table(
        pd.crosstab(
            rt['group'],
            rt['response'],
            dropna=False  # Include all categories even with zero counts
        ), shift_zeros=False
    ).test_ordinal_association(
        # col_scores=np.array([1, 2, 3, 4])
    )
)

null_mean   31.5
null_sd     5.350636442968966
pvalue      0.002044125135113983
statistic   48.0
zscore      3.0837453031745254

def plot_bliss(
    dataset: pd.DataFrame,
    single_1_name: str,
    single_2_name: str,
    combo_name: str,
    control_name: str = 'control',
    ln_scale: bool = False,
    backend_seaborn: bool = True,
    colormap: str = "viridis"
) -> None:
    """Plots tumor growth curves with 95% confidence intervals and Bliss independence line.

    Args:
        dataset (pd.DataFrame): DataFrame with columns 'day', 'group', 'animal', and 'tv' (tumor volume).
        single_1_name (str): Name of the first single-agent treatment group.
        single_2_name (str): Name of the second single-agent treatment group.
        combo_name (str): Name of the combination treatment group.
        control_name (str, optional): Name of the control group. Defaults to 'control'.
        ln_scale (bool, optional): If True, plot on natural log scale. Defaults to False.
        backend_seaborn (bool, optional): If True, plots with Seaborn, else uses Matplotlib.
        colormap (str, optional): Colormap to use for group colors. Defaults to "viridis".
    """

    # Filter data for relevant groups
    included_groups = [control_name, single_1_name, single_2_name, combo_name]
    dataset_list = dataset[dataset['group'].isin(included_groups)].copy()

    # Calculate summary statistics for each group and day
    summ_stats = summary_stats_tgi(dataset_list)

    # Calculate Bliss independence values
    def bliss_values(x: pd.DataFrame) -> pd.DataFrame:
        new_x = pd.DataFrame({
            'day':   [x['day'].unique()[0]],
            'group': ['Bliss'],
            'est':   [
                x[x['group'] == single_1_name]['est'].values[0] *
                x[x['group'] == single_2_name]['est'].values[0] /
                x[x['group'] == control_name]['est'].values[0]
            ]
        })
        return pd.concat([x, new_x], ignore_index=True)

    summ_stats = (
        summ_stats
        .groupby('day')
        .apply(bliss_values)
        .reset_index(drop=True)
    )

    # Prepare for plotting
    groups = [control_name, single_1_name, single_2_name, combo_name, 'Bliss']
    summ_stats['group'] = pd.Categorical(summ_stats['group'], categories=groups)

    # Plot
    if backend_seaborn:
        ax = sns.lineplot(
            data=TGI_dataset,
            x='day',
            y='tv',
            hue='group',
            style='group',
            dashes=False,
            markers=["o", "s", "^", "D"],
            estimator=gmean,
            errorbar=('ci', 95),
            n_boot=100,
            palette='viridis',
            err_style='band',
            lw=2,
            markersize=8)
        
        # plt.yscale('log')
        ax.set_xticks(TGI_dataset['day'].unique())
        ax.set_xticklabels(TGI_dataset['day'].unique(), fontsize=10)
        ax.set_xlabel("Days", fontsize=12)
        ax.set_ylabel("Tumor volume (mm³)", fontsize=12)
        ax.set_title("Tumor growth curves", fontsize=14)
        ax.grid(which='major', axis='y', linestyle='--')
        if ln_scale: ax.set_yscale("log")
    else:
        plot_tgi_stats(
            summ_stats[summ_stats['group'] != 'Bliss'],  # Plot only the non-Bliss groups
            ln_scale=ln_scale,
            colormap=colormap)

    # Add Bliss line
    bliss_data = summ_stats[summ_stats['group'] == 'Bliss']
    plt.plot(
        bliss_data['day'],
        bliss_data['est'],
        color='red',
        linestyle='dashed',
        linewidth=2,
        label='Bliss')

    # Show plot
    plt.legend(loc='upper left', fontsize=10);

plot_bliss(
    dataset=TGI_dataset,
    control_name='control',
    single_1_name='therapy_a',
    single_2_name='therapy_b',
    combo_name='therapy_combo_ab',
    backend_seaborn=False)

C:\Users\Sébastien\AppData\Local\Temp\ipykernel_4332\955682489.py:47: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
  .apply(bliss_values)

from patsy.contrasts import Treatment
from statsmodels.stats.contrast import ContrastResults

def bliss_synergy_test(
    dataset: pd.DataFrame,
    single_1_name: str,
    single_2_name: str,
    combo_name: str,
    time_points: List[int],
    control_name: str = 'control',
    baseline_adjusted: bool = False,
    baseline: int = 0,
) -> pd.DataFrame:
    """
    Tests for Bliss independence synergy between two single-agent treatments and their combination.

    Args:
        dataset (pd.DataFrame): TGI dataset with columns 'day', 'group', 'animal', and 'tv' (tumor volume).
        single_1_name (str): Name of the first single-agent treatment group.
        single_2_name (str): Name of the second single-agent treatment group.
        combo_name (str): Name of the combination treatment group.
        time_points (List[int]): List of time points (days) to analyze.
        control_name (str, optional): Name of the control group. Defaults to 'control'.
        baseline_adjusted (bool, optional): Whether to adjust for baseline tumor volume. Defaults to False.
        baseline (int, optional): Baseline day for adjustment. Required if `baseline_adjusted` is True. Defaults to None.

    Returns:
        pd.DataFrame:  A DataFrame with the following columns:
            - contrast: The comparison being made (e.g., "treatment_combo_ab vs. Bliss").
            - day: The time point (day) of the analysis.
            - synergy: exp(estimate), where 'estimate' is from the contrast.
            - lb95_synergy: Lower bound of the 95% confidence interval for synergy.
            - ub95_synergy: Upper bound of the 95% confidence interval for synergy.
            - pvalue: One-sided P value from the t-distribution.
            - df: Degrees of freedom.
            - baseline_adj: Whether the estimate was adjusted for baseline day.
    """

    # Filter data for relevant groups and time points (including baseline if adjusting)
    relevant_days = time_points + \
        [baseline] if baseline_adjusted else time_points
    groups = [control_name, single_1_name, single_2_name, combo_name]

    tv_data = dataset[
        dataset['group'].isin(groups) & dataset['day'].isin(relevant_days)
    ].copy()

    # Log-transform tumor volumes
    tv_data['log_tv'] = np.log(tv_data['tv'])

    # Baseline adjustment (if applicable)
    if baseline_adjusted:
        tv_data['log_tv_baseline'] = (
            tv_data
            .groupby('animal')
            .log_tv
            .transform(lambda x: x[tv_data['day'] == baseline].values[0])
        )  
        
        # Calculate adjusted lnTV (difference from baseline)
        tv_data['log_tv_adjusted'] = tv_data['log_tv'] - tv_data['log_tv_baseline']
        # Remove the baseline day
        tv_data = tv_data[tv_data['day'] != baseline]

    # Convert 'group' and 'day' to categorical variables with appropriate ordering
    tv_data['group'] = pd.Categorical(tv_data['group'], categories=groups, ordered=True)
    day_idx = sorted(tv_data['day'].unique())
    tv_data['day'] = pd.Categorical(
        tv_data['day'],
        categories=day_idx,
        ordered=True)

    # Perform analysis for each time point
    model_estimates = []
    response_variable = 'log_tv' if not baseline_adjusted else 'log_tv_adjusted'

    for day in day_idx:
        # Filter data for the current time point
        tv_data_day = tv_data[tv_data['day'] == day]

        # Group means (using the appropriate response variable based on baseline_adjusted)
        group_means = tv_data_day.groupby('group')[response_variable].mean()
        y_AB = group_means[combo_name]
        y_C = group_means[control_name]
        y_A = group_means[single_1_name]
        y_B = group_means[single_2_name]

        # Group sizes
        n_i = tv_data_day['group'].value_counts()

        # Calculate numerator
        numerator = (y_AB + y_C - y_A - y_B) * np.sqrt(np.sum(n_i) - 4)

        # Calculate denominator
        sum_sq_diff = 0
        for group, mean in group_means.items():
            group_data = tv_data_day[tv_data_day['group'] == group][response_variable]
            sum_sq_diff += np.sum((group_data - mean)**2)
        denominator = np.sqrt(sum_sq_diff * np.sum(1 / n_i))

        # Calculate t-statistic
        t_stat = numerator / denominator

        # Degrees of freedom
        df = np.sum(n_i) - 4

        # Calculate p-value (one-sided)
        p_value = t_dist.cdf(t_stat, df)

        # Calculate synergy estimate
        estimate = y_AB + y_C - y_A - y_B

        # Calculate confidence intervals
        alpha = 0.1  # For 90% confidence intervals
        t_critical = t_dist.ppf(1 - alpha / 2, df)
        se = (y_AB + y_C - y_A - y_B) / t_stat  # Standard error
        lb95_synergy = np.exp((estimate) - t_critical * se)
        ub95_synergy = np.exp((estimate) + t_critical * se)

        # Store results
        model_estimates.append({
            'contrast': f'{combo_name} vs. Bliss',
            'day': day,
            'synergy': np.exp(estimate),
            'lb95_synergy': lb95_synergy,
            'ub95_synergy': ub95_synergy,
            'pvalue': p_value,
            'df': df,
            'baseline_adj': baseline_adjusted,
        })

    return pd.DataFrame(model_estimates)

bliss_synergy_test(
    dataset=TGI_dataset,
    single_1_name='therapy_a',
    single_2_name='therapy_b',
    combo_name='therapy_combo_ab',
    time_points=[0, 3, 6, 9, 12, 15, 18],
)

# Generate another TGI_dataset with a substantial synergy
np.random.seed(12345)

# Simulate the vehicle (control) group data
vehicle_group = simulate_exp_growth(
    group_name="control",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=5,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Simulate data for the single-agent therapy A group
therapy_a_group = simulate_exp_growth(
    group_name="therapy_a",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=4,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Simulate data for the single-agent therapy B group
therapy_b_group = simulate_exp_growth(
    group_name="therapy_b",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=4,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Simulate data for the combination therapy A+B group
therapy_ab_group = simulate_exp_growth(
    group_name="therapy_combo_ab",
    duration=18,
    initial_tumor_volume=100,
    group_growth_rate_day_pct=1.1,
    growth_rate_cv=0.1,
    residual_cv=0.25,
    n_animals=15,
    max_tumor_volume=2000)

# Combine all group data into a single DataFrame
TGI_dataset_12345 = pd.concat(
    [
        vehicle_group,
        therapy_a_group,
        therapy_b_group,
        therapy_ab_group
    ],
    ignore_index=True)

plot_bliss(
    dataset=TGI_dataset_12345,
    single_1_name='therapy_a',
    single_2_name='therapy_b',
    combo_name='therapy_combo_ab',
    ln_scale=False,
    backend_seaborn=True)

bliss_synergy_test(
    dataset=TGI_dataset_12345,
    single_1_name='therapy_a',
    single_2_name='therapy_b',
    combo_name='therapy_combo_ab',
    time_points=[0, 3, 6, 9, 12, 15, 18],
    baseline_adjusted=True)

from scipy.stats import t, nct
from scipy.optimize import fsolve

# Parameters
df = 10  # Degrees of freedom
alpha = 0.05  # Significance level

# Calculate critical value
critical_value = t.ppf(1 - alpha, df)

# Generate x-values for plotting
x = np.linspace(t.ppf(0.001, df), t.ppf(0.999, df), 100)

# Generate y-values for central t-distribution
y_central = t.pdf(x, df)

# Function to calculate power
def calculate_power(delta):
    y_noncentral = nct.pdf(x, df, delta)
    power = y_noncentral[x > critical_value].sum() / y_noncentral.sum()  # Approximate power
    return power - target_power  # Adjust for fsolve

# Target power levels (which correspond to 1-beta)
target_powers = [0.2, 0.5, 0.8]

# Create subplots
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for i, target_power in enumerate(target_powers):
    # Find delta for target power
    delta_optimized = fsolve(calculate_power, 2)  # Start with an initial guess of 2

    # Generate y-values for non-central t-distribution with optimized delta
    y_noncentral_opt = nct.pdf(x, df, delta_optimized)

    # Plot the distributions
    axes[i].plot(x, y_central, label='Central t')
    axes[i].plot(x, y_noncentral_opt, label='Non-central t')

    # Shade the rejection region
    axes[i].fill_between(x[x > critical_value], 0, y_central[x > critical_value], alpha=0.5, color='mediumseagreen', label='Rejection Region')

    # Shade the beta area
    axes[i].fill_between(x[x < critical_value], 0, y_noncentral_opt[x < critical_value], alpha=0.5, color='orange', label='Beta')

    # Add vertical line for critical value
    axes[i].axvline(x=critical_value, color='black', linestyle='--', label='Critical value')

    # Add labels and legend
    axes[i].set_title(fr"$\beta=$ {1-target_power:.0%} | $\delta=$ {delta_optimized[0]:.2f}")
    axes[i].set_xlabel('t-value')
    axes[i].set_ylabel('Probability density')
    axes[i].legend()

# Adjust layout
plt.tight_layout();

from scipy.optimize import brentq

def pwr_gt2_grps(n, sigma, eff, two_tailed, n_groups, alpha, method):
    """
    Calculates the power for more than two groups.

    Args:
      n: Sample size per group.
      sigma: Residual standard deviation.
      eff: Effect size.
      two_tailed: Two-sided or one-sided test.
      n_groups: Number of groups.
      alpha: Significance level.
      method: "ANOVA" or "ANCOVA".

    Returns:
      Power of the test.
    """
    df_n = (n - 1) * n_groups - (method == "ANCOVA")
    se_n = (2 / n)**0.5 * sigma
    if two_tailed:
        alpha /= 2
    tcrit = t.ppf(1 - alpha, df_n)
    noncen = np.log(1 + eff) / se_n
    return 1 - nct.cdf(tcrit, df=df_n, nc=noncen) + nct.cdf(-tcrit, df=df_n, nc=noncen)

def pwr_2_grps(n1, n2, sigma1, sigma2, eff, two_tailed, alpha, method):
    """
    Calculates the power for two groups.

    Args:
      n1: Sample size for group 1.
      n2: Sample size for group 2.
      sigma1: Standard deviation for group 1.
      sigma2: Standard deviation for group 2.
      eff: Effect size.
      two_tailed: Two-sided or one-sided test.
      alpha: Significance level.
      method: "ANOVA" or "ANCOVA".

    Returns:
      Power of the test.
    """
    df_n = n1 + n2 - 2 - (method == "ANCOVA")
    se_n = (1 / n1 * sigma1**2 + 1 / n2 * sigma2**2)**0.5
    alpha = alpha / 2 if two_tailed else alpha
    tcrit = t.ppf(1 - alpha, df_n)
    noncen = np.log(1 + eff) / se_n
    return 1 - nct.cdf(tcrit, df=df_n, nc=noncen) + nct.cdf(-tcrit, df=df_n, nc=noncen)

def power_function(
        sigma: float | List[float] = .25,
        n: float | List[float] = None,
        eff: float = 0.5,
        power: float = 0.8,
        n_groups=2,
        alpha=0.05,
        two_tailed=True,
        method="ANCOVA") -> pd.DataFrame:
    """
    Calculates power or other parameters for a given set of inputs.

    Args:
      sigma: Residual standard deviation (or list for two groups).
      n: Sample size per group (or list for two groups).
      eff: Effect size.
      power: Desired power.
      n_groups: Number of group.
      alpha: Significance level.
      two_tailed: Two-sided test or one-sided (Defaults to True).
      method: "ANOVA" or "ANCOVA".

    Returns:
      A pandas DataFrame with calculated parameters.
    """
    pars = [
        sigma is None,
        n is None,
        eff is None,
        power is None]
    
    if sum(pars) != 1:
        return "Error: one and only one parameter should be NULL"
    
    par = pars.index(True)

    if par == 3 and n_groups > 2:  # Calculate 'power' for more than 2 groups
        df_n = (n - 1) * n_groups - (method == "ANCOVA")
        se_n = (2 / n)**0.5 * sigma
        pow = pwr_gt2_grps(
            n=n,
            sigma=sigma,
            eff=eff,
            two_tailed=two_tailed,
            n_groups=n_groups,
            alpha=alpha,
            method=method)
        return pd.DataFrame({
            "sigma": [sigma],
            "n_groups": [n_groups],
            "n": [n],
            "df": [df_n],
            "se": [se_n],
            "eff": [eff],
            "power": [pow],})

    elif par == 3 and n_groups == 2:  # Calculate 'power' for two groups
        n1, n2 = n[0], n[1] if isinstance(n, list) else n, n
        s1, s2 = sigma[0], sigma[1] if isinstance(sigma, list) else sigma, sigma
        df_n = (n1 + n2 - 2) - (method == "ANCOVA")
        se_n = (s1**2 / n1 + s2**2 / n2)**0.5
        pow = pwr_2_grps(
            n1=n1,
            n2=n2,
            sigma1=s1,
            sigma2=s2,
            eff=eff,
            two_tailed=two_tailed,
            alpha=alpha,
            method=method)
        return pd.DataFrame({
            "sigma1": [s1],
            "sigma2": [s2],
            "n1": [n1],
            "n2": [n2],
            "df": [df_n],
            "se": [se_n],
            "eff": [eff],
            "power": [pow],})

    elif par != 3 and n_groups > 2:  # Calculate other parameters
        if sigma is None:
            f = lambda sc: power - pwr_gt2_grps(
                n=n,
                sigma=sc,
                eff=eff,
                two_tailed=two_tailed,
                n_groups=n_groups,
                alpha=alpha,
                method=method)
            sc = brentq(f, 0.01, 20)  # looking up at values of sigma between 0.01 and 20
            df_n = (n - 1) * n_groups - (method == "ANCOVA")
            se_n = (2 / n)**0.5 * sc
            return pd.DataFrame({
                "sigma": [sc],
                "n_groups": [n_groups],
                "n": [n],
                "df": [df_n],
                "se": [se_n],
                "eff": [eff],
                "power": [power],})

        elif n is None:
            f = lambda ns: power - pwr_gt2_grps(
                n=ns,
                sigma=sigma,
                eff=eff,
                two_tailed=two_tailed,
                n_groups=n_groups,
                alpha=alpha,
                method=method)
            ns = brentq(f, 2, 400)  # looking up at values of n between 2 and 400
            df_n = (ns - 1) * n_groups - (method == "ANCOVA")
            se_n = (2 / ns)**0.5 * sigma
            return pd.DataFrame({
                "sigma": [sigma],
                "n_groups": [n_groups],
                "n": [ns],
                "df": [df_n],
                "se": [se_n],
                "eff": [eff],
                "power": [power],})

        else:  # eff is None:
            f = lambda ef: power - pwr_gt2_grps(
                n=n,
                sigma=sigma,
                eff=ef,
                two_tailed=two_tailed,
                n_groups=n_groups,
                alpha=alpha,
                method=method)
            ef = brentq(f, .1, 10)  # looking up at values of eff between 0.1 and 10
            df_n = (n - 1) * n_groups - (method == "ANCOVA")
            se_n = (2 / n)**0.5 * sigma
            return pd.DataFrame({
                "sigma": [sigma],
                "n_groups": [n_groups],
                "n": [n],
                "df": [df_n],
                "se": [se_n],
                "eff": [ef],
                "power": [power],})

    elif par != 3 and n_groups == 2:  # Calculate other parameters 2 groups
        if sigma is None:
            n1, n2 = (n[0], n[1]) if isinstance(n, list) else (n, n)
            f = lambda sc: power - pwr_2_grps(
                n1=n1,
                n2=n2,
                sigma1=sc,
                sigma2=sc,
                eff=eff,
                two_tailed=two_tailed,
                alpha=alpha,
                method=method)
            sc = brentq(f, 0.01, 20)  # looking up at values of sigma between 0.01 and 20
            df_n = (n1 + n2 - 2) - (method == "ANCOVA")
            se_n = (sc**2 / n1 + sc**2 / n2)**0.5
            return pd.DataFrame({
                "sigma1": [sc],
                "sigma2": [sc],
                "n1": [n1],
                "n2": [n2],
                "df": [df_n],
                "se": [se_n],
                "eff": [eff],
                "power": [power],})

        elif n is None:
            s1, s2 = (sigma[0], sigma[1]) if isinstance(sigma, list) else (sigma, sigma)
            f = lambda ns: power - pwr_2_grps(
                n1=ns,
                n2=ns,
                sigma1=s1,
                sigma2=s2,
                eff=eff,
                two_tailed=two_tailed,
                alpha=alpha,
                method=method)
            ns = brentq(f, 2, 400)  # looking up at values of n between 2 and 400
            df_n = (ns + ns - 2) - (method == "ANCOVA")
            se_n = (s1**2 / ns + s2**2 / ns)**0.5
            return pd.DataFrame({
                "sigma1": [s1],
                "sigma2": [s2],
                "n1": [ns],
                "n2": [ns],
                "df": [df_n],
                "se": [se_n],
                "eff": [eff],
                "power": [power],})

        else:  # eff is None:
            n1, n2 = (n[0], n[1]) if isinstance(n, list) else (n, n)
            s1, s2 = (sigma[0], sigma[1]) if isinstance(sigma, list) else (sigma, sigma)
            f = lambda ef: power - pwr_2_grps(
                n1=n1,
                n2=n2,
                sigma1=s1,
                sigma2=s2,
                eff=ef,
                two_tailed=two_tailed,
                alpha=alpha,
                method=method)
            ef = brentq(f, .1, 10)  # looking up at values of eff between 0.1 and 10
            df_n = (n1 + n2 - 2) - (method == "ANCOVA")
            se_n = (s1**2 / n1  + s2**2 / n2)**0.5
            return pd.DataFrame({
                "sigma1": [s1],
                "sigma2": [s2],
                "n1": [n1],
                "n2": [n2],
                "df": [df_n],
                "se": [se_n],
                "eff": [ef],
                "power": [power],})

power_function(
    sigma=.25,
    n=15,
    eff=None,
    power=.8,
    n_groups=4,
    alpha=.05,
    method='ANOVA',
    two_tailed=True
)

power_function(
    sigma=[.25, .35],
    n=[10, 15],
    eff=None,
    power=.8,
    n_groups=2,
    alpha=.05,
)

!python --version
print("-------------")
# List of packages we want to check the version
packages = ['numpy', 'pandas', 'matplotlib', 'seaborn', 'scipy', 'statsmodels', 'pingouin']

# Initialize an empty list to store the versions
versions = []

# Loop over the packages
for package in packages:
    # Get the version of the package
    output = !pip show {package} | findstr "Version"
    # If the output is not empty, get the version
    if output:
        version = output[0].split()[-1]
    else:
        version = 'Not installed'
    # Append the version to the list
    versions.append(version)

# Print the versions
for package, version in zip(packages, versions):
    print(f'{package}: {version}')

Python 3.12.7
-------------
numpy: 1.26.4
pandas: 2.2.3
matplotlib: 3.9.2
seaborn: 0.13.2
scipy: 1.14.1
statsmodels: 0.14.4
pingouin: 0.5.5

		tv
		count	mean	std	min	25%	50%	75%	max
group	day
control	0	15.0	102.337373	12.879400	73.966753	98.079914	102.927819	108.765041	125.245267
	3	15.0	135.822689	32.528511	83.542341	113.084926	127.282145	159.104481	200.072588
	6	15.0	189.766983	47.882196	120.659156	159.487281	183.848327	218.303025	279.768096
	9	15.0	301.584275	113.275488	120.438583	236.362885	282.724282	359.417484	563.845041
	12	15.0	530.467881	247.383477	205.588181	333.563703	456.265546	707.509967	970.984050
	15	15.0	1347.140701	670.145035	296.057302	870.395873	1286.366006	1795.646514	2831.206338
	18	15.0	3420.889383	2090.362199	463.809975	2017.002897	3228.041028	4630.789148	8345.411255
therapy_a	0	15.0	106.963038	28.635230	58.913855	91.820106	97.169935	113.085365	160.933938
	3	15.0	113.424589	45.696562	47.949666	84.560273	99.675973	131.490110	218.594997
	6	15.0	147.484120	57.812170	45.955902	118.466128	131.192080	187.781208	270.608753
	9	15.0	217.481620	99.489732	62.945368	155.360546	175.624492	272.227813	393.041179
	12	15.0	356.223597	178.620840	84.066759	240.887216	350.549399	423.646575	798.510798
	15	15.0	616.542902	287.643446	173.321385	412.393238	609.756321	768.338336	1300.752659
	18	15.0	1272.715957	616.831030	491.715154	830.020666	939.946099	1727.996591	2449.440952
therapy_b	0	15.0	101.622632	27.748183	68.683737	86.996987	94.462453	100.606739	173.368783
	3	15.0	116.502890	45.676643	51.471512	77.642927	111.360294	151.265746	190.062664
	6	15.0	159.943658	66.658128	71.677215	112.450580	152.966335	201.741230	285.400361
	9	15.0	223.773995	105.409381	71.919811	148.479905	211.324123	275.358286	430.777369
	12	15.0	320.154044	177.458595	85.357512	224.919608	283.700964	421.497751	704.788052
	15	15.0	559.365932	385.954352	115.675295	370.667577	443.944166	689.627500	1523.108747
	18	15.0	847.034379	613.333316	183.144309	508.833347	762.757727	939.605722	2468.363846
therapy_combo_ab	0	15.0	98.863617	21.913776	58.542840	80.452217	103.228822	112.582789	131.965342
	3	15.0	102.033695	34.218192	56.524456	81.264018	97.960340	112.704304	188.781830
	6	15.0	111.866578	58.485060	36.540022	75.760590	101.452698	128.681882	260.732431
	9	15.0	130.357569	70.031341	37.381473	74.024599	127.813497	172.560568	274.715958
	12	15.0	167.975680	109.856143	34.714108	84.479192	149.857299	208.723595	413.430485
	15	15.0	213.252799	123.938002	40.303314	115.875287	212.318566	265.859798	489.926716
	18	15.0	287.799543	203.230603	63.786859	150.290977	234.453907	333.979129	783.215058

group	control	therapy_a	therapy_b	therapy_combo_ab
day
0	0.0	-0.045200	0.006984	0.033944
3	0.0	0.164907	0.142243	0.248773
6	0.0	0.222815	0.157158	0.410506
9	0.0	0.278869	0.258005	0.567757
12	0.0	0.328473	0.396469	0.683344
15	0.0	0.542332	0.584775	0.841700
18	0.0	0.627958	0.752394	0.915870

	adjusted_geometric_mean
group	control	therapy_a	therapy_b	therapy_combo_ab
day
0	0.000000	0.000000	0.000000	0.000000
3	30.725044	2.022239	9.306823	0.925464
6	82.763423	32.416607	48.762187	2.790968
9	179.926451	92.035116	100.820564	15.480252
12	375.574628	209.007015	172.120106	40.514959
15	1068.437135	447.567110	344.807745	81.850178
18	2648.639697	1033.455049	578.431143	128.736380

	Coef.	Std.Err.	t	P>\|t\|	[0.025	0.975]
Intercept	4.586897	3.325829	1.379174	0.179159	-2.237141	11.410934
C(group, Treatment(reference='control'))[T.therapy_b]	-1.381330	0.267891	-5.156304	0.000020	-1.930998	-0.831662
log_tv_baseline	0.721260	0.718649	1.003633	0.324466	-0.753286	2.195806

	Coef.	Std.Err.	t	P>\|t\|	[0.025	0.975]
group[control]	4.586897	3.325829	1.379174	0.179159	-2.237141	11.410934
group[therapy_b]	3.205567	3.305763	0.969691	0.340808	-3.577298	9.988432
log_tv_baseline	0.721260	0.718649	1.003633	0.324466	-0.753286	2.195806

	day	animal	group	tv
0	0	control_1	control	104.825993
1	3	control_1	control	125.137266
2	6	control_1	control	155.847247
3	9	control_1	control	238.687772
4	12	control_1	control	335.119128
5	15	control_1	control	729.343194

group	control	therapy_a	therapy_b	therapy_combo_ab
day
0	101.536731	103.594835	98.736314	96.364433
3	132.261775	105.617074	108.043137	97.289897
6	184.300154	136.011443	147.498501	99.155402
9	281.463182	195.629951	199.556878	111.844686
12	477.111360	312.601850	270.856421	136.879392
15	1169.973866	551.161945	443.544059	178.214612
18	2750.176428	1137.049884	677.167458	225.100813

	contrast	day	tgi_delta_pct
1	therapy_b vs. control	6	33.293736
2	therapy_b vs. control	12	48.956823
3	therapy_b vs. control	18	77.538042

	contrast	day	tgi	lb95_tgi	ub95_tgi	pvalue	df	baseline_adj	correction
0	therapy_a vs. control	18	0.594287	NaN	0.399294	1.922899e-04	41.0	True	none
1	therapy_b vs. control	18	0.747207	NaN	0.625606	3.074967e-07	41.0	True	none

	contrast	day	tgi	lb95_tgi	ub95_tgi	pvalue	df	baseline_adj	correction
0	therapy_a vs. control	18	0.586554	NaN	0.378781	3.613441e-04	42.0	False	none
1	therapy_b vs. control	18	0.753773	NaN	0.630035	3.985567e-07	42.0	False	none

	contrast	day	tgi	lb95_tgi	ub95_tgi	pvalue	df	baseline_adj	correction
0	therapy_a vs. control	12	0.344803	NaN	0.056728	8.602930e-02	56.0	False	bonferroni
1	therapy_b vs. control	12	0.432299	NaN	0.182694	1.792068e-02	56.0	False	bonferroni
2	therapy_combo_ab vs. control	12	0.713108	NaN	0.586968	6.225085e-07	56.0	False	bonferroni

	day	animal	group	tv	log_tv	log_tv_baseline
6	18	control_1	control	2447.881533	7.802978	4.652302
13	18	control_2	control	1586.124262	7.369049	4.577525
20	18	control_3	control	2681.138968	7.893997	4.612758
27	18	control_4	control	463.809975	6.139475	4.634028
34	18	control_5	control	997.884643	6.905638	4.759049

	Coef.	Std.Err.	t	P>\|t\|	[0.025	0.975]
Intercept	7.863679	0.217587	36.140394	4.855426e-25	7.417972	8.309385
group[T.therapy_b]	-1.555499	0.307714	-5.055007	2.386084e-05	-2.185823	-0.925174

	Coef.	Std.Err.	t	P>\|t\|	[0.025	0.975]
Intercept	7.919420	0.195191	40.572714	3.267968e-43	7.528406	8.310434
group[T.therapy_a]	-0.883228	0.233706	-3.779235	3.835349e-04	-1.351396	-0.415060
group[T.therapy_b]	-1.401502	0.267171	-5.245708	2.468133e-06	-1.936710	-0.866294
group[T.therapy_combo_ab]	-2.502872	0.276547	-9.050442	1.494836e-12	-3.056862	-1.948882

	contrast	day	tgi	lb95_tgi	ub95_tgi	pvalue	df	baseline_adj	correction
0	therapy_a vs. control	0	-0.020270	NaN	-0.169660	5.965608e-01	56.0	False	none
1	therapy_b vs. control	0	0.027580	NaN	-0.114804	3.666954e-01	56.0	False	none
2	therapy_combo_ab vs. control	0	0.050940	NaN	-0.088023	2.624107e-01	56.0	False	none
3	therapy_a vs. control	3	0.201454	NaN	0.014365	3.963471e-02	56.0	False	none
4	therapy_b vs. control	3	0.183111	NaN	-0.008275	5.683480e-02	56.0	False	none
5	therapy_combo_ab vs. control	3	0.264414	NaN	0.092076	8.938337e-03	56.0	False	none
6	therapy_a vs. control	6	0.262011	NaN	0.048127	2.536545e-02	56.0	False	none
7	therapy_b vs. control	6	0.199683	NaN	-0.032264	7.441462e-02	56.0	False	none
8	therapy_combo_ab vs. control	6	0.461990	NaN	0.306064	7.359078e-05	56.0	False	none
9	therapy_a vs. control	9	0.304954	NaN	0.052699	2.719380e-02	56.0	False	none
10	therapy_b vs. control	9	0.291002	NaN	0.033684	3.424308e-02	56.0	False	none
11	therapy_combo_ab vs. control	9	0.602631	NaN	0.458414	3.159851e-06	56.0	False	none
12	therapy_a vs. control	12	0.344803	NaN	0.056728	2.867643e-02	56.0	False	none
13	therapy_b vs. control	12	0.432299	NaN	0.182694	5.973560e-03	56.0	False	none
14	therapy_combo_ab vs. control	12	0.713108	NaN	0.586968	2.075028e-07	56.0	False	none
15	therapy_a vs. control	15	0.528911	NaN	0.306588	9.585859e-04	56.0	False	none
16	therapy_b vs. control	15	0.620894	NaN	0.441981	4.888934e-05	56.0	False	none
17	therapy_combo_ab vs. control	15	0.847676	NaN	0.775790	2.263913e-11	56.0	False	none
18	therapy_a vs. control	18	0.586554	NaN	0.370586	4.394189e-04	56.0	False	none
19	therapy_b vs. control	18	0.753773	NaN	0.625154	3.658760e-07	56.0	False	none
20	therapy_combo_ab vs. control	18	0.918150	NaN	0.875395	2.636353e-14	56.0	False	none

	TGI_delta_arithmetic_mean
group	control	therapy_a	therapy_b	therapy_combo_ab
day
0	NaN	NaN	NaN	NaN
3	0.0	0.807033	0.555618	0.905329
6	0.0	0.536529	0.332937	0.851275
9	0.0	0.445318	0.386935	0.841935
12	0.0	0.417793	0.489568	0.838572
15	0.0	0.590634	0.632277	0.908107
18	0.0	0.648716	0.775380	0.943067

	contrast	coef	hr	se(coef)	z	pvalue	lb95_hr	up95_hr
0	therapy_a vs. control	-2.636960	0.071579	0.459606	-5.737436	9.612043e-09	0.029078	0.176198
1	therapy_b vs. control	-5.812376	0.002990	0.785134	-7.403039	1.331023e-13	0.000642	0.013932

	contrast	day	tgi	lb95_tgi	ub95_tgi	pvalue	df	baseline_adj	correction
0	therapy_a vs. control	3	0.217298	NaN	0.079788	3.642892e-03	55.0	True	dunnett
1	therapy_b vs. control	3	0.159973	NaN	0.012312	1.263017e-04	55.0	True	dunnett
2	therapy_combo_ab vs. control	3	0.224986	NaN	0.088368	2.331468e-15	55.0	True	dunnett
3	therapy_a vs. control	6	0.276538	NaN	0.100631	3.642892e-03	55.0	True	dunnett
4	therapy_b vs. control	6	0.177198	NaN	-0.022975	1.263017e-04	55.0	True	dunnett
5	therapy_combo_ab vs. control	6	0.433388	NaN	0.295142	2.331468e-15	55.0	True	dunnett
6	therapy_a vs. control	9	0.319469	NaN	0.102286	3.642892e-03	55.0	True	dunnett
7	therapy_b vs. control	9	0.269836	NaN	0.036678	1.263017e-04	55.0	True	dunnett
8	therapy_combo_ab vs. control	9	0.580168	NaN	0.445706	2.331468e-15	55.0	True	dunnett
9	therapy_a vs. control	12	0.358360	NaN	0.099632	3.642892e-03	55.0	True	dunnett
10	therapy_b vs. control	12	0.415513	NaN	0.179689	1.263017e-04	55.0	True	dunnett
11	therapy_combo_ab vs. control	12	0.697046	NaN	0.574437	2.331468e-15	55.0	True	dunnett
12	therapy_a vs. control	15	0.538075	NaN	0.333737	3.642892e-03	55.0	True	dunnett
13	therapy_b vs. control	15	0.610371	NaN	0.437908	1.263017e-04	55.0	True	dunnett
14	therapy_combo_ab vs. control	15	0.839677	NaN	0.768492	2.331468e-15	55.0	True	dunnett
15	therapy_a vs. control	18	0.593895	NaN	0.390426	3.642892e-03	55.0	True	dunnett
16	therapy_b vs. control	18	0.747548	NaN	0.620984	1.263017e-04	55.0	True	dunnett
17	therapy_combo_ab vs. control	18	0.914239	NaN	0.871108	2.331468e-15	55.0	True	dunnett

	contrast	day	tgi	lb95_tgi	ub95_tgi	pvalue	df	baseline_adj	correction
0	therapy_a vs. therapy_b	3	0.081024	0.234103	-0.102650	0.349802	27.0	True	none
1	therapy_a vs. therapy_b	6	0.127909	0.315878	-0.111707	0.257496	27.0	True	none
2	therapy_a vs. therapy_b	9	0.074951	0.322581	-0.263200	0.612064	27.0	True	none
3	therapy_a vs. therapy_b	12	-0.075091	0.252544	-0.546339	0.685970	27.0	True	none
4	therapy_a vs. therapy_b	15	-0.159271	0.225560	-0.735330	0.458732	27.0	True	none
5	therapy_a vs. therapy_b	18	-0.576189	-0.061831	-1.339706	0.025560	27.0	True	none

	contrast	time_window	tgr	lb95_tgr	ub95_tgr	pvalue
0	therapy_a vs. control	0-18	-0.046	NaN	-0.022972	3.131057e-04
1	therapy_b vs. control	0-18	-0.071	NaN	-0.047972	6.354270e-08
2	therapy_combo_ab vs. control	0-18	-0.132	NaN	-0.108972	5.629282e-23

	day	animal	group	tv
0	0	control_1	control	107.670789
1	3	control_1	control	213.533876
2	6	control_1	control	250.344526
3	9	control_1	control	486.383178
4	12	control_1	control	593.092860
5	15	control_1	control	1420.028888

	day	animal	group	tv
0	0	control_1	control	115.507965
1	3	control_1	control	330.419462
2	6	control_1	control	321.841189
3	9	control_1	control	386.864409
4	12	control_1	control	642.549477
5	15	control_1	control	1411.567508

	contrast	Gompertz_par	lb95	ub95	pvalue
0	therapy_a vs. control	0.010	NaN	0.016579	0.002998
1	therapy_b vs. control	0.016	NaN	0.022579	0.000006

	contrast	AUC_diff	lb95_AUC_diff	ub95_AUC_diff	pvalue
0	therapy_a vs. control	-7.286309	NaN	-3.828678	4.266910e-04
1	therapy_b vs. control	-7.996176	NaN	-4.538545	1.442333e-04
2	therapy_combo_ab vs. control	-15.740844	NaN	-12.283213	1.672120e-10

	contrast	chisq	pvalue
0	therapy_a vs. control	36.326328	1.668940e-09
1	therapy_b vs. control	42.811232	6.028503e-11

	contrast	day	synergy	lb95_synergy	ub95_synergy	pvalue	df	baseline_adj
0	therapy_combo_ab vs. Bliss	0	0.956588	0.788495	1.160515	0.351171	56	False
1	therapy_combo_ab vs. Bliss	3	1.127641	0.837315	1.518633	0.748753	56	False
2	therapy_combo_ab vs. Bliss	6	0.910917	0.635576	1.305541	0.333133	56	False
3	therapy_combo_ab vs. Bliss	9	0.806371	0.520427	1.249424	0.207284	56	False
4	therapy_combo_ab vs. Bliss	12	0.771307	0.460687	1.291362	0.201489	56	False
5	therapy_combo_ab vs. Bliss	15	0.852911	0.493712	1.473442	0.314173	56	False
6	therapy_combo_ab vs. Bliss	18	0.804011	0.443755	1.456737	0.270894	56	False

	contrast	day	synergy	lb95_synergy	ub95_synergy	pvalue	df	baseline_adj
0	therapy_combo_ab vs. Bliss	3	0.991736	0.795169	1.236895	0.475063	56	True
1	therapy_combo_ab vs. Bliss	6	1.171574	0.827840	1.658032	0.775551	56	True
2	therapy_combo_ab vs. Bliss	9	0.951826	0.642334	1.410438	0.417224	56	True
3	therapy_combo_ab vs. Bliss	12	0.669722	0.403424	1.111804	0.095637	56	True
4	therapy_combo_ab vs. Bliss	15	0.451375	0.269701	0.755428	0.006210	56	True
5	therapy_combo_ab vs. Bliss	18	0.339773	0.203861	0.566293	0.000414	56	True

Experimental design considerations and statistical analyses in preclinical tumor growth inhibition studies - Python tutorial¶

Introduction¶

Objectives¶

Toy sample generation¶

Analyzing tumor growth dynamics¶

Visualizing tumor growth curves¶

Customized plotting function¶

Versatile visualization with seaborn¶

Visualizing individual tumor growth trajectories¶

Experimental design¶

Overview of animal models in TGI studies¶

Endpoints in preclinical TGI studies¶

Quantifying treatment effects¶

Enhanced data analysis with Pandas¶

Measuring tumor volume in preclinical studies¶

Understanding and mitigating measurement errors in TGI studies¶

Handling discrepant measurements and missing data in TGI studies¶

Statistical analysis in preclinical cancer studies¶

Analyzing continuous tumor volume data at a single time point¶

How it works¶

Model-based estimation of the TGI¶

Practical example¶

TGI analysis for single time-point data in Python¶

The problem with multiple pairwise t-tests (j > 2)¶

Note on the regression model formula¶

Pre-calculating individual TGI-delta and running OLS without adjustment¶

Going beyond simple comparisons: ANCOVA and non-parametric methods¶

Analysis of covariance (ANCOVA)¶

Handling unequal variances¶

Non-parametric methods¶

Analyzing continuous data longitudinally¶

Hypothesis testing¶

TGI analysis over multiple time points¶

Why the degrees of freedom don't increase as expected¶

Analyzing longitudinal TGI data¶

Modeling tumor growth curves¶

The mathematics behing MEMs¶

Random effects: capturing individual animal variations¶

Between-animal vs. within-animal variation¶

Advantages of MEMs¶

Estimating TGR with MEMs with Python¶

Statistical distributions¶

The interaction effect¶

Examples¶

Modeling tumor growth with quadratic exponential models¶

Analyzing quadratic exponential tumor growth¶

Modeling long-term tumor dynamics with the Gompertzian model¶

Long-term TGI studies and the Gompertzian model¶

Summarizing tumor growth with AUC¶

AUC analysis with Python¶

Linear models: the workhorses of longitudinal TGI analysis¶

Time-to-event analysis in preclinical cancer studies¶

Kaplan-Meier and cumulative incidence¶

The Cox Proportional Hazards (CPH) model and beyond¶

Using Python for TTE analysis in TGI studies¶

Analysis of categorical endpoints in TGI studies¶

Analysis of synergistic effects¶

HSA¶

Bliss¶

Bliss independence synergy¶

Visualization¶

Inference¶

Powering studies¶

Statistical power¶

Calculating sample size for two-sample t-tests¶

Calculating sample size for ANOVA and ANCOVA¶

Conclusion¶

Session Information¶

Versatile visualization with `seaborn`¶