Programação em R: loops for – parte 2

Continuando o tema da semana anterior (se você não leu, leia antes de ler este post de hoje), hoje mostrarei como usar loops para analisar diversas variáveis automaticamente.

Por que precisaríamos fazer isso? Bom, pode ser que o objetivo do seu estudo é descobrir quais variáveis ambientais afetam uma determinada plantinha espécie. Por exemplo, pode querer avaliar como esta plantinha ou bichinho espécie é afetada por temperatura média do local em que vive, estrutura da vegetação, quantidade de floresta remanescente, etc. Neste caso você vai fazer uma única análise.

Mas pode acontecer que você quer saber como diferentes plantinhas espécies ou outras coisas variáveis são afetadas pela mesma variável explanatória. Por exemplo, no meu mestrado, eu queria saber como altura da vegetação, temperatura do ar, quantidade de gramíneas invasoras e outras variáveis são afetadas pela distância até a borda do remanescente florestal (ou savânico – trabalhei no cerrado). A variável explanatória é sempre a mesma, mas as variáveis-resposta mudam. Poderia eu fazer uma PCA ou outra análise multivariada? Poderia. Mas o meu interesse era em cada variável individualmente, não no seu conjunto, de modo que eu precisava fazer várias análises.

Outro exemplo: um conjunto de dados que ajudei a analisar consistia em várias espécies vegetais ao longo de um gradiente ambiental. Queríamos saber quais espécies respondem de forma positiva, negativa ou neutra ao gradiente – i.e. quais aumentam em abundância, diminuem em abundância, ou não respondem ao gradiente de forma alguma (“tô nem aíiiii… tô nem aíiii…).

Dá pra rodar a anaĺise manualmente espécie por espécie, e isso às vezes é bom para você conhecer melhor seus dados e conseguir entender melhor o que eles te mostram. Mas se você precisa repetir a análise porque decidiu fazer algo diferente, ou se você tem 142 espécies, isso fica algo bem tedioso de fazer. Ou se cada análise demora 5 minutos, o conjunto delas pode demorar o dia inteiro. Podemos ficar o dia inteiro no computador, ou podemos deixar ele rodando e ir dormir. Fiz muito isso. Só cuidado pra não dar curto :-)

Mas Pav, se eu fizer isso, não posso obter algumas correlações pelo acaso? – Ótima pergunta, jovem Padawan! E por isso mesmo que vou trabalhar com dados aleatórios, pra exemplificar essas benditas correlações pelo acaso.

Primeiro, vou definir uma semente (seed) para que todos os resultados sejam os mesmos. Definindo o seed, sempre que as simulações forem rodadas elas vão dar o mesmo resultado.

set.seed(28)

Vou simular uma variável aleatória com distribuição uniforme e 42 valores:

dados.x = runif(42,0,100)

Chamei de dados.x porque vai ser nossa variável x. O comando runif gera 42 valores, que variam de 0 a 100. Poderiam ser por exemplo a porcentagem de floresta ao redor do sítio amostral, ou a cobertura de dossel, ou qualquer outra coisa que varia de 0 a 100.

Agora vou gerar uma matriz de dados com 20 variáveis aleatórias com distribuição normal:

dados.y = matrix(rnorm(42*20, 20, 4),ncol=20)

O comando rnorm gera 42*20=840 valores, que seguem uma distribuição normal, com média 20 e desvio 4. O comando matrix coloca eles em uma matriz com 20 colunas. Cada coluna vai ser uma variável.

Vou juntar eles em um data.frame, porque nossos dados de campo costumam vir em data.frames:

dados = data.frame(dados.x, dados.y)

str(dados)

Eis a nossa estrutura de dados:

