```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

## 1. A Story About Precision and Uncertainty

A public health researcher conducted a community-based nutrition intervention in a rural district, teaching families about balanced diets and meal planning. After six months, she measured the average change in daily vegetable consumption among 80 participating families. The sample showed an average increase of 1.8 servings per day—an encouraging finding that she planned to present to the regional health committee to secure funding for expansion.

However, when the committee chair asked, "How confident are you that this 1.8 servings represents the true population effect?", she realized her point estimate told only part of the story. While 1.8 servings was her best estimate from the sample, the true population mean could plausibly be anywhere from 1.3 to 2.3 servings (or perhaps even wider). This range—the **confidence interval**—communicates the precision of her estimate and the degree of uncertainty inherent in sampling.

Understanding confidence intervals is fundamental to statistical inference. They provide a range of plausible values for population parameters, acknowledging that our sample-based estimates are subject to sampling variability. Rather than offering a single number that falsely implies certainty, confidence intervals honestly communicate what we can reasonably conclude from our data.

## 2. What are Confidence Intervals?

A **confidence interval** is a range of values that likely contains an unknown population parameter with a specified level of confidence. Most commonly, we construct 95% confidence intervals, meaning that if we repeated our study many times and calculated a confidence interval each time, approximately 95% of those intervals would contain the true population parameter.

### Key Concepts

**Confidence Level (1 - α):** The probability that the interval contains the true parameter when the procedure is repeated many times. Common choices are 90%, 95%, and 99%.

**Margin of Error:** The distance from the point estimate to the endpoints of the interval, determined by the standard error and critical value.

**Width of the Interval:** Reflects precision—narrower intervals indicate more precise estimates. Width depends on:
- Sample size (larger n → narrower intervals)
- Variability in the data (larger SD → wider intervals)
- Confidence level (higher confidence → wider intervals)

### Common Misconceptions

**INCORRECT interpretation:** "There is a 95% probability that the true parameter lies within this specific interval."

**CORRECT interpretation:** "If we repeated this study many times and calculated 95% CIs each time, approximately 95% of those intervals would contain the true parameter."

The true parameter is fixed (though unknown); the interval is what varies across repeated samples.

## 3. Mathematical Foundations

### Confidence Interval for a Mean

For a population mean μ with known variance, the confidence interval is:

$$
\bar{X} \pm z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}}
$$

When the population variance is unknown (the usual case), we use the t-distribution:

$$
\bar{X} \pm t_{\alpha/2, df} \cdot \frac{s}{\sqrt{n}}
$$

where:
- $\bar{X}$ is the sample mean
- $t_{\alpha/2, df}$ is the critical value from the t-distribution with df = n - 1
- $s$ is the sample standard deviation
- $n$ is the sample size

### Confidence Interval for a Proportion

For a population proportion p:

$$
\hat{p} \pm z_{\alpha/2} \cdot \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}
$$

where $\hat{p}$ is the sample proportion.

### Confidence Interval for the Difference Between Two Means

For independent samples:

$$
(\bar{X}_1 - \bar{X}_2) \pm t_{\alpha/2, df} \cdot \sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}
$$

### Confidence Interval for Regression Coefficients

For a regression coefficient $\beta_j$:

$$
\hat{\beta}_j \pm t_{\alpha/2, df} \cdot SE(\hat{\beta}_j)
$$

where $SE(\hat{\beta}_j)$ is the standard error of the coefficient.

## 4. Simulating Datasets for Different Scenarios

We'll create multiple datasets to demonstrate confidence interval construction across various statistical contexts.

```{r simulate_data}
set.seed(2025)

# Scenario 1: Single mean (Nutrition intervention - vegetable servings)
# Daily vegetable servings after intervention
vegetable_servings <- rnorm(80, mean = 1.8, sd = 0.9)

nutrition_data <- data.frame(
  servings = vegetable_servings
)

# Scenario 2: Two-group comparison (Sleep study)
# Sleep quality scores (0-100 scale) for two interventions
meditation_group <- rnorm(50, mean = 72, sd = 12)
exercise_group <- rnorm(50, mean = 68, sd = 12)

sleep_data <- data.frame(
  quality = c(meditation_group, exercise_group),
  intervention = factor(rep(c("Meditation", "Exercise"), each = 50))
)

# Scenario 3: Proportion (Screening program)
# COVID-19 vaccination uptake in two communities
n_urban <- 200
n_rural <- 180
vaccinated_urban <- rbinom(1, n_urban, 0.78)
vaccinated_rural <- rbinom(1, n_rural, 0.65)

# Scenario 4: Regression (Stress and productivity)
n <- 120
stress_level <- rnorm(n, mean = 50, sd = 15)
productivity <- 85 - 0.6 * stress_level + rnorm(n, mean = 0, sd = 8)

