Table of Contents

1. Introduction to Statistical Inference
   • 1.1 What is Statistical Inference?
   • 1.2 Population vs Sample
   • 1.3 Data Science Workflow
   • 1.4 Traditional R vs Tidyverse Approach
2. Probability Review
   • 2.1 Random Variables
   • 2.2 Probability Distributions
   • 2.3 Expected Value and Variance
3. Discrete Probability Distributions
   • 3.1 Binomial Distribution
   • 3.2 Poisson Distribution
4. Continuous Probability Distributions
   • 4.1 Normal Distribution
   • 4.2 Standard Normal Distribution

1. Introduction to Statistical Inference

1.1 What is Statistical Inference?

Statistical inference is the process of drawing conclusions about a population based on information obtained from a sample. It allows us to make probabilistic statements about population parameters using sample statistics.

The two main branches of statistical inference are:

• Estimation: Using sample data to estimate population parameters
• Hypothesis Testing: Using sample data to test claims about population parameters

1.2 Population vs Sample

Population

A population is the complete set of all items of interest in a statistical study. Population parameters are typically denoted by Greek letters:

• \(\mu\) (mu): Population mean
• \(\sigma\) (sigma): Population standard deviation
• \(p\): Population proportion

Sample

A sample is a subset of the population that is actually observed. Sample statistics are typically denoted by Roman letters:

• \(\bar{x}\) (x-bar): Sample mean
• \(s\): Sample standard deviation
• \(\hat{p}\) (p-hat): Sample proportion

Why Sample?

• Cost efficiency
• Time constraints
• Practical impossibility of measuring entire population
• Destructive testing scenarios

1.3 Data Science Workflow

The modern data science workflow for inferential statistics follows these steps:

1. Import: Bring data into R
2. Tidy: Organize data into tidy format
3. Transform: Manipulate data for analysis
4. Visualize: Create plots to understand data
5. Model: Apply statistical models
6. Communicate: Share results and insights

1.4 Traditional R vs Tidyverse Approach

Traditional R

# Traditional approach
asia <- gapminder[gapminder$$continent == "Asia", ]
mean(asia$$lifeExp)

Tidyverse Approach

# Tidyverse approach
library(dplyr)
gapminder %>% 
  filter(continent == "Asia") %>% 
  summarize(mean_exp = mean(lifeExp))

The tidyverse offers:

• More readable code with the pipe operator (%>%)
• Consistent grammar of data manipulation
• Integration with modern visualization tools

2. Probability Review

2.1 Random Variables

A random variable is a real-valued function defined on a sample space. For each outcome \(s\) in the sample space \(S\), a random variable \(X\) assigns a real number \(X(s)\).

Types of Random Variables

1. Discrete Random Variable: Can take a countable number of values
   • Example: Number of heads in 10 coin flips
   • Example: Number of defective items in a batch
2. Continuous Random Variable: Can take any value within an interval
   • Example: Height of students
   • Example: Time to complete an exam

2.2 Probability Distributions

Discrete Probability Distribution

For a discrete random variable \(X\), the probability distribution is defined by the probability mass function (PMF):

Properties:

1. \(p(x) \geq 0\) for all values of \(x\)
2. \[\sum p(x) = 1\]

Continuous Probability Distribution

For a continuous random variable \(X\), the probability distribution is defined by the probability density function (PDF) \(f(x)\).

For an interval \((a, b]\): \(P(a < X \leq b) = \int_a^b f(x)dx\)

Properties:

1. \(f(x) \geq 0\) for all \(x\)
2. \[\int_{-\infty}^{\infty} f(x)dx = 1\]

2.3 Expected Value and Variance

Expected Value (Mean)

The expected value represents the long-run average value of a random variable.

For discrete random variables: \(E(X) = \mu = \sum x \cdot p(x)\)

For continuous random variables: \(E(X) = \mu = \int_{-\infty}^{\infty} x \cdot f(x)dx\)

Variance

Variance measures the spread of a distribution around its mean.

For discrete random variables: \(Var(X) = \sigma^2 = E[(X - \mu)^2] = \sum (x - \mu)^2 \cdot p(x)\)

For continuous random variables: \(Var(X) = \sigma^2 = E[(X - \mu)^2] = \int_{-\infty}^{\infty} (x - \mu)^2 \cdot f(x)dx\)

Standard Deviation

The standard deviation is the square root of the variance: \(\sigma = \sqrt{Var(X)}\)

3. Discrete Probability Distributions

3.1 Binomial Distribution

The binomial distribution models the number of successes in a fixed number of independent trials.

Conditions for Binomial Distribution

1. Fixed number of trials (\(n\))
2. Two possible outcomes: Success or Failure
3. Constant probability of success (\(p\)) for each trial
4. Independent trials

Probability Mass Function

where:

• \(x\) = number of successes
• \(n\) = number of trials
• \(p\) = probability of success
• \(\binom{n}{x} = \frac{n!}{x!(n-x)!}\) (binomial coefficient)

Mean and Variance

• Mean: \(\mu = np\)
• Variance: \(\sigma^2 = np(1-p)\)

Example

If 40% of a class is female, what is the probability that exactly 6 of the first 10 students walking in will be female?

3.2 Poisson Distribution

The Poisson distribution models the number of events occurring in a fixed interval of time or space.

Conditions for Poisson Distribution

1. Events occur independently
2. Average rate of occurrence (\(\lambda\)) is constant
3. Two events cannot occur at exactly the same instant

Probability Mass Function

where:

• \(x\) = number of events
• \(\lambda\) = average rate of occurrence
• \(e\) = Euler’s number (≈ 2.71828)

Mean and Variance

• Mean: \(\mu = \lambda\)
• Variance: \(\sigma^2 = \lambda\)

Example

If there are an average of 3 striped trout per 100 yards in a stream, what is the probability of seeing exactly 5 striped trout in the next 100 yards?

4. Continuous Probability Distributions

4.1 Normal Distribution

The normal distribution is the most important continuous distribution in statistics, often called the “bell curve.”

Probability Density Function

where:

• \(\mu\) = mean
• \(\sigma\) = standard deviation
• \(\pi\) ≈ 3.14159
• \(e\) ≈ 2.71828

Properties

1. Symmetric about the mean
2. Mean = Median = Mode
3. Total area under the curve = 1
4. Approximately 68% of data falls within ±1σ of μ
5. Approximately 95% of data falls within ±2σ of μ
6. Approximately 99.7% of data falls within ±3σ of μ

Notation

\(X \sim N(\mu, \sigma^2)\) means “\(X\) follows a normal distribution with mean \(\mu\) and variance \(\sigma^2\)”

4.2 Standard Normal Distribution

The standard normal distribution is a special case of the normal distribution with \(\mu = 0\) and \(\sigma = 1\).

Z-Score Transformation

Any normal random variable \(X \sim N(\mu, \sigma^2)\) can be standardized: \(Z = \frac{X - \mu}{\sigma}\)

where \(Z \sim N(0, 1)\)

Using Z-Tables

To find probabilities for normal distributions:

1. Convert to z-score
2. Use standard normal table or R functions

Example

If gestational length \(X \sim N(39, 2^2)\) weeks, what percentage of gestations are less than 40 weeks?

Using the standard normal table: \(P(Z \leq 0.5) = 0.6915\)

Therefore, approximately 69.15% of gestations are less than 40 weeks.

R Implementation

# Using pnorm for normal probabilities
pnorm(40, mean = 39, sd = 2)  # Direct calculation
pnorm(0.5)                    # Using z-score

Part 2: Sampling and Point Estimation

5. Sampling and Sampling Distributions

5.1 Introduction to Sampling

Sampling is the process of selecting a subset of individuals from a population to estimate characteristics of the whole population.