‘data.frame’: 42 obs. of 21 variables:
$ dados.x: num 3.08 41 2.46 79.68 86.68 …
$ X1 : num 17.9 18.3 14.8 14.7 27.1 …
$ X2 : num 27.4 27.5 16 18.7 21.1 …
$ X3 : num 16.3 20.5 16.6 23.1 12.4 …
$ X4 : num 19.7 33.3 23.9 20 15 …
$ X5 : num 21.1 16.5 20.1 25.8 17.3 …
$ X6 : num 24.4 19.6 19.1 20 18.3 …
$ X7 : num 21.4 19.6 11.6 19.4 21.9 …
$ X8 : num 23 26.2 22.2 22.9 20.5 …
$ X9 : num 24.7 17 20.7 23.6 24.4 …
$ X10 : num 27.3 17.9 19.8 12.6 16.8 …
$ X11 : num 21 17.5 15.7 15.1 18.7 …
$ X12 : num 19 19.8 22.7 14 19.4 …
$ X13 : num 28.7 22.1 15.5 23.5 23.7 …
$ X14 : num 15 16.2 18.4 14.3 21.6 …
$ X15 : num 21.8 18.6 17.7 23.5 26.2 …
$ X16 : num 18.1 22.9 22.5 28.4 21.9 …
$ X17 : num 18.3 22.4 19.3 18.6 17.9 …
$ X18 : num 14.5 32.4 21.2 22.2 15.7 …
$ X19 : num 20.3 24.9 22.5 24.8 27.9 …
$ X20 : num 22 21.8 20.1 18.8 18.9 …
Temos 42 observações (réplicas) de 21 variáveis. A primeira se chama dados.x e as outras receberam nomes de X1 até X2.

Vou renomear a primeira para “Explanatoria”, porque sim:

names(dados)[1] = "Explanatoria"

Bom, e agora o loop… Vamos primeiro fazer um loop para plotar todos os gráficos. Usaremos o comando par para criar uma janela em que todos eles aparecerão:

foo.x = dados[,1]

par(mfrow=c(5,4), mar=c(2,2,2,2), oma=c(3,3,2,2))

for(i in 2:21) {
    foo.y = dados[,i]
    foo.name = names(dados)[i]
    plot(foo.y ~ foo.x, main=foo.name, xlab="", ylab="")
}

mtext(text="Variável explanatória", side=1, outer=T)
mtext(text="Variável resposta", side=2, outer=T)

loops_for2_1.png

Primeiro criei um objeto que não vai ser relevante no futuro, com a variável X, e chamei ele de foo.x. Os nomes foo, bar, foobar são nomes que são frequentemente dados para objetos que não têm muita função além de demonstrar alguma coisa ou serem usados para algo específico, tipo dentro de um loop.

A seguir criei um espaço para 20 gráfico (5 linhas, 4 colunas), com margens de 2 linhas entre os gráficos e margens externas de 3 linhas embaixo e à esquerda e 2 linhas em cima e à direita. Margens externas é o espaço que fica ao redor do conjunto dos 20 gráficos, neste exemplo.

Dentro do loop, comecei de 2 porque as variáveis resposta estão a partir da 2a coluna. Então, em cada rodada do loop, o R vai pegar a segunda, depois a terceira, depois a quarta… até a 21a linha e criar um objeto chamado foo.y com ela. Também vai criar um objeto chamado foo.name com o nome desta variável.

A seguir, faço o gráfico de foo.y em função de foo.x, e dou a ele o nome da variável. São gráficos exploratórios, para publicação eu faria algo um pouco mais bonito. Não vou dar nomes aos eixos X e Y deles, porque o nome do eixo Y é o nome do gráfico e o eixo X é o mesmo para todos.

Finalmente, uso os comandos mtext para colocar nomes nos eixos, abaixo e à esquerda de todos os gráficos. O primeiro coloca o texto “Variável explanatória” embaixo (side=1) na parte externa (outer=T) do gráfico, o segundo coloca “Variável resposta” à esquerda (side=2) na parte externa (outer=T) do gráfico. O comando mtext serve para colocar texto nas margens da figura – podem ser as internas, ao redor de cada gráfico, ou externas, ao redor do conjunto deles.

Vou agora ajustar modelos lineares, seguindo o mesmo esquema. Só que antes vou criar uma lista onde eles ficarão armazenados:

resultados.lm = list()

foo.x = dados[,1]

for(i in 2:21) {
    foo.y = dados[,i]
    foo.name = names(dados)[i]
    resultados.lm[[i-1]] = lm(foo.y ~ foo.x)
}

