Why does R t.test()'s confidence interval shift slightly w.r.t. alternative hypothesis?

Question

So I've noticed something extremely confusing about t.test, which is that the confidence interval it provides w.r.t a two-sided alt hypothesis seems to be inconsistent with what it provides for a one-sided test, my own CI function & excel's CONFIDENCE() (which are consistent with each other). In addition it appears to be inconsistent with my own two-tailed z test function, and all of these inconsistencies are pretty small (but bigger than floating point error)!

Here is my example: t.test---

sample = c(38, 30, 41, 28, 31)
t.test(sample) # CI = 26.65334, 40.54666

t.test(sample, alternative = 'less') # CI = -Inf, 38.93388
t.test(sample, alternative = 'greater') # CI = 28.26612,  Inf`

My own CI---

# tested 10/5/19
# xbar = sample mean, s = sample sigma, n = sample size
mean_conf_interval <- function(xbar, s, n, conf_lvl=0.95) {
  z_a2 <- qnorm((1-conf_lvl)/2, lower.tail=F)
  err_margin <- z_a2*s/sqrt(n)
  return(c(xbar-err_margin, xbar+err_margin))
}

mean_conf_interval(mean(sample),sd(sample),length(sample))
# CI = 28.26612,  38.93388

excel--- (rough translation to pseudo code)

E276:E280 = 38, 30, 41, 28, 31 # same sample as before
B277 = CONFIDENCE(0.05,  STDEV(E276:E280), COUNT(E276:E280)) # error margin
B278=AVERAGE(E276:E280) # xbar
# CI = B278-B277, B278+B277 = 28.69617167,  38.50382833

For the sake of brevity I won't go into comparison with my own two-sided-z-test function unless necessary. But what's going on? This is so strange that it actually looks like a bug...

Answer 1

You are using the standard normal distribution to calculate the 0.975 quantile as 1.96. You should be using the $t_{n-1}$ distribution, which will give a value somewhat larger than 1.96.

Try z_a2 <- qt((1-conf_lvl)/2,n-1,lower.tail=F) in your function and see what you get.

Edit: I have run the following code, getting the same confidence interval that I have in my comment.

mean_conf_interval <- function(xbar, s, n, conf_lvl=0.95) {
  z_a2 <- qnorm((1-conf_lvl)/2, lower.tail=F)
  err_margin <- z_a2*s/sqrt(n)
  return(c(xbar-err_margin, xbar+err_margin))
}

sample <- c(38, 30, 41, 28, 31)
xbar <- mean(sample)
s <- sd(sample)
n <- length(sample)
mean_conf_interval(xbar,s,n)