stress_data <- data.frame(
  stress = stress_level,
  productivity = productivity
)

# Scenario 5: Correlation (Social support and well-being)
social_support <- rnorm(100, mean = 60, sd = 15)
wellbeing <- 30 + 0.7 * social_support + rnorm(100, mean = 0, sd = 10)

support_data <- data.frame(
  support = social_support,
  wellbeing = wellbeing
)
```

## 5. Visualizing Confidence Intervals

Before calculating CIs, let's visualize our data to understand the distributions.

```{r visualize, fig.width=10, fig.height=8}
library(ggplot2)

# Plot 1: Distribution of vegetable servings
p1 <- ggplot(nutrition_data, aes(x = servings)) +
  geom_histogram(bins = 20, fill = "#3498DB", alpha = 0.7, color = "black") +
  geom_vline(aes(xintercept = mean(servings)), color = "#E74C3C", 
             linetype = "dashed", linewidth = 1) +
  labs(title = "Distribution of Daily Vegetable Servings",
       x = "Servings per Day", y = "Frequency") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# Plot 2: Sleep quality by intervention
p2 <- ggplot(sleep_data, aes(x = intervention, y = quality, fill = intervention)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 3, 
               fill = "red", color = "black") +
  labs(title = "Sleep Quality by Intervention Type",
       x = "Intervention", y = "Sleep Quality Score") +
  scale_fill_manual(values = c("#E74C3C", "#27AE60")) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "none"
  )

# Plot 3: Stress and productivity
p3 <- ggplot(stress_data, aes(x = stress, y = productivity)) +
  geom_point(alpha = 0.6, color = "#9B59B6") +
  geom_smooth(method = "lm", se = TRUE, color = "#E74C3C", fill = "#E74C3C", alpha = 0.2) +
  labs(title = "Stress Level vs. Productivity",
       x = "Stress Level", y = "Productivity Score") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# Plot 4: Social support and wellbeing
p4 <- ggplot(support_data, aes(x = support, y = wellbeing)) +
  geom_point(alpha = 0.6, color = "#F39C12") +
  geom_smooth(method = "lm", se = TRUE, color = "#E74C3C", fill = "#E74C3C", alpha = 0.2) +
  labs(title = "Social Support vs. Well-being",
       x = "Social Support Score", y = "Well-being Score") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# Arrange plots
library(gridExtra)
grid.arrange(p1, p2, p3, p4, ncol = 2)
```

## 6. Calculating Confidence Intervals in R

Now we calculate confidence intervals for each scenario using actual R output.

### Confidence Interval for a Single Mean

```{r ci_mean}
# Calculate mean and CI using t.test
t_result <- t.test(nutrition_data$servings, conf.level = 0.95)

# Extract values
mean_servings <- mean(nutrition_data$servings)
sd_servings <- sd(nutrition_data$servings)
n_servings <- length(nutrition_data$servings)
se_servings <- sd_servings / sqrt(n_servings)

# Calculate CI manually to verify
df <- n_servings - 1
t_critical <- qt(0.975, df)
margin_error <- t_critical * se_servings
ci_lower_manual <- mean_servings - margin_error
ci_upper_manual <- mean_servings + margin_error

# Display results
cat("Confidence Interval for Mean Vegetable Servings:\n")
cat("Sample size:", n_servings, "\n")
cat("Mean:", round(mean_servings, 3), "\n")
cat("SD:", round(sd_servings, 3), "\n")
cat("SE:", round(se_servings, 4), "\n")
cat("t critical value (df =", df, "):", round(t_critical, 3), "\n")
cat("Margin of error:", round(margin_error, 4), "\n")
cat("95% CI: [", round(t_result$conf.int[1], 3), ",", 
    round(t_result$conf.int[2], 3), "]\n")
cat("t-statistic:", round(t_result$statistic, 3), "\n")
cat("p-value:", format(t_result$p.value, scientific = TRUE, digits = 3), "\n")
```

### Confidence Interval for Difference Between Two Means

```{r ci_two_means}
# Independent samples t-test
t_sleep <- t.test(quality ~ intervention, data = sleep_data, 
                  var.equal = TRUE, conf.level = 0.95)

# Calculate statistics for each group
meditation_mean <- mean(sleep_data$quality[sleep_data$intervention == "Meditation"])
exercise_mean <- mean(sleep_data$quality[sleep_data$intervention == "Exercise"])
meditation_sd <- sd(sleep_data$quality[sleep_data$intervention == "Meditation"])
exercise_sd <- sd(sleep_data$quality[sleep_data$intervention == "Exercise"])
n_meditation <- sum(sleep_data$intervention == "Meditation")
n_exercise <- sum(sleep_data$intervention == "Exercise")

