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