diff --git a/_site.yml b/_site.yml index c4f899c126193faa17f055d8200ba7cf9a403f7c..b96c1ddae0b935c47eaf6a9979ec56de77764a90 100644 --- a/_site.yml +++ b/_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)" diff --git a/slides/08-aleatorizacao.Rnw b/slides/08-aleatorizacao.Rnw new file mode 100644 index 0000000000000000000000000000000000000000..e414869d435f108e2205eb1c26b32ee2809df429 --- /dev/null +++ b/slides/08-aleatorizacao.Rnw @@ -0,0 +1,181 @@ +%----------------------------------------------------------------------- + +\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} diff --git a/tutoriais/09-aleatorizacao.Rmd b/tutoriais/09-aleatorizacao.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..a0d4fe68eaee36445517f8ddc8663b91c5a56a34 --- /dev/null +++ b/tutoriais/09-aleatorizacao.Rmd @@ -0,0 +1,287 @@ +--- +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") +```