# 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