Revisa teste de aleatorização para DSBD.

scripts/ce089-04.R

 #-----------------------------------------------------------------------
-# 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)