A sintaxe é praticamente a mesma que para o gráfico. Só que agora, os resultados do modelo linear (lm) são guardados no objeto resultados.lm. Ou seja, o resultado da primeira iteração vai ser o primeiro elemento da lista resultados.lm, o resultado da segunda vai ser o segundo elemento, etc. Isso é definido pela parte resultados.lm[[i-1]]: o (i-1)-ésimo elemento desta lista vai ser o resultado do comando lm. Coloquei i-1 porque a contagem começa a partir de 2, então sem este -1 o primeiro elemento ia ficar vazio.

Podemos extrair coeficientes deste elemento. Tem diversas maneiras; eu vou usar uma mais complexa que também é mais flexível.

Como extrair informações de um único modelo? Digamos que eu quero saber o R-quadrado e o p (significância). Para ver isso em um modelo, eu usaria o comando summary:

summary(resultados.lm[[1]])

Entre outras coisas, ele mostra que o R2 múltiplo é de 0.01 e o p-valor associado à nossa variável explanatória, que recebeu aqui o nome de foo.x, é de 0.499. Mas e pra eu estrair isso?

Se eu olhar a estrutura do resultado do summary, terei algo assim:

str(summary(resultados.lm[[1]]))

List of 11
$ call : language lm(formula = foo.y ~ foo.x)
$ terms :Classes 'terms', 'formula' language foo.y ~ foo.x
.. ..- attr(*, "variables")= language list(foo.y, foo.x)
.. ..- attr(*, "factors")= int [1:2, 1] 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:2] "foo.y" "foo.x"
.. .. .. ..$ : chr "foo.x"
.. ..- attr(*, "term.labels")= chr "foo.x"
.. ..- attr(*, "order")= int 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
.. ..- attr(*, "predvars")= language list(foo.y, foo.x)
.. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
.. .. ..- attr(*, "names")= chr [1:2] "foo.y" "foo.x"
$ residuals : Named num [1:42] -2.56 -2.2 3.48 -4.98 -1.56 ...
..- attr(*, "names")= chr [1:42] "1" "2" "3" "4" ...
$ coefficients : num [1:2, 1:4] 20.1204 -0.0153 1.294 0.0223 15.5485 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "(Intercept)" "foo.x"
.. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
$ aliased : Named logi [1:2] FALSE FALSE
..- attr(*, "names")= chr [1:2] "(Intercept)" "foo.x"
$ sigma : num 4.14
$ df : int [1:3] 2 40 2
$ r.squared : num 0.0115
$ adj.r.squared: num -0.0132
$ fstatistic : Named num [1:3] 0.466 1 40
..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
$ cov.unscaled : num [1:2, 1:2] 9.75e-02 -1.46e-03 -1.46e-03 2.91e-05
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "(Intercept)" "foo.x"
.. ..$ : chr [1:2] "(Intercept)" "foo.x"
- attr(*, "class")= chr "summary.lm"

Este summary tem um moooonte de coisas. Inclusive uma parte contendo os coeficientes – onde está o p-valor – e também uma coisa chamada r.squared. Então eu posso pegar esses objetos diretamente:

summary(resultados.lm[[1]])$coefficients

Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.1204277 1.29404029 15.5485326 1.490934e-18
foo.x -0.0152642 0.02234891 -0.6829955 4.985467e-01

O p-valor que queremos está na segunda linha, quarta coluna:

summary(resultados.lm[[1]])$coefficients[2,4]

E temos aqui o nosso valor de 0.4985467.

Para o R2 é mais simples:

summary(resultados.lm[[1]])$r.squared

Aparece o valor de 0.01.

Parece complicado, e é um tanto complicadinho mesmo… Se soberem de formas mais simples de fazer isso, comentem! A vantagem de fazer assim é que consigo usar o mesmo esquema em modelos generalizados, modelos aditivos, modelos aditivos mistos… É uma forma bem flexível, embora provavelmente não a mais recomendada.

Com isso, posso fazer um loop para extrair essa informação para todos os modelos. Vou criar objetos (vetores) onde armazenar essas informações e fazer um loop para guardá-las:

pvalues = numeric(20)
r2s = numeric(20)

for(i in 1:20) {
    pvalues[i] = summary(resultados.lm[[i]])$coefficients[2,4]
    r2s[i] = summary(resultados.lm[[1]])$r.squared
}

O comando numeric(20) cria um vetor numérico com espaço para 20 números. E depois, dentro do loop (que vai de 1 a 20 porque temos 20 variáveis-resposta), definimos que o i-ésimo elemento de p-values vai ser o valor retirado do summary do i-ésimo modelo. Mesma coisa para o i-ésimo valor do r2s.

Olhando o objeto pvalues, tivemos o seguinte:

pvalues
[1] 0.49854667 0.60865682 0.92600464 0.50774634 0.04738943 0.90978881
[7] 0.15829905 0.77186511 0.88435453 0.77300637 0.34370308 0.24061792
[13] 0.29982024 0.97769126 0.51966803 0.98008657 0.07643105 0.94369858
[19] 0.42453532 0.28747774

Reparem que tivemos um p abaixo de 0.05 (no caso, de 0.047). Se definirmos como significativa a relação que tiver um p<0.05 (pra mim não faz muito sentido usar critérios arbitrários assim, a não ser que estejamos trabalhando com muitas variáveis resposta e usar algum critério acaba sendo necessário), teríamos encontrado uma relação significatiav pelo simples acaso. E isso, jovens Padawans, é porque não devemos coletar um monte de variáveis pra ver se dá alguma coisa – porque vai dar alguma coisa, mas essa alguma coisa não vai significar nada. Neste caso, se fossem dados reais, eu diria que provavelmente este resultado significativo veio pelo acaso. Mas se eu tivesse, digamos, 7 relações significativas, eu concluiria que provavelmente tem alguma coisa aí. Se tivesse uma relação com p muito baixo, da ordem de 0.00001, também concluiria que provavelmente tem alguma coisa aí.

Bom, e agora, vou fazer o gráfico mostrando as relações obtidas pelo modelo. Mas, mais do que isso, vou definir que a linha vai ser uma linha pontilhada se p>=0.05 e linha contínha se p<0.05. Isso é um jeito prático de diferenciar relações que foram mostradas significativas na análises das que não foram.

foo.x = dados[,1]

par(mfrow=c(5,4), mar=c(2,2,2,2), oma=c(3,3,2,2))

for(i in 2:21) {
    foo.y = dados[,i]
    foo.name = names(dados)[i]
    plot(foo.y ~ foo.x, main=foo.name, xlab="", ylab="")
    abline(resultados.lm[[i-1]], lty=ifelse(pvalues[i-1]<0.05,1,2),lwd=1.5)
}

mtext(text="Variável explanatória", side=1, outer=T)
mtext(text="Variável resposta", side=2, outer=T)

A diferença em relação ao primeiro gráfico é que adicionei o comando abline, que adiciona um linha reta. Ele pode adicionar este linha reta a partir do modelo – que sempre vai ser o modelo de número i-1 do objeto resultados.lm (porque a contagem começa de 2, como expliquei acima). O argumento lty define o tipo da linha. No caso, se o (i-1)-ésimo valor do objeto p-values for menor do que 0.05, a linha vai ser do tipo 1, ou contínua; caso contrário, vai ser do tipo 2, ou pontilhada. O argumento lwd=1.5 deixa a linha um pouco mais grossinha. Éis a figura, mostrando uma relação com p<0.05 em um dos 20 gráficos:

loops_for2_2.png

E é isso. Agora se você precisa avaliar como um monte de espécies respondem a algum gradiente ambiental, já sabe que pode fazer um loop para isso. Se esta abordagem é recomendada ou não, já é outra questão. No fundo, a questão é mais como você vai interpretar e discutir os resultados do que a estatística utilizada, e o objetivo deste post era dar mais um exemplo do uso de for. Semana que vem mostrarei como abrir vários arquivos simultaneamente, novamente usando for.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair /  Alterar )

Foto do Google

Você está comentando utilizando sua conta Google. Sair /  Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair /  Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair /  Alterar )

Conectando a %s