---
title: "RRR: Reproduzindo o Resampling com Rsampling"
author: "Paulo Inácio Prado"
date: "Outubro de 2015"
output:
  rmarkdown::html_vignette:
    fig_width: 5
    fig_height: 5
    fig_caption: true
vignette: >
  %\VignetteIndexEntry{Introdução ao Rsampling (PT-BR)}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, echo = FALSE}
knitr::opts_chunk$set(
    collapse=TRUE,
    comment = NA,
    prompt = TRUE
    )
set.seed(42)
```

## Apresentação

Este roteiro é uma introdução ao pacote **Rsampling**, que reproduz em R as funções
do programa *Resampling Stats*
(http://www.resample.com/).

Essas funções são usadas em um ciclo de trabalho que resume a lógica
dos testes de significância:

1. Defina uma estatística de interesse;
2. Defina a hipótese nula;
3. Obtenha a distribuição da estatística de interesse sob a hipótese nula;
4. Se a probabilidade da estatística de interesse observada ocorrer sob a hipótese nula
  é menor do que um valor crítico rejeite a hipótese nula.

A ideia do *Resampling Stats* é facilitar o entendimento dessa lógica,
fazendo o usuário executar cada um dos passos em um planilha,
com o auxílio de algumas macros. Um elemento muito efetivo para este
aprendizado é que a hipótese nula é simulada por aleatorização dos dados.
Isso também dá muita flexibilidade aos testes que podem ser feitos.
O manual do *Resampling Stats* é uma excelente introdução a esta metodologia,
e à lógica dos testes de significância [^2].

O objetivo do pacote **Rsampling**
é possibilitar este mesmo treinamento no R.
Assim, privilegiamos fidelidade à lógica
original e à didática em eventual detrimento de
desempenho computacional.

As seções após instruções de instalação
são exemplos de uso mais simples e comuns
do **Rsampling**. Consulte também as páginas
de ajuda do pacote para conhecer todas
as funcionalidades.


## Instalação

O Rsampling está no repositório oficial de pacotes do R (CRAN).
Para instalá-lo use

```{r installation CRAN, eval=FALSE}
install.packages("Rsampling")
```

Você pode também instalar versão de desenvolvimento,
que está no GitHub. Para isso você vai precisar da função `install_github` do pacote devtools:

```{r installation, eval=FALSE}
library(devtools)
install_github(repo = 'lageIBUSP/Rsampling')
```

Depois de instalar o pacote carregue-o em sua seção de R com 

```{r load library}
library(Rsampling)
```

## Embaralhando dentro de uma coluna para testar diferença entre grupos

O dataframe `embauba` tem os dados de presença
e ausência de lianas em embaúbas de dois morfotipos (brancas e vermelhas).

```{r inspecionando objeto embauba}
head(embauba)
summary(embauba)
```
Para mais detalhes sobre os dados e o estudo que os produziu consulte a
página de ajuda (`?embauba`).

### Hipótese do estudo

A hipótese deste estudo é
que as formigas removem lianas das embaúbas onde estão suas colônias.
A previsão é que embaúbas vermelhas seriam menos infestadas por lianas do que as
brancas, por abrigarem colônias de formigas mais frequentemente.
De fato, esta diferença é observada nas proporções de árvores
infestadas na amostra:

```{r proporcao de infestacao por morfo de embauba}
tapply(embauba$with.vines, embauba$morphotype, mean)
```
### Hipótese nula

A hipótese nula é de que as proporções de infestação são iguais
na população de onde vieram as amostras.
Sob esta hipótese, uma liana tem a mesma chance de estar em uma embaúba
branca ou vermelha.
Simulamos a hipótese nula
embaralhando as presenças de lianas entre plantas
na tabela de dados.

### Estatística de interesse

A cada simulação temos que calcular nossa
**estatística de interesse**, que é a
a diferença de infestação
entre os dois morfos.
Criamos uma função para isso:

```{r estatistica de interesse embaubas}
emb.ei <- function(dataframe){
    props <- tapply(dataframe$with.vines, dataframe$morphotype, mean)
    props[[1]] - props[[2]]
}
## Verificando
emb.ei(embauba)
```
### Distribuição da estatística sob a hipótese nula

Em seguida fazemos a simulação com a função
`Rsampling`:

```{r embaubas resampling, results="hide"}
emb.r <- Rsampling(type = "normal", dataframe = embauba,
                   statistics = emb.ei, cols = 2, ntrials = 1000)
