RRR: Reproduzindo o Resampling com Rsampling

In Rsampling: Ports the Workflow of "Resampling Stats" Add-in to R

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

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:

library(devtools)
install_github(repo = 'lageIBUSP/Rsampling')

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

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

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:

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:

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:

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:

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:

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.

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:

splot(azteca$extract.new, azteca$extract.old,
           groups.names=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:

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:

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:

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:

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:

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.

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:

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

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

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]:

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:

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:

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:

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

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:

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:

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:

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.

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

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

pielou.ei <- function(dataframe)
    brillouin( dataframe ) - brillouin( apply(dataframe,2,sum) )

Cujo valor é

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.

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.

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:

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

Rsampling documentation built on May 2, 2019, 2:09 a.m.