Cochran.Armitage <- function(n.i, r.i, x.i=F)
{
k <- length(n.i)
n <- sum(n.i)
r <- sum(r.i)
p.i <- r.i/n.i
p.m <- r/n
if (x.i == F) {
x.i <- 1:k
}
xx <- x.i-(x.m <- sum(n.i*x.i)/n)
b <- sum(n.i*(p.i-p.m)*xx)/sum(n.i*xx^2)
a <- p.m-b*x.m
xt <- b^2*sum(n.i*xx^2)/(p.m*(1-p.m))
xh <- n^2*(sum(r.i^2/n.i)-r^2/n)/r/(n-r)
xq <- xh-xt
res <- matrix(c(xt, xq, xh, 1, k-2, k-1, 1-pchisq(xt,1), 1-pchisq(xq,k-2), 1-pchisq(xh,k-1)), nrow=3)
colnames(res) <- c("Chi-sq.", "d.f.", "P")
rownames(res) <- c("Trend", "Quad.", "Homo.")
res
}
x.i <- c(10, 20, 30, 40, 50)
n.i <- c(30, 35, 47,21, 45)
r.i <- c(2, 4, 14, 13, 39)
Cochran.Armitage(n.i, r.i, x.i)
結果
Chi-sq. d.f. P
Trend 68.572732 1 1.110223e-16
Quad. 4.014484 3 2.599043e-01
Homo. 72.587216 4 6.439294e-15