```
**O que significa este comando?**

* `type = "normal"` escolhe uma randomização de todos os elementos
		(mais abaixo você verá outros tipos de randomização).
* `dataframe = embauba` indica a tabela com os dados
* `statistics = emb.ei` indica a função que calcula a(s)
	estatística(s) de interesse da tabela de dados.
* `cols = 2` indica que a randomização deve ser feita sobre a segunda
  coluna da tabela de dados.
* `ntrials = 1000` indica o número de repetições da simulação.

A distribuição das estatística de interesse
na simulação nem incluiu o valor observado:

```{r embaubas distribuicao nula, fig.cap="Distribuição das diferenças nas proporções de embaúbas brancas e vermelhas com lianas em 1000 simulações da hipótese nula de ausência de diferença nas populações amostradas. A linha vermelha indica a diferença observada. A região de aceitação da hipótese nula para 5% de significância está delimitada em cinza."}
dplot(emb.r, svalue = emb.ei(embauba), pside="Greater",
      main = "Distribuição da estatística de interesse sob H0",
      xlab = "Estatística de interesse")
```


### Decisão: rejeitamos a hipótese nula?

Seguindo o padrão das ciências biológicas,
adotamos o critério de rejeitar
a hipótese nula se a probabilidade
da estatística de interesse sob a hipótese nula
for menor que 5%.

No gráfico o que não está marcado em cinza são os 5%
mais extremos da distribuição da estatística sob a hipótese nula.
Então se a estatística observada estiver na região cinza não rejeitamos
a hipótese nula. Esta é a chamada \emph{região de aceitação} de H0.
Como o valor observado (linha vermelha) está fora da região de aceitação,
podemos rejeitar H0.
Você também pode verificar isso com um teste lógico no
R:

```{r embaubas teste} 
sum(emb.r >= emb.ei(embauba))/1000 < 0.05
```

**Conclusão:** rejeita-se a hipótese nula (p < 0,05).


## Embaralhando dentro de linhas para testar diferenças dentro de pares

O dataframe `azteca` tem o número de formigas *Azteca* sp
recrutadas por extratos aquosos de folhas novas e velhas de
embaúbas.

```{r inspecionando objeto azteca}
head(azteca)
summary(azteca)
```
Saiba mais sobre os dados em sua página de ajuda (`?azteca`).

### Hipótese do estudo

A hipótese do estudo é que
o recrutamento é mais intenso quando uma
folha nova é danificada.
A previsão para o experimento é
que o recrutamento por extrato de folhas novas
seja maior, o que ocorreu:

```{r pairplot azteca, fig.cap = "Número de formigas recrutadas por extratos de folhas novas e velhas de embaúbas. Os extratos foram aplicados em pares de folhas próximas em embaúbas que tinham colônias de formigas. As linhas ligam folhas do mesmo par experimental."}
splot(azteca$extract.new, azteca$extract.old,
           group.name=c("Folha nova","Folha velha"),
           ylab="N de formigas recrutadas",
           xlab="Tipo de extrato aplicado")
