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

Digamos que você precisa repetir um procedimento no computador umas cinco vezes – talvez rodar a mesma análise sobre cinco conjuntos de dados. O que você faz?

Provavelmente repete o procedimento cinco vezes e faz fazer outra coisa da vida, né? :-)

Mas digamos que você precisa repetir um procedimento vinte ou trinta vezes. E agora?

Talvez a tendência seja pensar “Que tédio! Mas vamos lá né.”, repetir o procedimento e depois, talvez já com um certo grau de irritação, ir fazer outra coisa da vida.

E se forem cem vezes? Dá pra passar um dia rodando análises de forma repetida, mas será este o melhor investimento do nosso tempo?

Neste post vou mostrar com usar loops (ou seja, procedimentos repetitivos) em R para automatizar uma tarefa. Especificamente, vou mostrar como usar a estrutura for para repetir um procedimento um determinado número de vezes. Darei dois exemplos: um gráfico de bolinhas que não serve pra nada a não ser demonstrar loops; e uma análise por permutações. Semana que vem darei mais dois exemplos: uma regressão com diversas variáveis-resposta no mesmo objeto; e uma análise de variáveis que estão em arquivos diferentes no computador.

(Este post assume uma familiaridade com o R. Se não a tens, tenho outros a respeito neste blog – tipo este e este – e tem muitos tutoriais por aí!)

Loops não são a única forma de repetir um procedimento. Funções do tipo apply muitas vezes são mais recomendadas – por exemplo, para calcular as médias das colunas em uma matriz, é bem mais simples usar apply. Mas loops são mais flexíveis e permitem fazer coisas mais complexas – que até poderiam ser feitas usando apply, mas talvez de um jeito mais complicado. Digo isso porque já vi muitas críticas ao uso de loops em R, e muita coisa é mais simples de fazer de outro jeito – mas nem tudo, jovens Padawans, nem tudo.

Mas o que são loops? Exemplo com bolinhas aletórias

Loops, ou repetições, ou iterações, são precedimentos que se repetem várias (e várias e várias (e várias)) vezes, e podem retornar algum valor, por exemplo o resultado de um teste. Para começar com um exemplo artificial, vamos fazer gráficos de bolinhas aleatórias!

(Recomendo usarem o R, não o R Studio, para essa prática, porque às vezes o R Studio mostra os gráficos e outras coisas apenas depois que o loop termina, e não é o que quero aqui).

Para fazer um gráfico de bolinhas aleatórias, posso fazer o seguinte procedimento:

x = runif(15)
y = runif(15)
plot (y~x)

loopsFor_fig0.png

Criei 15 valores de x e 15 valores de Y e plotei uns contra os outros; o comando runif cria valores de uma distribuição aleatória uniforme, que no caso variam de 0 a 1. O gráfico de vocês deve ter ficado diferentes, afinal, são valores aleatórios!

Poderia ser feito mais simples também:

plot((runif(15)) ~ (runif(15)) )

Mas vou deixar do jeito mais complexo porque esse jeito pode ser expandido com mais facilidade.

E se eu quiser repetir isso várias vezes? Por exemplo, 100 vezes?

Posso fazer assim:

for(i in 1:100) {
    x = runif(15)
    y = runif(15)
    plot (y~x)
}

Rodando este código, vemos uma série de gráficos passando na velocidade da luz. rs

Reparem que o que está entre { } ficou alinhado mais à frente do que está fora do { } (se o wordpress não desconfigurar tudo, rs). Isso é muito útil, pois mostra qual é a parte que será repetida 100 vezes.

A expressão 1:100 cria uma sequência: 1, 2, 3, …, 99, 100. Para cada valor desta sequência, será criado um objeto, que eu chamei de i (porque sim), que irá assumir os valores de 1, 2, 3, …, 99, 100. Ou seja, é como se fosse feito o seguinte:

i = 1
x = runif(15)
y = runif(15)
plot(y~x)
i = 2
x = runif(15)
y = runif(15)
plot(y~x)
i = 3
x = runif(15)
y = runif(15)
plot(y~x)
...
i = 100
x = runif(15)
y = runif(15)
plot(y~x)

Mas Pav, pra que serve este i? Só pra complicar a vida?

