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)

# 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")
Critical value for a Type 1 error of 0.05 =  0.5201484 

b) Determine the critical value for alpha = 0.05. (mu != 0)

# 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")
Critical values for a Type 1 error of 0.05 =  -0.619795 0.619795 

c) Function for generating random samples and testing significance

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"))
Null Hypothesis is rejected in  55  out of  1000  cases ( 5.5 %) 
[1] 55
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"))
Null Hypothesis is rejected in  49  out of  1000  cases ( 4.9 %) 
[1] 49

d) Calculating the power for mu > 0.

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") 
Power of the test for (a) is = 0.09198572 if true mean is 0.1 

e) Calculating the power for mu != 0.

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") 
Power of the test for (b) is = 0.06153262 if true mean is 0.1 

f) Calculating the power for mu > 0 (var = 10).

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") 
Power of the test for (b) is = 0.06119084 if true mean is 0.1, var = 10

g) Calculating the power for mu != 0 (var = 10).

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") 
Power of the test for (b) is = 0.0511463 if true mean is 0.1, var = 10 

h) Power analysis using the function in (c).

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"))
Null Hypothesis is rejected in  83  out of  1000  cases ( 8.3 %) 
[1] 83
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"))
Null Hypothesis is rejected in  65  out of  1000  cases ( 6.5 %) 
[1] 65
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)

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")
Likelihood Ratio Test Statistic =  3.585091 
pval <- pchisq(like_ratio, 1, lower.tail = FALSE)

cat("P-value = ", pval, "\n")
P-value =  0.05830021