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