RRRR: Reproduzindo o Resampling com Rsampling para Regressões

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)

Exemplos de regressão

O dataframe rhyzophora tem medidas de árvores de mangue em solos lodosos mais e menos instáveis.

> head(rhyzophora)
  soil.instability canopy.trunk     root n.roots
1             high     937.6087 19.78200     117
2             high    1957.0434 23.66775     152
3             high    1403.0969 20.72400     106
4             high     785.6869  6.73530      91
5             high    1431.3668 15.70000     161
6             high    1208.7816 16.97170     125
> summary(rhyzophora)
 soil.instability  canopy.trunk         root            n.roots      
 high  :12        Min.   : 279.9   Min.   : 0.5212   Min.   : 19.00  
 medium:12        1st Qu.: 839.7   1st Qu.: 5.9357   1st Qu.: 54.50  
                  Median :1067.8   Median :14.7619   Median : 88.50  
                  Mean   :1179.9   Mean   :12.7651   Mean   : 87.25  
                  3rd Qu.:1408.9   3rd Qu.:17.3053   3rd Qu.:110.25  
                  Max.   :3675.4   Max.   :40.5531   Max.   :172.00

Saiba mais sobre os dados em sua página de ajuda (?rhyzophora).

Hipótese do estudo

A hipótese é que árvores em solos mais instáveis investem mais em estruturas de sustentação. Uma previsão é que a relação entre o torque da árvore e o investimento em raízes de sustentação deve ser diferente nos dois tipos de solo. Para representar o torque foi usada a razão entre a a área da copa e do tronco. O investimento em raízes foi expresso em número de raízes de sustentação e a área coberta por elas.

Os dados sugerem uma relação positiva entre a variável de torque e o número de raízes. Também parece que os pontos das árvores amostradas nos dois tipos de solo separam-se, sugerindo uma relação diferente:

> plot(n.roots ~ canopy.trunk, data=rhyzophora, type="n",
+      xlab="área copa / área tronco", ylab="número de raízes")
> points(n.roots ~ canopy.trunk, data=rhyzophora,
+        subset=soil.instability=="medium")
> points(n.roots ~ canopy.trunk, data=rhyzophora,
+        subset=soil.instability=="high", pch=19)
> legend("topright", c("Média","Alta"), title="Instabilidade do solo", pch=c(1,19))

Relação entre o número de raízes de sustentação e razão área da copa / área do tronco em árvores de mangue em solos mais e menos instáveis.

Embaralhando linhas dentro de estratos

Hipótese nula

Para ilustrar randomizações restritas a estratos vamos testar a hipótese nula mais básica de que que não há relação em nenhum dos dois tipos de solos. Simulamos isso embaralhando os valores da variável de torque entre árvores de cada tipo de solo.

Estatística de interesse

Temos uma estatística de interesse para cada solo, que são as inclinações das regressões lineares:

> rhyz.ei <- function(dataframe){
+     m1 <- lm(n.roots ~ canopy.trunk, data=dataframe,
+              subset=soil.instability=="medium")
+     m2 <- lm(n.roots ~ canopy.trunk, data=dataframe,
+              subset=soil.instability=="high")
+     c(med = coef(m1)[[2]],
+       high = coef(m2)[[2]])
+ }
> ## Valore observados
> rhyz.ei(rhyzophora)
       med       high 
0.01813085 0.05570843

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

Simulamos a hipótese nula de ausência de relação embaralhando os valores da variável de torque entre árvores do mesmo tipo de solo:

> rhyz.r <- Rsampling(type = "normal_rand", dataframe = rhyzophora,
+                     statistics = rhyz.ei, stratum = rhyzophora$soil.instability,
+                         cols = 2, ntrials = 1000)

O argumento stratum = rhyzophora$soil.instability, que indica que o embaralhamento da coluna 2 deve ser feito dentro de cada tipo de solo.

Como há mais de uma estatística de interesse, a função Rsampling retorna uma matriz em cada linha é uma estatística e as colunas são as repetições

> rhyz.r[,1:3]
            [,1]        [,2]        [,3]
med  -0.01012567 -0.01486993 0.006809489
high  0.04071075 -0.01222508 0.047974435

Valores iguais ou maiores que as inclinações observadas parecem bem raros na distribuição de valores sob a hipótese nula:

