Error bars in plots

A common irritation in making plots with R is the difficulty one has with superimposing error bars, or more precise, confidence intervals (c.i.'s) on plots. It can be actually quite simple to do, although you keep forgetting how if you don't often use it (like I do). It is most convenient to define the function superpose.eb which accepts the x and y coordinates of the points around which the confidence intervals should be marked, and a lower bound and upper bound which indicate how much lower/higher than y the lower/higher extreme of the c.i. reaches. By default a symmetric c.i. is assumed, so that the upper bound parameter is optional:

superpose.eb <-
function (x, y, ebl, ebu = ebl, length = 0.08, ...)
    arrows(x, y + ebu, x, y - ebl, angle = 90, code = 3,
    length = length, ...)

Its use is demonstrated in a barplot example (the RT data were kindly provided by Mariëtte Huizinga):

RT = matrix(c(814, 967, 500, 574, 424, 472, 394, 446), 2, 4) # data
colnames(RT) = c("7", "11", "15", "21")
rownames(RT) = c("repetition", "alternation")
eblb = matrix(c(14,21,12,18,12,18,13,19),2,4) # 1.96 * s.d. of data

x.abscis <- barplot(RT, beside=TRUE, col=0:1, ylim=c(0,1200),
    main="RT as a function of Age with 95%-confidance bars",
    xlab="Age (yrs)")
superpose.eb(x.abscis, RT, eblb, col="orange", lwd=2)

The function barplot returns the abscissa at which the bars are plotted, and these are used as the x coordinates passed to superpose.eb. The result lookes like this:

Farrel Buchinsky said... How does your function definition compare to the built in function, "errbar"

The function errbar is not realy built-in; its a function from the Hmisc package... I don't know all the differences, but one difference appears to be that in errbar you always have to fully compute the confidence interval by yourself, as in set.seed(1)
x <- 1:10
y <- x + rnorm(10)
delta <- runif(10)
errbar( x, y, y + delta, y - delta ) # you have to add and subtract delta from y by yourself
while superpose.eb can be called providing only means and associated s.d.'s (for asymptotic 68% c.i.'s), which is more intuitive for the people that ask me "how-to?" things like this (they apparently find it hard to think about how to calculate the locations for the lines to be drawn -- don't ask me why!). A more important difference is that errbar generates a plot if you do not set the add=TRUE parameter. Furthermore, superpose.eb has much much much less code than errbar (even if you would strip the optional plotting). Last but not least: setting col="orange" actually works in superpose.eb, as opposed to errbar.

Random.org - True Random Number Service

You probably know that computer generated numbers are not truely random, not even those generated by R. The random numbers generated by R are usualy very useful for all kinds of statistical purposes. However, they are not suitable in certain cases. For example for .... eh... Well, for example for conducting 'Psi' experiment (e.g., testing someones clairvoyance abilities). Another use is strong cryptography. If you wish to have true random numbers, you have to resort to other devices than computers. Or you can get some numbers from random.org, which are generated by constantly sampling the atmosphere. Other online random number generators are LavaRnd and

As an example, if you are using MSIE click here for 100 recent random integers ranging from 1 to 100, otherwise click here.

The block transpose of a matrix

Sometimes it's neat to be able to transpose the subblocks of a matrix: The "block transpose" of an n·k × m matrix can be defined as

Y =

⌈ Y1 ⌉ | ⋮ | ⌊ Yk ⌋

→

YBn

⌈ Y1' ⌉ | ⋮ | ⌊ Yk' ⌋

where the Y_i's are n × m matrices. The following routine performs this function in R:

# the function bt computes the block transpose
bt <-
function(y, n) {
  if(nrow(y)%%n!=0)
     stop(n,' is not a valid block size for this matrix.')
  m = ncol(y)
  k = nrow(y) %/% n
  matrix( aperm(array(t(y), c(m,n,k)), c(2,1,3)), , n, byrow=TRUE)
}

An example:

> idx = expand.grid(1:2,1:3,letters[1:4])
> m = matrix(paste(idx[,3],idx[,2],idx[,1],sep=''),nc=2,byrow=TRUE)
> m
      [,1]  [,2]
 [1,] "a11" "a12"
 [2,] "a21" "a22"
 [3,] "a31" "a32"
 [4,] "b11" "b12"
 [5,] "b21" "b22"
 [6,] "b31" "b32"
 [7,] "c11" "c12"
 [8,] "c21" "c22"
 [9,] "c31" "c32"
[10,] "d11" "d12"
[11,] "d21" "d22"
[12,] "d31" "d32"
> bt(m,3)
     [,1]  [,2]  [,3]
[1,] "a11" "a21" "a31"
[2,] "a12" "a22" "a32"
[3,] "b11" "b21" "b31"
[4,] "b12" "b22" "b32"
[5,] "c11" "c21" "c31"
[6,] "c12" "c22" "c32"
[7,] "d11" "d21" "d31"

How to multiply a matrix by a Kronecker product more or less efficiently

Suppose you would like to calculate the product

(A ⊗ B)   U' where U is an m ·n×k matrix, A is m×m, and B is n×n. Calculating the Kronecker product (A⊗ B) first and then multiplying it by U is feasible for small values of m and n, but computationally inefficient. For large values of m or n however, it is not only computationally inefficient, it is also not feasible memory wise: suppose m=151 and n=100, then (A⊗ B) represents a 15,100 × 15,100 matrix. If one of them is of type double, this requires more than 1739MB or 1.7GB of storage. We better do it as follows: If A is an identity matrix, it is easy to see that every m×k sub-block of U' is multiplied by B, so we better multiply each block separately. We can do a similar thing for A by rearanging the elements of U using the Commutation matrix operation (see Magnus & Neudecker, 1999). (One could also reason that each row u_i of U is the vec of a matrix U_i and use the relation vec(BU_iA') = (A⊗ B)vec(U_i).)

To calculate Y = (A⊗ B)U', we use the fact that (A⊗ B)=(A⊗ I_m)(I_n⊗ B), and compute (A⊗ I_m)Û and (I_n⊗ B)U' separately:

Û = B %*% matrix(t(U), m)   # m x (n . k)
Ŭ = A %*% matrix(t(Û), n)   # n x (k . m)
Y = matrix(t(Ŭ), k)         # k x (m . n)

That's it. The following simple and non-optimized function implements this:

function(A, B, U)
{
  k = nrow(U)
  if(missing(A)) A = diag(ncol(U)%/%ncol(B))
  if(missing(B)) B = diag(ncol(U)%/%ncol(A))
  U = B %*% matrix(t(U), ncol(B))
  U = A %*% matrix(t(U), ncol(A))
  U = matrix(t(U), k)
  U
}

It accepts A and B, one of which may be missing (in which case it is taken equal to the identity matrix of suitable order), and U. It returns Y'.

The following routine generalizes this to the product (A₁ &otimes A₂ ⊗ ··· ⊗ A_n)U'. It is used as kronecker.U(A1,A2, …, An, t(U)). It is not time efficient because the implementation is recursive.

function(A, B, U)

{

kronecker.U <- function(...)
{ 
 M = list(...);
 if(length(M)==2){
  return(M[[1]] %*% M[[2]]);
 }
 else{
  n = length(M);
  E = M[[n]];
  A = M[[n-1]];
  B = M[1:(n-2)];
  nce = ncol(E);
  nca = ncol(A);
  nra = nrow(A);
  ncb = prod(sapply(B,ncol));
  if(ncb*nca!=nrow(E))
   stop("non-conformable arguments");
  V = matrix( t(A%*%matrix(E, nca)), ncb);
  cl = match.call();
  cl[[n+1]] = V;
  cl[2:(n-1)] = B;
  cl[[n]] = NULL;
  W = eval(cl);
  Z = matrix(t(matrix(W,,nra)),,nce);
  return(Z);
 }
}