```

### Hipótese nula

A hipótese nula é de que o recrutamento provocado pelos estratos
é o mesmo. Note que para controlar outras fontes de variação o
experimento foi pareado.
Então para simular a hipótese nula temos que
embaralhar o número de formigas recrutadas **dentro** de cada par de
folhas.

### Estatística de interesse

A cada simulação temos que calcular nossa
**estatística de interesse**, que é a
média da diferença das folhas de cada par.
Uma função para isso:

```{r estatistica de interesse azteca}
azt.ei <- function(dataframe){
    diferencas <- with(dataframe, extract.new - extract.old)
    mean(diferencas)
}
## Valor observado
azt.ei(azteca)
```

No experimento o extrato de folhas novas recrutou em média
`r round(azt.ei(azteca),1)` formigas que o extrato de folha velha, em cada par.

### Distribuição da estatística sob a hipótese nula

Como os pares são linhas em nosso dataframe,
simulamos a hipótese nula embaralhando os valores
dentro de cada linha:

```{r azteca resampling, results="hide"}
azt.r <- Rsampling(type = "within_rows", dataframe = azteca,
                   statistics = azt.ei, cols = 2:3, ntrials = 1000)
```

Mudamos o argumento `type = "within_rows"`, para indicar que
os valores devem ser embaralhados dentro das linhas.
O argumento `cols = 2:3` indica as colunas do dataframe
que têm as contagens.

Uma diferença igual ou maior que a observada foi muito rara
na distribuição da estatística de interesse:

```{r azteca distribuicao nula, fig.cap="Distribuição das diferenças do número de formigas recrutadas por extratos de folhas novas e velhas de embaúba em pares experimentais, em 1000 simulações da hipótese nula de ausência de diferença. A linha vermelha indica a diferença observada. A região de aceitação da hipótese nula para 5% de significância está delimitada em cinza."}
dplot(azt.r, svalue = azt.ei(azteca), pside="Greater",
      main = "Distribuição da estatística de interesse sob H0",
      xlab = "Estatística de interesse")
```

### Decisão: rejeitamos a hipótese nula?

Novamente o gráfico mostra que o valor observado da estatística está fora da região de aceitação da hipótese nula sob nosso critério de significância (5% de chance de erro).
O mesmo resultado é verificado com o teste lógico:

```{r azteca teste} 
sum(azt.r >= azt.ei(azteca))/1000 < 0.05
```

**Conclusão:** rejeita-se a hipótese nula (p<0,05).

#### Coda: testes unicaudais e bicaudais

Até agora testamos hipóteses de que um valor **igual ou maior** que o observado
pode ser gerado pela hipótese nula. É um teste **unicaudal** ou **unidirecional**, como
seria também o teste de que um valor igual ou menor pode ser gerado pela hipótese nula.
Nos testes unicaudais a região de aceitação é toda a distribuição nula exceto seus 5% mais extremos.

Mas pode interessar o teste de que há diferenças, sem especificar sua direção. Por exemplo,
o conhecimento prévio poderia indicar a hipótese de que extratos de folhas jovens e velhas devem recrutar
números diferentes de formigas, mas sem a expectativa de qual extrato recrutaria mais. Este é um caso
de teste **bicaudal**, quando a região de aceitação é o centro da distribuição nula, exceto seus
2,5% mais extremos de cada lado:

```{r azteca distribuicao nula bicaudal, fig.cap="Distribuição das diferenças do número de formigas recrutadas por extratos de folhas novas e velhas de embaúba em pares experimentais, em 1000 simulações da hipótese nula de ausência de diferença. A região de aceitação da hipótese nula para 5% de significância para teste bicaudal está delimitada em cinza."}
dplot(azt.r, svalue = azt.ei(azteca), pside="Two sided",
      main = "Teste bicaudal",
      xlab = "Estatística de interesse")