Why Sample?

• Census limitations: Measuring entire populations is often impractical
• Cost efficiency: Samples are less expensive to collect
• Time constraints: Samples provide quicker results
• Destructive testing: Some measurements destroy the item being tested

Types of Sampling

1. Random Sampling: Every member of the population has an equal chance of selection
2. Stratified Sampling: Population divided into strata, then random samples from each
3. Cluster Sampling: Population divided into clusters, some clusters randomly selected
4. Systematic Sampling: Select every kth element from a list

5.2 Sampling With and Without Replacement

Sampling With Replacement

• Each selected item is returned to the population before the next draw
• Same item can be selected multiple times
• Observations are independent

Sampling Without Replacement

• Selected items are not returned to the population
• Same item cannot be selected twice
• Observations are not independent

R Implementation

# Sampling with replacement
sample(1:10, size = 5, replace = TRUE)

# Sampling without replacement
sample(1:10, size = 5, replace = FALSE)

5.3 Sampling Distributions

A sampling distribution is the probability distribution of a statistic obtained from all possible samples of a specific size from a population.

Key Concepts

1. Sample Statistic: A value calculated from sample data (e.g., \(\bar{x}\), \(\hat{p}\))
2. Population Parameter: The true value in the population (e.g., \(\mu\), \(p\))
3. Sampling Variability: Different samples yield different statistics

The Bowl Example

The classic bowl activity demonstrates sampling concepts:

• Bowl contains red and white balls (unknown proportion)
• Take samples of 50 balls using a “shovel”
• Calculate proportion of red balls in each sample
• Observe variation in sample proportions

library(moderndive)
library(tidyverse)

# Virtual bowl sampling
bowl %>% 
  rep_sample_n(size = 50, reps = 1000) %>%
  group_by(replicate) %>%
  summarize(prop_red = mean(color == "red"))

5.4 Central Limit Theorem

The Central Limit Theorem (CLT) is fundamental to statistical inference:

Statement of CLT

For a random sample of size \(n\) from a population with mean \(\mu\) and standard deviation \(\sigma\):

1. The sampling distribution of \(\bar{X}\) has mean \(\mu_{\bar{X}} = \mu\)
2. The sampling distribution of \(\bar{X}\) has standard deviation \(\sigma_{\bar{X}} = \frac{\sigma}{\sqrt{n}}\)
3. As \(n\) increases, the sampling distribution of \(\bar{X}\) approaches a normal distribution

Mathematical Expression

\(\bar{X} \sim N\left(\mu, \frac{\sigma^2}{n}\right)\) for large \(n\)

Implications

1. Sample means are normally distributed for large samples
2. Larger samples produce more precise estimates
3. We can make probability statements about sample means

5.5 Standard Error

The standard error (SE) measures the variability of a sample statistic.

Standard Error of the Mean

If \(\sigma\) is unknown, we estimate it with \(s\): \(SE(\bar{X}) = \frac{s}{\sqrt{n}}\)

Standard Error of a Proportion

If \(p\) is unknown, we estimate it with \(\hat{p}\): \(SE(\hat{p}) = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\)

Effect of Sample Size

As \(n\) increases:

• Standard error decreases
• Estimates become more precise
• Sampling distribution becomes narrower

6. Point Estimation

6.1 Introduction to Point Estimation

A point estimate is a single value used to estimate a population parameter.

Common Point Estimates

• Sample mean (\(\bar{x}\)) estimates population mean (\(\mu\))
• Sample proportion (\(\hat{p}\)) estimates population proportion (\(p\))
• Sample variance (\(s^2\)) estimates population variance (\(\sigma^2\))

Properties of Good Estimators

1. Unbiased: \(E(\hat{\theta}) = \theta\)
2. Consistent: \(\hat{\theta}\) converges to \(\theta\) as \(n \to \infty\)
3. Efficient: Minimum variance among unbiased estimators

6.2 Method of Moments

The method of moments equates sample moments to population moments to derive estimators.

First Moment (Mean)

Population: \(E(X) = \mu\) Sample: \(\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i\)

Second Moment (Variance)

Population: \(E(X^2) = \sigma^2 + \mu^2\) Sample: \(\frac{1}{n}\sum_{i=1}^n X_i^2\)

Example: Poisson Distribution

For \(X \sim \text{Poisson}(\lambda)\):

• \[E(X) = \lambda\]
• Method of moments estimator: \(\hat{\lambda} = \bar{X}\)

6.3 Maximum Likelihood Estimation (MLE)

Maximum likelihood estimation finds parameter values that maximize the probability of observing the sample data.

Likelihood Function

For independent observations \(x_1, x_2, ..., x_n\): \(L(\theta) = \prod_{i=1}^n f(x_i\|\theta)\)

Log-Likelihood Function

Finding MLE

1. Write the likelihood function
2. Take the logarithm
3. Differentiate with respect to \(\theta\)
4. Set derivative equal to zero
5. Solve for \(\theta\)

Example: Normal Distribution

For \(X_1, ..., X_n \sim N(\mu, \sigma^2)\):

Log-likelihood: \(\ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi) - \frac{n}{2}\log(\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n(x_i - \mu)^2\)

MLEs:

• \[\hat{\mu}_{MLE} = \bar{x}\]
• \[\hat{\sigma}^2_{MLE} = \frac{1}{n}\sum_{i=1}^n(x_i - \bar{x})^2\]

6.4 Examples of Point Estimation

Example 1: Binomial Distribution

For \(Y \sim \text{Binomial}(n, p)\):

• MLE: \(\hat{p} = \frac{Y}{n}\)

Example 2: Exponential Distribution

For \(X_1, ..., X_n \sim \text{Exponential}(\lambda)\):

• MLE: \(\hat{\lambda} = \frac{1}{\bar{x}}\)

Example 3: Poisson Distribution

For \(X_1, ..., X_n \sim \text{Poisson}(\mu)\):

• MLE: \(\hat{\mu} = \bar{x}\)

R Implementation

# Sample data
set.seed(123)
data <- rpois(100, lambda = 5)

# Method of moments estimate
mom_estimate <- mean(data)

# Maximum likelihood estimate (same for Poisson)
mle_estimate <- mean(data)

print(paste("Estimate:", mle_estimate))

6.5 Sampling Distribution of Estimators

The sampling distribution of an estimator shows how the estimator varies across all possible samples.

For Sample Mean

If \(X_1, ..., X_n \sim N(\mu, \sigma^2)\): \(\bar{X} \sim N\left(\mu, \frac{\sigma^2}{n}\right)\)

For Sample Proportion

If \(n\) is large: \(\hat{p} \sim N\left(p, \frac{p(1-p)}{n}\right)\)

Bootstrap Distribution

When theoretical distributions are unknown, we can use bootstrap methods:

library(infer)

# Generate bootstrap distribution
boot_dist <- data %>%
  specify(response = variable) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean")

# Visualize
boot_dist %>% visualize()

6.6 Properties of Estimators

Bias

The bias of an estimator \(\hat{\theta}\) is: \(\text{Bias}(\hat{\theta}) = E(\hat{\theta}) - \theta\)

An estimator is unbiased if \(\text{Bias}(\hat{\theta}) = 0\).

Mean Squared Error (MSE)

Consistency

An estimator \(\hat{\theta}*n\) is consistent if: \(\lim*{n \to \infty} P(\|\hat{\theta}_n - \theta\| > \epsilon) = 0\) for any \(\epsilon > 0\)

Efficiency

Among unbiased estimators, the one with minimum variance is called efficient.

The relative efficiency of two estimators is: \(\text{RE}(\hat{\theta}_1, \hat{\theta}_2) = \frac{\text{Var}(\hat{\theta}_2)}{\text{Var}(\hat{\theta}_1)}\)

