--- author: "Jinhyun Ju" output: html_document --- Quantitative Genomics and Genetics 2016 ====== Homework 4 R-code Key ------ -- 11 March 2016 -- Author: Jin Hyun Ju (jj328@cornell.edu) ### Problem 2 #### a) Determine the critical value for alpha = 0.05. (mu > 0) ```{r, comment = NA, fig.show = 'hold', fig.height = 4, fig.width = 8, fig.align='center'} # we are only interested in the positive half of the distribution. critical_a <- qnorm(0.05, mean = 0, sd=1 / sqrt(10), lower.tail = FALSE) cat("Critical value for a Type 1 error of 0.05 = ", critical_a, "\n") ``` #### b) Determine the critical value for alpha = 0.05. (mu != 0) ```{r, comment = NA, fig.show = 'hold', fig.height = 4, fig.width = 8, fig.align='center'} # In this case we are interested in the both sides of the distribution. critical_b <- qnorm(0.05 / 2, mean = 0, sd=1 / sqrt(10)) cat("Critical values for a Type 1 error of 0.05 = ", critical_b, -critical_b, "\n") ``` #### c) Function for generating random samples and testing significance ```{r, comment = NA, fig.show = 'hold', fig.height = 4, fig.width = 8, fig.align='center'} sampling_function<- function(number_of_samples = 1000, sample_size = 100, population_mean = 0, population_var = 1, c_a = qnorm(0.05, 0, 1/sqrt(100)), test_direction = c("one_sided","two_sided")){ sample_means <- c() population_sd <- sqrt(population_var) for ( i in 1:number_of_samples ){ sample_means[i] <- mean(rnorm(sample_size, mean = population_mean, sd = population_sd)) } if (test_direction == "one_sided"){ test_sample <- sum(sample_means >= c_a) } else { c_a <- abs(c_a) test_sample <- sum(sample_means >= c_a | sample_means <= -c_a) } cat("Null Hypothesis is rejected in ", test_sample," out of ",number_of_samples," cases", "(",100 * test_sample/number_of_samples ,"%)", "\n") return(test_sample) } sampling_function(number_of_samples = 1000, sample_size = 10, population_mean = 0, population_var = 1, c_a = qnorm(0.05, 0, 1/sqrt(10), lower.tail = FALSE), test_direction = c("one_sided")) sampling_function(number_of_samples = 1000, sample_size = 10, population_mean = 0, population_var = 1, c_a = abs(qnorm(0.05/2, 0, 1/sqrt(10))), test_direction = c("two_sided")) ``` #### d) Calculating the power for mu > 0. ```{r, comment = NA} type_2_error_d <- pnorm(critical_a, mean = 0.1, sd = 1 / sqrt(10)) power_of_test_d <- 1 - type_2_error_d cat("Power of the test for (a) is =", power_of_test_d, "if true mean is 0.1 \n") ``` #### e) Calculating the power for mu != 0. ```{r, comment = NA} type_2_error_e <- pnorm(abs(critical_b), mean = 0.1, sd = 1 / sqrt(10)) - pnorm(-abs(critical_b), mean = 0.1, sd = 1 / sqrt(10)) power_of_test_e <- 1 - type_2_error_e cat("Power of the test for (b) is =", power_of_test_e, "if true mean is 0.1 \n") ``` #### f) Calculating the power for mu > 0 (var = 10). ```{r, comment = NA} critical_f <- qnorm(0.05, mean = 0, sd=sqrt(10/10), lower.tail = FALSE) type_2_error_f <- pnorm(critical_f, mean = 0.1, sd = sqrt(10/10)) power_of_test_f <- 1 - type_2_error_f cat("Power of the test for (b) is =", power_of_test_f, "if true mean is 0.1, var = 10\n") ``` #### g) Calculating the power for mu != 0 (var = 10). ```{r, comment = NA} critical_g <- qnorm(0.05 / 2, mean = 0, sd=sqrt(10 / 10)) type_2_error_g <- pnorm(abs(critical_g), mean = 0.1, sd = sqrt(10 / 10)) - pnorm(-abs(critical_g), mean = 0.1, sd = sqrt(10 / 10)) power_of_test_g <- 1 - type_2_error_g cat("Power of the test for (b) is =", power_of_test_g, "if true mean is 0.1, var = 10 \n") ``` #### h) Power analysis using the function in (c). ```{r, comment = NA} sampling_function(number_of_samples = 1000, sample_size = 10, population_mean = 0.1, population_var = 1, c_a = qnorm(0.05, 0, 1/sqrt(10), lower.tail = FALSE), test_direction = c("one_sided")) sampling_function(number_of_samples = 1000, sample_size = 10, population_mean = 0.1, population_var = 1, c_a = qnorm(0.05/2, 0, 1/sqrt(10), lower.tail = FALSE), test_direction = c("two_sided")) likelihood_normal <- function(input_mean, input_sd, x){ temp = 1 for ( i in x){ temp = temp * (1/ sqrt( 2 * pi * input_sd^2)) * exp( - ((i - input_mean) ^ 2) / (2 * input_sd^2)) } return(temp) } ``` #### i) ```{r, comment = NA} likelihood_normal <- function(input_mean, input_sd, x){ temp = 1 for ( i in x){ temp = temp * (1/ sqrt( 2 * pi * input_sd^2)) * exp( - ((i - input_mean) ^ 2) / (2 * input_sd^2)) } return(temp) } x <- c(1.19, 0.33, 0.19, -1.89, 0.49, 0.08, 0.98, 2.92, 1.31, 1.39) sd_mle <- sqrt(mean((x-mean(x))^2)) like_null <- likelihood_normal(input_mean = 0, input_sd = sd_mle, x = x) like_alt <- likelihood_normal(input_mean = mean(x), input_sd = sd_mle, x = x) like_ratio <- -2 * log(like_null / like_alt) cat("Likelihood Ratio Test Statistic = ", like_ratio, "\n") pval <- pchisq(like_ratio, 1, lower.tail = FALSE) cat("P-value = ", pval, "\n") ```