Conclui arquivo com exemplos de métodos intensivos.

scripts/metodos-intensivos.Rmd

@@ -19,18 +19,44 @@ opts_chunk$set(cache = TRUE,
                message = FALSE)
 ```
 
+Esse documento contém exemplos de aplicação de métodos
+computacionalmente intensivos em problemas de inferência estatística.
+Tais problemas são tipicamente testes de hipóteses, estimação intervalar
+e estudo de propriedades de estimadores.  Os métodos exemplificados são:
+
+  * Testes de aleatorização
+  * Métodos Jackknife.
+  * Métodos bootstrap.
+  * Métodos Monte Carlo.
+
+As sessões a seguir estão dividas em exemplos.
+
 # Testes de aleatorização
 
 ## Uma senhora toma chá
 
+Visite
+
+  * <https://en.wikipedia.org/wiki/Lady_tasting_tea>.
+  * Série de vídeos.
+    * <https://www.youtube.com/watch?v=vYVr50hjFbQ>.
+    * <https://www.youtube.com/watch?v=xh20btybjp4>.
+    * <https://www.youtube.com/watch?v=JMDycaD-YwY>.
+    * <https://www.youtube.com/watch?v=dQT6T3Zo-oA>.
+    * <https://www.youtube.com/watch?v=9sgQMHtLPNU>.
+  * <https://www.youtube.com/watch?v=lgs7d5saFFc>.
+
 ```{r}
 #-----------------------------------------------------------------------
-# Uma senhora toma chá.
+# O curioso caso da senhora que toma chá.
 
 # Quantidade de maneiras de gerar dois grupos de 4 xícaras usando 8.
 n_poss <- choose(8, 4)
 n_poss
 
+# E se fossem 12 xícaras, 6 de cada?
+choose(12, 6)
+
 # De 4 xícaras selecionadas para um grupo, X representa o número de
 # acertos.
 n_acertar <- sapply(4:0,
@@ -40,15 +66,18 @@ n_acertar <- sapply(4:0,
                     })
 n_acertar
 
-# Pr(Acerto total) = 1/70 < 5%.
-cumsum(n_acertar)/n_poss
+# Pr(Acerto total) = Pr(X = 4) = 1/70 < 5%.
+cbind(X = 4:0, "Pr(X = x)" = n_acertar/n_poss)
 
 # Comprovando por simulação.
+# Uma particular configuração das xícaras.
 v <- rep(0:1, each = 4)
+
+# Com qual frequência, ao acaso, a mesma configuração ocorre?
 mean(replicate(n = 100000,
                expr = all(sample(v) == v)))
 
-# Essa forma de fazer a conta é mais realista mas dá na mesma.
+# Essa forma de fazer é mais realista mas o resultado é o memso.
 mean(replicate(n = 100000,
                expr = all(sample(v) == sample(v))))
 ```
@@ -59,10 +88,11 @@ mean(replicate(n = 100000,
 #-----------------------------------------------------------------------
 # 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)
+# Comprimentos de mandíbulas de cães pré-históricos.
+m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112) # machos.
+f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111) # fêmeas.
 
+#  Distribuição acumulada empírica.
 plot(ecdf(m), xlim = range(c(m, f)), col = "cyan")
 lines(ecdf(f), col = "magenta")
 rug(m, col = "cyan")
@@ -72,11 +102,11 @@ rug(f, col = "magenta")
 d <- mean(m) - mean(f)
 d
 
-# Todos as combinações possíveis.
+# Todos as combinações possíveis entre 20 objetos tomados 10 a 10.
 choose(n = 20, k = 10)
 
-#--------------------------------------------
-# Com todas as combinações possíveis (exaustão).
+#-----------------------------------------------------------------------
+# Solução com todas as permutações possíveis (exaustão).
 
 # Para construir todas as combinações possíveis.
 k <- combn(x = 1:20, m = 10)
@@ -85,33 +115,34 @@ dim(k)
 # Vetor com os valores dos dois grupos concatenados.
 mf <- c(m, f)
 
-# Vetor cheio de zeros.
+# Vetor de zeros.
 g <- integer(20)
 
-# Calcula a diferença para todo arranjo possível.
+# Calcula a diferença para cada uma das permutações.
 D <- apply(k,
            MARGIN = 2,
            FUN = function(i) {
-               # Alguns elementos do vetor convertidos para 1.
+               # Alguns zeros são convertidos para 1 criando 2 grupos.
                g[i] <- 1L
                # Cálculo da estatística de teste.
                -diff(tapply(mf, g, FUN = mean))
            })
 
-# Histograma da distribuição da estatística sib H_0.
+# Histograma da distribuição da estatística sob H_0.
 hist(D, col = "gray50")
 rug(D)
 abline(v = d, col = 2)
 
+# Distribuição acumulada empírica da estatística sob H_0.
 plot(ecdf(D), cex = 0)
 rug(D)
 abline(v = d, col = 2)
 
-# P-valor do teste.
+# P-valor do teste (bilateral).
 2 * sum(D >= d)/length(D)
 
-#--------------------------------------------
-# Com simulação (não necessáriamente exaustivo).
+#-----------------------------------------------------------------------
+# Com reamostragem sem reposição (não necessáriamente exaustivo).
 
 # Variáveis que indentifica os grupos.
 g <- rep(1:2, each = 10)
@@ -136,7 +167,7 @@ plot(ecdf(D), cex = 0)
 rug(D)
 abline(v = d, col = 2)
 
-# P-valor do teste.
+# P-valor do teste (bilateral).
 2 * sum(D >= d)/length(D)
 ```
 