mean_diff <- meditation_mean - exercise_mean

# Display results
cat("Confidence Interval for Difference in Sleep Quality:\n")
cat("Meditation group: n =", n_meditation, ", M =", round(meditation_mean, 2), 
    ", SD =", round(meditation_sd, 2), "\n")
cat("Exercise group: n =", n_exercise, ", M =", round(exercise_mean, 2), 
    ", SD =", round(exercise_sd, 2), "\n")
cat("Mean difference:", round(mean_diff, 3), "\n")
cat("95% CI for difference: [", round(t_sleep$conf.int[1], 3), ",", 
    round(t_sleep$conf.int[2], 3), "]\n")
cat("t-statistic:", round(t_sleep$statistic, 3), "\n")
cat("df:", round(t_sleep$parameter, 1), "\n")
cat("p-value:", round(t_sleep$p.value, 4), "\n")
```

### Confidence Interval for a Proportion

```{r ci_proportion}
# Urban vaccination rate
prop_test_urban <- prop.test(vaccinated_urban, n_urban, conf.level = 0.95)
p_urban <- vaccinated_urban / n_urban
se_urban <- sqrt(p_urban * (1 - p_urban) / n_urban)

# Rural vaccination rate
prop_test_rural <- prop.test(vaccinated_rural, n_rural, conf.level = 0.95)
p_rural <- vaccinated_rural / n_rural
se_rural <- sqrt(p_rural * (1 - p_rural) / n_rural)

# Difference in proportions
prop_diff_test <- prop.test(c(vaccinated_urban, vaccinated_rural), 
                             c(n_urban, n_rural), conf.level = 0.95)

cat("Confidence Intervals for Vaccination Rates:\n\n")

cat("Urban community:\n")
cat("Vaccinated:", vaccinated_urban, "out of", n_urban, "\n")
cat("Proportion:", round(p_urban, 4), "\n")
cat("SE:", round(se_urban, 4), "\n")
cat("95% CI: [", round(prop_test_urban$conf.int[1], 4), ",", 
    round(prop_test_urban$conf.int[2], 4), "]\n\n")

cat("Rural community:\n")
cat("Vaccinated:", vaccinated_rural, "out of", n_rural, "\n")
cat("Proportion:", round(p_rural, 4), "\n")
cat("SE:", round(se_rural, 4), "\n")
cat("95% CI: [", round(prop_test_rural$conf.int[1], 4), ",", 
    round(prop_test_rural$conf.int[2], 4), "]\n\n")

cat("Difference in proportions:\n")
cat("Urban - Rural =", round(p_urban - p_rural, 4), "\n")
cat("95% CI for difference: [", round(prop_diff_test$conf.int[1], 4), ",", 
    round(prop_diff_test$conf.int[2], 4), "]\n")
cat("χ² statistic:", round(prop_diff_test$statistic, 3), "\n")
cat("p-value:", format(prop_diff_test$p.value, scientific = TRUE, digits = 3), "\n")
```

### Confidence Intervals for Regression Coefficients

```{r ci_regression}
# Fit linear regression model
stress_model <- lm(productivity ~ stress, data = stress_data)
model_summary <- summary(stress_model)

# Extract coefficients
coeffs <- coef(stress_model)
beta_0 <- coeffs[1]
beta_1 <- coeffs[2]

# Get confidence intervals
ci_coeffs <- confint(stress_model, level = 0.95)

# Extract standard errors
se_beta_0 <- model_summary$coefficients[1, "Std. Error"]
se_beta_1 <- model_summary$coefficients[2, "Std. Error"]

# Extract t-values and p-values
t_beta_0 <- model_summary$coefficients[1, "t value"]
t_beta_1 <- model_summary$coefficients[2, "t value"]
p_beta_0 <- model_summary$coefficients[1, "Pr(>|t|)"]
p_beta_1 <- model_summary$coefficients[2, "Pr(>|t|)"]

# R-squared
r_squared <- model_summary$r.squared
adj_r_squared <- model_summary$adj.r.squared

cat("Regression Model: Productivity ~ Stress\n\n")

cat("Intercept (β₀):\n")
cat("Estimate:", round(beta_0, 3), "\n")
cat("SE:", round(se_beta_0, 3), "\n")
cat("95% CI: [", round(ci_coeffs[1, 1], 3), ",", round(ci_coeffs[1, 2], 3), "]\n")
cat("t-value:", round(t_beta_0, 3), "\n")
cat("p-value:", format(p_beta_0, scientific = TRUE, digits = 3), "\n\n")

cat("Stress coefficient (β₁):\n")
cat("Estimate:", round(beta_1, 4), "\n")
cat("SE:", round(se_beta_1, 4), "\n")
cat("95% CI: [", round(ci_coeffs[2, 1], 4), ",", round(ci_coeffs[2, 2], 4), "]\n")
cat("t-value:", round(t_beta_1, 3), "\n")
cat("p-value:", format(p_beta_1, scientific = TRUE, digits = 3), "\n\n")