```

## Aleatorização com reposição

O dataframe `peucetia` tem os dados de um experimento de escolha de substrato
por aranhas do gênero *Peucetia*.
Vinte e sete aranhas foram mantidas em placas de Petri
cobertas com dois substratos (folhas com e sem tricomas glandulosos).
Em seis inspeções a cada placa registrou-se
se cada aranha estava sobre as folhas com tricomas.

```{r inspecionando objeto peucetia}
head(peucetia)
```
Saiba mais sobre os dados em sua página de ajuda (`?peucetia`).

### Hipótese do estudo

A hipótese do estudo é que
as aranhas preferem caçar em plantas
com pelos glandulosos, onde a captura
de presas é mais fácil.
A previsão para o experimento é
que as aranhas devem estar a maior parte
do tempo nas folhas com tricomas.
De fato, a maioria das aranhas esteve
nas folhas com tricomas em 4 ou mais inspeções:

```{r barplot peucetia, fig.cap = "Número de inspeções em que as 27 aranhas foram registradas em folhas com tricomas, em um experimento de preferência por substratos."}
## Número de inspeções em que estava em folha com tricomas
n.insp <- apply(peucetia, 1, sum)
barplot(table(factor(n.insp, levels=0:6)),
        xlab="N de inspeções em que estava na folha com tricoma",
        ylab="N de aranhas")

```

### Hipótese nula

A hipótese nula é de que não há preferência.
Como metade das placas estavam cobertas com cada
tipo de folha, a expectativa nula
é que as aranhas estivessem
na área coberta por folhas com tricomas em metade das inspeções,
em média.
Esta expectativa tem a premissa que cada inspeção
é um evento independente.

### Estatística de interesse

A cada simulação temos que calcular nossa
**estatística de interesse**, que é a
média do número de inspeções em que as aranhas estavam sobre folhas com tricomas.
Uma função para isso:

```{r estatistica de interesse peucetia}
peu.ei <- function(dataframe){
    mean(apply(dataframe, 1, sum))
}
## Valor observado
peu.ei(peucetia)
```

As aranhas foram registradas em média `r
round(peu.ei(peucetia),2)`
das 6 inspeções na área coberta por folhas com tricomas.

### Distribuição da estatística sob a hipótese nula

Para simular nossa hipótese nula criamos um
*data frame* com a mesma estrutura, em que cada
aranha esteja metade das inspeções nas folhas com tricomas

```{r peucetia H0}
peu.H0 <- matrix( rep(c(TRUE,FALSE), each = 3),
                 nrow = nrow(peucetia), ncol = ncol(peucetia), byrow=TRUE)
## Converte em data.frame
peu.H0 <- data.frame(peu.H0)
## verificando
head(peu.H0)
```
E agora simulamos a hipótese nula amostrando
com reposição cada linha [^3]:

```{r peucetia resampling, results="hide"}
peu.r <- Rsampling(type = "within_rows", dataframe = peu.H0,
                   statistics = peu.ei, ntrials = 1000, replace=TRUE)
```

O argumento `replace = TRUE`, indica amostragem com reposição.
No caso isso equivale a sortear uma posição independente
para cada aranha a cada inspeção. A probabilidade da aranha estar
na folha com tricomas é 0,5 a cada sorteio.

Uma média igual ou maior que a observada não ocorreu
na distribuição da estatística de interesse simulada:


```{r peucetia distribuicao nula, fig.cap="Distribuição do número médio de inspeções em que as aranhas estavam em folhas com tricomas, em 1000 simulações da hipótese nula de ausência de preferência por substrato. A linha vermelha indica a média observada. A região de aceitação da hipótese nula para 5% de significância está delimitada em cinza."}
dplot(peu.r, svalue = peu.ei(peucetia), pside="Greater",
      main = "Distribuição da estatística de interesse sob H0",
      xlab = "Estatística de interesse")