> par(mfrow=c(1,2))
> dplot(rhyz.r[1,], svalue=rhyz.ei(rhyzophora)[1], pside="Greater",
+       main="Média instabilidade", xlab="Inclinações sob H0")
> dplot(rhyz.r[2,], svalue=rhyz.ei(rhyzophora)[2], pside="Greater",
+       main="Alta instabilidade", xlab="Inclinações sob H0")

Distribuição das inclinações da regressão linear do número de raízes em função da razão das áreas da copa e tronco, em 1000 simulações da hipótese nula de ausência de relação. As linhas vermelhas indicam as inclinações observadas. A região de aceitação da hipótese nula a 5% está em cinza. Em laranja o número de valores da distribuição nula maiores que os observados.

> par(mfrow=c(1,1))

Decisão: rejeitamos a hipótese nula?

As inclinações observadas para os dois grupos estão fora da região de aceitação da hipótese nula unicaudal ¹ a 5% de significância. Podemos verificar isso com um teste lógico aplicado a cada estatística de interesse:

> sum(rhyz.r[1,] >= rhyz.ei(rhyzophora)[1])/1000 < 0.05
[1] TRUE
> sum(rhyz.r[2,] >= rhyz.ei(rhyzophora)[2])/1000 < 0.05
[1] TRUE

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

Comparação das inclinações

A hipótese principal do estudo é que a relação entre torque e sustentação é diferente nos dois tipos de solo. Supondo que a relação linear existe, ela pode diferir quanto à inclinação ou intercepto.

Hipótese nula

Começamos testando a hipótese nula de que a inclinação das regressões lineares não difere entre solos.

Estatística de interesse

A estatística de interesse é a diferença entre as inclinações, que parece pequena:

> rhyz.ei2 <- function(dataframe){
+     m1 <- lm(n.roots ~ canopy.trunk, data=dataframe,
+              subset=soil.instability=="medium")
+     m2 <- lm(n.roots ~ canopy.trunk, data=dataframe,
+              subset=soil.instability=="high")
+     coef(m1)[[2]] - coef(m2)[[2]]
+ }
> ## Valores observados
> rhyz.ei2(rhyzophora)
[1] -0.03757759

Simulação da hipótese nula

Simulamos a nova hipótese nula embaralhando as árvores entre os tipos de solos (primeira coluna da tabela de dados):

> rhyz.r2 <- Rsampling(type = "normal_rand", dataframe = rhyzophora,
+                     statistics = rhyz.ei2,
+                         cols = 1, ntrials = 1000)

Decisão: rejeitamos a hipótese nula?

Neste caso não podemos descartar a hipótese nula:

> sum(rhyz.r2 > rhyz.ei2(rhyzophora))/1000 < 0.05
[1] FALSE

Comparando interceptos

Decidimos aceitar a hipótese nula de que as inclinações são iguais. A interpretação biológica disso é que nos dois tipos de solo o número de raízes de sustentação segue a mesma relação de proporcionalidade com a variável de torque.

Este fator de proporcionalidade é a inclinação da regressão linear aplicada a todas as árvores, que estimamos ajustando a regressão:

> lm(n.roots ~ canopy.trunk, data=rhyzophora)

Call:
lm(formula = n.roots ~ canopy.trunk, data = rhyzophora)

Coefficients:
 (Intercept)  canopy.trunk  
    66.30197       0.01775

Ou seja, a cada aumento de 100 unidades da variável de torque em média 1.8 raízes são adicionadas.

Note que esta proporcionalidade se mantém se adicionarmos qualquer constante. Por isso o modelo linear é expresso por

E[Y] = α + βX

Em que E[Y] é o valor esperado da resposta (número de raízes), β é a inclinação ou fator de proporcionalidade, e X a variável preditora (torque). O intercepto α não altera a proporcionalidade, apenas desloca a reta mais para cima ou mais para baixo.

Ou seja, retas com a mesma inclinação mas interceptos diferentes são paralelas. No nosso caso isso expressaria que árvores com mesmo valor da razão copa/troco sempre têm mais raízes em um dos tipos de solo.

Hipótese nula