Part 3: Confidence Intervals - Theory

7. Introduction to Confidence Intervals

7.1 What is a Confidence Interval?

A confidence interval (CI) is a range of values that likely contains the true population parameter with a specified level of confidence.

Key Components

1. Point estimate: Best single estimate of the parameter
2. Margin of error: How far the interval extends from the point estimate
3. Confidence level: Probability that the method produces an interval containing the true parameter

General Form

Interpretation

A 95% confidence interval means that if we repeated the sampling process many times and calculated a CI each time, approximately 95% of these intervals would contain the true population parameter.

7.2 The Traditional (Standard Error) Method

The traditional method for constructing confidence intervals relies on:

1. The sampling distribution of the statistic
2. The standard error
3. A critical value from the appropriate distribution

General Formula

where:

• \(\hat{\theta}\) = point estimate
• \(z_{\alpha/2}\) = critical value
• \(SE(\hat{\theta})\) = standard error
• \(\alpha\) = significance level (1 - confidence level)

Common Critical Values

| Confidence Level | \(\alpha\) | \(z_{\alpha/2}\) | | —————- | ——– | ————– | | 90% | 0.10 | 1.645 | | 95% | 0.05 | 1.960 | | 99% | 0.01 | 2.576 |

7.3 Confidence Interval for Population Mean

Case 1: Known Population Standard Deviation (\(\sigma\))

When \(\sigma\) is known and either:

• The population is normally distributed, or
• The sample size is large (\(n \geq 30\))

Case 2: Unknown Population Standard Deviation

When \(\sigma\) is unknown, we use the \(t\)-distribution:

where:

• \(t_{\alpha/2, n-1}\) = critical value from \(t\)-distribution with \(n-1\) degrees of freedom
• \(s\) = sample standard deviation

Conditions for Using the t-Distribution

1. Random sample
2. Population is approximately normal, or \(n \geq 30\)
3. Independence of observations

Example

Customer support response times: \(\bar{x} = 57\) minutes, \(s = 15\) minutes, \(n = 200\)

For a 95% CI: \(SE = \frac{15}{\sqrt{200}} = 1.06\) \(\text{Margin of Error} = 1.96 \times 1.06 = 2.08\) \(\text{CI} = 57 \pm 2.08 = [54.92, 59.08]\)

We are 95% confident that the true mean response time is between 54.92 and 59.08 minutes.

R Implementation

# Using t.test function
data <- c(12, 16, 14, 18, 20, 10, 15, 14, 16, 13)
t.test(data, conf.level = 0.95)$$conf.int

# Manual calculation
x_bar <- mean(data)
s <- sd(data)
n <- length(data)
se <- s / sqrt(n)
t_critical <- qt(0.975, df = n - 1)
margin_error <- t_critical * se
ci_lower <- x_bar - margin_error
ci_upper <- x_bar + margin_error

7.4 Confidence Interval for Population Proportion

For a population proportion \(p\), we use the sample proportion \(\hat{p}\):

Conditions

1. Random sample
2. Independence
3. Success-failure condition: \(n\hat{p} \geq 10\) and \(n(1-\hat{p}) \geq 10\)

Example

Customer satisfaction survey: 200 out of 250 customers are satisfied

\(\hat{p} = \frac{200}{250} = 0.80\) \(SE = \sqrt{\frac{0.80 \times 0.20}{250}} = 0.0253\) \(\text{Margin of Error} = 1.96 \times 0.0253 = 0.0496\) \(\text{CI} = 0.80 \pm 0.0496 = [0.750, 0.850]\)

We are 95% confident that between 75% and 85% of all customers are satisfied.

R Implementation

# Using prop.test
prop.test(x = 200, n = 250, conf.level = 0.95)$$conf.int

# Manual calculation
p_hat <- 200/250
se <- sqrt(p_hat * (1 - p_hat) / 250)
margin_error <- 1.96 * se
ci_lower <- p_hat - margin_error
ci_upper <- p_hat + margin_error

7.5 Margin of Error and Width of Confidence Intervals

Margin of Error Components

Factors Affecting CI Width

1. Sample size: Larger \(n\) → narrower CI
2. Confidence level: Higher confidence → wider CI
3. Population variability: Larger \(\sigma\) → wider CI

Relationship with Sample Size

For means: \(\text{Width} \propto \frac{1}{\sqrt{n}}\)

Quadrupling the sample size cuts the CI width in half.

7.6 Sample Size Determination

For Population Mean

To achieve a margin of error \(E\) with confidence level \((1-\alpha)\):

For Population Proportion

If \(p\) is unknown, use \(p = 0.5\) for the most conservative estimate.

Example

To estimate proportion of students who prefer games within 4% margin of error with 95% confidence:

7.7 Interpreting Confidence Intervals

Correct Interpretations

1. “We are 95% confident that the true population mean lies between [lower, upper]”
2. “If we repeated this procedure many times, 95% of the resulting intervals would contain the true parameter”
3. “The method used produces intervals that contain the true parameter 95% of the time”

Incorrect Interpretations

1. ❌ “There is a 95% probability that the true mean is in this interval”
2. ❌ “95% of the data falls within this interval”
3. ❌ “95% of sample means fall within this interval”

Key Points

• The parameter is fixed; the interval is random
• Confidence is in the method, not in a specific interval
• Once calculated, the interval either contains the parameter or it doesn’t

7.8 Finite Population Correction

When sampling from a finite population without replacement, and the sample size is more than 5% of the population:

The corrected standard error becomes: \(SE_{corrected} = SE \times \text{FPC}\)

Example

Population of 1,200 students, sample of 400: \(\text{FPC} = \sqrt{\frac{1200-400}{1200-1}} = 0.816\)

The standard error is reduced by approximately 18.4%.

7.9 Confidence Intervals in R

Using Base R

# For means
t.test(data, conf.level = 0.95)$$conf.int

# For proportions
prop.test(x = successes, n = total, conf.level = 0.95)$$conf.int

Using tidyverse and infer

library(infer)

# For means
data %>%
  specify(response = variable) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean") %>%
  get_confidence_interval(level = 0.95, type = "percentile")

# For proportions
data %>%
  specify(response = variable, success = "yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_confidence_interval(level = 0.95, type = "percentile")

7.10 Visual Representation

Confidence intervals can be visualized using:

library(ggplot2)

# Create data frame for plotting
ci_data <- data.frame(
  estimate = c(x_bar),
  lower = c(ci_lower),
  upper = c(ci_upper)
)

# Plot with error bars
ggplot(ci_data, aes(x = "Sample", y = estimate)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  labs(y = "Estimate", title = "95% Confidence Interval")

Part 4: Confidence Intervals - Advanced & Bootstrap Methods

8. Two-Sample Confidence Intervals

8.1 Confidence Interval for Difference in Means

Independent Samples

When comparing means from two independent populations:

where: \(SE(\bar{x}_1 - \bar{x}_2) = \sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}\)

Degrees of Freedom (Welch’s approximation)

Example

Comparing incomes between Cleveland and Sacramento:

• Cleveland: \(\bar{x}_1 = 45,000\), \(s_1 = 15,000\), \(n_1 = 212\)
• Sacramento: \(\bar{x}_2 = 48,000\), \(s_2 = 18,000\), \(n_2 = 175\)

# Using t.test
t.test(cleveland_income, sacramento_income, conf.level = 0.95)$$conf.int

# Manual calculation
x1_bar <- 45000; s1 <- 15000; n1 <- 212
x2_bar <- 48000; s2 <- 18000; n2 <- 175
se_diff <- sqrt(s1^2/n1 + s2^2/n2)
df <- ((s1^2/n1 + s2^2/n2)^2) / ((s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1))
t_critical <- qt(0.975, df)
margin_error <- t_critical * se_diff
ci_lower <- (x1_bar - x2_bar) - margin_error
ci_upper <- (x1_bar - x2_bar) + margin_error

Paired Samples

For paired data (e.g., before-after measurements):

1. Calculate differences: \(d_i = x_{1i} - x_{2i}\)
2. Find mean difference: \(\bar{d}\)
3. Find standard deviation of differences: \(s_d\)

Example: Weight Loss Program

Before-after weights for 10 participants:

before <- c(65, 60, 63, 58, 70, 72, 80, 81, 83, 69)
after <- c(55, 52, 60, 55, 60, 66, 70, 66, 63, 60)
differences <- after - before

# Using t.test for paired data
t.test(after, before, paired = TRUE, conf.level = 0.95)$$conf.int

# Manual calculation
d_bar <- mean(differences)
s_d <- sd(differences)
n <- length(differences)
se_d <- s_d / sqrt(n)
t_critical <- qt(0.975, df = n - 1)
margin_error <- t_critical * se_d
ci_lower <- d_bar - margin_error
ci_upper <- d_bar + margin_error

8.2 Confidence Interval for Difference in Proportions

For two independent proportions:

where: \(SE(\hat{p}_1 - \hat{p}_2) = \sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1} + \frac{\hat{p}_2(1-\hat{p}_2)}{n_2}}\)