@@ -149,12 +180,12 @@ abline(v = d, col = 2)
 # N = 5 para um par de medidas.
 x <- c(4.1, 8.3, 2.9, 10.8, 9.5)
 y <- c(3.7, 5.1, 1.0, 7.7, 8.9)
-
 cbind(x, y)
 
 plot(x, y)
 
-# Estatística de teste na amostra original.
+# Estatística de teste na amostra original: coeficiente de correlação de
+# Pearson.
 r0 <- cor(x, y, method = "pearson")
 r0
 
@@ -166,6 +197,11 @@ head(X)
 # As 120 correlações obtidas para cada arranjo.
 r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
 
+#  Gráfico da distribuição da estatística de teste.
+plot(ecdf(r), cex = 0)
+rug(r)
+abline(v = r0, col = 2)
+
 # P-valor do teste.
 2 * sum(r >= r0)/length(r)
 ```
@@ -181,25 +217,38 @@ r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
 #-----------------------------------------------------------------------
 # Índice de Moran para medir dependência espacial.
 
-# Coordenadas dos eventos em uma malha regular.
+# Coordenadas dos eventos em uma malha regular 8 x 8.
 x <- 1:8
 y <- 1:8
 
-# Para criar a matriz de pesos.
-ind <- expand.grid(i = 1:length(x), j = 1:length(y))
+# Construção da matriz de pesos que determina a vizinhança entre
+# observações.
+ind <- expand.grid(i = 1:length(x),
+                   j = 1:length(y))
+#  Função que determina o peso entre duas localizações na malha.
 f <- function(i, j) {
     u <- min(3, sum(abs(ind[i, ] - ind[j, ])))
-    c(0, 1, sqrt(1/2), 0)[u + 1]
+    w <- c(0, 1, sqrt(1/2), 0)[u + 1]
+    return(w)
 }
-w <- matrix(0, nrow(ind), nrow(ind))
+#  Cria os pesos, matriz (8^2) x (8^2) = 64 x 64.
+w <- matrix(0, nrow = nrow(ind), ncol = nrow(ind))
 for (i in 1:nrow(ind)) {
     for (j in 1:nrow(ind)) {
         w[i, j] <- f(i, j)
     }
 }
+#  Normaliza.
 w <- w/sum(w)
 
-image(w)
+#  Gráfico. Valores claros indicam maior peso entre observações.
+image(w, asp = 1, col = gray.colors(100))
+
+#  Lógica do índica de Moran: correlação entre valores observados e
+#  média dos vizinhos. Exemplo com valores simulados.
+xx <- rnorm(64)
+cor(cbind("Valores observados" = xx,
+          "Média dos vizinhos" = as.vector(xx %*% w)))
 
 # Índice de Moran.
 moran <- function(x, weights) {
@@ -211,7 +260,7 @@ moran <- function(x, weights) {
     as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
 }
 
-# Teste de permutação com gráfico.
+# Teste de permutação com saída gráfica.
 ppt <- function(z, w, N = 10000, stat, ...) {
     # Índice de Moran por reamostragem.
     sim <- replicate(N,
@@ -238,6 +287,9 @@ z <- matrix(rexp(length(x) * length(y),
             length(x))
 image(log(z), main = "Com dependência")
 
+cor(cbind("Valores observados" = as.vector(z),
+          "Média dos vizinhos" = as.vector(as.vector(z) %*% w)))
+
 # Índice de Moran com dados originais.
 stat <- moran(z, w)
 stat
@@ -250,6 +302,9 @@ ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
 z <- matrix(rnorm(length(x) * length(y), 0, 1/2), length(x))
 image(z, main = "Sem dependência")
 
+cor(cbind("Valores observados" = as.vector(z),
+          "Média dos vizinhos" = as.vector(as.vector(z) %*% w)))
+
 # Índice de Moran com dados originais.
 stat <- moran(z, w)
 stat
@@ -266,9 +321,12 @@ ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
 #-----------------------------------------------------------------------
 # Aplicação de Jackknife em modelos de regressão não linear.
 
-library(latticeExtra)
 library(car)
 library(alr3)
+library(latticeExtra)
+
+# Documentação do conjunto de dados.
+# help(segreg)
 
 # Curva "palpite".
 start <- list(th0 = 75, th1 = 0.5, th2 = 50)
@@ -277,7 +335,7 @@ xyplot(C ~ Temp, data = segreg) +
                       0 * (x < th2), lty = 2),
           data = start)
 
-# Ajuste.
+# Ajuste do modelo não linear aos dados (modelo platô-linear).
 n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
               0 * (Temp < th2),
           data = segreg,
@@ -286,7 +344,7 @@ n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
 # Estimativas e medidas de ajuste.
 summary(n0)
 
-# Observados e preditos.
+# Valores observados curva ajustada.
 xyplot(C ~ Temp, data = segreg) +
     layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
                       0 * (x < th2)),
@@ -295,11 +353,10 @@ xyplot(C ~ Temp, data = segreg) +
 #-----------------------------------------------------------------------
 # Estimativas Jackknife para os parâmetros.
 
-
 theta <- coef(n0) # Estimativas com todos os dados.
 n <- nrow(segreg) # Número de observações.
 
-# Matriz vazia para armazenar os pseudo-valores.
+# Matriz vazia para armazenar os pseudo-valores durante iterações.
 thetaj <- matrix(0, nrow = n, ncol = length(theta))
 colnames(thetaj) <- names(theta)
 
@@ -314,16 +371,19 @@ for (i in 1:n) {
     thetaj[i, ] <- n * theta - (n - 1) * coef(n1)
 }
 
-# Os primeiros pseudo valores.
+# Alguns pseudo-valores.
 head(thetaj)
+tail(thetaj)
 
-# Estimativas pontuais.
+# Estimativas pontuais de Jackknife.
 jk <- colMeans(thetaj)
-cbind(MLE = theta, Jack = jk, Vício = theta - jk)
+rbind(MLE = theta,
+      Jackknife = jk,
+      Vício = theta - jk)
 
 # Erros padrões.
-cbind(MLE = summary(n0)$coefficients[, "Std. Error"],
-      Jack = apply(thetaj, MARGIN = 2, sd)/sqrt(n))
+rbind(MLE = summary(n0)$coefficients[, "Std. Error"],
+      Jackknife = apply(thetaj, MARGIN = 2, sd)/sqrt(n))
 
 # Pacote com função pronta para Jackknife para modelos não lineares.
 library(nlstools)
@@ -340,6 +400,8 @@ summary(j0)
 #-----------------------------------------------------------------------
 # Inferência para a DL 50.
 
+# Dados no pacote labestData: https://github.com/pet-estatistica/labestData.
+
 # library(labestData)
 # data(PaulaTb3.12)
 # str(PaulaTb3.12)
@@ -349,40 +411,50 @@ summary(j0)
 # dput(PaulaTb3.12)
 
 PaulaTb3.12 <-
-structure(list(vol = c(3.7, 3.5, 1.25, 0.75, 0.8, 0.7, 0.6, 1.1,
-0.9, 0.9, 0.8, 0.55, 0.6, 1.4, 0.75, 2.3, 3.2, 0.85, 1.7, 1.8,
-0.4, 0.95, 1.35, 1.5, 1.6, 0.6, 1.8, 0.95, 1.9, 1.6, 2.7, 2.35,
-1.1, 1.1, 1.2, 0.8, 0.95, 0.75, 1.3), razao = c(0.825, 1.09,
-2.5, 1.5, 3.2, 3.5, 0.75, 1.7, 0.75, 0.45, 0.57, 2.75, 3, 2.33,
-3.75, 1.64, 1.6, 1.415, 1.06, 1.8, 2, 1.36, 1.35, 1.36, 1.78,
-1.5, 1.5, 1.9, 0.95, 0.4, 0.75, 0.03, 1.83, 2.2, 2, 3.33, 1.9,
-1.9, 1.625), resp = c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
-1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1,
-1, 1, 0, 0, 1)), .Names = c("vol", "razao", "resp"), row.names = c(NA,
-39L), class = "data.frame")
+    structure(list(
+        vol = c(3.7, 3.5, 1.25, 0.75, 0.8, 0.7, 0.6, 1.1, 0.9, 0.9, 0.8,
+                0.55, 0.6, 1.4, 0.75, 2.3, 3.2, 0.85, 1.7, 1.8, 0.4,
+                0.95, 1.35, 1.5, 1.6, 0.6, 1.8, 0.95, 1.9, 1.6, 2.7,
+                2.35, 1.1, 1.1, 1.2, 0.8, 0.95, 0.75, 1.3),
+        razao = c(0.825, 1.09, 2.5, 1.5, 3.2, 3.5, 0.75, 1.7, 0.75,
+                  0.45, 0.57, 2.75, 3, 2.33, 3.75, 1.64, 1.6, 1.415,
+                  1.06, 1.8, 2, 1.36, 1.35, 1.36, 1.78, 1.5, 1.5, 1.9,
+                  0.95, 0.4, 0.75, 0.03, 1.83, 2.2, 2, 3.33, 1.9, 1.9,
+                  1.625),
+        resp = c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1,
+                 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1,
+                 0, 0, 1)),
+        .Names = c("vol", "razao", "resp"),
+        row.names = c(NA, 39L),
+        class = "data.frame")
 
 layout(1)
-plot(resp ~ vol, data = PaulaTb3.12)
 
+#  Visualização dos dados.
+xyplot(resp ~ vol,
+       data = PaulaTb3.12,
+       groups = resp)
+
+#  Ajuste do modelo considerando resposta binária.
 m0 <- glm(resp ~ vol,
           data = PaulaTb3.12,
           family = binomial)
 summary(m0)
 
-# DL_50.
+# Estimativa da DL_50 a partir do modelo ajustado.
 dl <- -coef(m0)[1]/coef(m0)[2]
 
-# str(m0$family)$link
-
+# Curva ajustada.
 plot(resp ~ vol, data = PaulaTb3.12)
 curve(m0$family$linkinv(coef(m0)[1] + coef(m0)[2] * x),
       add = TRUE)
 abline(v = dl, h = 0.5, lty = 2)
 
+#  Número de observações e vetor de índices.
 n <- nrow(PaulaTb3.12)
 i <- 1:n
 
-# Estimativas por Jackknife.
+# Estimativas por Jackknife: pseudo-valores.
 DL <- sapply(i,
              FUN = function(i) {
                  m0 <- glm(resp ~ vol,
@@ -394,24 +466,30 @@ DL <- sapply(i,
                  n * dl - (n - 1) * DL
              })
 
+#  Estimativa pontual e erro-padrão por Jackknife.
 m <- mean(DL)
 e <- sd(DL)/sqrt(n)
 
+#  Comparando com as estimativas de Wald (usa método Delta).
 cbind(Est = rbind(MLE = dl,
-                  Jack = m),
+                  Jackknife = m),
       StdErr = rbind(MLE = car::deltaMethod(m0, g = "-(Intercept)/vol")["SE"],
                      Jack = e))
 
+# Distribuição acumulada empírica.
 plot(ecdf(DL))
 rug(DL)
+abline(v = c(m, dl), lty = c(1, 2), col = 2)
 
+#  Densidade empírica.
 plot(density(DL))
 rug(DL)
+abline(v = c(m, dl), lty = c(1, 2), col = 2)
 
 # NOTE: alguma pista sobre porque a distribuição fica bimodal?
 ```
 