Nossa hipótese nula é que os interceptos das regressões lineares não diferem entre os tipos de solo. Se isso é verdade a regressão linear ajustada a todos os dados deve prever bem os valores da resposta. Se não for verdade os pontos de um tipo de solo tenderão a ficar abaixo da reta, enquanto os do outro tipo de solo tenderão a ficar acima.

Já ajustamos essa regressão acima, e podemos adicionar a reta ao gráfico:

> plot(n.roots ~ canopy.trunk, data=rhyzophora, type="n",
+      xlab="área copa / área tronco", ylab="número de raízes")
> points(n.roots ~ canopy.trunk, data=rhyzophora,
+        subset=soil.instability=="medium")
> points(n.roots ~ canopy.trunk, data=rhyzophora,
+        subset=soil.instability=="high", pch=19)
> abline(lm(n.roots ~ canopy.trunk, data=rhyzophora))
> legend("topright", c("Média","Alta"), title="Instabilidade do solo", pch=c(1,19))

Relação entre o número de raízes de sustentação e razão área da copa / área do tronco em árvores de mangue em solos mais e menos instáveis. A reta é a regressão linear ajustada a todos os pontos.

Parece que de fato esta regressão subestima o número de raízes das árvores amostradas no solo mais instável e faz o oposto para as árvores do solo menos instável. Isso faz com que os resíduos desta regressão sejam positivos para árvores do solo instável e negativos para as outras.

Estatística de interesse

Nossa estatística de interesse é a diferença da médias dos resíduos das árvores em cada tipo de solo. Os resíduos são calculados da regressão aplicada a todos os dados:

> rhyz.ei3 <- function(dataframe){
+     m1 <- lm(n.roots ~ canopy.trunk, data=dataframe)
+     res.media <- tapply(resid(m1), dataframe$soil.instability, mean)
+     res.media[[1]] - res.media[[2]]
+ }
> ## Valores observados
> rhyz.ei3(rhyzophora)
[1] 51.60582

Simulação da hipótese nula

Simulamos a nova hipótese nula do mesmo jeito: embaralhando as árvores entre os tipos de solos (primeira coluna da tabela de dados).

> rhyz.r3 <- Rsampling(type = "normal_rand", dataframe = rhyzophora,
+                     statistics = rhyz.ei3,
+                         cols = 1, ntrials = 1000)

Decisão: rejeitamos a hipótese nula?

Descartamos a hipótese nula:

> sum(rhyz.r3 > rhyz.ei3(rhyzophora))/1000 < 0.05
[1] TRUE

Portanto há um intercepto para cada tipo de solo. Podemos estimá-los incluindo o efeito de solo no ajuste da regressão ²:

> (rhyz.ancova <- lm(n.roots ~ soil.instability + canopy.trunk  -1,
+                    data=rhyzophora))

Call:
lm(formula = n.roots ~ soil.instability + canopy.trunk - 1, data = rhyzophora)

Coefficients:
  soil.instabilityhigh  soil.instabilitymedium            canopy.trunk  
              83.25788                29.10476                 0.02633

E adicionamos as retas ao gráfico:

> cfs <- coef(rhyz.ancova)
> plot(n.roots ~ canopy.trunk, data=rhyzophora, type="n",
+      xlab="área copa / área tronco", ylab="número de raízes")
> points(n.roots ~ canopy.trunk, data=rhyzophora,
+        subset=soil.instability=="medium", col="blue")
> points(n.roots ~ canopy.trunk, data=rhyzophora,
+        subset=soil.instability=="high", col="red")
> abline(cfs[1],cfs[3], col="red")
> abline(cfs[2],cfs[3], col="blue")
> legend("topright", c("Média","Alta"), title="Instabilidade do solo", col=c("blue", "red"))

Relação entre o número de raízes de sustentação e razão área da copa / área do tronco em árvores de mangue em solos mais e menos instáveis. As retas são regressões lineares de mesma inclinação mas interceptos diferentes para cada tipo de solo.

Como não faz sentido neste caso esperar que o número de raízes diminua com a variável de torque fizemos um teste unicaudal.↩︎
Detalhe técnico: Acrescentamos o termo -1 na fórmula da regressão para indicar ao R que queremos as estimativas de cada intercepto. Caso contrário teríamos a estimativa de um intercepto e da diferença dele em relação ao outro.↩︎

- Instalação
- Exemplos de regressão