Logistic Regression: Biological Applications

This page covers logistic regression applied to microbiome and clinical data. The core method is described in the main logistic regression page.

A recurring challenge in biological applications is that the number of predictors (taxa, genes, metabolites) far exceeds the number of samples. Standard logistic regression breaks down in this setting. LASSO regularisation handles this by shrinking most coefficients to zero, effectively selecting a sparse set of predictors that best discriminate the groups.

---

Microbiome Disease Classification with LASSO

LASSO logistic regression is well suited to microbiome data: it performs variable selection automatically, handles compositional data after CLR transformation, and produces a sparse, interpretable model.

library(phyloseq)
library(microbiome)
library(glmnet)
library(pROC)

data(GlobalPatterns)
gp <- GlobalPatterns

# 1. Define binary outcome: human vs. non-human samples
sample_data(gp)$is_human <- ifelse(
  sample_data(gp)$SampleType %in% c("Feces", "Skin", "Tongue"), 1, 0)

# 2. Filter rare taxa
gp_filtered <- filter_taxa(gp,
  function(x) sum(x > 3) > (0.1 * length(x)), TRUE)

# 3. CLR transformation
gp_clr <- microbiome::transform(gp_filtered, "clr")
otu_clr <- t(as(otu_table(gp_clr), "matrix"))
outcome <- sample_data(gp_clr)$is_human

# 4. LASSO logistic regression with cross-validated lambda
set.seed(123)
cv_fit <- cv.glmnet(otu_clr, outcome,
                    family  = "binomial",
                    alpha   = 1,          # alpha = 1 for LASSO
                    nfolds  = 10)

plot(cv_fit)
# Left dashed line: lambda.min (minimum CV error)
# Right dashed line: lambda.1se (most regularised model within 1 SE)

# 5. Discriminant taxa at lambda.min
coef_mat <- coef(cv_fit, s = "lambda.min")
selected  <- rownames(coef_mat)[coef_mat[, 1] != 0]
selected  <- selected[selected != "(Intercept)"]
cat("Selected taxa:", length(selected), "\n")
print(coef_mat[selected, ])

Positive coefficients correspond to taxa enriched in human samples; negative coefficients correspond to taxa enriched in environmental samples. The magnitude reflects the strength of association on the log-odds scale.

# 6. Model performance
pred_probs <- predict(cv_fit, newx = otu_clr,
                      s = "lambda.min", type = "response")

roc_obj <- roc(outcome, pred_probs[, 1])
plot(roc_obj, print.auc = TRUE, col = "steelblue", lwd = 2,
     main = "LASSO: Human vs. Environmental Microbiome")
auc(roc_obj)

Training AUC is optimistic

The AUC above is computed on the same data used to fit the model. For an honest performance estimate, use a held-out test set or nested cross-validation.

---

Clinical Classification Workflow

A complete workflow for a clinical dataset with multiple continuous and categorical predictors:

library(caret)
library(pROC)
library(car)

# --- Data preparation ---
# Replace with your own data
data   <- your_clinical_data
outcome_var <- "disease"

# Check outcome balance
table(data[[outcome_var]])
# Severe imbalance (< 1:5) may require resampling strategies

# --- Train / test split ---
set.seed(123)
train_idx  <- createDataPartition(data[[outcome_var]], p = 0.7, list = FALSE)
train_data <- data[train_idx, ]
test_data  <- data[-train_idx, ]

# --- Fit model ---
model <- glm(disease ~ age + bmi + smoking_status + crp,
             data   = train_data,
             family = binomial)

# Check for warnings: convergence failure or complete separation
summary(model)

# Multicollinearity
vif(model)   # Values above 5 warrant attention

# Odds ratios with confidence intervals
exp(cbind(OR = coef(model), confint(model)))

# --- Cross-validation on training set ---
set.seed(123)
cv_result <- train(
  factor(disease) ~ age + bmi + smoking_status + crp,
  data      = train_data,
  method    = "glm",
  family    = binomial,
  trControl = trainControl(
    method          = "cv",
    number          = 10,
    classProbs      = TRUE,
    summaryFunction = twoClassSummary),
  metric = "ROC")

print(cv_result)

# --- Test set evaluation ---
test_probs <- predict(model, newdata = test_data, type = "response")
test_class <- ifelse(test_probs > 0.5, 1, 0)

confusionMatrix(factor(test_class),
                factor(test_data$disease),
                positive = "1")

roc_test <- roc(test_data$disease, test_probs)
plot(roc_test, print.auc = TRUE, col = "steelblue", lwd = 2)

---

Handling Imbalanced Data

When one class is rare (for example, a disease with 5% prevalence), a model that always predicts "healthy" achieves 95% accuracy but has no clinical value. Several strategies help:

Adjust the classification threshold. Rather than 0.5, choose a threshold that optimises sensitivity and specificity for your application. The ROC curve shows all possible tradeoffs:

# Find threshold maximising Youden's J (sensitivity + specificity - 1)
coords(roc_obj, "best", ret = c("threshold", "sensitivity", "specificity"))

Adjust prior probabilities. If the training set prevalence does not match the expected real-world prevalence, correct for this:

# Set prior reflecting true population prevalence
model_weighted <- glm(disease ~ .,
                      data    = train_data,
                      family  = binomial,
                      weights = ifelse(train_data$disease == 1,
                                       10,   # upweight rare class
                                       1))

Use LASSO with class weighting for high-dimensional imbalanced data:

weights <- ifelse(outcome == 1,
                  sum(outcome == 0) / sum(outcome == 1),
                  1)

cv_fit <- cv.glmnet(X, outcome,
                    family  = "binomial",
                    alpha   = 1,
                    weights = weights)

---

Complete Separation

Complete separation occurs when a predictor (or combination of predictors) perfectly predicts the outcome in the training data. The log-odds coefficient for that predictor becomes infinite and glm() issues a warning about fitted probabilities of 0 or 1.

This is common with small datasets or highly informative biomarkers. Firth's penalised likelihood provides a stable solution:

library(logistf)

model_firth <- logistf(disease ~ biomarker + age, data = df)
summary(model_firth)
exp(coef(model_firth))   # Odds ratios