Subtracting columns in a matrix created by user-defined function

Question

I'm having some trouble subtracting two columns generated by my user-defined function in R.

My data:

dat <- data.frame(matrix(rnorm(1000), nrow=50))  

My function:

CI <- function(Int, dat){
  alpha = (1-Int) / 2
  z <- qnorm(1-alpha)
  N <- sum(is.na(dat)==FALSE)
  av <- colMeans(dat, na.rm=TRUE)                
  me <- z*(N/sqrt(N))
  lower <- av - me
  upper <- av + me

  return(cbind(lower, upper))
}

The function is trying to manually determine the upper and lower part of a Confidence Interval for a column of values (part of a bigger scheme). What I then want to do is subtract the upper from the lower CI to get the range. This is my attemp:

CI(Int=0.99, dat=dat)
CIdat <- CI(Int=0.99, dat=dat)
up <-as.numeric(CIdat[,2])
low <-as.numeric(CIdat[,1])
up+low
up-low

The up + low works, but the up - low gives a uniform result:

> up+low
     [1] -0.24418497  0.26253152 -0.35225948  0.27564574 -0.10129327 -0.14313671  0.11122208  0.30469330 -0.09498985
    [10] -0.05658091 -0.19336386 -0.31854479 -0.05433397 -0.03539858  0.07648641 -0.46213684  0.66736452  0.42824967
    [19]  0.29344554 -0.06257925
> up-low
     [1] 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097
    [13] 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097 162.9097

Am I missing something obvious?

Answer 1

Have a look at up and low. The values are far too large for a standard-normal distribution.

mean(up)      mean(low)
[1] 81.47071  [1] -81.43904

You probably have a thinking error somewhere in your CI() function. If I understand it correctly, you have a multidimensional standard-normal distribution dat and want to get the span of the confidence interval for each dimension (column) of dat.

I would suggest doing it as follows:

CI_span <- function(Int, dat) {

    alpha <- (1 - Int) / 2
    p_range <- c(alpha, 1 - alpha)
    apply(dat, 2, function(iColDat) as.numeric(diff(quantile(iColDat, probs = p_range))))

}

This returns the spans of each columns CI. R has a built in function for computing empirical quantiles called quantile.

Hth, D