STA 506 Assignment 9: Fundamentals of Hypothesis Testing and Power Analysis

1 Introduction

For this assignment, we will be using various statistical methods to compare two drugs and their effects on Low-Density Lipoprotein (LDL): the current standard drug, which lowers LDL by 25 mg/dL with a standard deviation of 15 mg/dL, and a drug a pharmeceutical company develops called “CholestFix,” which a 3 month clinical trial of 5 participants determined that this drug lowers LDL by 29 mg/dL on average. We will assume that this drug also has a standard deviation of 15 mg/dL, and that both LDL changes follow a normal distribution.

2 Part A

Based on the results of the clinical trial, researchers believe that CholestFix is the more effective drug. For this part, we will formally test that hypothesis with a z-test. The reason we use a z-test is because our data follows a normal distribution and our “population” (the standard drug effect) standard deviation is known. The z-test will be conducted at significance level \(\alpha\) = 0.05.

The null hypothesis, \(H_0\), is that the mean LDL decrease from the two drugs is equal. In other words, there is no difference in the two drugs’ effects. the alternative hypothesis, \(H_a\), is that there is a difference in the two drugs effects, and that CholestFix is the more effective drug. In other words, the average LDL decrease from CholestFix is greater than the LDL decrease from the standard drug. To put into symbols, we have:

\[H_0:\mu_C = \mu_S\] \[H_a: \mu_C \gt \mu_S\]

Another way we can write the hypotheses:

\[H_0: \mu_C - \mu_S = 0\] \[H_a: \mu_C - \mu_S \gt 0\]

Since we have a directional hypothesis and our test uses the normal distribution, our test will be a one-tailed test using the right tail. The test statistic formula will be:

\[z = \frac{\bar{x} - \mu_0}{\sigma / \sqrt{n}}\]

Our test is as follows:

## Perform Z-Test
mu_c <- 29 # mean of CholestFix
mu_s <- 25 # mean of standard drug 
sd <- 15 # standard deviation of both drugs
n <- 5 # number of participants in clinical trial


z <- (mu_c - mu_s)/(sd/sqrt(n))

pval <- 1 - pnorm(z) # Probability in right tail

test <- c(z, pval) # Z score and p-value
test

[1] 0.5962848 0.2754925

The z score for the test is 0.5962 and our p-value is 0.2755. Since our p-value is greater than 0.05, we will fail to reject our null hypothesis, and we can conclude that there is not enough evidence to suggest that there is a difference between the effectiveness of the standard drug and the effectiveness of CholestFix, and we can attribute our results to randomness.

3 Part B

For this part, we are interested in calculating the probability of a true improvement given some new numbers. We are now interested in seeing if we can get a true improvement given a new sample size of 50 and an effect size of \(\delta\) = 4 mg/dL. We are assuming no changes to either drug for this part.

It looks like we’re interested in calculating the post-hoc power of the test. To do so, we will need several things: first is the noncentrality parameter, \(\lambda\). The non-centrality parameter is: \[\lambda = \frac{\mu - \mu_0}{\sigma / \sqrt{n}}\] which equals: \[\lambda = \frac{\delta}{\sigma / \sqrt{n}}\]

Pluggin in what we want, we get:

\[\lambda = \frac{4}{15 / \sqrt{50}}\] which gives us \[\lambda = 1.8856 \ldots\].

delta <- 4 # Effect Size
sd <- 15 # Standard Deviation
n <- 50 # New Sample Size

lambda <- delta/(sd/sqrt(n))
lambda

[1] 1.885618

Then we need our critical value corresponding to the significance level:

crit.alpha <- qnorm(0.95, 0, 1) # Crit Value of 95% of Normal Distribution
crit.alpha

[1] 1.644854

Now we can use the formula to find the power of the test:

## Calculating Power for Effect Size of 4 mg/dL

# Power calculation
power <- 1 - pnorm(crit.alpha - lambda)
power

[1] 0.5951312

The probability that we detect a true improvement, or the power of our test, is 0.5951. This means that we can correctly reject the null hypothesis when the mean is 29 mg/dL about 59.51% of the time.

4 Part C

For this part, we are interested in calculating the a-priori minimum sample size to detect an effect size of \(\delta\) = 4 mg/dL, with power of 0.8 and significance level of 0.05.

The formula for estimating minimal sample size is \[n \approx \left( \frac{(z_{1-\alpha} + z_{1-\beta}) \cdot \sigma}{\delta} \right)^2\]

We already have our z score for 1 - \(\alpha\), which is 1.645. Now we must get our z score corresponding to 1-\(\beta\):

beta <- 0.2 # 1 - 0.8

crit.beta <- qnorm(1 - beta) # crit z value for 1 - Beta
crit.beta

[1] 0.8416212

Plugging in what we have now, we can estimate our necessary sample size:

## Estimating Sample Size for Certain Conditions

# Sample Size Estimate Formula
sample.size <- (((crit.alpha + crit.beta)*sd)/delta)^2 
sample.size

[1] 86.94221

Our estimated a-priori minimum sample size to achieve an effect size of 4 mg/dL with power of 0.8 and significance level 0.05 is 87.

What’s interesting is how small the original sample size and the new theoretical sample size is compared to the minimum sample size needed for to get the power we desire. If we use the original sample size of 5, the power is only 0.1472. It seems that increasing sample size increases the power.

4.1 Part D

For this part, we will examine how power changes as our sample size changes. To do this, we will create a function that uses different values of the sample size and calculates the corresponding power. We will then plot those values and see what happens. We will use the same standard deviation, significance level, and effect size from before.

## Check Various Sample Sizes and Power

sd <- 15
delta <- 4
alpha = 0.05

# Crit Value: Same as Before
crit.alpha <- qnorm(1- alpha)

# Sequence of Sample Sizes
samp.sizes <- seq(5, 200, by = 5)

# Function to Compute Power
val.power <- sapply(samp.sizes, function(n) {
  lambda <- delta/(sd / sqrt(n))  # non-centrality parameter
  1 - pnorm(crit.alpha - lambda) # calculating power
})

# Plot Power vs Sample Size
plot(samp.sizes, val.power, type = "l", ylim = c(0,1), main = "Plot of Power vs Sample Size for Effect Size of 4 mg/dL", xlab = "Sample Size", ylab = "Power")

# Plot Reference Line for Power = 0.8
abline(h = 0.8, lty = 2)

We can see that we have a somewhat steady increase in power until about 150 or so, when the increase starts to level off. We can also see that our sample size calculation from part C holds up, as the graph crosses the 0.8 mark when the sample size is around at around 85-87.