1james.stein <- function(y, group)
2{
3  s <- !(is.na(y)|is.na(group))
4  y <- y[s];
5  group <- as.character(group[s])
6  ## as.char -> unused levels OK
7  k <- length(unique(group))
8  if(k<3)
9    stop("must have >=3 groups")
10
11  stats <- function(w) {
12    bar <- mean(w)
13    ss  <- sum((w-bar)^2)
14    n <- length(w)
15    ##if(n<2)
16    ##  stop("a group has n<2")
17
18    c(n=length(w), mean=bar, ss=ss, var=ss/n/(n-1))
19  }
20
21  Z <- stats(y)
22  st <- tapply(y, group, FUN=stats)
23  nams <- names(st)
24  z <- matrix(unlist(st),ncol=4,byrow=TRUE)
25  ssb <- stats(z[,2])["ss"]
26  shrink <- 1 - (k-3)*z[,4]/ssb
27  shrink[z[,1]==1] <- 0
28  shrink <- pmin(pmax(shrink,0),1)
29  list(n=z[,1], mean=z[,2],
30       shrunk.mean=structure(Z["mean"]*(1-shrink)+shrink*z[,2], names=nams),
31       shrink=shrink)
32}
33