Statistical Analysis

Topics: statistics, t-test, ANOVA, chi-square, correlation, regression

Statistical analysis is at the heart of medical research and clinical decision-making. It allows us to determine whether observed patterns in patient data are due to chance or represent meaningful differences or associations.

Descriptive Statistics

Descriptive statistics summarize your dataset and provide the foundation for deeper analysis.

In medicine, these metrics describe patient demographics, laboratory results, and treatment outcomes — helping you quickly understand what’s typical or unusual in your data.

🧪 Common Measures

  • Measures of central tendency:

    • Mean (e.g., average blood pressure)
    • Median (e.g., middle value of patient age)
    • Mode (e.g., most common diagnosis)
  • Measures of variability:

    • Range (e.g., youngest to oldest patient)
    • Standard deviation (e.g., how spread out glucose levels are)
  • Frequencies and proportions:

    • Counts and percentages (e.g., 60% of patients are female)

Let’s define a function to describe clinical data.

describe_data <- function(x, column_name) {
  # Generate comprehensive descriptive statistics.
  mode_val <- as.numeric(names(sort(table(x), decreasing = TRUE)[1]))
  iqr_val  <- IQR(x, na.rm = TRUE)

  cat(sprintf("\n=== Descriptive Statistics for %s ===\n", column_name))
  cat(sprintf("Count:              %d\n",    sum(!is.na(x))))
  cat(sprintf("Mean:               %.2f\n",  mean(x,   na.rm = TRUE)))
  cat(sprintf("Median:             %.2f\n",  median(x, na.rm = TRUE)))
  cat(sprintf("Mode:               %.2f\n",  mode_val))
  cat(sprintf("Standard Deviation: %.2f\n",  sd(x,     na.rm = TRUE)))
  cat(sprintf("Range:              %.2f\n",  diff(range(x, na.rm = TRUE))))
  cat(sprintf("IQR:                %.2f\n",  iqr_val))
}

Example Use

describe_data(df$Glucose,        "Fasting Glucose (mg/dL)")
describe_data(df$Age,            "Patient Age (years)")

Outputs:

=== Descriptive Statistics for Fasting Glucose (mg/dL) ===
Count:              100
Mean:               137.46
Median:             140.50
Mode:               96.00
Standard Deviation: 36.40
Range:              130.00
IQR:                63.50
=== Descriptive Statistics for Patient Age (years) ===
Count:              100
Mean:               46.90
Median:             47.50
Mode:               53.00
Standard Deviation: 16.91
Range:              57.00
IQR:                25.50

Univariate Tests

t-tests

One-sample t-test

Check if the average fasting glucose differs significantly from the normal fasting level of 100 mg/dL.

population_mean <- 100
result <- t.test(df$Glucose, mu = population_mean)

cat("\n=== One-Sample T-Test ===\n")
cat(sprintf("Testing if mean glucose differs from %d mg/dL\n", population_mean))
cat(sprintf("Sample mean:   %.2f\n",  mean(df$Glucose)))
cat(sprintf("T-statistic:   %.4f\n",  result$statistic))
cat(sprintf("P-value:       %.4f\n",  result$p.value))

Output

=== One-Sample T-Test ===
Testing if mean glucose differs from 100 mg/dL
Sample mean:   137.46
T-statistic:   10.2905
P-value:       0.0000

💡 Clinical meaning: A significant p-value (< 0.05) suggests that your patient group’s glucose levels are not within normal limits.

Two-sample t-test

Compare blood pressure between diabetic and hypertensive patients.

bp_diabetic     <- df$Blood_Pressure[df$Diagnosis == "Diabetes"]
bp_hypertensive <- df$Blood_Pressure[df$Diagnosis == "Hypertension"]

result <- t.test(bp_diabetic, bp_hypertensive)

cat("\n=== Two-Sample T-Test ===\n")
cat(sprintf("T-statistic: %.4f\n", result$statistic))
cat(sprintf("P-value:     %.4f\n", result$p.value))

Output

=== Two-Sample T-Test ===
T-statistic: -1.1210
P-value:     0.2662

💡 Clinical meaning: Tests if there’s a significant difference in blood pressure between the two groups.

Paired t-test example (before/after data)

Evaluate if blood pressure medication effectively reduces systolic BP.

set.seed(42)
bp_before <- sample(df$Blood_Pressure, 30)
bp_after  <- bp_before - rnorm(30, mean = 5, sd = 10)   # Simulated reduction

result <- t.test(bp_before, bp_after, paired = TRUE)

cat("\n=== Paired T-Test ===\n")
cat(sprintf("Mean before: %.2f\n", mean(bp_before)))
cat(sprintf("Mean after:  %.2f\n", mean(bp_after)))
cat(sprintf("P-value:     %.4f\n", result$p.value))

Output

=== Paired t-test ===
Mean before: 139.30
Mean after:  135.05
P-value:     0.0328

💡 Clinical meaning: Determines whether treatment significantly lowered blood pressure.

ANOVA (Analysis of Variance)

One-way Analysis of Variance (ANOVA)