cat("Model fit:\n")
cat("R² =", round(r_squared, 4), "\n")
cat("Adjusted R² =", round(adj_r_squared, 4), "\n")
```

### Confidence Interval for Correlation Coefficient

```{r ci_correlation}
# Pearson correlation test
cor_test <- cor.test(support_data$support, support_data$wellbeing, 
                     conf.level = 0.95)

r_value <- cor_test$estimate
r_lower <- cor_test$conf.int[1]
r_upper <- cor_test$conf.int[2]

cat("Correlation: Social Support and Well-being\n\n")
cat("Pearson's r:", round(r_value, 4), "\n")
cat("95% CI: [", round(r_lower, 4), ",", round(r_upper, 4), "]\n")
cat("t-statistic:", round(cor_test$statistic, 3), "\n")
cat("df:", cor_test$parameter, "\n")
cat("p-value:", format(cor_test$p.value, scientific = TRUE, digits = 3), "\n")
```

## 7. Interpreting Confidence Intervals

Let's interpret each confidence interval in practical context.

### Vegetable Servings (95% CI: [`r round(t_result$conf.int[1], 2)`, `r round(t_result$conf.int[2], 2)`])

We are 95% confident that the true mean increase in daily vegetable servings for the population is between `r round(t_result$conf.int[1], 2)` and `r round(t_result$conf.int[2], 2)` servings. Our point estimate of `r round(mean_servings, 2)` servings represents our best guess, but the interval acknowledges sampling uncertainty.

**Practical interpretation:** Even at the lower bound (`r round(t_result$conf.int[1], 2)` servings), the intervention appears beneficial. The interval does not include zero, consistent with the significant t-test (p `r ifelse(t_result$p.value < 0.001, "< .001", paste("=", round(t_result$p.value, 3)))`).

### Sleep Quality Difference (95% CI: [`r round(t_sleep$conf.int[1], 2)`, `r round(t_sleep$conf.int[2], 2)`])

The difference in sleep quality between meditation and exercise groups is estimated at `r round(mean_diff, 2)` points (Meditation - Exercise), with a 95% CI of [`r round(t_sleep$conf.int[1], 2)`, `r round(t_sleep$conf.int[2], 2)`].

**Practical interpretation:** `r ifelse(t_sleep$conf.int[1] > 0, "Since the entire interval is positive, we can be confident that meditation produces higher sleep quality than exercise. The effect ranges from as small as", ifelse(t_sleep$conf.int[2] < 0, "Since the entire interval is negative, we can be confident that exercise produces higher sleep quality. The effect ranges from as small as", "Since the interval includes zero, we cannot confidently conclude that one intervention is superior. The true difference could range from"))` `r round(abs(t_sleep$conf.int[1]), 2)` points to as large as `r round(abs(t_sleep$conf.int[2]), 2)` points.

### Vaccination Rates (Urban: [`r round(prop_test_urban$conf.int[1], 3)`, `r round(prop_test_urban$conf.int[2], 3)`]; Rural: [`r round(prop_test_rural$conf.int[1], 3)`, `r round(prop_test_rural$conf.int[2], 3)`])

The urban vaccination rate is estimated at `r round(p_urban * 100, 1)`% (95% CI: [`r round(prop_test_urban$conf.int[1] * 100, 1)`%, `r round(prop_test_urban$conf.int[2] * 100, 1)`%]), while the rural rate is `r round(p_rural * 100, 1)`% (95% CI: [`r round(prop_test_rural$conf.int[1] * 100, 1)`%, `r round(prop_test_rural$conf.int[2] * 100, 1)`%]).

**Practical interpretation:** The intervals do not overlap substantially, suggesting a meaningful difference. The difference in proportions (`r round((p_urban - p_rural) * 100, 1)`%, 95% CI: [`r round(prop_diff_test$conf.int[1] * 100, 1)`%, `r round(prop_diff_test$conf.int[2] * 100, 1)`%]) indicates urban communities have higher vaccination rates.

### Stress-Productivity Relationship (β₁: [`r round(ci_coeffs[2, 1], 3)`, `r round(ci_coeffs[2, 2], 3)`])

For each 1-point increase in stress level, productivity decreases by `r round(abs(beta_1), 3)` points (95% CI: [`r round(ci_coeffs[2, 1], 3)`, `r round(ci_coeffs[2, 2], 3)`]).

**Practical interpretation:** The entire confidence interval is negative, indicating that stress reliably predicts lower productivity. The effect could be as small as a `r round(abs(ci_coeffs[2, 2]), 3)`-point decrease or as large as a `r round(abs(ci_coeffs[2, 1]), 3)`-point decrease per unit of stress.

### Social Support-Wellbeing Correlation (r: [`r round(r_lower, 3)`, `r round(r_upper, 3)`])

The correlation between social support and well-being is `r round(r_value, 3)` (95% CI: [`r round(r_lower, 3)`, `r round(r_upper, 3)`]).

**Practical interpretation:** We are 95% confident the true population correlation is between `r round(r_lower, 3)` and `r round(r_upper, 3)`, indicating a `r ifelse(r_lower > 0.5, "strong", ifelse(r_lower > 0.3, "moderate to strong", "moderate"))` positive relationship between social support and well-being.

## 8. Visualizing Confidence Intervals

Let's create visualizations that display confidence intervals clearly.

```{r visualize_ci, fig.width=10, fig.height=8}
# Create summary data for visualization
ci_data <- data.frame(
  Parameter = c("Vegetable Servings\n(Mean)",
                "Sleep Quality Diff\n(Med - Exe)",
                "Urban Vaccination\n(Proportion)",
                "Rural Vaccination\n(Proportion)",
                "Stress Effect\n(β₁)"),
  Estimate = c(mean_servings, mean_diff, p_urban, p_rural, beta_1),
  Lower = c(t_result$conf.int[1], t_sleep$conf.int[1], 
            prop_test_urban$conf.int[1], prop_test_rural$conf.int[1],
            ci_coeffs[2, 1]),
  Upper = c(t_result$conf.int[2], t_sleep$conf.int[2],
            prop_test_urban$conf.int[2], prop_test_rural$conf.int[2],
            ci_coeffs[2, 2]),
  Type = c("Mean", "Difference", "Proportion", "Proportion", "Slope")
)

