Adiciona conteúdo sobre testes de aleatorização.

_site.yml

@@ -45,6 +45,10 @@ navbar:
         href: slides/07-jackknife.pdf
       - text: "Jackknife (tutorial)"
         href: tutoriais/08-jackknife.html
+      - text: "Testes de aleatorização (slides)"
+        href: slides/08-aleatorizacao.pdf
+      - text: "Testes de aleatorização (html)"
+        href: tutoriais/09-aleatorizacao.html
       - text: "----------"
       - text: "Arquivos complementares"
       - text: "GNA Uniformes (2015)"

slides/08-aleatorizacao.Rnw

+%-----------------------------------------------------------------------
+
+\documentclass[serif, professionalfont, usenames, dvipsnames]{beamer}
+\usepackage[T1]{fontenc}
+
+% ATTENTION: preamble.tex contains all style definitions.
+\input{config/preamble.tex}
+% \usepackage[backend=bibtex, style=authoryear]{biblatex}
+\addbibresource{config/refs.bib}
+
+<<include = FALSE>>=
+source("config/setup.R")
+@
+
+%-----------------------------------------------------------------------
+
+\title{Testes de aleatorização}
+\subtitle{Fundamentos e aplicações}
+\date{\small{ \Sexpr{sprintf('Atualizado em %s', Sys.Date())}}}
+
+%-----------------------------------------------------------------------
+
+\begin{document}
+
+{\setbeamertemplate{footline}{}
+  \frame{\titlepage} %--------------------------------------------------
+}
+
+\begin{frame}{}
+
+  {\large Justificativas}
+
+  \begin{itemize}
+  \item Métodos computacionalmente intensivos para inferência
+    estatística são usados quando as abordagens tradicionais não são
+    adequadas.
+    \begin{itemize}
+    \item Resultados assintóticos em pequenas amostras.
+    \item Violação de pressupostos.
+    \item Não existência de mecanísmos de inferência específicos.
+    \end{itemize}
+  \item Tais métodos se baseiam em reamostragem e/ou simulação.
+  \item Podem ser aplicados em muitos contextos.
+  \end{itemize}
+
+  {\large Objetivos}
+
+  \begin{itemize}
+    \item TODO
+  \end{itemize}
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{O funcionamento dos testes de hipótese}
+
+Para recapitular, a lógica dos testes de hipótese frequentistas é:
+
+\begin{enumerate}
+\item Definir a \textbf{hipótese nula} ($H_0$)e hipótese alternativa ($H_a$).
+\item Definir uma \textbf{estatística de teste} apropriada para avaliar
+  a hipótese em questão.
+\item Estabelecer a \textbf{região crítica} para tomar decisão sobre a
+  hipóteses estabelecidas.
+\item Calcular a estatística de teste a partir dos dados observados.
+\item Verificar se o valor observado da estística de teste está região
+  de aceitação da hipótese nula e tomar a decisão correspondente.
+\end{enumerate}
+
+A região crítica é baseada na \textbf{distribuição amostral} da
+estatística de teste sob a hipótese nula.  Pode envolver o
+estabelecimento de suposições.
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{Testes de Aleatorização}
+
+  \begin{itemize}
+  \item Abordagem baseada em permutação das observações
+    (\textit{permutation tests}).
+  \item São considerados testes livre de distribuição.
+  \item Faz suposições sobre o processo gerador dos dados.
+  \item Duas formas de cálculo da estatística de teste:
+    \begin{itemize}
+    \item Exaustiva: no conjunto de todos os arranjos possíveis
+      $\rightarrow$ distribuição amostral exata.
+    \item Por reamostragem: amostra do conjunto completo de arranjos com
+      reamostragem sem reposição.
+    \end{itemize}
+  \item IMPORTANTE: Sob a hipótese nula os dados são \hi{permutáveis}.
+  \item Esse é o principal conceito e requisito dos testes de
+    aleatorização.
+  \end{itemize}
+
+  % \myurl{http://genomicsclass.github.io/book/pages/permutation_tests.html}
+  % \myurl{http://www.columbia.edu/~ahd2125/post/2015/6/8/}
+  % \myurl{https://thomasleeper.com/Rcourse/Tutorials/permutationtests.html}
+  % \myurl{https://www.ohbmbrainmappingblog.com/blog/a-brief-overview-of-permutation-testing-with-examples}
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}{Uma senhora toma chá}
+
+  \begin{itemize}
+  \item Aconteceu com R. A. Fisher e Muriel Bristol.
+  \item Fisher descreve em seu livro em 1935.
+  \item A senhora \hi{declarou} saber discriminar bebida conforme a
+    ordem em que chá e leite eram adicionados à xícara.
+  \item $H_0$: a senhora \hi{não sabe distinguir} entre as misturas da
+    bebida, i.g. a classifição feita por ela é aleatória.
+  \item Experimento: 8 xícaras, 4 de cada tipo servidas aleatoriamente.
+  \item Resposta: a classificação de 4 xícaras de um tipo.
+  \end{itemize}
+
+\end{frame}
+
+%-----------------------------------------------------------------------
+\begin{frame}
+
+\begin{block}{Perguntas}
+  \begin{itemize}
+  \item Quantos arranjos entre as xícaras são possíveis?
+  \item Qual a chance da senhora acertar todas por mero acaso (sob $H_0$)?
+  \item Qual a chance de acertar 3 em 4?
+  \item Qual a região crítica?
+  \end{itemize}
+\end{block}
+
+\pause
+
+\begin{block}{Respostas}
+  \begin{itemize}
+  \item $\binom{8}{4} = \frac{8!}{4!(8-4)!} = 70$.
+  \item É $1/70$ pois só existe uma forma correta no universo das $70$
+    possibilidades.
+  \item ``Arranjos de 3 corretos em 4 selecionados'' $\times$
+    ``arranjos de 1 errado em 4 selecionados'':
+    $\binom{4}{3} \cdot \binom{4}{1} = 16$, então $16/70 \approx 0.23$.
+  \item Ao nível de 5\%, a hipótese nula será rejeitada apenas se a
+    senhora acertar as 4 xícaras pois 1/70 $\approx$ 0.14 $<$ 0.05.
+  \end{itemize}
+
+\end{block}
+\end{frame}
+
+%-----------------------------------------------------------------------
+{
+  \usebackgroundtemplate{\includegraphics[height=\paperheight, width=\paperwidth]{../img/looking-ahead.jpg}}
+  % \setbeamersize{text margin left=30mm}
+
+  \begin{frame}[b]{}
+
+    \hspace*{0.5\linewidth}
+    \begin{minipage}[t]{0.5\linewidth}
+
+      \hi{Próxima aula}
+      \begin{itemize}
+      \item Mais sobre Testes de aleatorização.
+      \end{itemize}
+
+      \hi{Avisos}
+      \begin{itemize}
+      \item Sabatina estará disponível a partir de Qua.
+      \end{itemize}
+
+      \vspace{3em}
+    \end{minipage}
+
+\end{frame}
+}
+
+%-----------------------------------------------------------------------
+% \begin{frame}[t, fragile, allowframebreaks]
+%   \frametitle{Referências bibliográficas}
+%
+%   \printbibliography[heading=none]
+% \end{frame}
+
+%-----------------------------------------------------------------------
+\end{document}

tutoriais/09-aleatorizacao.Rmd

+---
+title: "Aplicações de testes de aleatorização"
+author: Prof. Walmes M. Zeviani
+date: '`r Sys.Date()`'
+bibliography: ../config/Refs.bib
+csl: ../config/ABNT-UFPR-2011-Mendeley.csl
+---
+
+# 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}
+#-----------------------------------------------------------------------
+# 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,
+                    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) = 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 é mais realista mas o resultado é o memso.
+mean(replicate(n = 100000,
+               expr = all(sample(v) == sample(v))))
+```
+
+# Teste para a diferença de médias
+
+```{r}
+#-----------------------------------------------------------------------
+# Exemplo com teste para a diferença de médias.
+
+# 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")
+rug(f, col = "magenta")
+
+# Diferença de média observada.
+d <- mean(m) - mean(f)
+d
+
+# Todos as combinações possíveis entre 20 objetos tomados 10 a 10.
+choose(n = 20, k = 10)
+
+#-----------------------------------------------------------------------
+# 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)
+dim(k)
+
+# Vetor com os valores dos dois grupos concatenados.
+mf <- c(m, f)
+
+# Vetor de zeros.
+g <- integer(20)
+
+# Calcula a diferença para cada uma das permutações.
+D <- apply(k,
+           MARGIN = 2,
+           FUN = function(i) {
+               # 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 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 (bilateral).
+2 * sum(D >= d)/length(D)
+
+#-----------------------------------------------------------------------
+# Com reamostragem sem reposição (não necessáriamente exaustivo).
+
+# Variáveis que indentifica os grupos.
+g <- rep(1:2, each = 10)
+
+# Apenas para conferir.
+cbind(g, mf)
+
+# Replicando a diferença para grupos formados por aleatorização.
+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")
+rug(D)
+abline(v = d, col = 2)
+
+plot(ecdf(D), cex = 0)
+rug(D)
+abline(v = d, col = 2)
+
+# P-valor do teste (bilateral).
+2 * sum(D >= d)/length(D)
+```
+
+# Teste para a correlação
+
+```{r}
+#-----------------------------------------------------------------------
+# Teste de aleatorização para a correlação.
+
+# 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: coeficiente de correlação de
+# Pearson.
+r0 <- cor(x, y, method = "pearson")
+r0
+
+# Todas as permutações possiveis: 5! = 120 do vetor x.
+X <- gtools::permutations(n = length(x), r = length(x), v = x)
+str(X)
+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)
+```
+
+# Índice de Moran
+
+  * <http://gis.stackexchange.com/questions/161887/significance-test-for-morans-i-using-monte-carlo-simulation>.
+  * <https://en.wikipedia.org/wiki/Moran's_I>.
+  * <http://rstudio-pubs-static.s3.amazonaws.com/268058_dd37b98f25a4496b9f0a7eb2fcf892cd.html>.
+  * <http://rspatial.org/analysis/rst/3-spauto.html>.
+
+```{r}
+#-----------------------------------------------------------------------
+# Índice de Moran para medir dependência espacial.
+
+# Coordenadas dos eventos em uma malha regular 8 x 8.
+x <- 1:8
+y <- 1:8
+
+# 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, ])))
+    w <- c(0, 1, sqrt(1/2), 0)[u + 1]
+    return(w)
+}
+#  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)
+
+#  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) {
+    # Tamanho da amostra.
+    n <- length(x)
+    # Valores normalizados.
+    z <- as.vector((x - mean(x))/sd(x))
+    # Índice de Moran.
+    as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
+}
+
+# Teste de permutação com saída gráfica.
+ppt <- function(z, w, N = 10000, stat, ...) {
+    # Índice de Moran por reamostragem.
+    sim <- replicate(N,
+                     moran(sample(z), w))
+    # Determina o p-valor.
+    p.value <- mean((all <- c(stat, sim)) >= stat)
+    # Histograma da distribuição empírica sob H_0.
+    hist(sim,
+         sub = paste("p =", round(p.value, 4)),
+         xlim = range(all),
+         ...)
+    abline(v = stat, col = "#903030", lty = 3, lwd = 2)
+    return(p.value)
+}
+
+# Observações simuladas.
+set.seed(17)
+par(mfrow = c(2, 3))
+
+# Dados com dependência espacial ---------------------------------------
+# Indução de autocorrelação por meio de uma tendência.
+z <- matrix(rexp(length(x) * length(y),
+                 outer(x, y^2)),
+            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
+
+hist(z)
+ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
+
+# Dados sem dependência espacial ---------------------------------------
+# Geração de de um conjunto de dados sob hipótese nula.
+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
+
+hist(z)
+ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
+```