Bom, também ele na verdade é muito importante. E ele poderia receber qualquer nome; eu gosto de usar letras – i, j, k, l, m… Acho que isso é relativamente comum entre quem programa. E a utilidade disso? Bom, digamos que queremos salvar estas figuras, ou seja, ter 100 arquivos de figuras aleatórias. Como fazer isso?

Nós poderíamos chamá-los de “fig_runif_01.png”, “fig_runif_02.png”… Mas escrever 100 nomes assim vai ser tedioso. Então podemos usar a função paste:

i = 1
paste("fig_runif_",i,".png", sep="")

O resultado é “fig_runif_1.png”. O argumento sep=”” serve para não termos espaços entre as coisas.

Mas é melhor ainda fazermos assim:

i = 1
paste("fig_runif_",formatC(i, digits=2, flag=0, format="d"),".png",sep="")

Este argumento, formatC, faz com que sempre tenhamos um número com ao menos três dígitos (digits=2 – ao menos dois dígitos além do número) – ou seja, 001 ao invés de 1 e 015 ao invés de 15 – isso é importante para que arquivos fiquem na ordem certa no computador; flag=0 significa que teremos “007” ao invés de ” 7″ (sem este argumento, temos espaços no lugar dos zeros, e já pensou como seria se James Bond fosse o Agente Espaço Espaço Sete?), format=”d” significa que estamos lidando com números inteiros.

Vamos criar uma pasta onde essas figuras vão ficar e salvar as figuras lá dentro, senão seu computador pode ficar todo bagunçado (ou mais bagunçado, rs):

dir.create("praticaR_loops_for")
setwd("praticaR_loops_for")

E vamos fazer um loop para criar essas figuras todas aí:

for(i in 1:100) {
    x = runif(15)
    y = runif(15)
    cores = sample(1:10, 1)
    nome <- paste("fig_runif_", formatC(i, digits=2, flag=0, format="d"), ".png", sep="")

    png(filename=nome, height=20, width=20, res=300, unit="cm")
    plot(y ~ x, col=cores, pch=16)
    dev.off()
}

Demora um pouquinho, então podemos incluir um contador pra você saber quanto já foi e quanto falta:

for(i in 1:100) {
    x = runif(15)
    y = runif(15)
    cores = sample(1:10, 1)
    nome <- paste("fig_runif_", formatC(i, digits=2, flag=0, format="d"), ".png", sep="")

    png(filename=nome, height=20, width=20, res=300, unit="cm")
    plot(y ~ x, col=cores, pch=16, cex=4)
    dev.off()

    print(i)
}

O que este comando todo significa? Bom, vamos por partes

Desta vez criei o objeto “cores”, é um número aleatório entre 1 e 10 para cada iteração (ou rodada) do loop. Apenas pra deixar os gráficos mais coloridos e bunitenhos.

Também criei o objeto “nome” – e é aí que o objeto i, que assume valores de 1 a 100, mostra sua importância pela primeira vez! Sem este objeto, os arquivos teriam todos o mesmo nome. Mas com este objeto, temos esta lindeza aqui:

loopsFor_fig1.png

Depois usei o comando png. Este comando faz com que o gráfico, ao invés de ser plotado na tela, vai ser salvo diretamente em um arquivo PNG – que é o formato que eu mais gosto pra gráficos, por ser quase tão leve quanto um GIF e tão flexível quanto um JPG. Não usem JPG para gráficos, pois vai perder qualidade; JPG é bom pra fotos, que ficariam pesadas demais em outro formato. Os argumento do comando png incluem: filename – o nome do arquivo; height e width – a altura e altura da imagem a ser criada; res – a resolução, em DPI; e unit – a unidade em que a altura e largura são medidas, no caso usei cm.

Dentro do plot, é a mesma coisa que antes, mas adicionei o argumento col=cores para que as bolinhas fiquem coloridas, pch=16 para serem bolinhas preenchidas, e cex=4, que deixa as bolinhas quatro vezes maiores do que antes (cex vem de character expansion).

Depois tem o comando dev.off(). Este comando é muito importante. Ele fala pro R “beleza, terminamos com esse arquivo PNG aqui, fecha ele e bora pro próximo!”. Se não usarmos o dev.off, teremos um monte de arquivos sem nenhuma imagem neles. O dev.off é tipo um comando para salvar aquilo que foi feito.

Finalmente, print(i) é para mostrar na tela quantas iterações ou rodadas já passaram. Às vezes fazemos análises em que cada iteração demora cinco minutos ou mais – nestes casos, ter este contador pode facilitar bastante a nossa vida.