```

### Decisão: rejeitamos a hipótese nula?

Novamente temos um teste unicaudal, e o valor observado da estatística de interesse não
está na região de aceitação da hipótese nula (5%).
Confirmamos com o teste lógico do nosso critério de significância:

```{r peucetia teste} 
sum(peu.r >= peu.ei(peucetia))/1000 < 0.05
```

**Conclusão:** rejeita-se a hipótese nula (p < 0,05).

## Uma hipótese nula mais realista?

No exemplo anterior simulamos a hipótese
nula sorteando uma posição para cada aranha a cada
inspeção. A premissa é que a posição da aranha
em uma inspeção não afeta sua posição nas outras,
ou seja, que as inspeções são
eventos independentes.

Mas e se há uma correlação temporal na posição das aranhas?
Isso pode acontecer com aranhas que se movem a uma
frequência menor que o intervalo entre as
inspeções. Se isso é verdade, registros seguidos em um
tipo de folha podem indicar apenas tendência a ficar no
mesmo lugar, e não preferência. Nesse caso a hipótese
nula deve manter o número de inspeções em cada tipo
de folha, alterando apenas o tipo.

### Hipótese nula

A proporção das inspeções em que as aranhas
permanecem em um dos substratos não depende do
tipo de substrato (folha com ou sem tricomas).

Portanto a hipótese nula é sobre a independência entre número
de inspeções e tipo de substrato. Simulamos este cenário
embaralhando número de ocasiões entre substratos,
para cada aranha. Para isso vamos criar um *data frame*
com número de inspeções em cada substrato:

```{r peucetia n de inspeções}
## N de inspeções em folha com tricoma
tric <- apply(peucetia, 1, sum)
## N de inspeções em folha lisa
lisa <- apply(peucetia, 1, function(x) sum(x==0))
## Monta o data frame
peu.H0b <- data.frame(tric=tric, lisa = lisa)
## Primeiras linhas
head(peu.H0b)
```


### Estatística de interesse

Uma mesma estatística de interesse pode ser aplicada a
diferentes hipóteses nulas. Então mantemos a mesma do
exemplo anterior: média do número de inspeções em que
as aranhas foram registradas nas folhas com tricomas.

Mas como o *data frame* que será aleatorizado
mudou, criamos uma nova função no R para
calcular a estatística de interesse

```{r peucetia statistics 2}
peu.ei2 <- function(dataframe) mean(dataframe$tric)
## Verificando
peu.ei2(peu.H0b)
```


### Distribuição da estatística de interesse sob a hipótese nula

Simulamos a hipótese nula embaralhando as linhas
do *data frame* com número de inspeções por substrato:

```{r peucetia resampling 2, results="hide"}
peu.r2 <- Rsampling(type = "within_rows", dataframe = peu.H0b,
                   statistics = peu.ei2, ntrials = 1000)
```

A distribuição nula mudou bastante de forma, comparada com a seção anterior.
Mas uma média igual ou maior que a observada ainda foi muito rara:

```{r peucetia distribuicao nula 2, fig.cap="Distribuição do número médio de inspeções em que as aranhas estavam em folhas com tricomas, em 1000 simulações da hipótese nula de ausência de preferência por substrato, considerando tendência das aranhas permanecerem onde estão. A linha vermelha indica a média observada."}
dplot(peu.r2, svalue = peu.ei2(peu.H0b), pside="Greater",
      main = "Distribuição da estatística de interesse sob H0",
      xlab = "Estatística de interesse")
