Determining the n sample size to reject null at alpha =0.01 in R script

user_14390594 0
Viewed 60

Create a simulated dataset of 100 observations, where x is a random normal variable with mean 0 and standard deviation 1, and y = 0.1 + 2 * X + e, where epsilon is also a random normal error with mean 0 and sd 1.

set.seed(1)
# simulate a data set of 100 observations 
x <- rnorm(100)
y.1 <- 0.1 + 2*x + rnorm(100)

Now extract the first 5 observations.

y1.FirstFive <- (y.1[1:5]) # extract first 5 observations from y 
x.FirstFive <- (x[1:5]) # extract first 5 observations from x 


y1.FirstFive # extracted 5 observations from y1
[1] -1.7732743  0.5094025 -2.4821789  3.4485904  0.1044309
x.FirstFive # extracted 5 observations from x 
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

Assuming the mean and sd of the sample that you calculated from the first five observations would not change, what is the minimum total number of additional observations you would need to be able to conclude that the true mean of the population is different from 0 at the p = 0.01 confidence level?

alpha <- 0.01
mu <- 0

for (i in 5:2000) {
# Recalculate the standard error and CI
    stand_err <- Sd_y1 / sqrt(i)
    ci <- sample_mean_y1 + c(qt(alpha/2, i-1), qt(1-alpha/2, i-1))*stand_err
    if (ci[2] < mu)
        break     # condition met, exit loop
} 
i
[1] 2000

Here, I wrote a loop that iteratively increases n from the initial n=5 to n=2000, uses pt to find the p value (given a fixed y-bar and sd), and stops when p < 0.01. However I keep getting the wrong output. Such that, the output is always the number of the maximum range that I give (here, it is 2000) instead of giving me the specific minimum n sample in order to reject the null that mu_y = 0 at the p=0.01 level. Any suggestions as to how to fix the code?

  • additional info: the sd of y1.FirstFive = 2.3 and mean of y1.FirstFive = -0.04
statistics r
2 Answers

Assuming:

Sd_y1 = sd(y1.FirstFive)
sample_mean_y1 = mean(y1.FirstFive)
sample_mean_y1
[1] -0.03860587

As pointed out by @jblood94, you need to go for larger sample size.

You don't need a for loop for this, most of your functions are vectorized, so something like this:

n = 5:30000
stand_err = Sd_y1 / sqrt(n)
ub = sample_mean_y1 + qt(1-alpha/2, n-1)*stand_err
n[min(which(ub<0))]
[1] 23889
user_122584590

It's because n > 2000.

set.seed(1)
x <- rnorm(100)
y.1 <- 0.1 + 2*x + rnorm(100)

Sd_y1 <- sd(y.1[1:5])
sample_mean_y1 <- mean(y.1[1:5])
alpha <- 0.01
sgn <- 2*(sample_mean_y1 > 0) - 1

f <- function(n) qt(alpha/2, n - 1)*Sd_y1 + sgn*sample_mean_y1*sqrt(n)

upper <- 2
while (f(upper) < 0) upper <- upper*2
(n <- ceiling(uniroot(f, lower = upper/2, upper = upper, tol = 0.5)$root))
#> [1] 23889
user_94634890
Related