# Plot 1: Forest plot of CIs
p_forest <- ggplot(ci_data, aes(x = Parameter, y = Estimate, color = Type)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.3, linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  labs(title = "Confidence Intervals Across Different Analyses",
       subtitle = "Points show estimates; error bars show 95% CIs",
       x = "Parameter", y = "Estimate Value") +
  scale_color_manual(values = c("#E74C3C", "#3498DB", "#27AE60", "#F39C12")) +
  coord_flip() +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "right"
  )

# Plot 2: Regression with CI
p_reg_ci <- ggplot(stress_data, aes(x = stress, y = productivity)) +
  geom_point(alpha = 0.5, color = "#3498DB") +
  geom_smooth(method = "lm", se = TRUE, color = "#E74C3C", 
              fill = "#E74C3C", alpha = 0.2) +
  labs(title = "Stress-Productivity Relationship with 95% CI",
       subtitle = paste0("β₁ = ", round(beta_1, 3), 
                        ", 95% CI [", round(ci_coeffs[2, 1], 3), 
                        ", ", round(ci_coeffs[2, 2], 3), "]"),
       x = "Stress Level", y = "Productivity Score") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

# Plot 3: Means with error bars (sleep study)
sleep_summary <- data.frame(
  Intervention = c("Meditation", "Exercise"),
  Mean = c(meditation_mean, exercise_mean),
  SE = c(meditation_sd / sqrt(n_meditation), exercise_sd / sqrt(n_exercise))
)
sleep_summary$CI_Lower <- sleep_summary$Mean - 1.96 * sleep_summary$SE
sleep_summary$CI_Upper <- sleep_summary$Mean + 1.96 * sleep_summary$SE

p_means_ci <- ggplot(sleep_summary, aes(x = Intervention, y = Mean, fill = Intervention)) +
  geom_bar(stat = "identity", alpha = 0.7, width = 0.6) +
  geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), 
                width = 0.2, linewidth = 1) +
  labs(title = "Sleep Quality by Intervention with 95% CI",
       x = "Intervention", y = "Mean Sleep Quality Score") +
  scale_fill_manual(values = c("#27AE60", "#E74C3C")) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "none"
  )

# Plot 4: Proportions with CIs
prop_summary <- data.frame(
  Community = c("Urban", "Rural"),
  Proportion = c(p_urban, p_rural),
  Lower = c(prop_test_urban$conf.int[1], prop_test_rural$conf.int[1]),
  Upper = c(prop_test_urban$conf.int[2], prop_test_rural$conf.int[2])
)

p_prop_ci <- ggplot(prop_summary, aes(x = Community, y = Proportion, fill = Community)) +
  geom_bar(stat = "identity", alpha = 0.7, width = 0.6) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), 
                width = 0.2, linewidth = 1) +
  labs(title = "Vaccination Rates by Community with 95% CI",
       x = "Community Type", y = "Vaccination Rate") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_fill_manual(values = c("#9B59B6", "#F39C12")) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "none"
  )