Testes por permutações

Agora que já sabemos como um loop funciona, vamos criar um teste por permutações. A ideia de um teste por permutações é simularmos o que poderia acontecer na hipótese nula, e com base nisso calcularmos a significância, ou seja, a probabilidade de que um resultado observado nos nossos dados poderia ter sido observado se a hipótese nula for verdadeira.

Vamos criar dois conjuntos de observações, cada um com 42 pontos. Os dois conjuntos vão ter distribuição normal; o primeiro terá média 5 e desvio-padrão 1, o segundo terá média 6 e desvio-padrão 2:

dados1 = rnorm(42, 5, 1)
dados2 = rnorm(42, 6, 1)

O comando rnorm cria valores que vêm de uma distribuição normal. O primeiro argumento é o número de valores a serem criados, o segundo é a média e o terceiro é o desvio-padrão.

Vamos ver como é a distribuição destes pontos?

par(mfrow=c(2,1), mar=c(3,3,2,2))
hist(dados1, breaks=11)
hist(dados2, breaks=11)

loopsFor_fig2.png

O comando par especifica como vai ser o gráfico. No caso, o argumento mfrow=c(2,1) significa que teremos dois gráficos, distribuiídos em duas linhas e uma coluna. Sempre linhas antes de colunas. O argumento mar=c(3,3,2,2) significa que as margens (ou seja, espaço ao redor de cada gráfico) inferior, à esquerda, superior e à direta terão respectivamente 3, 3, 2 e 2 unidades. Sempre baixo-esquerda-cima-direta.

O comando hist cria um histograma, e breaks define quantas divisões teremos nele.

Um problema é que o eixo X nos dois gráficos é diferente, dificultando a comparação. Para resolver isso, podemos fazer assim:

dados.junto = c(dados1, dados2)
amplitude = range(dados.junto)

Juntei os dois conjuntos de dados em um único objeto, e usei o comando range para ver de quanto a quanto estes valores variam. No meu caso, é de 2.5 até 8.6, aproximadamente.

Refazendo o gráfico:

par(mfrow=c(2,1), mar=c(3,3,2,2))
hist(dados1, breaks=11, xlim=amplitude)
hist(dados2, breaks=11, xlim=amplitude)

loopsFor_fig3.png

Interessante como dados criados a partir de uma distribuição normal nem sempre se assemelham a uma distribuição normal né? :-)

E pelos gráficos parece que a segunda amostra realmente está mais deslocada para à direita (tem mais valores maiores). Se estes fossem dados reais, iríamos querer talvez saber se a diferença é significativa ou se esta diferença apareceu pelo acaso.

Só que esta pergunta – se duas amostras são diferentes – na verdade não faz muito sentido. É claro que elas são diferentes – em algum aspecto! Os valores não são exatamente os mesmos. Mas os valores não serem exatamente os mesmos não é relevante. O relevante é aquilo que está por detrás dos valores.

Uma pergunta mais informativa seria: as médias das duas amostras são signficativamente diferentes? Ou seja, qual a probabilidade delas terem vindo da mesma população estatística, ou de duas populações com a mesma média, e a diferença observada entre as médias se deve ao acaso?

É esta pergunta que um teste t de Student nos responderia. E é esta pergunta que responderemos usando um teste de permutações.

Primeiro, calculamos a diferença entre as duas médias:

mean(dados1) - mean(dados2)

A minha diferença foi de -1.17; a de vocês deve ter sido diferente, afinal, estamos trabalhando com números aleatórios.

Vamos armazenar esta diferença:

dif.real = mean(dados1) - mean(dados2)

Chamei de “dif.real” porque é a diferença observada com nossos dados “reais”, sem aleatorizar.

Como simular a hipótese nula? Bom, podemos juntar os nossos dados e pegar aleatoriamente 42 deles para representarem a amostra 1 e os 42 restantes para representar a amostra 2.

dados.junto = c(dados1, dados2)
dados.junto.perm = sample(dados.junto)
dados1.perm = dados.junto.perm[1:42]
dados2.perm = dados.junto.perm[43:84]

Novamente juntei os dados – sei que já tinha feito isso antes, mas acho mais prático quando tudo que for relevante para o loop estiver próximo. Depois aleatorizei a ordem deles usando o comando sample. Peguei os primeiros 42 valores para corresponderem ao dados1 e os restantes como representando o conjunto 2.

