drug-discovery-prompts/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-statistical-modeling/references/troubleshooting.md

title, task, lineage_type, upstream_source, upstream_sha, imported_at, prompt_class, upstream_changes, author, validated

title	task	lineage_type	upstream_source	upstream_sha	imported_at	prompt_class	upstream_changes	author	validated
Troubleshooting Guide		import	https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-statistical-modeling/references/troubleshooting.md	e2520a96	2026-06-26	prompt	accepted	upstream	false

Troubleshooting Guide

Common statistical modeling issues and solutions.

Convergence Issues

Problem: Model doesn't converge

Symptoms:

Warning: Maximum iterations reached
ConvergenceWarning: Maximum Likelihood optimization failed to converge

Solutions:

1. Increase max iterations:

# Logistic regression
model = smf.logit('y ~ x + z', data=df).fit(disp=0, maxiter=500)

# Ordinal logit
model = OrderedModel(y, X, distr='logit').fit(method='bfgs', maxiter=500)

# Cox model (lifelines handles this automatically)

2. Scale predictors:

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = pd.DataFrame(scaler.fit_transform(X), columns=X.columns)
model = OrderedModel(y, X_scaled, distr='logit').fit(method='bfgs')

3. Try different optimization method:

# Try Nelder-Mead instead of BFGS
model = OrderedModel(y, X, distr='logit').fit(method='nm')

4. Check for separation (see below)

---

Separation Problems

Problem: Perfect or quasi-complete separation

Symptoms:

• Very large coefficients (>10)
• Very large standard errors
• Warning: "Perfect separation detected"

Check for separation:

# Crosstab of predictor vs outcome
for pred in ['exposure', 'treatment']:
    print(f"\n{pred} by outcome:")
    print(pd.crosstab(df[pred], df['outcome']))

# Look for cells with 0 counts - that's separation

Solutions:

1. Remove problematic predictor:

# If one predictor causes separation, exclude it
model = smf.logit('outcome ~ age + sex', data=df).fit(disp=0)  # Exclude 'exposure'

2. Use Firth logistic regression (penalized likelihood):

# Requires logistf package (not standard)
# Alternative: use Ridge penalty in sklearn

from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(penalty='l2', C=1.0, solver='lbfgs')
lr.fit(X, y)

3. Collapse categories:

# If categorical predictor has sparse levels, combine them
df['stage_collapsed'] = df['stage'].replace({'I': 'Early', 'II': 'Early',
                                               'III': 'Late', 'IV': 'Late'})

4. Increase sample size (if possible)

---

Multicollinearity

Problem: Predictors highly correlated

Symptoms:

• Large standard errors
• VIF > 10
• Coefficients change dramatically when adding/removing predictors

Check for multicollinearity:

from statsmodels.stats.outliers_influence import variance_inflation_factor

X = df[['age', 'bmi', 'weight', 'height']].copy()
X['const'] = 1

for i, col in enumerate(X.columns[:-1]):
    vif = variance_inflation_factor(X.values, i)
    print(f"{col}: VIF = {vif:.2f}")
    if vif > 10:
        print(f"  ⚠️  High multicollinearity")

Solutions:

1. Remove correlated predictors:

# Check correlation matrix
corr_matrix = df[['age', 'bmi', 'weight', 'height']].corr()
print(corr_matrix)

# Remove one of highly correlated pairs (r > 0.8)
# E.g., remove weight if weight and BMI are r=0.9
model = smf.ols('outcome ~ age + bmi', data=df).fit()  # Exclude weight

2. Use Ridge regression (L2 regularization):

from sklearn.linear_model import Ridge

ridge = Ridge(alpha=1.0)
ridge.fit(X, y)

3. Principal Component Analysis:

from sklearn.decomposition import PCA

pca = PCA(n_components=3)
X_pca = pca.fit_transform(X)
# Use principal components as predictors

---

Heteroscedasticity

Problem: Non-constant variance of residuals

Symptoms:

• Breusch-Pagan test p < 0.05
• Residual plot shows funnel shape

Check:

from statsmodels.stats.diagnostic import het_breuschpagan

residuals = model.resid
bp_stat, bp_p, _, _ = het_breuschpagan(residuals, model.model.exog)