# Arrange plots
grid.arrange(p_forest, p_reg_ci, p_means_ci, p_prop_ci, ncol = 2)
```

## 9. Factors Affecting Confidence Interval Width

Let's demonstrate how sample size, variability, and confidence level affect CI width.

```{r ci_factors, fig.width=10, fig.height=6}
# Simulate data with different sample sizes
set.seed(2025)
sample_sizes <- c(20, 50, 100, 200, 500)
true_mean <- 50
true_sd <- 15

ci_widths <- data.frame()

for (n in sample_sizes) {
  # Generate sample
  sample_data <- rnorm(n, mean = true_mean, sd = true_sd)
  
  # Calculate CI
  t_res <- t.test(sample_data, conf.level = 0.95)
  
  ci_widths <- rbind(ci_widths, data.frame(
    n = n,
    Mean = mean(sample_data),
    Lower = t_res$conf.int[1],
    Upper = t_res$conf.int[2],
    Width = t_res$conf.int[2] - t_res$conf.int[1]
  ))
}

# Plot: CI width vs sample size
p_width <- ggplot(ci_widths, aes(x = n, y = Width)) +
  geom_line(color = "#3498DB", linewidth = 1.2) +
  geom_point(color = "#3498DB", size = 3) +
  labs(title = "Effect of Sample Size on Confidence Interval Width",
       subtitle = "95% CI for population mean",
       x = "Sample Size (n)", y = "CI Width") +
  scale_x_continuous(breaks = sample_sizes) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

# Effect of confidence level
conf_levels <- c(0.80, 0.90, 0.95, 0.99)
sample_fixed <- rnorm(100, mean = true_mean, sd = true_sd)

ci_conf_levels <- data.frame()

for (conf in conf_levels) {
  t_res <- t.test(sample_fixed, conf.level = conf)
  
  ci_conf_levels <- rbind(ci_conf_levels, data.frame(
    Confidence = paste0(conf * 100, "%"),
    Conf_Numeric = conf,
    Lower = t_res$conf.int[1],
    Upper = t_res$conf.int[2],
    Width = t_res$conf.int[2] - t_res$conf.int[1]
  ))
}

# Plot: CI width vs confidence level
p_conf <- ggplot(ci_conf_levels, aes(x = Conf_Numeric, y = Width)) +
  geom_line(color = "#E74C3C", linewidth = 1.2) +
  geom_point(color = "#E74C3C", size = 3) +
  labs(title = "Effect of Confidence Level on CI Width",
       subtitle = "Fixed sample size (n = 100)",
       x = "Confidence Level", y = "CI Width") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "gray90", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    panel.border = element_blank(),
    axis.line.x.top = element_blank(),
    axis.line.y.right = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

grid.arrange(p_width, p_conf, ncol = 2)

