Bootstrap confidence interval in R using 'replicate' and 'quantile' I have around one thousand measurements (numbers). All these measurements are my observations and I have calculated a 95% confidence interval for the mean and for the variance by using the normal formulas (without software).
Then, I have used the replicate function in R with around one hundred thousand simulations for my measurements (observations) and the parameter "replace" has been set to "true".
After that, I used the apply function with the three parameters: my measurements as the data, 2 as the margin and mean as the function.
Then to get a 95% confidence interval this way, I used the quantile function containing the variable for the apply function I used, and then 0.025 and 0.975 combined as the second parameter for the quantile function.
In that way, I got almost exactly the same 95% confidens interval as calculated with the normal formula (without software).
So now I wanted to do exactly the same thing for the variance, i.e. use replicate, apply and quantile to get a 95% confidence interval for the variance. So I just changed the third parameter "mean" to "var" in the apply function. I then noticed that the outputted 95% confidence interval for the variance (from the quantile function) is a little bit different (it is wider) than the one I calculated by the normal formula without software.
So my question is:
Did I use the replicate, apply and quantile function correctly for the variance confidence interval? I know I did use the functions right for the mean, since I got almost exactly the same result there as calculated by normal formula.
If I did use the functions correctly for the variance, why is the confidence interval a little bit different? Is it because as sample size increases, there might be more observations far away from the mean resulting in a bigger variance?
 A: From your description, I am not sure exactly what you are doing.
But maybe I can give you some helpful ideas.
1) For normal data, the standard CI for the population mean $\mu$ uses
Student's t distribution which is symmetrical. The 95% CI is of the
form $\bar X \pm t^*S/\sqrt{n},$ where $t^*$ cuts 2.5% from the upper
tail of $\mathsf{T}(\nu= n-1)$.
By contrast, the standard CI for the population variance $\sigma^2$ uses
the chi-squared distribution which is not symmetrical. The 95% CI is of
the form $\left(\frac{(n-1)S^2}{U}, \frac{(n-2)S^2}{L}\right),$ where
$L$ and $U$ cut 2.5% of the probability from the lower and upper tails,
respectively, of $\mathsf{Chisq}(\nu = n-1).$
Unless your bootstrapping procedure corrects for the bias in the case
of the variance (with its skewed distribition), you cannot expect an accurate result.
2) It is not exactly 'fair' to compare CIs (using Student's t and
chi-squared distributions) based on the assumption data are normal
with results of nonparametric bootstrap CIs. The assumption that
data are normal provides 'information' for the t and chi-squared
intervals that is not used in nonparametric bootstrapping.
3) In making bootstrap CIs for variation, I have found it better
to find CIs for $\sigma$ rather than $\sigma^2,$ possibly because the former
has the same units as the data.  Also, for scale parameters, I have found
that it often works better to bootstrap ratios rather than differences.
4) Finally, along lines of @HansEngler's suggestion, I will generate $n = 1000$ observations from $\mathsf{Norm}(\mu = 100, \sigma=15),$ find the chi-sqared CI for $\sigma$ and compare it with
a bias-corrected nonparametric bootstrap. 
First, here are my fake data and the corresponding
standard CI for $\sigma,$  which turns out to be $(9.55, 10.43).$
[I have provided set.seed statements, so that you can replicate
the exact simulations I have used.]
set.seed(1234)
n = 1000;  mu = 100;  sg = 10;  x = rnorm(n, mu, sg)
a = mean(x); s = sd(x);  a;  s
## 99.73403    # sample mean
## 9.973377    # sample SD
# CI for sg
UL = qchisq(c(.975,.025), 999); UL
## 1088.487  913.301
sqrt((n-1)*var(x)/UL)
##  9.554619 10.430810  # CI for sg (includes 9.973, as it must)

Now for the (bias corrected) nonparametric bootstrap CI: I will bootstrap the ratio
$R = S/\sigma.$ If I knew the distribution of $R,$ then I could find
$L$ and $U$ with $P(L \le R = S/\sigma \le U) = 0.95$ so that 
$P(S/U \le \sigma \le S/L) = .95$ and a 95%
CI for $\sigma$ would be of the form $(S/U,\, S/L).$ By bootstrapping
$R,\,$ I can estimate $U$ and $L.$  I use
the observed $s = 9.97$ temporarily as a proxy for unknown $\sigma.$
Suffixes .re indicate bootstrapped quantities. [On this site
with relatively few experienced R users,
I try to use only the most fundamental R functions; I will leave it
to you to write more elegant R code, which is obviously possible.]  
set.seed(1235)
B = 10^5;  r = numeric(B)
for(i in 1:B) {
  s.re = sd(sample(x,n,repl=T))
  r[i] = s.re/s }
L.re = quantile(r, .025);  U.re = quantile(r, .975)
c(s/U.re, s/L.re)  
##    97.5%      2.5% 
## 9.537866 10.462202 

So the nonparametric bootstrap CI is $(9.54, 10.46),$ which
is not a bad match for the standard chi-squared CI $(9.55, 10.43)$
obtained above. [The large sample size ($n$ = 1000) has mostly
obviated my comment in (2) above, which still stands for smaller $n$.
Even so, it is possible that the information that the data are normal
accounts for the fact that the normal-based CI is a bit shorter.]
Note: If you do a parametric bootsrap CI, you may get
closer to the result for the chi-squared CI. In parametric
bootstrapping, one re-samples from a parametric distribution (here normal)
with parameters suggested by the data (rather than from the data themselves). 
In particular, in the code above one would substitute the code 
s.re = sd(rnorm(n, a, s)) for the line with sample. With this change
and bootstrap seed 1066, I got the 95% parametric bootstrap CI $(9.55, 10.43).$
