diff --git a/scripts/ce089-04.R b/scripts/ce089-04.R
index 80ccfcbb41e6c767e3e836e45621dad12f1e89bb..bed0879f3db6c79e9a5bbc359102d03cdc1ab88a 100644
--- a/scripts/ce089-04.R
+++ b/scripts/ce089-04.R
@@ -1,10 +1,52 @@
#-----------------------------------------------------------------------
-# Teste de aleatórização.
+# Teste de aleatorização.
+
+# Uma senhora toma chá.
+#
+# A senhora (Muriel Bristol) declarou saber distringuir entre bedidas
+# segundo a ordem com que o chá ou leite foram colocados na xícara. O
+# experimento consistiu em servir 8 xícaras da bebida, 4 em cada order
+# (chá/leite e leite/chá). No final a senhora deveria indicar quais as
+# xícaras que estão em cada um dos dois grupos. DETALHE: ela sabe
+# distinguir mas não sabe classificar.
+
+# Quantidade de maneiras de gerar dois grupos de 4 xícaras usando 8.
+n_poss <- choose(8, 4)
+n_poss
+
+# De 4 xícaras selecionadas para um grupo, X representa o número de
+# acertos.
+n_acertar <- sapply(4:0,
+ FUN = function(i) {
+ # (formas de acertar x em 4) * (formas de errar x em 4).
+ choose(4, i) * choose(4, 4 - i)
+ })
+n_acertar
+
+# Pr(Acerto total) = 1/70 < 5%.
+cumsum(n_acertar)/n_poss
+
+# Comprovando por simulação.
+v <- rep(0:1, each = 4)
+mean(replicate(n = 100000,
+ expr = all(sample(v) == v)))
+
+# Essa forma de fazer a conta é mais realista mas dá na mesma.
+mean(replicate(n = 100000,
+ expr = all(sample(v) == sample(v))))
+
+#-----------------------------------------------------------------------
+# Exemplo com teste para a diferença de médias.
# Comprimentos de crânios de cães pré históricos.
m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112)
f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111)
+plot(ecdf(m), xlim = range(c(m, f)), col = "cyan")
+lines(ecdf(f), col = "magenta")
+rug(m, col = "cyan")
+rug(f, col = "magenta")
+
# Diferença de média.
d <- mean(m) - mean(f)
d
@@ -12,20 +54,54 @@ d
# Aplicando um teste t.
t.test(x = m, y = f, var.equal = TRUE)
+# Todos as combinações possíveis.
+choose(n = 20, k = 10)
+
+#--------------------------------------------
+# Com todas as combinações possíveis (exaustão).
+
+# Para construir todas.
+k <- combn(x = 1:20, m = 10)
+dim(k)
+
# Vetor com os valores dos dois grupos.
mf <- c(m, f)
+g <- integer(20)
+D <- apply(k,
+ MARGIN = 2,
+ FUN = function(i) {
+ g[i] <- 1L
+ -diff(tapply(mf, g, FUN = mean))
+ })
+
+hist(D, col = "gray50")
+rug(D)
+abline(v = d, col = 2)
+
+plot(ecdf(D), cex = 0)
+rug(D)
+abline(v = d, col = 2)
+
+# P-valor do teste.
+2 * sum(D >= d)/length(D)
+
+#--------------------------------------------
+# Com simulação (não necessáriamente exaustivo).
+
# Variáveis que indentifica os grupos.
g <- rep(1:2, each = 10)
-cbind(mf, g)
+# Apenas para conferir.
+cbind(g, mf)
# Replicando a diferença para grupos formados por aleatorização.
-D <- replicate(999, {
- y <- sample(mf, size = length(mf), replace = FALSE)
- -diff(tapply(y, g, FUN = mean))
+D <- replicate(9999, {
+ gg <- sample(g)
+ -diff(tapply(mf, gg, FUN = mean))
})
+# Tem que juntar a estatística observada com as simuladas.
D <- c(D, d)
hist(D, col = "gray50")
@@ -53,38 +129,48 @@ plot(y ~ as.numeric(g))
fobs <- anova(lm(y ~ g))[1, "F value"]
fsim <- replicate(9999, {
- y <- sample(y, size = 12, replace = FALSE)
- anova(lm(y ~ g))[1, "F value"]
+ gg <- sample(g)
+ anova(lm(y ~ gg))[1, "F value"]
})
+# Junta a estatística observada com as obtidas por simulação.
fsim <- c(fsim, fobs)
# P-valor do teste.
sum(fsim >= fobs)/length(fsim)
+# (maneiras de agrupar 4 em 12) * (maneiras de agrupar 4 em 8).
choose(12, 4) * choose(8, 4)
#-----------------------------------------------------------------------
# Teste de aleatorização para a correlação.
+# N = 5 para um par de medidas.
x <- c(4, 8, 2, 10, 9)
y <- c(3, 5, 1, 7, 8)
-r0 <- cor(x, y)
+cbind(x, y)
+
+# r0 <- cor(x, y)
r0 <- cor(x, y, method = "spearman")
+r0
library(gtools)
-# Todas as permutações possiveis = 5! = 120.
-Y <- permutations(n = length(y), r = length(y), v = y)
-str(Y)
+# Todas as permutações possiveis = 5! = 120 do vetor x.
+X <- permutations(n = length(x), r = length(x), v = x)
+str(X)
+head(X)
-r <- apply(Y, MARGIN = 1, FUN = cor, y = x)
+# As 120 correlações obtidas para cada arranjo.
+r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
# P-valor do teste.
sum(r >= r0)/length(r)
#-----------------------------------------------------------------------
+# Exemplo da matriz de distância geográfica e similaridade entre
+# espécies.
# s <- scan()
# dput(s)
@@ -96,6 +182,8 @@ sim <- matrix(NA,
nrow = length(cont),
dimnames = list(cont, cont))
sim[lower.tri(sim)] <- s
+
+# Coeficientes de associação entre espécies de tesourinha.
sim
# j <- scan()
@@ -104,9 +192,11 @@ j <- c(1, 2, 1, 2, 3, 2, 1, 1, 2, 3, 4, 3, 2, 3, 4, 5, 4, 3, 1, 2, 3, 2,
1, 4, 3, 5, 4, 1)
jum <- sim
jum[lower.tri(jum)] <- j
+
+# Número de saltos entre continentes conforme movimento das placas.
jum
-# Correlação entre as distância (fere a suposição de independência).
+# Correlação entre os valores (fere a suposição de independência).
k <- cor(s, j)
# Empilha os valores de similaridade.
@@ -139,6 +229,12 @@ hist(K)
abline(v = k, col = 2)
plot(density(K))
+rug(K)
+abline(v = k, col = 2)
+
+plot(ecdf(K))
+rug(K)
+abline(v = k, col = 2)
# Mantel Matrix Randomization Test.
RSiteSearch("mantel randomization")
@@ -150,7 +246,8 @@ jumd
library(ade4)
-mr <- mantel.randtest(jumd, simd, nrepet = 999)
+# Passar valores negativos porque a hipótese H_a é "greater".
+mr <- mantel.randtest(-jumd, simd, nrepet = 999)
mr
plot(mr)