https://github.com/cran/Hmisc
Tip revision: 4edd006bd91028b616abe3761e300f0f4d5d4ad7 authored by Charles Dupont on 26 December 2008, 00:00:00 UTC
version 3.8-0
version 3.8-0
Tip revision: 4edd006
t.test.cluster.s
t.test.cluster <- function(y, cluster, group, conf.int=.95)
{
## See:
## Donner A, Birkett N, Buck C, Am J Epi 114:906-914, 1981.
## Donner A, Klar N, J Clin Epi 49:435-439, 1996.
## Hsieh FY, Stat in Med 8:1195-1201, 1988.
group <- as.factor(group)
cluster <- as.factor(cluster)
s <- !(is.na(y)|is.na(cluster)|is.na(group))
y <- y[s];
cluster <- cluster[s];
group <- group[s]
n <- length(y)
if(n<2)
stop("n<2")
gr <- levels(group)
if(length(gr)!=2)
stop("must have exactly two treatment groups")
n <- table(group)
nc <- tapply(cluster, group, function(x)length(unique(x)))
bar <- tapply(y, group, mean)
u <- oldUnclass(group)
y1 <- y[u==1];
y2 <- y[u==2]
c1 <- factor(cluster[u==1]);
c2 <- factor(cluster[u==2]) #factor rids unused lev
b1 <- tapply(y1, c1, mean);
b2 <- tapply(y2, c2, mean)
m1 <- table(c1);
m2 <- table(c2)
if(any(names(m1)!=names(b1)))
stop("logic error 1")
if(any(names(m2)!=names(b2)))
stop("logic error 2")
if(any(m2 < 2))
stop(paste('The following clusters contain only one observation:',
paste(names(m2[m2 < 2]), collapse=' ')))
M1 <- mean(y1);
M2 <- mean(y2)
ssc1 <- sum(m1*((b1-M1)^2));
ssc2 <- sum(m2*((b2-M2)^2))
if(nc[1]!=length(m1))
stop("logic error 3")
if(nc[2]!=length(m2))
stop("logic error 4")
df.msc <- sum(nc)-2
msc <- (ssc1+ssc2)/df.msc
v1 <- tapply(y1,c1,var);
v2 <- tapply(y2,c2,var)
ssw1 <- sum((m1-1)*v1);
ssw2 <- sum((m2-1)*v2)
df.mse <- sum(n)-sum(nc)
mse <- (ssw1+ssw2)/df.mse
na <- (sum(n)-(sum(m1^2)/n[1]+sum(m2^2)/n[2]))/(sum(nc)-1)
rho <- (msc-mse)/(msc+(na-1)*mse)
r <- max(rho, 0)
C1 <- sum(m1*(1+(m1-1)*r))/n[1]
C2 <- sum(m2*(1+(m2-1)*r))/n[2]
v <- mse*(C1/n[1]+C2/n[2])
v.unadj <- mse*(1/n[1]+1/n[2])
de <- v/v.unadj
dif <- diff(bar)
se <- sqrt(v)
zcrit <- qnorm((1+conf.int)/2)
cl <- c(dif-zcrit*se, dif+zcrit*se)
z <- dif/se
P <- 2*pnorm(-abs(z))
stats <-
matrix(NA, nrow=20, ncol=2,
dimnames=list(c("N","Clusters","Mean",
"SS among clusters within groups",
"SS within clusters within groups",
"MS among clusters within groups","d.f.",
"MS within clusters within groups","d.f.",
"Na","Intracluster correlation",
"Variance Correction Factor","Variance of effect",
"Variance without cluster adjustment","Design Effect",
"Effect (Difference in Means)",
"S.E. of Effect",paste(format(conf.int),"Confidence limits"),
"Z Statistic","2-sided P Value"), gr))
stats[1,] <- n
stats[2,] <- nc
stats[3,] <- bar
stats[4,] <- c(ssc1, ssc2)
stats[5,] <- c(ssw1, ssw2)
stats[6,1] <- msc
stats[7,1] <- df.msc
stats[8,1] <- mse
stats[9,1] <- df.mse
stats[10,1] <- na
stats[11,1] <- rho
stats[12,] <- c(C1, C2)
stats[13,1] <- v
stats[14,1] <- v.unadj
stats[15,1] <- de
stats[16,1] <- dif
stats[17,1] <- se
stats[18,] <- cl
stats[19,1] <- z
stats[20,1] <- P
attr(stats,'class') <- "t.test.cluster"
stats
}
print.t.test.cluster <- function(x, digits, ...)
{
## if(!missing(digits)).Options$digits <- digits 6Aug00
if(!missing(digits)) {
oldopt <- options(digits=digits)
on.exit(options(oldopt))
}
cstats <- t(apply(x,1,format))
## cstats <- format(x)
attr(cstats,'class') <- NULL
cstats[is.na(x)] <- ""
invisible(print(cstats, quote=FALSE))
}