ANOVA is an extension to the t-test, where there can be more than 2 groups being compared. Here we compare mean blood pressure across different diagnosis groups.

anova_result <- aov(Blood_Pressure ~ Diagnosis, data = df)
summary_result <- summary(anova_result)

f_stat  <- summary_result[[1]]$`F value`[1]
p_value <- summary_result[[1]]$`Pr(>F)`[1]

cat("\n=== One-Way ANOVA ===\n")
cat(sprintf("F-statistic: %.4f, P-value: %.4f\n", f_stat, p_value))

Output:

=== One-Way ANOVA ===
F-statistic: 1.8568, P-value: 0.1617

💡 Clinical meaning: Tests if blood pressure differs across diagnosis groups.

Post-hoc analysis is required when ANOVA is significant (p-value < 0.05).

if (p_value < 0.05) {
  cat("\nPost-hoc pairwise comparisons (Tukey HSD):\n")
  print(TukeyHSD(anova_result))
}

A p-value threshold of 0.05 is the conventional standard. Going above (e.g., 0.10) increases the risk of:

  • **Type I errors **— falsely rejecting a true null hypothesis
  • Results that are less likely to replicate in future studies
  • Conclusions that reviewers may view as tentative or inconclusive
Chi-Square Tests

Chi-square test of independence

Check if Diagnosis is associated with a binary risk factor (e.g., high blood pressure threshold).

df$High_BP <- ifelse(df$Blood_Pressure >= 140, "High", "Normal")

contingency <- table(df$Diagnosis, df$High_BP)
result      <- chisq.test(contingency)

cat("\n=== Chi-Square Test of Independence ===\n")
print(contingency)
cat(sprintf("P-value: %.4f\n", result$p.value))

Output:

=== Chi-Square Test of Independence ===
              
               High Normal
  Diabetes       18     18
  Healthy         8     15
  Hypertension   27     14
P-value: 0.0520

💡 Clinical meaning: A significant association suggests high blood pressure may relate to diagnosis.

Chi-square goodness of fit test

Test if diagnosis distribution follows expected population proportions.

observed            <- table(df$Diagnosis)
expected_proportions <- c(Diabetes = 0.33, Healthy = 0.34, Hypertension = 0.33)
expected            <- expected_proportions * nrow(df)

result <- chisq.test(observed, p = expected_proportions)
cat("\n=== Chi-Square Goodness of Fit ===\n")
cat(sprintf("P-value: %.4f\n", result$p.value))

Output:

=== Chi-Square Goodness of Fit ===
P-value: 0.0558

Correlation and Regression Analysis

Correlation

💡 Clinical meaning: Identify linear relationships (e.g., Age positively correlating with blood pressure).

cat("\n=== Correlation Analysis ===\n")
numeric_cols  <- df |> select(Age, Glucose, Blood_Pressure)
corr_matrix   <- cor(numeric_cols, use = "complete.obs")
print(round(corr_matrix, 3))

Output:

=== Correlation Analysis ===
                  Age Glucose Blood_Pressure
Age             1.000  -0.160         -0.103
Glucose        -0.160   1.000          0.195
Blood_Pressure -0.103   0.195          1.000
Simple Linear Regression

Predict Glucose from Age.

model   <- lm(Glucose ~ Age, data = df)
summary_m <- summary(model)

slope     <- coef(model)["Age"]
r_squared <- summary_m$r.squared
p_value   <- summary_m$coefficients["Age", "Pr(>|t|)"]

cat("\n=== Linear Regression ===\n")
cat(sprintf("Slope: %.2f, R²: %.3f, P-value: %.4f\n", slope, r_squared, p_value))

Output:

=== Linear Regression ===
Slope: -0.34, R²: 0.026, P-value: 0.1125

💡 Clinical meaning: Each 1-year increase in age predicts an average change of -0.34 mg/dL in glucose.

Logistic Regression

Predict presence of high blood pressure (Yes/No) using Age.

Load packages

library(dplyr)

Prepare data

df <- df |>
  mutate(Hypertensive = as.integer(Blood_Pressure >= 140))

set.seed(42)
train_idx <- sample(seq_len(nrow(df)), size = floor(0.7 * nrow(df)))
train_df  <- df[train_idx, ]
test_df   <- df[-train_idx, ]

Fit logistic regression model

log_model <- glm(Hypertensive ~ Age + Glucose, data = train_df, family = binomial)

Predict and evaluate

pred_prob  <- predict(log_model, newdata = test_df, type = "response")
pred_class <- ifelse(pred_prob > 0.5, 1, 0)

accuracy <- mean(pred_class == test_df$Hypertensive)

cat("\n=== Logistic Regression: Predicting Hypertension ===\n")
cat(sprintf("Accuracy: %.3f\n", accuracy))

cat("\nConfusion Matrix:\n")
print(table(Predicted = pred_class, Actual = test_df$Hypertensive))

Output:

=== Logistic Regression: Predicting Hypertension ===
Accuracy: 0.500

Confusion Matrix:
         Actual
Predicted  0  1
        0  5  5
        1 10 10