5
$\begingroup$

I am using R and plain English to express my question. Let us say I have a "true"/made up population, which is normally distributed with a mean of 500000 and a standard deviation of 13000:

mean <- 500000
sd <- 13000
population <- rnorm(10000000, mean=mean, sd=sd)

I can then repeatedly sample from this population using different samples sizes ranging from 3 to 150. For each sample, I also calculate the 95% confidence interval using a t distribution. I think this is safe as at the beginning my sample size is pretty low (i.e. 3). This is not really exploited here but I think this would be correct?

results <- NULL

for (sample_size_ in 3:150) {
    
    for (sample_number_ in 1:10) {
        
        sample_ <- sample(population, sample_size_)
        ci_test <- t.test(sample_, conf.level=0.95)

        df <- data.frame(
            sample_size = sample_size_
            , low_ci = ci_test$conf.int[[1]] 
            , up_ci = ci_test$conf.int[[2]]
            , mean = mean(sample_)
        )

        results <- rbind(results, df)
        
    }
}

I can then plot the data frame, which nicely shows how the sample means gets closer to the population mean (red dotted line) as the sample size increases:

ggplot(results, aes(y = mean, x = sample_size)) +
    geom_point() +
    geom_hline(yintercept=mean, linetype="dotted", color = "red", size=3)

This is the output (the actual code uses a theme):

enter image description here

I guess this is straightforward but how can I determine these blue lines expressing the 95% theoretical uncertainty of the sample mean being inside the blue area?

enter image description here

Is there a formula and does it depend on the knowledge of that the true distribution is normal? I guess it should not, as the sample means will be normally distributed according to the CLT?

Thanks.

PS:

This is the current/final code:

mu_ <- 500000
sd_ <- 13000
z_value <- qnorm(.975)
max_sample_size = 500
repeats = 2

results <- NULL
upper_and_lower <- NULL

for (sample_size_ in 3:max_sample_size) {
    
    for (sample_number_ in 1:repeats) {
        
        sample_ <- rnorm(sample_size_, mean=mu_, sd=sd_)
        ci_test <- t.test(sample_, conf.level=0.95)

        df <- data.frame(
            sample_size = sample_size_
            , low_ci = ci_test$conf.int[[1]] 
            , up_ci = ci_test$conf.int[[2]]
            , mean = mean(sample_)
        )

        results <- rbind(results, df)
        
    }
    
    df <- data.frame(
        sample_size = sample_size_
        , upper = mu_ + (z_value * (sd_/sqrt(sample_size_)))
        , lower = mu_ - (z_value * (sd_/sqrt(sample_size_)))
    )
    upper_and_lower <- rbind(upper_and_lower, df)
}

ggplot() +
    geom_point(data=results, aes(y = mean, x = sample_size)) +
    geom_line(data=upper_and_lower, aes(y = upper, x = sample_size), color = "blue", size=2) +
    geom_line(data=upper_and_lower, aes(y = lower, x = sample_size), color = "blue", size=2)

Resulting graph with theoretical 95% uncertainty band:

enter image description here

$\endgroup$
8
  • $\begingroup$ Even though you 'know' $\sigma = 13000$ (from simulation, you're pretending all you know for each t CI is $n, \bar X_n, S_n.$ Then 95% CIs from t.test in R will be about $50000 \pm 2 S_n/\sqrt{n}$ for $n > 30.$ For smooth curves, maybe use upper curve $50000 + 2(13000)/\sqrt{n}$ from $n = 30$ to $150.$ [Lower curve with $-.$] $\endgroup$
    – BruceET
    Commented Dec 27, 2020 at 17:03
  • 1
    $\begingroup$ Unrelated to the question, watch out when you define variables with names used by functions. When I do simulations like this and define mean and standard deviation parameters, I call them “mu” and “s” (should be “sigma”, but I’m content to call it “s”). $\endgroup$
    – Dave
    Commented Dec 27, 2020 at 17:41
  • 1
    $\begingroup$ The blues lines that you imagine are not confidence intervals. Confidence intervals are a function of the data/sample mean $\bar{X}$, data/sample deviation $S$ and sample size $n$ and a constant $c$ which depends on the desired confidence level$$\bar{X}_n\pm c \frac{S_n}{\sqrt{n}}$$Your plot demonstrates how the distribution of $\bar{X}_n$ has a variance that scales with $n^{-0.5}$, and you could make boundaries that contain some $x\%$ of the observations, but that is not the same concept as confidence intervals $\endgroup$ Commented Dec 29, 2020 at 16:28
  • 1
    $\begingroup$ @cs0815 The confidence interval is $\bar{X}_n\pm c \frac{S_n}{\sqrt{n}}$. It depends on the sample/observation and you can not plot it as values that are for some given $n$ fixed. (But maybe you are not looking for a confidence interval, and instead you look for a prediction interval or tolerance interval?) $\endgroup$ Commented Dec 29, 2020 at 17:16
  • 1
    $\begingroup$ @cs0815 $c$ is a parameter that depends on the percentage for the confidence region. It will also depend on $n$. See for more information the link that I referenced in the comment. $\endgroup$ Commented Jan 1, 2021 at 17:16

1 Answer 1

6
$\begingroup$

Nice experiment. The blue lines will be at $\mu \pm z_{\alpha/2} \sigma/\sqrt{n}$ where $\alpha = 0.05$ and $\alpha \mapsto z_\alpha$ is the upper quantile function of the standard normal and $n$ is sample_size_. In this case, this is exact since you are generating from the normal distribution and the sample mean will have the distribution $N(\mu, \sigma^2/n)$. In general, it holds approximately based on the CLT for large $n$.

If you are trying to interpret the CI, that would be a bit different. Putting a CI around each of the points vertically, for each sample_size_, roughly 95% of those intervals will cross the dotted red line (you cover the true mean $\approx$ 95% of the time, hence the coverage probability is 0.95). Since your replication size is small (10), this is not that accurate. Try increasing sample_number_ to say 1000 to see it better. What you will observe is that although the length of the CIs will be smaller as you increase $n$, still $\approx$ 95% of them keep covering the red line no matter what $n$ is.


PS. I am assuming that instead of these two lines:

population <- rnorm(10000000, mean=mean, sd=sd)
sample_ <- sample(population, sample_size_)

you would do something like this:

sample_ <- rnorm(sample_size_, mean=mean, sd=sd)

That is, sample from the normal distribution, not a finite population drawn from the normal distribution as you are doing here, although the results will be close in your case (unless your sample size starts approaching 10000000).

$\endgroup$
6
  • $\begingroup$ Thanks. Reg PS yes. I just pretended that I have a fixed population. $\endgroup$
    – cs0815
    Commented Dec 27, 2020 at 19:15
  • $\begingroup$ Is the standard deviation in your formula the population one? I guess so? $\endgroup$
    – cs0815
    Commented Dec 27, 2020 at 19:18
  • $\begingroup$ @cs0815, yes, $\mu$ and $\sigma$ are the population mean and s.d. $\endgroup$
    – passerby51
    Commented Dec 27, 2020 at 19:19
  • 1
    $\begingroup$ @cs0815, no problem. Yes, it looks good to me. $\endgroup$
    – passerby51
    Commented Dec 27, 2020 at 21:52
  • 1
    $\begingroup$ Thanks for your help $\endgroup$
    – cs0815
    Commented Dec 27, 2020 at 21:54

Not the answer you're looking for? Browse other questions tagged or ask your own question.