--- title: "Hypothesis Testing" author: "JJB + Course" date: "9/21/2018" output: html_document: toc: true toc_float: collapse: false --- # Testing Frameworks ## Example: Types of Sampling Distribution ```{r ex-distribution-overlay} library("ggplot2") ggplot(data.frame(x = c(-5, 5)), aes(x)) + stat_function(fun = dnorm, aes(colour = "Normal")) + stat_function(fun = dt, args = list(df = 1), aes(colour = "t, df = 1")) + stat_function(fun = dt, args = list(df = 2), aes(colour = "t, df = 2")) + stat_function(fun = dt, args = list(df = 3), aes(colour = "t, df = 3")) + stat_function(fun = dt, args = list(df = 5), aes(colour = "t, df = 5")) + stat_function(fun = dt, args = list(df = 10), aes(colour = "t, df = 10")) + stat_function(fun = dt, args = list(df = 20), aes(colour = "t, df = 20")) + labs(title = "PDFs for Normal and t Distributions", y = "f(x)", x = "x", colour = "Distributions") + theme_bw() + theme(axis.text = element_text(size = 12), axis.title = element_text(size = 13), plot.title = element_text(size = 16)) ``` ## Example: Different Sample Groups From the graph, we think the two groups might be significantly different. How can we check though? ```{r initial-height-dis} # Simulate data set.seed(183) n = 300 tdata = data.frame(Sex = rep(c("Male", "Female"), each = n), Height = c(rnorm(n, 6.2), rnorm(n, mean=5))) # Compute mean statistics library("dplyr") data_means = tdata %>% group_by(Sex) %>% summarise(height_mean = mean(Height)) # Create a density plot of two groups library("ggplot2") ggplot(tdata, aes(x = Height, fill = Sex, color = Sex)) + geom_density(alpha = 0.5) + geom_vline(data = data_means, aes(xintercept = height_mean, color = Sex), linetype = "dashed", size=1) + theme_bw() + labs(title = "Comparison of Height Distributions between Sex", y = "Density of Height Distribution") ``` Generate a boxplot as an alternative ```{r alternative-plot} ggplot(tdata, aes(x = Sex, y = Height, fill = Sex, color = Sex)) + geom_boxplot(alpha = 0.5) + theme_bw() + labs(title = "Comparison of Height Distributions between Sex") ``` Would this graph be a good indicator of difference? ```{r height-graph-v2} set.seed(183) n = 300 tdata = data.frame(Sex = rep(c("Male", "Female"), each = n), Height = c(rnorm(n, 5.2), rnorm(n, mean=5))) library("dplyr") data_means = tdata %>% group_by(Sex) %>% summarise(height_mean = mean(Height)) library("ggplot2") ggplot(tdata, aes(x = Height, fill = Sex, color = Sex)) + geom_density(alpha = 0.5) + geom_vline(data = data_means, aes(xintercept = height_mean, color = Sex), linetype = "dashed", size=1) + theme_bw() + labs(title = "Comparison of Height Distributions between Sex", y = "Density of Height Distribution") ``` ## Example: Critical Regions for Hypotheses (two vs. one-sided) Two-sided example ```{r two-sided-vis-t-test} library("ggplot2") # Determine values for viewing t-distribution df = 20 min_val = -5 max_val = 5 graph_bounds = data.frame(x = c(min_val, max_val)) # Pick out significance levels alpha = 0.05 lower_tail = alpha / 2 upper_tail = 1 - alpha / 2 # Calculate critical values lower_tail_crit = round(qt(lower_tail, df = df), digits = 3) upper_tail_crit = round(qt(upper_tail, df = df), digits = 3) # Graph the distribution with critical region g2 = ggplot(graph_bounds, aes(x = x)) + stat_function(fun = dt, args = list(df = df)) + stat_function(fun = dt, args = list(df = df), xlim = c(min_val, lower_tail_crit), geom = "area", fill = "blue") + stat_function(fun = dt, args = list(df = df), xlim = c(max_val, upper_tail_crit), geom = "area", fill = "blue") + geom_vline(xintercept = c(lower_tail_crit, upper_tail_crit), linetype = "dashed", size=0.5) + scale_x_continuous(breaks = c(lower_tail_crit, 0, upper_tail_crit)) + labs(title = "Two-Sided t-Test", subtitle = paste0("with alpha = ", alpha, "/2 and df = ", df), x = "t-Value", y = "Density") + theme(legend.position = "none", axis.text = element_text(size = 20)) + theme_bw() g2 ``` One-sided example ```{r one-sided-vis-t-test} # Pick out significance levels upper_tail = 1 - alpha upper_tail_crit = round(qt(upper_tail, df = df), digits = 3) # Graph the distribution with critical region ggplot(graph_bounds, aes(x = x)) + stat_function(fun = dt, args = list(df = df)) + stat_function(fun = dt, args = list(df = df), xlim = c(max_val, upper_tail_crit), geom = "area", fill = "blue") + geom_vline(xintercept = upper_tail_crit, linetype = "dashed", size=0.5) + scale_x_continuous(breaks = c(0, upper_tail_crit)) + labs(title = "One-Sided t-Test", subtitle = paste0("with alpha = ", alpha, " and df = ", df), x = "t-Value", y = "Density") + theme(legend.position = "none", axis.text = element_text(size = 20)) + theme_bw() ``` ## Example: Generating Sample Data ```{r} # Set seed for reproducibility set.seed(881) # Generate data n = 10 x1 = round(rnorm(n), 1) x2 = round(rnorm(n) + 1, 1) ``` ## Example: Calculating a t-statistic ```{r free-hand-t-calculation} # Compute means of each group x1_mu = mean(x1) x2_mu = mean(x2) # Compute length and degrees of freedom n1 = length(x1) n2 = length(x2) ndf = n1 + n2 - 2 # Calculate pooled variance s2 = ((n1 - 1) * var(x1) + (n2 - 1) * var(x2)) / ndf # Compute the t-statistic tstat = (mean(x1) - mean(x2)) / sqrt(s2 * (1 / n1 + 1 / n2)) ``` ## Example: Floating Point Stability ```{r numerics-problematic} # Numerics are problematic 0.10 + 0.05 == 0.15 # Allow for tolerance with numerics via an epsilon neighborhood all.equal(0.10 + 0.05, 0.15) # Lack of output stability thoughâ€¦ all.equal(0.12, 0.19) # Check for whether it is true isTRUE(all.equal(0.12, 0.19)) ``` ## Example: Calculate P-value Two-sided hypothesis p-value ```{r calc-p-val-two-sided} two_sided_p = 2 * (1 - pt(abs(tstat), ndf)) two_sided_p ``` Upper hypothesis p-value ```{r calc-p-val-one-sided} one_sided_p = 1 - pt(tstat, ndf) one_sided_p ``` ```{r alt-calc-pval-one-sided} # Alternatively, start from the right (upper) tail # instead of the left (lower) tail. pt(tstat, ndf, lower.tail = FALSE) ``` Two-sided example ```{r two-sided-vis-t-test-stat} library("ggplot2") # Determine values for viewing t-distribution df = 18 min_val = -5 max_val = 5 graph_bounds = data.frame(x = c(min_val, max_val)) # Pick out significance levels alpha = 0.05 lower_tail = alpha / 2 upper_tail = 1 - alpha / 2 # Calculate critical values lower_tail_crit = round(qt(lower_tail, df = df), digits = 3) upper_tail_crit = round(qt(upper_tail, df = df), digits = 3) # Graph the distribution with critical region ggplot(graph_bounds, aes(x = x)) + stat_function(fun = dt, args = list(df = df)) + stat_function(fun = dt, args = list(df = df), xlim = c(min_val, lower_tail_crit), geom = "area", fill = "blue", alpha = 0.75) + stat_function(fun = dt, args = list(df = df), xlim = c(min_val, qt(two_sided_p, df = df)), geom = "area", fill = "orange", alpha = 0.75) + stat_function(fun = dt, args = list(df = df), xlim = c(max_val, upper_tail_crit), geom = "area", fill = "blue") + geom_vline(xintercept = c(lower_tail_crit, upper_tail_crit), linetype = "dashed", size=0.5) + geom_vline(xintercept = qt(two_sided_p, df = df), linetype="dotted", size=0.75, color = "orange") + scale_x_continuous(breaks = c(lower_tail_crit, 0, upper_tail_crit)) + labs(title = "Two-Sided t-Test", subtitle = paste0("with alpha = ", alpha, "/2 and df = ", df), x = "t-Value", y = "Density") + theme(legend.position = "none", axis.text = element_text(size = 20)) + theme_bw() ``` One-sided example ```{r one-sided-vis-t-test-stat} # Pick out significance levels upper_tail = 1 - alpha upper_tail_crit = round(qt(upper_tail, df = df), digits = 3) # Graph the distribution with critical region ggplot(graph_bounds, aes(x = x)) + stat_function(fun = dt, args = list(df = df)) + stat_function(fun = dt, args = list(df = df), xlim = c(max_val, upper_tail_crit), geom = "area", fill = "blue") + stat_function(fun = dt, args = list(df = df), xlim = c(max_val, qt(1-one_sided_p, df = df)), geom = "area", fill = "orange", alpha = 0.75) + geom_vline(xintercept = upper_tail_crit, linetype = "dashed", size=0.5) + geom_vline(xintercept = qt(1-one_sided_p, df = df), linetype="dotted", size=0.75, color = "orange") + scale_x_continuous(breaks = c(0, upper_tail_crit)) + labs(title = "One-Sided t-Test", subtitle = paste0("with alpha = ", alpha, " and df = ", df), x = "t-Value", y = "Density") + theme(legend.position = "none", axis.text = element_text(size = 20)) + theme_bw() ``` ## Example: Calculate Critical t-value ```{r calc-t-val-two-sided} # Significance Level alpha = 0.05 # Critical value for two-sided test qt(1 - alpha/2, ndf) ``` ```{r calc-t-val-onesided} # Significance Level alpha = 0.05 # Critical value for one-sided test qt(1 - alpha, ndf) ``` ## Example: Significant values for a t-Distribution table Create a t-distribution table ```{r build-t-table} alpha = c(0.4, 0.25, 0.1, 0.05, 0.025, 0.01, 0.005, 0.0005) df = c(1:30, 40, 60, Inf) t_table = outer(df, alpha, function(df, alpha) qt(1 - alpha, df = df)) rownames(t_table) = df colnames(t_table) = alpha knitr::kable(t_table) ``` ## Example: t-test Implemented To begin, we start by implementing the algorithm as if it just another function. After we have a stable implementation, then we can being adding in additional features. ```{r my-test-imp} my_ttest = function(x1, x2, test = c("two-sided", "lower", "upper"), alpha = 0.05) { # Force `test` to hold a pre-defined value of either: # "two-sided", "lower", or "upper" test = match.arg(test) # Compute length and degrees of freedom n1 = length(x1) n2 = length(x2) ndf = n1 + n2 - 2 # Calculate t-statistic s2 = ((n1 - 1) * var(x1) + (n2 - 1) * var(x2)) / ndf tstat = (mean(x1) - mean(x2)) / sqrt(s2 * (1 / n1 + 1 / n2)) # Compute tail probability tail_prob = switch(test, "two-sided" = 2 * (1 - pt(abs(tstat), ndf)), "lower" = pt(tstat, ndf), "upper" = 1 - pt(tstat, ndf)) # Format and return results results = list(tstat = tstat, df = ndf, reject = tail_prob < alpha, prob = tail_prob) return(results) } ``` Verify output of unpooled function against base R implementation ```{r simulate-data} # Simulate some data set.seed(881) n = 10 x1 = round(rnorm(n) , 1) x2 = round(rnorm(n) + 1 , 1) test_result = my_ttest(x1, x2) test_result # Check against built in all.equal(test_result[-3], t.test(x1,x2, var.equal = TRUE)[1:3], check.attributes = FALSE) ```