## fisher.test
X <- rbind(c( 7,  29),
           c(13, 208))
fisher.test(X)

## By hand
M <- sum(X[,1])
N <- sum(X)
n <- sum(X[1,])
x <- 0:20
Prob <- choose(M,x)*choose(N-M,n-x)/choose(N,n)
names(Prob) <- x
barplot(Prob, border="white", las=1)
sum(Prob[Prob <= Prob["7"]])

## Approximate interval
SE <- sqrt(sum(1/X))
OR <- X[1,1]*X[2,2]/(X[1,2]*X[2,1])
log(OR) + qnorm(c(0.025, 0.975))*SE      ## CI for log(OR)
exp(log(OR) + qnorm(c(0.025, 0.975))*SE) ## CI for OR

## Confidence interval agreement
Y <- rbind(c(4475, 65),
           c(1597, 31))
fisher.test(Y)$conf.int
SE <- sqrt(sum(1/Y))
OR <- Y[1,1]*Y[2,2]/(Y[1,2]*Y[2,1])
exp(log(OR) + qnorm(c(0.025, 0.975))*SE)

## chisq.test
chisq.test(X)
chisq.test(X, correct=FALSE)

E <- outer(apply(X, 1, sum), apply(X, 2, sum))/sum(X) ## Expected cell counts
sum((X-E)^2/E)
1-pchisq(sum((X-E)^2/E), 1)

x <- seq(0, 10, 0.01)
plot(x, dchisq(x, 1), type="l")
abline(v=sum((X-E)^2/E), col="red")

## Test agreement
fisher.test(Y)$p.value
chisq.test(Y)
chisq.test(Y, correct=FALSE)

## Party identification data
Z <- as.table(rbind(c(762, 327, 468), c(484, 239, 477)))
dimnames(Z) <- list(gender = c("Female","Male"),
                    party = c("Democrat","Independent", "Republican"))
Z

## chisq.test: Two-by-three table
chisq.test(Z, correct=FALSE)
E <- outer(apply(Z, 1, sum), apply(Z, 2, sum))/sum(Z)
1-pchisq(sum((Z-E)^2/E), 2)

## fisher.test: Two-by-three table
fisher.test(Z)

## Odds ratios: Two-by-three table
fisher.test(Z[,-2])
fisher.test(Z[,-3])

## Harmine data
x <- c(16, 21, 22)     ## Control
y <- c(22, 23, 33, 37) ## 100 micrograms of harmine
X <- sum(x)
Y <- sum(y)

## poisson.test: Two sample
poisson.test(c(Y, X), c(4, 3))

## Comparison: t-test
t.test(x, y)

## poisson.test: One sample
poisson.test(Y, 4)
poisson.test(X, 3)

## Equivalent to
f.u <- function(lam) ppois(X, 3*lam) - 0.025
uniroot(f.u, c(20,30))$root
f.l <- function(lam) 1-ppois(X-1, 3*lam) - 0.025
uniroot(f.l, c(1, 20))$root