# Display numerical results
cat("Effect of Sample Size (95% CI):\n")
print(ci_widths)
cat("\n\nEffect of Confidence Level (n = 100):\n")
print(ci_conf_levels[, c("Confidence", "Lower", "Upper", "Width")])
```

**Key findings:**

1. **Sample size:** As n increases from `r min(sample_sizes)` to `r max(sample_sizes)`, CI width decreases from `r round(max(ci_widths$Width), 2)` to `r round(min(ci_widths$Width), 2)` (a `r round((1 - min(ci_widths$Width)/max(ci_widths$Width)) * 100, 1)`% reduction).

2. **Confidence level:** Increasing confidence from 80% to 99% widens the CI from `r round(min(ci_conf_levels$Width), 2)` to `r round(max(ci_conf_levels$Width), 2)` (a `r round((max(ci_conf_levels$Width)/min(ci_conf_levels$Width) - 1) * 100, 1)`% increase).

## 10. Reporting Confidence Intervals in APA Style

Here's how to report each analysis with confidence intervals in APA format:

**Single Mean (Nutrition Intervention):**

> The nutrition intervention led to an average increase of `r round(mean_servings, 2)` servings of vegetables per day (*SD* = `r round(sd_servings, 2)`), which was significantly different from zero, *t*(`r df`) = `r round(t_result$statistic, 2)`, *p* `r ifelse(t_result$p.value < 0.001, "< .001", paste("=", round(t_result$p.value, 3)))`, 95% CI [`r round(t_result$conf.int[1], 2)`, `r round(t_result$conf.int[2], 2)`].

**Two-Group Comparison (Sleep Quality):**

> Independent samples t-test revealed that meditation (*M* = `r round(meditation_mean, 2)`, *SD* = `r round(meditation_sd, 2)`) produced `r ifelse(mean_diff > 0, "higher", "lower")` sleep quality scores than exercise (*M* = `r round(exercise_mean, 2)`, *SD* = `r round(exercise_sd, 2)`), *t*(`r round(t_sleep$parameter, 1)`) = `r round(t_sleep$statistic, 2)`, *p* = `r round(t_sleep$p.value, 3)`, 95% CI [`r round(t_sleep$conf.int[1], 2)`, `r round(t_sleep$conf.int[2], 2)`]. The mean difference of `r round(abs(mean_diff), 2)` points suggests a `r ifelse(abs(mean_diff) < 3, "small", ifelse(abs(mean_diff) < 6, "moderate", "large"))` practical effect.

**Proportion (Vaccination Rates):**

> Vaccination rates differed significantly between urban (`r round(p_urban * 100, 1)`%, 95% CI [`r round(prop_test_urban$conf.int[1] * 100, 1)`%, `r round(prop_test_urban$conf.int[2] * 100, 1)`%]) and rural (`r round(p_rural * 100, 1)`%, 95% CI [`r round(prop_test_rural$conf.int[1] * 100, 1)`%, `r round(prop_test_rural$conf.int[2] * 100, 1)`%]) communities, χ²(1) = `r round(prop_diff_test$statistic, 2)`, *p* `r ifelse(prop_diff_test$p.value < 0.001, "< .001", paste("=", round(prop_diff_test$p.value, 3)))`. The difference in rates was `r round((p_urban - p_rural) * 100, 1)` percentage points, 95% CI [`r round(prop_diff_test$conf.int[1] * 100, 1)`%, `r round(prop_diff_test$conf.int[2] * 100, 1)`%].

**Regression Coefficient (Stress-Productivity):**

> Stress level significantly predicted productivity, *b* = `r round(beta_1, 3)`, *SE* = `r round(se_beta_1, 3)`, *t*(`r model_summary$df[2]`) = `r round(t_beta_1, 2)`, *p* `r ifelse(p_beta_1 < 0.001, "< .001", paste("=", round(p_beta_1, 3)))`, 95% CI [`r round(ci_coeffs[2, 1], 3)`, `r round(ci_coeffs[2, 2], 3)`]. For each one-point increase in stress, productivity decreased by `r round(abs(beta_1), 3)` points. The model explained `r round(r_squared * 100, 1)`% of the variance in productivity (*R*² = `r round(r_squared, 3)`).

**Correlation:**

> Social support and well-being were significantly positively correlated, *r*(`r cor_test$parameter`) = `r round(r_value, 3)`, *p* `r ifelse(cor_test$p.value < 0.001, "< .001", paste("=", round(cor_test$p.value, 3)))`, 95% CI [`r round(r_lower, 3)`, `r round(r_upper, 3)`]. This indicates a `r ifelse(abs(r_value) < 0.3, "weak", ifelse(abs(r_value) < 0.5, "moderate", "strong"))` positive relationship.

## 11. Practical Considerations

### Choosing the Confidence Level

**95% is conventional** but not mandatory:

- **90% CI:** Wider tolerance for error, narrower intervals (exploratory research)
- **95% CI:** Standard in most fields (balances precision and confidence)
- **99% CI:** More conservative, wider intervals (high-stakes decisions, clinical trials)

**Trade-off:** Higher confidence → wider intervals → less precision

### Sample Size Planning

To achieve a target CI width:

$$
n = \left(\frac{2 \cdot z_{\alpha/2} \cdot \sigma}{W}\right)^2
$$

where W is the desired total width.

**Example:** To estimate a mean within ±2 units with 95% confidence when σ = 10:

```{r sample_size_ci}
# Parameters
z_value <- qnorm(0.975)  # 95% confidence
sigma <- 10
desired_width <- 4  # Total width (±2 on each side)

# Calculate required n
n_required <- ceiling((2 * z_value * sigma / desired_width)^2)

cat("Sample size needed for CI width of", desired_width, "units:\n")
cat("n =", n_required, "\n")
```

You would need n = `r n_required` participants.

### Confidence Intervals vs. P-values

**Confidence intervals are superior** to p-values alone because they:

1. **Provide magnitude information:** Not just "significant" but "how large?"
2. **Show precision:** Narrow CIs indicate precise estimates
3. **Aid interpretation:** Can assess practical significance
4. **Remain informative when non-significant:** A wide CI suggests low power, not necessarily no effect

**Example:** A non-significant result with CI [-0.5, 8.2] suggests uncertainty, whereas [-0.1, 0.3] suggests a genuinely small effect.

### Common Misinterpretations

**WRONG:** "There's a 95% probability the true parameter is in this interval."

**RIGHT:** "95% of intervals constructed this way will contain the true parameter."

**WRONG:** "Values outside the CI are impossible."

**RIGHT:** "Values outside the CI are less plausible given the data."

**WRONG:** "Overlapping CIs mean no significant difference."

**RIGHT:** "Use formal tests for comparisons; CI overlap is not equivalent to hypothesis testing."

### Reporting Best Practices

1. **Always report CIs** with point estimates
2. **Specify the confidence level** (usually 95%)
3. **Include units** and context
4. **Report actual limits**, not just "significant"
5. **Interpret practically**, not just statistically
6. **Use consistent precision** (2-3 decimal places typical)

## 12. Advanced Topics

### Bootstrap Confidence Intervals

For non-normal distributions or complex statistics:

```{r bootstrap_ci}
# Bootstrap CI for median
set.seed(2025)
n_boot <- 10000
boot_medians <- numeric(n_boot)