Example

Comparing satisfaction rates between two products:

• Product A: 180 satisfied out of 200 (\(\hat{p}_1 = 0.90\))
• Product B: 150 satisfied out of 200 (\(\hat{p}_2 = 0.75\))

# Using prop.test
prop.test(x = c(180, 150), n = c(200, 200), conf.level = 0.95)$$conf.int

# Manual calculation
p1_hat <- 180/200; n1 <- 200
p2_hat <- 150/200; n2 <- 200
se_diff <- sqrt(p1_hat*(1-p1_hat)/n1 + p2_hat*(1-p2_hat)/n2)
z_critical <- qnorm(0.975)
margin_error <- z_critical * se_diff
ci_lower <- (p1_hat - p2_hat) - margin_error
ci_upper <- (p1_hat - p2_hat) + margin_error

9. Bootstrap Methods for Confidence Intervals

9.1 Introduction to Bootstrap

The bootstrap is a resampling method that estimates the sampling distribution by repeatedly sampling with replacement from the original data.

Key Concepts

1. Resampling: Drawing samples from the original data
2. With replacement: Same observation can be selected multiple times
3. Bootstrap distribution: Distribution of statistics from resampled data

Advantages

• No distributional assumptions required
• Works for complex statistics
• Provides empirical sampling distribution

9.2 Bootstrap Procedure

1. Start with original sample of size \(n\)
2. Draw a bootstrap sample of size \(n\) with replacement
3. Calculate statistic on bootstrap sample
4. Repeat steps 2-3 many times (typically 1000+)
5. Use bootstrap distribution to construct CI

R Implementation

library(infer)

# Bootstrap distribution for mean
boot_dist <- data %>%
  specify(response = variable) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean")

# Visualize bootstrap distribution
boot_dist %>% visualize()

9.3 Bootstrap Confidence Intervals

Percentile Method

Uses percentiles of the bootstrap distribution:

For a 95% CI, use the 2.5th and 97.5th percentiles:

# Percentile CI
percentile_ci <- boot_dist %>%
  get_confidence_interval(level = 0.95, type = "percentile")

Standard Error Method

Uses bootstrap standard error with normal approximation:

where \(SE_{boot}\) is the standard deviation of the bootstrap distribution.

# Standard error CI
se_ci <- boot_dist %>%
  get_confidence_interval(level = 0.95, type = "se", 
                         point_estimate = mean(original_data))

9.4 Bootstrap for Different Statistics

For Proportions

boot_prop <- data %>%
  specify(response = variable, success = "yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop")

For Difference in Means

boot_diff_means <- data %>%
  specify(response ~ group) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "diff in means", order = c("A", "B"))

For Correlation

boot_cor <- data %>%
  specify(x ~ y) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "correlation")

9.5 Visualizing Bootstrap Results

# Bootstrap distribution with CI
boot_dist %>%
  visualize() +
  shade_confidence_interval(endpoints = percentile_ci)

# Comparison of methods
library(gridExtra)
p1 <- boot_dist %>%
  visualize() +
  shade_ci(endpoints = percentile_ci, color = "blue", fill = "lightblue") +
  ggtitle("Percentile Method")

p2 <- boot_dist %>%
  visualize() +
  shade_ci(endpoints = se_ci, color = "red", fill = "pink") +
  ggtitle("Standard Error Method")

grid.arrange(p1, p2, ncol = 2)

10. Sample Size Determination

10.1 Sample Size for Estimating a Mean

To achieve margin of error \(E\) with confidence level \((1-\alpha)\):

When \(\sigma\) is Unknown

1. Use estimate from pilot study
2. Use range/4 as rough estimate
3. Use conservative estimate

Example

To estimate mean income within \(500 with 95% confidence, assuming\)\sigma = \(4,500\):

10.2 Sample Size for Estimating a Proportion

Conservative Approach

When \(p\) is unknown, use \(p = 0.5\):

Example

To estimate proportion within 3% with 95% confidence:

10.3 Sample Size for Comparing Two Means

where:

• \(\beta\) = probability of Type II error
• \(\mu_1 - \mu_2\) = minimum detectable difference

10.4 Effect of Confidence Level and Margin of Error

Relationship with sample size:

• Halving margin of error → quadruples sample size
• Increasing confidence level → increases sample size

Trade-offs:

• Precision vs. cost
• Confidence vs. feasibility

10.5 Finite Population Correction

For finite populations:

where \(N\) is population size.

10.6 R Implementation for Sample Size

# For proportions
library(pwr)
pwr.p.test(h = ES.h(p1 = 0.5, p2 = 0.47), 
           sig.level = 0.05, 
           power = 0.80)

# Custom function for mean
sample_size_mean <- function(sigma, margin_error, conf_level = 0.95) {
  z <- qnorm((1 + conf_level) / 2)
  n <- ceiling((z * sigma / margin_error)^2)
  return(n)
}

# Custom function for proportion
sample_size_prop <- function(p, margin_error, conf_level = 0.95) {
  z <- qnorm((1 + conf_level) / 2)
  if (is.null(p)) p <- 0.5  # Conservative estimate
  n <- ceiling((z^2 * p * (1-p)) / margin_error^2)
  return(n)
}

10.7 Practical Considerations

1. Pilot studies: Small preliminary studies to estimate parameters
2. Resource constraints: Budget, time, personnel
3. Non-response: Plan for 10-30% non-response rate
4. Stratification: May require different sample sizes per stratum
5. Cluster sampling: Requires larger samples due to design effect

Part 5: Hypothesis Testing Framework

11. Introduction to Hypothesis Testing

11.1 What is Hypothesis Testing?

Hypothesis testing is a formal procedure for using sample data to evaluate claims about population parameters. It provides a framework for making decisions under uncertainty.

Key Components

1. Null hypothesis (H₀): The claim being tested (usually status quo)
2. Alternative hypothesis (H₁ or Hₐ): The competing claim
3. Test statistic: A value calculated from sample data
4. p-value: Probability of observing data as extreme as ours, assuming H₀ is true
5. Decision rule: Criteria for rejecting or failing to reject H₀

The Criminal Trial Analogy

| Criminal Trial | Hypothesis Test | | —————————- | —————————– | | Defendant presumed innocent | H₀ presumed true | | Prosecution must prove guilt | Must show evidence against H₀ | | “Beyond reasonable doubt” | Significance level α | | Guilty verdict | Reject H₀ | | Not guilty verdict | Fail to reject H₀ |