if bp_p < 0.05:
    print("⚠️  Heteroscedasticity detected")

Solutions:

1. Robust standard errors:

model_robust = model.get_robustcov_results(cov_type='HC3')
print(model_robust.summary())

2. Log transformation (if outcome is right-skewed):

df['outcome_log'] = np.log(df['outcome'] + 1)
model = smf.ols('outcome_log ~ x + z', data=df).fit()

3. Weighted Least Squares:

# Weight by inverse variance
residuals_sq = model.resid ** 2
weights = 1 / residuals_sq
wls_model = sm.WLS(y, X, weights=weights).fit()

---

Non-Normality of Residuals

Problem: Residuals not normally distributed

Symptoms:

• Shapiro-Wilk test p < 0.05
• Q-Q plot deviates from line

Check:

from scipy import stats as scipy_stats

residuals = model.resid
sw_stat, sw_p = scipy_stats.shapiro(residuals)

if sw_p < 0.05:
    print("⚠️  Residuals not normally distributed")

Solutions:

1. Transform outcome:

# Log transformation
df['outcome_log'] = np.log(df['outcome'] + 1)

# Square root transformation
df['outcome_sqrt'] = np.sqrt(df['outcome'])

# Box-Cox transformation
from scipy.stats import boxcox
df['outcome_bc'], lambda_param = boxcox(df['outcome'] + 1)

2. Use robust regression:

from statsmodels.robust.robust_linear_model import RLM

rlm_model = RLM.from_formula('outcome ~ x + z', data=df).fit()

3. Use non-parametric methods:

# Bootstrap confidence intervals instead of t-tests

Note: For large samples (n > 30), non-normality is less critical due to Central Limit Theorem.

---

Missing Data

Problem: Missing values in predictors or outcome

Check:

# Count missing values
missing = df.isnull().sum()
print("\nMissing values:")
print(missing[missing > 0])

# Percentage missing
pct_missing = (missing / len(df) * 100)
print("\nPercentage missing:")
print(pct_missing[pct_missing > 0])

Solutions:

1. Complete case analysis (delete rows with missing):

df_complete = df.dropna(subset=['outcome', 'x', 'z'])
model = smf.ols('outcome ~ x + z', data=df_complete).fit()

2. Mean/median imputation (simple):

df['age'].fillna(df['age'].mean(), inplace=True)

3. Multiple imputation (best practice):

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

imputer = IterativeImputer(random_state=42)
df_imputed = pd.DataFrame(imputer.fit_transform(df),
                          columns=df.columns)

4. Missing indicator method:

# Create indicator for missingness
df['age_missing'] = df['age'].isnull().astype(int)
df['age'].fillna(df['age'].mean(), inplace=True)

# Include indicator in model
model = smf.ols('outcome ~ age + age_missing + x', data=df).fit()

---

Proportional Hazards Violation

Problem: PH assumption violated in Cox model

Check:

# Test PH assumption
results = cph.check_assumptions(df, p_value_threshold=0.05, show_plots=False)

if len(results) > 0:
    print(f"⚠️  PH violated for: {results}")

Solutions:

1. Stratify by problematic variable:

# Don't estimate HR for treatment, but adjust for it
cph_strat = CoxPHFitter()
cph_strat.fit(df, duration_col='time', event_col='event',
              strata=['treatment'])

2. Time-varying coefficients:

# Allow coefficient to change over time (advanced)
# Interact predictor with time
df['treatment_time'] = df['treatment'] * df['time']
cph.fit(df[['time', 'event', 'treatment', 'treatment_time', 'age']],
        duration_col='time', event_col='event')

3. Use parametric survival model:

from lifelines import WeibullAFTFitter

# Accelerated failure time model (no PH assumption)
wf = WeibullAFTFitter()
wf.fit(df, duration_col='time', event_col='event')

---

Proportional Odds Violation

Problem: PO assumption violated in ordinal logit

Check:

# Fit binary logits at each cutpoint, compare coefficients
# See ordinal_logistic.md for full test

Solutions:

1. Partial proportional odds model:

# Allow some predictors to vary across cutpoints (requires mord package)

2. Multinomial logistic regression:

from sklearn.linear_model import LogisticRegression

# Treat outcome as nominal (loses ordering information)
model = LogisticRegression(multi_class='multinomial', solver='lbfgs')
model.fit(X, y)

3. Adjacent category logit (alternative ordinal model)

---

Small Sample Size

Problem: Too few observations per predictor

Rule of thumb:

• Linear regression: ≥20 observations per predictor
• Logistic regression: ≥10 events per predictor
• Ordinal logit: ≥10 observations per outcome level per predictor
• Cox regression: ≥10 events per predictor

Check:

n = len(df)
n_predictors = len(['x', 'z', 'age', 'sex'])  # Your predictors

print(f"Observations per predictor: {n / n_predictors:.1f}")

# For logistic
n_events = df['outcome'].sum()
print(f"Events per predictor: {n_events / n_predictors:.1f}")

if n_events / n_predictors < 10:
    print("⚠️  Small sample size - results may be unreliable")

Solutions:

1. Reduce number of predictors:

# Only include most important predictors
model = smf.logit('outcome ~ exposure + age', data=df).fit(disp=0)

2. Use penalized regression:

from sklearn.linear_model import LogisticRegression

# Ridge penalty helps with small samples
lr = LogisticRegression(penalty='l2', C=1.0)
lr.fit(X, y)

3. Exact logistic regression (for very small samples):

# Requires R or specialized packages

---

Outliers and Influential Points

Problem: Outliers affecting model fit

Check:

from statsmodels.stats.outliers_influence import OLSInfluence

influence = model.get_influence()

# Cook's distance
cooks_d = influence.cooks_distance[0]
influential = cooks_d > 4 / len(df)

print(f"Influential points: {influential.sum()}")
print(f"Indices: {df.index[influential].tolist()}")

Solutions:

1. Remove outliers (if justified):

# Remove points with Cook's distance > 4/n
df_clean = df[~influential]
model_clean = smf.ols('outcome ~ x + z', data=df_clean).fit()

2. Robust regression:

from statsmodels.robust.robust_linear_model import RLM

# Downweights outliers automatically
rlm_model = RLM.from_formula('outcome ~ x + z', data=df).fit()

3. Winsorize extreme values:

from scipy.stats.mstats import winsorize

# Cap extreme values at 5th and 95th percentiles
df['outcome_wins'] = winsorize(df['outcome'], limits=[0.05, 0.05])

---

Model Selection Uncertainty

Problem: Unsure which predictors to include

Solutions:

1. Forward selection:

# Start with null model, add predictors one by one
# Keep if p < 0.05 or AIC improves

2. Backward elimination:

# Start with full model, remove predictors one by one
# Remove if p > 0.10 or AIC improves

3. LASSO for variable selection:

from sklearn.linear_model import LogisticRegressionCV

# LASSO automatically selects variables
lasso = LogisticRegressionCV(penalty='l1', solver='saga', cv=5)
lasso.fit(X, y)

# Non-zero coefficients are selected
selected = X.columns[lasso.coef_[0] != 0]
print(f"Selected variables: {selected.tolist()}")

4. Use domain knowledge:

# Always include clinically important confounders
# Age, sex are usually important in biomedical studies

---

Package-Specific Issues

statsmodels singular matrix error

Problem:

LinAlgError: singular matrix

Cause: Perfect multicollinearity (one predictor is linear combination of others)

Solution:

# Check correlation matrix
corr = X.corr()
print(corr)

# Remove one of perfectly correlated predictors
# Or use pd.get_dummies(..., drop_first=True)

lifelines convergence warning

Problem:

ConvergenceWarning: Newton-Raphson failed to converge

Solutions:

# 1. Check for separation
# 2. Scale predictors
# 3. Use robust=True
cph.fit(df, duration_col='time', event_col='event', robust=True)

---

Quick Diagnostic Checklist

Before finalizing analysis:

• Check for missing data
• Check variable distributions (outliers, skewness)
• Check for multicollinearity (VIF < 10)
• Check model convergence
• Check sample size adequacy
• Run residual diagnostics
• Test model assumptions
• Compare alternative models
• Perform sensitivity analyses
• Interpret results in context