for (i in 1:n_boot) {
  boot_sample <- sample(nutrition_data$servings, replace = TRUE)
  boot_medians[i] <- median(boot_sample)
}

# Calculate percentile CI
boot_ci_lower <- quantile(boot_medians, 0.025)
boot_ci_upper <- quantile(boot_medians, 0.975)
observed_median <- median(nutrition_data$servings)

cat("Bootstrap 95% CI for Median Vegetable Servings:\n")
cat("Observed median:", round(observed_median, 3), "\n")
cat("Bootstrap CI: [", round(boot_ci_lower, 3), ",", round(boot_ci_upper, 3), "]\n")
cat("Based on", n_boot, "bootstrap samples\n")
```

### Bonferroni Correction for Multiple Comparisons

When constructing multiple CIs, adjust confidence level:

```{r bonferroni}
# For 5 comparisons, maintain family-wise error rate of 0.05
n_comparisons <- 5
alpha_adjusted <- 0.05 / n_comparisons
conf_level_adjusted <- 1 - alpha_adjusted

cat("Multiple Comparisons Adjustment:\n")
cat("Number of comparisons:", n_comparisons, "\n")
cat("Unadjusted α:", 0.05, "\n")
cat("Bonferroni-adjusted α:", round(alpha_adjusted, 4), "\n")
cat("Adjusted confidence level:", round(conf_level_adjusted * 100, 2), "%\n")

# Example: CI for vegetable servings with Bonferroni correction
t_bonf <- t.test(nutrition_data$servings, conf.level = conf_level_adjusted)

cat("\nOriginal 95% CI: [", round(t_result$conf.int[1], 3), ",", 
    round(t_result$conf.int[2], 3), "]\n")
cat("Bonferroni", round(conf_level_adjusted * 100, 2), "% CI: [", 
    round(t_bonf$conf.int[1], 3), ",", round(t_bonf$conf.int[2], 3), "]\n")
cat("Width increase:", round((t_bonf$conf.int[2] - t_bonf$conf.int[1]) - 
                              (t_result$conf.int[2] - t_result$conf.int[1]), 3), "\n")
```

## 13. Conclusion

Confidence intervals are indispensable tools for quantifying uncertainty in statistical estimates. They transform point estimates into ranges that communicate both our best guess and the precision of that estimate. Unlike p-values, which only indicate whether an effect exists, confidence intervals reveal the magnitude and plausibility of different effect sizes.

**Key takeaways:**

1. **CIs quantify uncertainty:** They provide a range of plausible values for population parameters
2. **Interpretation requires care:** The confidence level refers to the procedure, not individual intervals
3. **Width matters:** Narrower intervals indicate more precise estimates
4. **Always report CIs:** They provide essential context beyond p-values
5. **Consider practical significance:** Statistical significance (CI excludes zero) doesn't guarantee meaningful impact

In our scenarios, confidence intervals revealed:

- Nutrition intervention increased vegetable servings by `r round(mean_servings, 2)` per day (95% CI: [`r round(t_result$conf.int[1], 2)`, `r round(t_result$conf.int[2], 2)`])
- Meditation improved sleep quality by `r round(mean_diff, 2)` points over exercise (95% CI: [`r round(t_sleep$conf.int[1], 2)`, `r round(t_sleep$conf.int[2], 2)`])
- Urban vaccination rate was `r round((p_urban - p_rural) * 100, 1)`% higher than rural (95% CI: [`r round(prop_diff_test$conf.int[1] * 100, 1)`%, `r round(prop_diff_test$conf.int[2] * 100, 1)`%])
- Each stress point reduced productivity by `r round(abs(beta_1), 3)` points (95% CI: [`r round(ci_coeffs[2, 1], 3)`, `r round(ci_coeffs[2, 2], 3)`])
- Social support and well-being correlated at r = `r round(r_value, 3)` (95% CI: [`r round(r_lower, 3)`, `r round(r_upper, 3)`])

By embracing confidence intervals, we move beyond simplistic "significant/not significant" dichotomies to nuanced, informative quantification of uncertainty—the hallmark of rigorous, transparent scientific communication.