11.2 Null and Alternative Hypotheses

Null Hypothesis (H₀)

• Statement of “no effect” or “no difference”
• Represents the status quo
• Must contain equality (=, ≤, or ≥)
• What we assume true until proven otherwise

Alternative Hypothesis (H₁ or Hₐ)

• Statement we’re trying to find evidence for
• Represents the research claim
• Determines if test is one-tailed or two-tailed

Types of Tests

1. Two-tailed test: H₁: μ ≠ μ₀
   • Tests for difference in either direction
   • Example: “The mean is different from 50”
2. Right-tailed test: H₁: μ > μ₀
   • Tests if parameter is greater
   • Example: “The mean is greater than 50”
3. Left-tailed test: H₁: μ < μ₀
   • Tests if parameter is less
   • Example: “The mean is less than 50”

Examples

Claim: Average response time is 5 minutes
H₀: μ = 5
H₁: μ ≠ 5 (two-tailed)

Claim: New drug lowers blood pressure
H₀: μ₁ ≥ μ₂ (new drug mean ≥ placebo mean)
H₁: μ₁ < μ₂ (new drug mean < placebo mean)

Claim: More than 60% support the proposal
H₀: p ≤ 0.60
H₁: p > 0.60

11.3 Type I and Type II Errors

Error Types

Type I Error (False Positive)

• Rejecting H₀ when it’s actually true
• Probability = α (significance level)
• Example: Convicting an innocent person

Type II Error (False Negative)

• Failing to reject H₀ when it’s actually false
• Probability = β
• Example: Acquitting a guilty person

Power of a Test

• Power = 1 - β
• Probability of correctly rejecting false H₀
• Affected by:
  • Sample size (larger n → higher power)
  • Effect size (larger effect → higher power)
  • Significance level α (larger α → higher power)
  • Population variability (smaller σ → higher power)

11.4 Significance Level and p-values

Significance Level (α)

• Threshold for rejecting H₀
• Pre-specified before testing
• Common values: 0.05, 0.01, 0.10
• Represents acceptable Type I error rate

p-value

• Probability of observing data as extreme as ours (or more extreme), assuming H₀ is true
• Smaller p-value = stronger evidence against H₀
• Compare to α to make decision

Decision Rule

• If p-value ≤ α: Reject H₀
• If p-value > α: Fail to reject H₀

Interpreting p-values

• p < 0.01: Very strong evidence against H₀
• 0.01 ≤ p < 0.05: Strong evidence against H₀
• 0.05 ≤ p < 0.10: Weak evidence against H₀
• p ≥ 0.10: Little or no evidence against H₀

11.5 Test Statistics

General Form

Common Test Statistics

1. Z-statistic (known σ or large samples): \(Z = \frac{\bar{X} - \mu_0}{\sigma/\sqrt{n}}\)
2. t-statistic (unknown σ): \(t = \frac{\bar{X} - \mu_0}{s/\sqrt{n}}\)
3. Z-statistic for proportions: \(Z = \frac{\hat{p} - p_0}{\sqrt{p_0(1-p_0)/n}}\)

Critical Values

Values that separate rejection and non-rejection regions:

• Two-tailed: ±z_{α/2} or ±t_{α/2,df}
• Right-tailed: z_α or t_{α,df}
• Left-tailed: -z_α or -t_{α,df}

11.6 Steps in Hypothesis Testing

1. State the hypotheses
   • Identify H₀ and H₁
   • Determine if one-tailed or two-tailed
2. Set significance level
   • Choose α (typically 0.05)
3. Check conditions
   • Random sampling
   • Independence
   • Normality (if required)
4. Calculate test statistic
   • Use appropriate formula
5. Find p-value
   • Use distribution tables or software
6. Make decision
   • Compare p-value to α
7. State conclusion
   • In context of the problem

Example: Testing Mean Response Time

A company claims average response time is 5 minutes. Sample of 100 calls shows \(\bar{x} = 5.5\) minutes, \(s = 2\) minutes.

# 1. Hypotheses
# H₀: μ = 5
# H₁: μ ≠ 5

# 2. Significance level
alpha <- 0.05

# 3. Check conditions (assumed met)

# 4. Calculate test statistic
x_bar <- 5.5
mu_0 <- 5
s <- 2
n <- 100
t_stat <- (x_bar - mu_0) / (s / sqrt(n))

# 5. Find p-value
p_value <- 2 * pt(abs(t_stat), df = n - 1, lower.tail = FALSE)

# 6. Decision
if (p_value <= alpha) {
  decision <- "Reject H₀"
} else {
  decision <- "Fail to reject H₀"
}

# 7. Conclusion
cat("t-statistic:", t_stat, "\n")
cat("p-value:", p_value, "\n")
cat("Decision:", decision, "\n")

11.7 Relationship Between Confidence Intervals and Hypothesis Tests

For two-tailed tests at significance level α:

• H₀ is rejected if and only if the (1-α)100% CI doesn’t contain the hypothesized value
• The CI contains all values that would not be rejected by the hypothesis test

Example

95% CI for μ: [52, 58]

• H₀: μ = 50 → Reject (50 not in CI)
• H₀: μ = 55 → Fail to reject (55 in CI)

11.8 Statistical vs. Practical Significance

Statistical Significance

• p-value < α
• Evidence against H₀
• Affected by sample size

Practical Significance

• Size of the effect
• Real-world importance
• Context-dependent

Key Point

Large samples can make tiny differences statistically significant, but they may not be practically important.

11.9 Common Misconceptions

What p-values are NOT

1. ❌ Probability that H₀ is true
2. ❌ Probability that H₁ is true
3. ❌ Probability of making a Type I error (that’s α)
4. ❌ Measure of effect size

What “fail to reject H₀” does NOT mean

1. ❌ H₀ is proven true
2. ❌ H₀ is accepted
3. ❌ There is no effect

11.10 Using R for Hypothesis Testing

Traditional Methods

# One-sample t-test
t.test(data, mu = hypothesized_value, alternative = "two.sided")

# Two-sample t-test
t.test(group1, group2, alternative = "two.sided")

# Paired t-test
t.test(x, y, paired = TRUE, alternative = "two.sided")

# Proportion test
prop.test(x = successes, n = total, p = hypothesized_prop,
          alternative = "two.sided")

Modern Approach with infer

library(infer)

# Generate null distribution
null_dist <- data %>%
  specify(response = variable) %>%
  hypothesize(null = "point", mu = hypothesized_value) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean")

# Calculate p-value
p_value <- null_dist %>%
  get_p_value(obs_stat = observed_stat, direction = "two-sided")

# Visualize
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = observed_stat, direction = "two-sided")

11.11 Reporting Results

A complete hypothesis test report should include:

1. Statement of hypotheses
2. Significance level used
3. Test statistic value
4. p-value
5. Decision (reject/fail to reject)
6. Conclusion in context

Example Report

“We tested whether the mean response time differs from 5 minutes (H₀: μ = 5 vs. H₁: μ ≠ 5) at the α = 0.05 significance level. Based on a sample of 100 calls with mean 5.5 minutes and standard deviation 2 minutes, the test statistic was t = 2.5 with p-value = 0.014. Since p < 0.05, we reject H₀ and conclude there is sufficient evidence that the mean response time differs from 5 minutes.”

Part 6: One-Sample Tests

12. One-Sample Tests

12.1 One-Sample Z-Test

The z-test is used when:

• Population standard deviation (σ) is known
• Sample size is large (n ≥ 30)
• Population is normally distributed (or large sample)

Test Statistic

Decision Rule