```

### Decisão: rejeitamos a hipótese nula?

O valor observado da estatística de interesse não está na região de aceitação.
Aplicando nosso critério de significância:

```{r peucetia teste 2} 
sum(peu.r2 >= peu.ei(peucetia))/1000 < 0.05
```

**Conclusão:** rejeita-se a hipótese nula (p < 0,05).


## Zeros estruturais

Em alguns conjuntos de dados há observações com frequência zero
que são considerados impossíveis de ocorrer ou de se observar. Por exemplo,
o *dataframe* `pielou` tem o número de registros de dez
espécies de pulgões em doze espécies de plantas do gênero *Solidago*.

```{r pielou inspecionando objeto}
pielou
```

Para saber mais sobre este conjunto de dados consulte a página de ajuda (`?pielou`).
Há várias ocorrências com frequência zero.
Vamos simular uma hipótese nula supondo que
essas frequências são estruturais, ou seja, que indicam
associações inseto-planta que não podem ocorrer.

### Hipótese do estudo

Nossa hipótese de estudo é que há ou houve partilha de recursos entre as espécies
de pulgões. Neste caso, as associações observadas devem ter resultado em redução da
sobreposição de nichos dos insetos.

### Hipótese nula

Nossa hipótese nula é que a sobreposição de nicho não difere da esperada
caso as ocorrências dos pulgões nas plantas sejam independentes.

### Estatística de interesse

Esses dados foram usados para exemplificar um cálculo de amplitude e sobreposição
de nichos. A expressão proposta pela autora
para a sobreposição média de nichos é a diferença entre
a índice de Brillouin de todos os valores e da soma das colunas da
tabela. O índice de Brillouin é uma medida de diversidade em uma coleção de valores $x_i$:

$$H = \frac{1}{N} \log N! \ - \ \frac{1}{N} \sum \log x_i !$$

Onde $N = \sum x_i$. Vamos criar uma função para fazer este cálculo

```{r pielou indice de brillouin}
brillouin <- function(x, base=10) {
    N <- sum(x)
    lfactorial(N)/(log(base)*N)  -  sum(lfactorial(x)/log(base))/N
}
```
E então criamos um função para o cálculos da estatística de interesse

```{r pielou estatistica brillouin}
pielou.ei <- function(dataframe)
    brillouin( dataframe ) - brillouin( apply(dataframe,2,sum) )
```
Cujo valor é

```{r pielou estatistica de interesse}
pielou.ei(pielou)
```

### Distribuição da estatística de interesse sob a hipótese nula

Para simular nossa hipótese nula, embaralhamos os números de ocorrências registradas
de cada espécie de pulgão entre as plantas. Com isso mantemos criamos uma situação em que
as hierarquias de preferência de cada espécie de pulgão são mantidas, mas tornam-se independentes.
Além disso, usamos a opção `fix.zeroes = TRUE` para indicar que os valores zero não devem
ser embaralhados.

```{r , results="hide"}
pielou.r1 <- Rsampling(type = "within_rows", dataframe = pielou,
                   statistics = pielou.ei, ntrials = 1000, fix.zeroes = TRUE)
```

O valor observado é maior do que a maioria dos valores na distribuição nula. Como nossa
hipótese é unicaudal (sobreposição observada menor do que o esperado pelo acaso),
o valor observado está na região de aceitação da hipótese nula.

```{r pielou nula 2, fig.cap="Distribuição da sobreposição média de uso de plantas hospedeiras por espécies de pulgões, em 1000 simulações da hipótese nula de independência das espécies de inseto pelas plantas. As plantas sem ocorrência observadas dos pulgões foram consideradas não disponíveis (zeros estruturais). A linha vermelha indica a média observada."}
dplot(pielou.r1, svalue = pielou.ei(pielou), pside="Lesser",
      main = "Distribuição da estatística de interesse sob H0",
      xlab = "Estatística de interesse", xlim=c(0.3,0.6))
```

### Decisão: rejeitamos a hipótese nula?

O valor observado da estatística de interesse está na região de aceitação.
Aplicando nosso critério de significância:

```{r  teste 2} 
sum(pielou.r1 <= pielou.ei(pielou))/1000 < 0.05
```

**Conclusão:** não se rejeita a hipótese nula (p > 0,05).


[^2]: Statistics.com LCC. 2009. Resampling Stats Add-in for Excel User’s Guide.
http://www.resample.com/content/software/excel/userguide/RSXLHelp.pdf

[^3]: Há maneiras mais otimizadas de fazer isso, mas esta reproduz a lógica de sorteios de uma urna do *Resampling Stats*