Instalar com Codex ou Claude Copie este prompt, cole no Codex, Claude ou outro assistente e deixe que ele revise a página da skill e instale para você.
Um comando direto ignora o prompt de revisão. Verifique a origem antes de executá-lo.
✅ Read full question → Actual outcome = treatment response (PR vs non-PR)
Always check data columns first: print(df.columns.tolist())
import pandas as pd
import numpy as np
# Load data
df = pd.read_csv('data.csv')
# Check structureprint(f"Observations: {len(df)}")
print(f"Variables: {len(df.columns)}")
print(f"Missing: {df.isnull().sum().sum()}")
# Detect variable typesfor col in df.columns:
n_unique = df[col].nunique()
if n_unique == 2:
print(f"{col}: binary")
elif n_unique <= 10and df[col].dtype == 'object':
print(f"{col}: categorical ({n_unique} levels)")
elif df[col].dtype in ['float64', 'int64']:
print(f"{col}: continuous (mean={df[col].mean():.2f})")
Phase 1: Model Fitting
Goal: Fit appropriate model based on outcome type.
Linear Regression
import statsmodels.formula.api as smf
# R-style formula (recommended)
model = smf.ols('outcome ~ predictor1 + predictor2 + age', data=df).fit()
# Resultsprint(f"R-squared: {model.rsquared:.4f}")
print(f"AIC: {model.aic:.2f}")
print(model.summary())
Logistic Regression
# Fit model
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)
# Odds ratios
ors = np.exp(model.params)
ci = np.exp(model.conf_int())
for var in ['exposure', 'age', 'sex_M']:
print(f"{var}: OR={ors[var]:.4f}, CI=({ci.loc[var, 0]:.4f}, {ci.loc[var, 1]:.4f})")
Ordinal Logistic Regression
from statsmodels.miscmodels.ordinal_model import OrderedModel
# Prepare ordered outcome
severity_order = ['Mild', 'Moderate', 'Severe']
df['severity'] = pd.Categorical(df['severity'], categories=severity_order, ordered=True)
y = df['severity'].cat.codes
# Fit model
X = pd.get_dummies(df[['exposure', 'age', 'sex']], drop_first=True, dtype=float)
model = OrderedModel(y, X, distr='logit').fit(method='bfgs', disp=0)
# Odds ratios
ors = np.exp(model.params[:len(X.columns)])
print(f"Exposure OR: {ors[0]:.4f}")
Cox Proportional Hazards
from lifelines import CoxPHFitter
# Fit model
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age']],
duration_col='time', event_col='event')
# Hazard ratiosprint(f"HR (treatment): {cph.hazard_ratios_['treatment']:.4f}")
print(f"Concordance: {cph.concordance_index_:.4f}")
See references/ for detailed examples of each model type.
When data has multiple features (genes, miRNAs, metabolites, etc.), there are TWO approaches:
Question: "What is the F-statistic comparing [feature] expression across groups?"
DECISION TREE:
│
├─ Does question specify "the F-statistic" (singular)?
│ │
│ ├─ YES, singular → Likely asking for SPECIFIC FEATURE(S) F-statistic
│ │ │
│ │ ├─ Are there thousands of features (genes, miRNAs)?
│ │ │ YES → Per-feature approach (Method B below)
│ │ │
│ │ └─ Is there one feature of interest?
│ │ YES → Single feature ANOVA (Method A below)
│ │
│ └─ NO, asks about "all features" or "genes" (plural)?
│ YES → Aggregate approach or per-feature summary
│
└─ When unsure: Calculate PER-FEATURE and report summary statistics
Method A: Aggregate ANOVA (all features combined)
Use when: Testing overall expression differences across all features
Result: Single F-statistic representing global effect
# Flatten all features across all samples per group
groups_agg = []
for celltype in ['CD4', 'CD8', 'CD14']:
samples = df[df['celltype'] == celltype]
# Flatten: all features × all samples in this group
all_values = expression_matrix.loc[:, samples.index].values.flatten()
groups_agg.append(all_values)
f_stat_agg, p_value = stats.f_oneway(*groups_agg)
print(f"Aggregate F-statistic: {f_stat_agg:.4f}")
# Result: Very large F-statistic (e.g., 153.8)
Use when: Testing EACH feature individually (most common in genomics)
Result: Distribution of F-statistics (one per feature)
# Calculate F-statistic FOR EACH FEATURE separately
per_feature_f_stats = []
for feature in expression_matrix.index: # For each gene/miRNA/metabolite
groups = []
for celltype in ['CD4', 'CD8', 'CD14']:
samples = df[df['celltype'] == celltype]
# Get expression of THIS feature in THIS cell type
values = expression_matrix.loc[feature, samples.index].values
groups.append(values)
f_stat, _ = stats.f_oneway(*groups)
ifnot np.isnan(f_stat):
per_feature_f_stats.append((feature, f_stat))
# Summary statistics
f_values = [f for _, f in per_feature_f_stats]
print(f"Per-feature F-statistics:")
print(f" Median: {np.median(f_values):.4f}")
print(f" Mean: {np.mean(f_values):.4f}")
print(f" Range: [{np.min(f_values):.4f}, {np.max(f_values):.4f}]")
# Find features in specific range (e.g., 0.76-0.78)
target_features = [(name, f) for name, f in per_feature_f_stats
if0.76 <= f <= 0.78]
if target_features:
print(f"Features with F ∈ [0.76, 0.78]: {len(target_features)}")
for name, f_val in target_features:
print(f" {name}: F = {f_val:.6f}")
Key Differences:
Aspect
Method A (Aggregate)
Method B (Per-feature)
Interpretation
Overall expression difference
Feature-specific differences
Result
1 F-statistic
N F-statistics (N = # features)
Typical value
Very large (e.g., 153.8)
Small to large (e.g., 0.1 to 100+)
Use case
Global effect size
Gene/biomarker discovery
Common in
Rarely used
Genomics, proteomics, metabolomics ⭐
Real-world example (BixBench bix-36-q1):
Question: "What is the F-statistic comparing miRNA expression across immune cell types?"
Expected: 0.76-0.78
Method A (aggregate): 153.836 ❌ WRONG
Method B (per-miRNA): Found 2 miRNAs with F ∈ [0.76, 0.78] ✅ CORRECT
Default assumption for gene expression data: Use Method B (per-feature)
Phase 2: Model Diagnostics
Goal: Check model assumptions and fit quality.
OLS Diagnostics
from scipy import stats as scipy_stats
from statsmodels.stats.diagnostic import het_breuschpagan
# Residual normality
residuals = model.resid
sw_stat, sw_p = scipy_stats.shapiro(residuals)
print(f"Shapiro-Wilk: p={sw_p:.4f} (normal if p>0.05)")
# Heteroscedasticity
bp_stat, bp_p, _, _ = het_breuschpagan(residuals, model.model.exog)
print(f"Breusch-Pagan: p={bp_p:.4f} (homoscedastic if p>0.05)")
# VIF (multicollinearity)from statsmodels.stats.outliers_influence import variance_inflation_factor
X = model.model.exog
for i inrange(1, X.shape[1]): # Skip intercept
vif = variance_inflation_factor(X, i)
print(f"{model.model.exog_names[i]}: VIF={vif:.2f}")
Proportional Hazards Test
# Test PH assumption for Cox model
results = cph.check_assumptions(df, p_value_threshold=0.05, show_plots=False)
iflen(results) == 0:
print("✅ Proportional hazards assumption met")
else:
print(f"⚠️ PH violated for: {results}")
See references/troubleshooting.md for common diagnostic issues.
Phase 3: Interpretation
Goal: Generate publication-quality summary.
Odds Ratio Interpretation
definterpret_odds_ratio(or_val, ci_lower, ci_upper, p_value):
"""Interpret odds ratio with clinical meaning."""if or_val > 1:
pct_increase = (or_val - 1) * 100
direction = f"{pct_increase:.1f}% increase in odds"else:
pct_decrease = (1 - or_val) * 100
direction = f"{pct_decrease:.1f}% decrease in odds"
sig = "significant"if p_value < 0.05else"not significant"
ci_contains_null = ci_lower <= 1 <= ci_upper
returnf"{direction} (OR={or_val:.4f}, 95% CI [{ci_lower:.4f}, {ci_upper:.4f}], p={p_value:.6f}, {sig})"
Common BixBench Patterns
Pattern 1: Odds Ratio from Ordinal Regression
Question: "What is the odds ratio of disease severity associated with exposure?"
Solution:
Identify ordinal outcome (mild/moderate/severe)
Fit ordinal logistic regression (proportional odds model)
Extract OR = exp(coefficient for exposure)
Report with CI and p-value
Pattern 2: Percentage Reduction in Odds
Question: "What is the percentage reduction in OR after adjusting for confounders?"
Question: "What is the odds ratio for the interaction between A and B?"
Solution:
# Fit model with interaction
model = smf.logit('outcome ~ A * B + age', data=df).fit(disp=0)
# Interaction OR
interaction_coef = model.params['A:B']
interaction_or = np.exp(interaction_coef)
print(f"Interaction OR: {interaction_or:.4f}")
Pattern 4: Survival Analysis
Question: "What is the hazard ratio for treatment?"
Solution:
Load survival data (time, event, covariates)
Fit Cox proportional hazards model
Extract HR = exp(coefficient)
Report with CI and concordance index
Pattern 5: Multi-feature ANOVA (Gene Expression)
Question: "What is the F-statistic comparing miRNA expression across cell types?"
Solution:
Identify that data has multiple features (genes/miRNAs)
Use per-feature ANOVA (NOT aggregate)
Calculate F-statistic for EACH feature separately
If question asks for "the F-statistic" (singular):
Check if specific features match expected range
Report those feature(s) F-statistics
If question asks for summary: report median/mean/distribution
Critical: For gene expression data, default to per-feature ANOVA. Aggregate ANOVA gives F-statistics ~200× larger and is rarely correct.
See references/bixbench_patterns.md for 15+ question patterns.
Statsmodels vs Scikit-learn
Use Case
Library
Reason
Inference (p-values, CIs, ORs)
statsmodels
Full statistical output
Prediction (accuracy, AUC)
scikit-learn
Better prediction tools
Mixed-effects models
statsmodels
Only option
Regularization (LASSO, Ridge)
scikit-learn
Better optimization
Survival analysis
lifelines
Specialized library
General rule: Use statsmodels for BixBench questions (they ask for p-values, ORs, HRs).
Confounders addressed: Adjusted analyses if applicable
Input Validation
This skill accepts requests that match the documented purpose of tooluniverse-statistical-modeling and include enough context to complete the workflow safely.
Do not continue the workflow when the request is out of scope, missing a critical input, or would require unsupported assumptions. Instead respond:
tooluniverse-statistical-modeling only handles its documented workflow. Please provide the missing required inputs or switch to a more suitable skill.