• Two-tailed: Reject H₀ if |Z| > z_{α/2}
• Right-tailed: Reject H₀ if Z > z_α
• Left-tailed: Reject H₀ if Z < -z_α

Example: Rice Consumption

Testing if average rice consumption differs from 3.9 kg/year, with known σ = 4.5:

# Given data
x_bar <- 4.2     # Sample mean
mu_0 <- 3.9      # Hypothesized mean
sigma <- 4.5     # Known population SD
n <- 150         # Sample size
alpha <- 0.05    # Significance level

# Calculate z-statistic
z_stat <- (x_bar - mu_0) / (sigma / sqrt(n))

# Calculate p-value (two-tailed)
p_value <- 2 * (1 - pnorm(abs(z_stat)))

# Critical values
z_critical <- qnorm(1 - alpha/2)

# Decision
if (abs(z_stat) > z_critical) {
  print("Reject H₀")
} else {
  print("Fail to reject H₀")
}

# Using R's built-in function (requires BSDA package)
library(BSDA)
z.test(x = sample_data, mu = 3.9, sigma.x = 4.5, 
       alternative = "two.sided")

12.2 One-Sample t-Test

The t-test is used when:

• Population standard deviation is unknown
• Sample standard deviation (s) is used instead
• Population is approximately normal (or large sample)

Test Statistic

Degrees of Freedom

df = n - 1

Example: Call Center Response Time

Testing if mean response time equals 120 minutes:

# Sample data
response_times <- c(115, 125, 118, 122, 130, 119, 121, 117, 124, 116)
mu_0 <- 120
alpha <- 0.05

# Using t.test function
result <- t.test(response_times, mu = mu_0, alternative = "two.sided")
print(result)

# Manual calculation
x_bar <- mean(response_times)
s <- sd(response_times)
n <- length(response_times)
se <- s / sqrt(n)
t_stat <- (x_bar - mu_0) / se
df <- n - 1
p_value <- 2 * pt(-abs(t_stat), df)

# Confidence interval
t_critical <- qt(1 - alpha/2, df)
margin_error <- t_critical * se
ci_lower <- x_bar - margin_error
ci_upper <- x_bar + margin_error

cat("t-statistic:", t_stat, "\n")
cat("p-value:", p_value, "\n")
cat("95% CI: [", ci_lower, ",", ci_upper, "]\n")

Using infer Package

library(infer)

# Calculate observed statistic
obs_stat <- response_times %>%
  specify(response = response_times) %>%
  calculate(stat = "t", mu = mu_0)

# Generate null distribution
null_dist <- response_times %>%
  specify(response = response_times) %>%
  hypothesize(null = "point", mu = mu_0) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "t", mu = mu_0)

# Get p-value
p_value <- null_dist %>%
  get_p_value(obs_stat = obs_stat, direction = "two-sided")

# Visualize
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = obs_stat, direction = "two-sided")

12.3 One-Sample Test for Proportion

Used to test claims about population proportions.

Test Statistic

Conditions

1. Random sample
2. Independence
3. Success-failure condition: np₀ ≥ 10 and n(1-p₀) ≥ 10

Example: Customer Satisfaction

Testing if satisfaction rate equals 80%:

# Given data
n <- 100         # Sample size
x <- 73          # Number of successes
p_0 <- 0.80      # Hypothesized proportion
alpha <- 0.05

# Calculate sample proportion
p_hat <- x / n

# Calculate test statistic
z_stat <- (p_hat - p_0) / sqrt(p_0 * (1 - p_0) / n)

# Calculate p-value (two-tailed)
p_value <- 2 * pnorm(-abs(z_stat))

# Using prop.test
result <- prop.test(x = x, n = n, p = p_0, 
                    alternative = "two.sided", 
                    correct = FALSE)

# Using infer
satisfaction_data <- data.frame(
  satisfied = c(rep("yes", x), rep("no", n - x))
)

# Observed statistic
obs_stat <- satisfaction_data %>%
  specify(response = satisfied, success = "yes") %>%
  calculate(stat = "prop")

# Null distribution
null_dist <- satisfaction_data %>%
  specify(response = satisfied, success = "yes") %>%
  hypothesize(null = "point", p = p_0) %>%
  generate(reps = 1000, type = "draw") %>%
  calculate(stat = "prop")

# P-value
p_value <- null_dist %>%
  get_p_value(obs_stat = obs_stat, direction = "two-sided")

12.4 One-Tailed Tests

Right-Tailed Test (Greater Than)

H₀: μ ≤ μ₀ H₁: μ > μ₀

# Right-tailed t-test
t.test(data, mu = mu_0, alternative = "greater")

# Right-tailed proportion test
prop.test(x, n, p = p_0, alternative = "greater")

Left-Tailed Test (Less Than)

H₀: μ ≥ μ₀ H₁: μ < μ₀

# Left-tailed t-test
t.test(data, mu = mu_0, alternative = "less")

# Left-tailed proportion test
prop.test(x, n, p = p_0, alternative = "less")

12.5 Chi-Square Test for Variance

Testing claims about population variance σ²:

Test Statistic

Degrees of Freedom

df = n - 1

Example

# Sample data
sample_data <- c(10, 12, 9, 11, 13, 10, 12, 11, 9, 10)
sigma_0_squared <- 4  # Hypothesized variance
alpha <- 0.05

# Calculate test statistic
n <- length(sample_data)
s_squared <- var(sample_data)
chi_squared <- (n - 1) * s_squared / sigma_0_squared

# Calculate p-value (two-tailed)
p_value_lower <- pchisq(chi_squared, df = n - 1)
p_value_upper <- 1 - pchisq(chi_squared, df = n - 1)
p_value <- 2 * min(p_value_lower, p_value_upper)

# Critical values
chi_lower <- qchisq(alpha/2, df = n - 1)
chi_upper <- qchisq(1 - alpha/2, df = n - 1)

12.6 Sign Test (Non-parametric)

Used when normality assumption is violated.

Procedure

1. Count values above and below hypothesized median
2. Use binomial distribution with p = 0.5

# Sign test for median
library(BSDA)
SIGN.test(data, md = hypothesized_median, 
          alternative = "two.sided")

# Manual implementation
above <- sum(data > hypothesized_median)
below <- sum(data < hypothesized_median)
n_eff <- above + below
p_value <- 2 * pbinom(min(above, below), n_eff, 0.5)

12.7 Practical Examples

Example 1: Age at First Marriage

Testing if mean age at first marriage > 23 years:

# Load data
age_data <- read.csv("age_at_marriage.csv")
mu_0 <- 23

# Check normality
hist(age_data$$age, main = "Distribution of Age at First Marriage")
qqnorm(age_data$$age)
qqline(age_data$$age)

# Perform t-test
result <- t.test(age_data$$age, mu = mu_0, alternative = "greater")
print(result)

# Using infer
obs_stat <- age_data %>%
  specify(response = age) %>%
  calculate(stat = "mean")

null_dist <- age_data %>%
  specify(response = age) %>%
  hypothesize(null = "point", mu = mu_0) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean")

p_value <- null_dist %>%
  get_p_value(obs_stat = obs_stat, direction = "greater")

# Visualize
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = obs_stat, direction = "greater")

Example 2: Proportion of Female Students

Testing if proportion of females = 0.50:

# Survey data
total_students <- 500
female_students <- 237

# Traditional approach
prop.test(x = female_students, n = total_students, 
          p = 0.50, alternative = "two.sided")

# Using infer
gender_data <- data.frame(
  gender = c(rep("female", female_students), 
             rep("male", total_students - female_students))
)

obs_stat <- gender_data %>%
  specify(response = gender, success = "female") %>%
  calculate(stat = "prop")