-## Uma situação problemática
+## Condição problemática para uso do método Jackknife
 
 ```{r}
 #-----------------------------------------------------------------------
@@ -419,6 +497,7 @@ rug(DL)
 
 # Ordena o vetor de valores.
 x <- sort(precip)
+x
 
 # Calcula a mediana com os dados originais.
 M <- median(x)
@@ -453,6 +532,9 @@ length(x)
 library(alr3)
 str(turk0)
 
+# Documentação do conjunto de dados.
+# help(turk0)
+
 # Gráfico dos valores observados.
 plot(Gain ~ A, data = turk0)
 
@@ -481,14 +563,16 @@ cbind(coef(n0), coef(n1))
 #-----------------------------------------------------------------------
 # Criando funções que fazem reamostragem dos dados e ajustam o modelo.
 
+# Número de observações e índices.
 n <- nrow(turk0)
 i <- 1:n
 
-myboot0 <- function(formula, start) {
+#  Função que obtém estimativas bootstrap.
+myboot <- function(formula, start) {
     n0 <- nls(formula = formula,
               data = turk0[sample(i, n, replace = TRUE), ],
               start = start)
-    coef(n0)
+    return(coef(n0))
 }
 
 #-----------------------------------------------------------------------
@@ -496,18 +580,28 @@ myboot0 <- function(formula, start) {
 
 # B = 999 reamostragens.
 b0 <- replicate(999,
-                myboot0(Gain ~ Int + Ass * A/(Mei + A), start))
+                myboot(Gain ~ Int + Ass * A/(Mei + A), start))
 b1 <- replicate(999,
-                myboot0(Gain ~ Int + Ass * (1 - 2^(-A/Mei)), start))
+                myboot(Gain ~ Int + Ass * (1 - 2^(-A/Mei)), start))
 
 # Conjunto de 1000 estimativas bootstrap.
 b0 <- cbind(b0, coef(n0))
 b1 <- cbind(b1, coef(n1))
 str(b0)
 
+#  Estimativas originais.
+mle_est <- rbind(n0 = coef(n0),
+                 n0 = coef(n1))
+
+# Estimativas bootstrap.
+bot_est <- rbind(n0 = rowMeans(b0),
+                 n1 = rowMeans(b1))
+
+mle_est
+bot_est
+
 # Calcula o vício absoluto (bootstrap - original).
-rbind(n0 = rowMeans(b0) - coef(n0),
-      n1 = rowMeans(b1) - coef(n1))
+bot_est - mle_est
 
 # Variância do estimador.
 rbind(n0 = apply(b0, MARGIN = 1, var),
@@ -531,8 +625,12 @@ rbind(n0 = apply(b0, 1, FUN = eqm),
 
 bts <- cbind(rep(1:2, each = ncol(b0)),
              as.data.frame(rbind(t(b0), t(b1))))
+
+# Matriz de diagramas de dispersão das estimativas nas amostras
+# bootstrap.
 splom(bts[, -1], groups = bts[, 1])
 
+# Cada uma das curvas que foram ajustas às amostras bootstrap.
 par(mfrow = c(1, 2))
 plot(Gain ~ A, data = turk0, type = "n")
 apply(b0,
@@ -567,8 +665,12 @@ p1 <- apply(b1, MARGIN = 2,
                 b[1] + b[2] * (1 - 2^(-0.2/b[3]))
             })
 
-cbind(Est = rbind(mean(p0), mean(p1)),
-      StdErr = rbind(sd(p0), sd(p1)))
+# Estimativa e erro-padrão para o valor predito em x = 0.4.
+est <- cbind(rbind(mean(p0), mean(p1)),
+             rbind(sd(p0), sd(p1)))
+dimnames(est) <- list(c("n0", "n1"),
+                      c("Est", "StdErr"))
+est
 ```
 
 ## Usando o pacote `boot`
@@ -578,7 +680,6 @@ cbind(Est = rbind(mean(p0), mean(p1)),
 # Usando o pacote boot.
 
 library(boot)
-library(latticeExtra)
 
 # dput(as.list(coef(n1)))
 start <- list(Int = 622.958054146589,
@@ -661,8 +762,9 @@ confint(b0, type = "bca")
 #-----------------------------------------------------------------------
 # Teste para independência de processo pontual.
 
-plot(NULL, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1)
-lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0))
+# layout(1)
+# plot(x = NULL, y = NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1)
+# lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0))
 # xy <- locator(n = 20, type = "p", pch = 19)
 # dput(lapply(xy, round, digits = 3))
 
@@ -835,9 +937,6 @@ system.time(m2 <- replicate(1000, rss(m = m)))
 # Distribuição da média na amostra por conjuntos ordenados imperfeito.
 system.time(m3 <- replicate(1000, rssur(m = m, rho = 0.75)))
 
-library(lattice)
-library(latticeExtra)
-
 densityplot(~m1 + m2 + m3,
             auto.key = TRUE,
             layout = c(NA, 1)) +