if (FALSE) {
library("HelpersMG")
# Symmetry of Lepidochelys olivacea scutes
table <- t(data.frame(SriLanka=c(200, 157), AfricaAtl=c(19, 12),
Guyana=c(8, 6), Suriname=c(162, 88),
MexicoPac1984=c(42, 34), MexicoPac2014Dead=c(8, 9),
MexicoPac2014Alive=c(13, 12),
row.names =c("Symmetric", "Asymmetric")))
table
contingencyTable.compare(table)
table <- t(data.frame(SriLanka=c(200, 157), AfricaAtl=c(19, 12), Guyana=c(8, 6),
Suriname=c(162, 88), MexicoPac1984=c(42, 34),
MexicoPac2014Dead=c(8, 9),
MexicoPac2014Alive=c(13, 12), Lepidochelys.kempii=c(99, 1),
row.names =c("Symmetric", "Asymmetric")))
table
contingencyTable.compare(table)
# Conformity to a model
table <- matrix(c(33, 12, 25, 75), ncol = 2, byrow = TRUE)
probs <- c(0.5, 0.5)
contingencyTable.compare(table, probs=probs)
# Conformity to a model
table <- matrix(c(33, 12), ncol = 2, byrow = TRUE)
probs <- c(0.5, 0.5)
contingencyTable.compare(table, probs=probs)
# Conformity to a model
table <- matrix(c(33, 12, 8, 25, 75, 9), ncol = 3, byrow = TRUE)
probs <- c(0.8, 0.1, 0.1)
contingencyTable.compare(table, probs=probs)
# Comparison of chisq.test() and this function
table <- matrix(c(NA, NA, 25, 75), ncol = 2, byrow = TRUE)
pv <- NULL
aw <- NULL
par(new=FALSE)
n <- 100
for (GroupA in 0:n) {
table[1, 1] <- GroupA
table[1, 2] <- n-GroupA
pv <- c(pv, chisq.test(table)$p.value)
aw <- c(aw, contingencyTable.compare(table, criterion="BIC")[1])
}
x <- 0:n
y <- pv
y2 <- aw
plot(x=x, y=y, type="l", bty="n", las=1, xlab="Number of type P in Group B", ylab="Probability",
main="", lwd=2)
lines(x=x, y=y2, type="l", col="red", lwd=2)
# w-value
(l1 <- x[which(aw>0.05)[1]])
(l2 <- rev(x)[which(rev(aw)>0.05)[1]])
aw[l1]
pv[l1]
aw[l2+2]
pv[l2+2]
# p-value
l1 <- which(pv>0.05)[1]
l2 <- max(which(pv>0.05))
aw[l1]
pv[l1]
aw[l2]
pv[l2]
y[which(y2>0.05)[1]]
y[which(rev(y2)>0.05)[1]]
par(xpd=TRUE)
text(x=25, y=1.15, labels="Group A: 25 type P / 100", pos=1)
segments(x0=25, y0=0, x1=25, y1=1, lty=3)
# plot(1, 1)
v1 <- c(expression(italic("p")*"-value"), expression("after "*chi^2*"-test"))
v2 <- c(expression(italic("w")*"-value for A"), expression("and B identical models"))
legend("topright", legend=c(v1, v2),
y.intersp = 1,
col=c("black", "black", "red", "red"), bty="n", lty=c(1, 0, 1, 0))
segments(x0=0, x1=n, y0=0.05, y1=0.05, lty=2)
text(x=101, y=0.05, labels = "0.05", pos=4)
}
Run the code above in your browser using DataLab