null_dist <- gender_data %>%
  specify(response = gender, success = "female") %>%
  hypothesize(null = "point", p = 0.50) %>%
  generate(reps = 1000, type = "draw") %>%
  calculate(stat = "prop")

p_value <- null_dist %>%
  get_p_value(obs_stat = obs_stat, direction = "two-sided")

12.8 Effect Size for One-Sample Tests

Cohen’s d for Means

Interpretation:

• Small effect: d ≈ 0.2
• Medium effect: d ≈ 0.5
• Large effect: d ≈ 0.8

# Calculate Cohen's d
cohens_d <- abs(mean(data) - mu_0) / sd(data)

Effect Size for Proportions

# Calculate h
h <- 2 * asin(sqrt(p_hat)) - 2 * asin(sqrt(p_0))

12.9 Power Analysis

For One-Sample t-Test

library(pwr)

# Power calculation
pwr.t.test(n = NULL, d = 0.5, sig.level = 0.05, 
           power = 0.80, type = "one.sample")

# Sample size calculation
pwr.t.test(d = 0.5, sig.level = 0.05, power = 0.80, 
           type = "one.sample")

For One-Sample Proportion Test

# Power calculation
pwr.p.test(h = ES.h(p1 = 0.75, p2 = 0.80), 
           sig.level = 0.05, power = 0.80)

12.10 Reporting Results

Template for Reporting

“A one-sample [test type] was conducted to determine whether [parameter] differed from [hypothesized value]. The sample [statistic] was [value] ([measure of spread]). The test was [significant/not significant], t(df) = [test statistic], p = [p-value]. [Effect size measure] indicated a [small/medium/large] effect size. We [reject/fail to reject] the null hypothesis and conclude that [contextual conclusion].”

Example Report

“A one-sample t-test was conducted to determine whether the mean response time differed from 120 minutes. The sample mean was 118.7 minutes (SD = 4.2). The test was not significant, t(9) = -0.98, p = 0.353. Cohen’s d = 0.31 indicated a small effect size. We fail to reject the null hypothesis and conclude that there is insufficient evidence to suggest the mean response time differs from 120 minutes.”

Part 7: Two-Sample Tests

13. Two-Sample Tests

13.1 Independent Samples t-Test

The independent samples t-test compares means from two unrelated groups.

Assumptions

1. Independent random samples
2. Approximately normal distributions (or large samples)
3. Equal variances (classical t-test) or unequal variances (Welch’s t-test)

Test Statistic (Equal Variances)

where pooled standard deviation: \(s_p = \sqrt{\frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1 + n_2 - 2}}\)

Test Statistic (Unequal Variances - Welch’s)

Degrees of Freedom (Welch’s)

Example: Income Comparison

Comparing average income between Cleveland and Sacramento:

# Load data
income_data <- read.delim("cleSac.txt")
cleveland <- income_data$$income[income_data$$metro_area == "Cleveland_OH"]
sacramento <- income_data$$income[income_data$$metro_area == "Sacramento_CA"]

# Check assumptions
# Normality
par(mfrow = c(1, 2))
hist(cleveland, main = "Cleveland Income")
hist(sacramento, main = "Sacramento Income")

# Equal variances
var.test(cleveland, sacramento)

# Perform t-test (Welch's by default)
t_test_result <- t.test(cleveland, sacramento, 
                        alternative = "two.sided",
                        var.equal = FALSE)
print(t_test_result)

# Manual calculation
x1_bar <- mean(cleveland)
x2_bar <- mean(sacramento)
s1 <- sd(cleveland)
s2 <- sd(sacramento)
n1 <- length(cleveland)
n2 <- length(sacramento)

t_stat <- (x1_bar - x2_bar) / sqrt(s1^2/n1 + s2^2/n2)
df <- ((s1^2/n1 + s2^2/n2)^2) / 
      ((s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1))
p_value <- 2 * pt(-abs(t_stat), df)

# Using infer
library(infer)
obs_diff <- income_data %>%
  specify(income ~ metro_area) %>%
  calculate(stat = "diff in means", 
            order = c("Cleveland_OH", "Sacramento_CA"))

null_dist <- income_data %>%
  specify(income ~ metro_area) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", 
            order = c("Cleveland_OH", "Sacramento_CA"))

p_value <- null_dist %>%
  get_p_value(obs_stat = obs_diff, direction = "two-sided")

# Visualize
null_dist %>%
  visualize() +
  shade_p_value(obs_stat = obs_diff, direction = "two-sided")

13.2 Paired Samples t-Test

Used when comparing two related measurements (e.g., before-after, matched pairs).

Test Statistic

where:

• \(\bar{d}\) = mean of differences
• \(s_d\) = standard deviation of differences
• \(\mu_d\) = hypothesized mean difference (usually 0)

Example: Weight Loss Program

# Before and after weights
before <- c(65, 60, 63, 58, 70, 72, 80, 81, 83, 69)
after <- c(55, 52, 60, 55, 60, 66, 70, 66, 63, 60)

# Calculate differences
differences <- after - before

# Check normality of differences
hist(differences, main = "Distribution of Weight Differences")
qqnorm(differences)
qqline(differences)

# Paired t-test
paired_result <- t.test(after, before, 
                        paired = TRUE, 
                        alternative = "less")
print(paired_result)

# Manual calculation
d_bar <- mean(differences)
s_d <- sd(differences)
n <- length(differences)
t_stat <- d_bar / (s_d / sqrt(n))
p_value <- pt(t_stat, df = n - 1)

# Using infer
weight_data <- data.frame(
  id = rep(1:10, 2),
  weight = c(before, after),
  time = rep(c("before", "after"), each = 10)
)

obs_diff <- weight_data %>%
  specify(weight ~ time) %>%
  calculate(stat = "diff in means", 
            order = c("after", "before"))

null_dist <- weight_data %>%
  specify(weight ~ time) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", 
            order = c("after", "before"))

p_value <- null_dist %>%
  get_p_value(obs_stat = obs_diff, direction = "less")

Example: Water Quality Testing

Comparing zinc concentration in surface vs. bottom water:

# Load and prepare data
zinc_data <- read.csv("zinc_tidy.csv")
zinc_diff <- zinc_data %>%
  group_by(loc_id) %>%
  summarize(pair_diff = diff(concentration))

# Test if surface water has lower concentration
t.test(zinc_diff$$pair_diff, 
       mu = 0, 
       alternative = "less")

# Visual comparison
library(ggplot2)
ggplot(zinc_data, aes(x = depth, y = concentration)) +
  geom_boxplot() +
  geom_line(aes(group = loc_id), alpha = 0.3) +
  geom_point(aes(color = depth), size = 3) +
  labs(title = "Zinc Concentration by Depth",
       x = "Water Depth",
       y = "Zinc Concentration")

13.3 Two-Sample Test for Proportions

Compares proportions between two independent groups.

Test Statistic

where pooled proportion: \(\hat{p} = \frac{x_1 + x_2}{n_1 + n_2}\)

Example: Gender and Weight Loss

Testing if weight loss rates differ by gender:

# Data
male_success <- 45
male_total <- 100
female_success <- 62
female_total <- 100

# Two-proportion test
prop.test(x = c(male_success, female_success),
          n = c(male_total, female_total),
          alternative = "two.sided",
          correct = FALSE)

# Manual calculation
p1_hat <- male_success / male_total
p2_hat <- female_success / female_total
p_pool <- (male_success + female_success) / (male_total + female_total)
se_diff <- sqrt(p_pool * (1 - p_pool) * (1/male_total + 1/female_total))
z_stat <- (p1_hat - p2_hat) / se_diff
p_value <- 2 * pnorm(-abs(z_stat))

