Solution to Homework 6

Say we have the following: we have data \(X_1,..,X_n \sim N(0, 1)\) and \(Y_1,..,Y_m \sim N(\mu, 1)\). We wish to test

\[H_0:\mu=0\text{ vs. }H_a:\mu>0\]

In R this can be done with the command

t.test(x, y, alternative = "less")$p.value

For example

x <- rnorm(10)
y <- rnorm(13, 0)
t.test(x, y, alternative = "less")$p.value

## [1] 0.2919114

y <- rnorm(13, 2)
t.test(x, y, alternative = "less")$p.value

## [1] 0.0004251926

Problem 1

Verify that this test has the correct type I error for the case \(n=m=10\), \(\alpha=0.01, 0.05, 0.10\).

B <- 10000
pvals <- rep(0, B)
for(i in 1:B) {
  x <- rnorm(10)
  y <- rnorm(10)
  pvals[i] <- t.test(x, y, alternative = "less")$p.value
}
c(sum(pvals<0.01), sum(pvals<0.05), sum(pvals<0.1))/B

## [1] 0.0098 0.0510 0.0998

Problem 2

What is the power of this test if \(n=m=10\), \(\mu=1.1,\alpha=0.05\)?

pwr <- function(n, m, mu, alpha=0.05) {
  B <- 10000
  pvals <- rep(0, B)
  for(i in 1:B) {
    x <- rnorm(n)
    y <- rnorm(m, mu)
    pvals[i] <- t.test(x, y, alternative = "less")$p.value
  }
  round(100*sum(pvals<alpha)/B, 1)  
}
pwr(10, 10, 1.1, 0.05)

## [1] 76.3

Problem 3

How large does \(\mu\) have to be so that the power is \(95\%\) for the case \(n=m=10\) if the test is done with \(\alpha=0.05\)?

mu <- 1.1
repeat {
  mu <- mu+0.1
  if(pwr(10, 10, mu)>95) break
}
c(mu, pwr(10, 10, mu))

## [1]  1.6 96.4

Problem 4

How large does \(\mu\) have to be so that the power is \(95\%\) for the case \(n=m=10\) if the test is done with \(\alpha=0.01\)?

mu <- 1.6
repeat {
  mu <- mu+0.1
  if(pwr(10, 10, mu, 0.01)>95) break
}
c(mu, pwr(10, 10, mu, 0.01))

## [1]  2.0 96.2

Problem 5

What sample size \(n=m\) is needed so the test has a power of \(95\%\) if \(\mu=0.5\) and \(\alpha=0.05\)?

n <- 10
repeat {
  n <- n+10
  if(pwr(n, n, 0.5)>95) break
}
n <- n-10
repeat {
  n <- n+1
  if(pwr(n, n, 0.5)>95) break
}
c(n, pwr(n, n, 0.5))

## [1] 88.0 94.8

Problem 6

If \(n=10\), is there an \(m\) so that the test has a power of \(95\%\) if \(\mu=0.5\) and \(\alpha=0.05\)?

m <- 10*1:20
p <- 0*m
for(i in seq_along(m))
  p[i] <- pwr(10, m[i], 0.5)
plot(m, p)

So there does not seem to be such an m.