Adicionei “.perm” no nome dos objetos para indicar que são dados que foram sujeitos a permutação, ou seja, cujos números foram colocados em ordem aleatória.

Podemos ver como ficaria o gráfico na hipótese nula:

par(mfrow=c(2,1), mar=c(3,3,2,2))
hist(dados1.perm, breaks=11, xlim=amplitude)
hist(dados2.perm, breaks=11, xlim=amplitude)

loopsFor_fig4.png

E a diferença entre as médias?

mean(dados1.perm) - mean(dados2.perm)

A diferença é de 0.12 nos meus dados – bem pequena.

Em um teste por permutações, repetiríamos este procedimento um grande número de vezes. Tipo 5000 vezes. Só que a ideia de um teste por permutações é que a diferença real também poderia ser observada por acaso – então incluímos ela entre os 5000 valores. De modo que teremos a diferença real e 4999 diferenças observadas. Por que 5000 e não 10000? Pode ser 10000 sem nenhum problema; é até melhor. Mas tem uns estudos que mostram que 5000 costuma ser suficiente para uma significância a 1%, ou algo assim. E por que 5000 e não 5010, por exemplo? Sei lá, pra ficar bonitinho, uai. rs

Vamos entar criar um objeto, que chamarei de “Nperm”, que vai dizer quantas vezes o procedimento será repetido.

Nperm = 5000

E agora precisamos criar um objeto onde as diferenças calculadas com os dados aleatórios serão armazenadas. Chamarei ele de “dif.perm”:

dif.perm = numeric(Nperm)

Este comando diz que o objeto “dif.perm” é um vetor numérico, com espaço para Nperm (ou seja, 5000) números.

O primeiro valor dele será a diferença real:

dif.perm[1] = dif.real

E agora fazemos o loop!

for (i in 2:Nperm) {
    dados.junto.perm = sample(dados.junto)
    dados1.perm = dados.junto.perm[1:42]
    dados2.perm = dados.junto.perm[43:84]
    dif.perm[i] = mean(dados1.perm) - mean(dados2.perm)
}

Repararam que comecei a contagem de 2? É porque o primeiro valor é a diferença real.

Dentro do loop, faço a aleatorização, calculo a diferença entre as médias, e guardo ela no vetor “dif.perm”. Cada diferença é guardada como um valor deste loop. Não usei o comando print(i) porque o procedimento é ridiculamente rápido.

E aqui o objeto i é mais importante ainda do que antes. Na primeira rodada, i=2 – então o segundo número vetor “dif.perm” vai ser a primeira diferença simulada na hipótese nula. Depois, i=3 – então a próxima diferença calculada vai para o terceiro número deste vetor. E isso até chegar a i=5000, quando o vetor inteiro terá sido preenchido. Sem o i, não teríamos como armazenar estes valores.

Para visualizar o resultado, podemos fazer um histograma dos valores simulados e colocar a diferença real nele:

hist(dif.perm)
abline(v=dif.real, lwd=2, col="red")

loopsFor_fig5.png

O comando abline adiciona uma linha ao gráfico; abline(v=dif.real) adiciona uma linha no valor da diferença real. lwd=2 define uma linha mais grossa, col=”red” faz com que a linha seja vermelha.

Reparem que a diferença observada está bem na extremidade do gráfico – ou seja, observar uma diferença como essa se a hipótese nula for correta seria bem difícil.

Podemos calcular a significância como a proporção de diferenças absolutas que foram maiores ou iguais à diferença absoluta observada:

signif = sum(abs(dif.perm) >= abs(dif.real)) / Nperm

No caso, a significância foi de 0.0002, ou seja, 1 em 5000 – o valor mais significativo possível de ser observado com 4999 permutações junto com o valor real. A probabilidade da hipótese nula estar correta é realmente muito, mas muito pequena.

Por hoje é só, pessoal!

3 pensamentos sobre “Programação em R: loops for – parte 1

  1. Pingback: Programação em R: loops for – parte 2 – Mais Um Blog de Ecologia e Estatística

  2. Essa dica do formato PNG foi muito útil! Resolve muita dor de cabeça na hr de modificar a resolução e o formato dos gráficos para se enquadrar nos temos das revistas…

    Curtir

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