# Using infer
gender_data <- data.frame(
  gender = c(rep("male", male_total), rep("female", female_total)),
  success = c(rep(c("yes", "no"), c(male_success, male_total - male_success)),
              rep(c("yes", "no"), c(female_success, female_total - female_success)))
)

obs_diff <- gender_data %>%
  specify(success ~ gender, success = "yes") %>%
  calculate(stat = "diff in props", 
            order = c("male", "female"))

null_dist <- gender_data %>%
  specify(success ~ gender, success = "yes") %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in props", 
            order = c("male", "female"))

p_value <- null_dist %>%
  get_p_value(obs_stat = obs_diff, direction = "two-sided")

13.4 Mann-Whitney U Test (Wilcoxon Rank-Sum)

Non-parametric alternative to independent samples t-test.

When to Use

• Data not normally distributed
• Ordinal data
• Small sample sizes
• Presence of outliers

# Mann-Whitney U test
wilcox.test(cleveland, sacramento, 
            alternative = "two.sided")

# With continuity correction
wilcox.test(cleveland, sacramento, 
            alternative = "two.sided", 
            correct = TRUE)

# Visual comparison using violin plots
library(ggplot2)
ggplot(income_data, aes(x = metro_area, y = income)) +
  geom_violin(fill = "lightblue") +
  geom_boxplot(width = 0.1) +
  labs(title = "Income Distribution by City",
       x = "Metropolitan Area",
       y = "Income")

13.5 Wilcoxon Signed-Rank Test

Non-parametric alternative to paired t-test.

# Wilcoxon signed-rank test
wilcox.test(after, before, 
            paired = TRUE, 
            alternative = "less")

# Manual implementation
differences <- after - before
ranks <- rank(abs(differences))
signed_ranks <- ranks * sign(differences)
W <- sum(signed_ranks[signed_ranks > 0])

13.6 Chi-Square Test for Independence

Tests whether two categorical variables are independent.

Test Statistic

where:

• O = observed frequency
• E = expected frequency

Example: Species and Petal Width

# Prepare data
iris_subset <- iris[iris$$Species != "virginica", ]
iris_subset$$Species <- factor(iris_subset$$Species)
iris_subset$$Petal.Width.Cat <- cut(iris_subset$$Petal.Width, 
                                   breaks = quantile(iris_subset$$Petal.Width, 
                                                    probs = c(0, 0.5, 1)), 
                                   include.lowest = TRUE,
                                   labels = c("Below", "Above"))

# Create contingency table
cont_table <- table(iris_subset$$Petal.Width.Cat, iris_subset$$Species)
print(cont_table)

# Chi-square test
chi_test <- chisq.test(cont_table)
print(chi_test)

# Expected values
print(chi_test$$expected)

# Residuals
print(chi_test$$residuals)

# Visualization
library(ggplot2)
ggplot(iris_subset, aes(x = Species, fill = Petal.Width.Cat)) +
  geom_bar(position = "dodge") +
  labs(title = "Distribution of Petal Width by Species",
       x = "Species",
       y = "Count",
       fill = "Petal Width")

13.7 Effect Size for Two-Sample Tests

Cohen’s d for Independent Samples

Cohen’s d for Paired Samples

# Effect size for independent samples
library(effsize)
cohen.d(cleveland, sacramento)

# Effect size for paired samples
cohen.d(after, before, paired = TRUE)

# Manual calculation
d_independent <- abs(mean(cleveland) - mean(sacramento)) / 
                 sqrt(((length(cleveland) - 1) * var(cleveland) + 
                       (length(sacramento) - 1) * var(sacramento)) / 
                      (length(cleveland) + length(sacramento) - 2))

d_paired <- abs(mean(differences)) / sd(differences)

13.8 Power Analysis for Two-Sample Tests

library(pwr)

# Power for independent t-test
pwr.t2n.test(n1 = 100, n2 = 120, d = 0.5, 
             sig.level = 0.05, 
             alternative = "two.sided")

# Sample size for paired t-test
pwr.t.test(d = 0.3, sig.level = 0.05, 
           power = 0.80, type = "paired")

# Power for two proportions
pwr.2p.test(h = ES.h(p1 = 0.45, p2 = 0.60), 
            sig.level = 0.05, power = 0.80)

13.9 Comprehensive Example: Medicine Effectiveness

# Load data
medicine_data <- read.csv("medicine_survey.txt")

# Research question: Do females have less weight than males before treatment?

# Subset data
female_weight <- medicine_data$$KgBefore[medicine_data$$Gender == "Female"]
male_weight <- medicine_data$$KgBefore[medicine_data$$Gender == "Male"]

# Exploratory analysis
library(ggplot2)
ggplot(medicine_data, aes(x = Gender, y = KgBefore, fill = Gender)) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.5) +
  labs(title = "Weight Distribution by Gender",
       x = "Gender",
       y = "Weight Before Treatment (kg)")

# Check assumptions
shapiro.test(female_weight)
shapiro.test(male_weight)
var.test(female_weight, male_weight)

# Perform t-test
result <- t.test(female_weight, male_weight, 
                 alternative = "less",
                 var.equal = FALSE)
print(result)

# Effect size
cohen.d(female_weight, male_weight)

# Bootstrap confidence interval
library(infer)
boot_dist <- medicine_data %>%
  specify(KgBefore ~ Gender) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "diff in means", 
            order = c("Female", "Male"))

ci_boot <- boot_dist %>%
  get_confidence_interval(level = 0.95, type = "percentile")

# Report results
cat("Mean weight (Female):", mean(female_weight), "kg\n")
cat("Mean weight (Male):", mean(male_weight), "kg\n")
cat("Difference:", mean(female_weight) - mean(male_weight), "kg\n")
cat("95% CI (traditional):", result$$conf.int, "\n")
cat("95% CI (bootstrap):", ci_boot$$lower_ci, "to", ci_boot$$upper_ci, "\n")
cat("p-value:", result$$p.value, "\n")

13.10 Reporting Two-Sample Test Results

Template

“A [test type] was conducted to compare [variable] between [group 1] and [group 2]. There was a [significant/non-significant] difference in [variable] between [group 1] (M = [mean1], SD = [sd1]) and [group 2] (M = [mean2], SD = [sd2]); t(df) = [test statistic], p = [p-value]. The effect size was [small/medium/large] (d = [effect size]). These results suggest that [contextual conclusion].”

Example Report

“An independent samples t-test was conducted to compare weight before treatment between females and males. There was a significant difference in weight between females (M = 65.4, SD = 8.2) and males (M = 72.1, SD = 9.5); t(198) = -5.33, p < .001. The effect size was medium (d = 0.75). These results suggest that females had significantly lower weight than males before treatment.”

Part 8: R Implementation & Case Studies

14. R Implementation: Traditional vs Simulation Methods

14.1 Traditional Methods Overview

Traditional statistical methods rely on mathematical theory and distributional assumptions.

Key Functions in Base R

# t-tests
t.test()           # One-sample, two-sample, paired t-tests
var.test()         # F-test for equality of variances

# Proportion tests
prop.test()        # One-sample, two-sample proportion tests
binom.test()       # Exact binomial test

# Non-parametric tests
wilcox.test()      # Wilcoxon rank-sum/signed-rank tests
chisq.test()       # Chi-square test

# ANOVA
aov()              # Analysis of variance
summary()          # ANOVA table

Example: Traditional t-test workflow

# Load data
data <- read.csv("sample_data.csv")

# Check assumptions
shapiro.test(data$$response)  # Normality
var.test(group1, group2)     # Equal variances

# Perform test
result <- t.test(response ~ group, data = data)

# Extract results
result$$statistic  # Test statistic
result$$p.value    # p-value
result$$conf.int   # Confidence interval