Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- roteiro:meta_chuva [2012/05/28 12:14] – criada adalardo
+++ roteiro:meta_chuva [2024/01/09 18:18] (atual) – edição externa 127.0.0.1
@@ Linha 1: / Linha 1: @@
 ====== Chuva de Propágulos ======
-{{:ecovirt:chuvaprop.jpg?150  | Ênio e Beto andando na chuva http://desenhos.kids.sapo.pt/a-chuva.htm}}
+{{:roteiro:chuvaprop.jpg?150 | Ênio e Beto andando na chuva http://desenhos.kids.sapo.pt/a-chuva.htm}}
-Neste exercício, iremos simular a dinâmica de uma metapopulação onde as probabilidades de extinção e de colonização de manchas pela nossa //espécie virtual// são constantes (modelo de chuva de propágulos). Neste modelo há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia. Relembrando, a variação da fração de manchas ocupadas no tempo é descrito pela seguinte fórmula geral:
+Neste exercício, iremos simular a dinâmica de uma metapopulação onde as probabilidades de extinção e de colonização de manchas são constantes (modelo de chuva de propágulos). Neste modelo há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia. Relembrando, a variação da fração de manchas ocupadas no tempo é descrito pela seguinte fórmula geral:
-$$(df)/(dt) = I - E $$
+$$\frac{df}{dt}=I - E $$
-, onde I é a taxa de entrada de migrantes, e E a taxa de saída. A partir dele podemos definir um modelo simples para a dinâmica de ocupação de manchas que formam uma metapopulação:
+ onde $I$ é a taxa de entrada de migrantes, e $E$ a taxa de saída. A partir dele podemos definir um modelo simples para a dinâmica de ocupação de manchas que formam uma metapopulação:
-$$(df)/(dt) = p_i(1 - f)- p_e f $$
+$$\frac{df}{dt}=p_i(1 - f)-p_e f $$
-onde $$p_i$$ é a probabilidade de imigração ou colonização, $$p_e$$ é a probabilidade de extinção e //f// é a fração de manchas ocupadas (número de manchas ocupadas / número total de manchas). Esses, portanto são os parâmetros do nosso modelo.
+onde $p_i$ é a taxa de imigração ou colonização, $ p_e$ é a taxa de extinção e $ f $ é a fração de manchas ocupadas (número de manchas ocupadas / número total de manchas).
+Em primeiro lugar, vamos estabelecer a probabilidade de colonização de manchas vazias (pi), a probabilidade de extinção em manchas ocupadas (pe) e a fração inicial de manchas ocupadas (fi) como 30%, 15% e 40%, respectivamente.
+<code>
+pi=0.3
+pe=0.15
+fi=0.4
+</code>
-Vamos imaginar agora um cenário de uma paisagem virtual com um certo número de manchas de habitat que podem ou não estar ocupadas pela nossa espécie. Por praticidade vamos definir nossa paisagem um número de colunas ($$n_(col)=20$$) e de linhas ($$n_(row)=20$$), totalizando 400 manchas. Vamos definir também uma proporção de manchas ocupadas no tempo inicial (ocupância no $$t_0 = 0.25$$, ou 100 manchas ocupadas).  A nossa espécie virtual tem uma probabilidade de imigração de 10% ($$p_i= 0.1$$) e uma probabilidade de extinção de 5% ($$p_e=0.05$$).
+Agora que já temos os parâmetros do nosso modelo (pi, pe e fi) vamos criar nossa paisagem virtual. Essa paisagem deve ser constituída por manchas de habitat que podem estar ocupadas ou não. Para tanto, criaremos uma matriz de 10 linhas e 10 colunas, sendo que cada célula dessa matriz representará uma mancha. Serão, portanto, 100 manchas.
+<code>
+paisag=array(0,dim=c(10,10,1))
+paisag
+</code>
-<box green 70% | Nossa pegunta>
+Muito bem, mas temos somente a matriz no tempo inicial. Para acompanhar a dinâmica de ocupação de manchas é preciso criarmos uma terceira dimensão: o tempo. Por enquanto vamos criar 10 passos no tempo.
-Bom, até aqui sabemos que nossa paisagem tem uma certa capacidade em acomodar a espécies, definida pelo número de manchas de habitat disponível e que 1/4 dessa manchas estão ocupadas no momento. Definimos também que a espécie tem uma probabilidade de colonização maior que a de extinção nas manchas.
+<code>
- Nosso pergunta é: qual o destino da população, nessa paisagem caso as condições ambientais não se modifiquem e as taxas de extinção e migração permaneçam constantes?
+paisag=array(0,dim=c(10,10,10))
-</box>
+paisag
+</code>
+Ótimo! Agora temos nossa paisagem em 10 tempos diferentes.
-===== Fazendo as iterações =====
+Finalmente, antes de começarmos a brincadeira, precisamos definir quais manchas estarão ocupadas no tempo inicial. Vamos combinar que quando a mancha está ocupada ela recebe o valor 1 e quando está vazia recebe o valor 0. Do jeito que a matriz está, repleta de zeros, nenhuma mancha está ocupada. Então vamos preencher ao acaso algumas manchas, usando um fi de 40%:
+<code>
+fi
+-fi
+rep(1, fi*100)
+ocor<-rep(1,fi*100)
+nocor<-rep(0,(1-fi)*100)
+nocor
+c(ocor,nocor)
+sample(c(ocor,nocor))
+paisag[,,1]<-sample(c(ocor,nocor))
+paisag[,,1]
+image(paisag[,,1], col=c("white", "green"))
+</code>
+Agora sim, estamos prontos para a simulação.
-A  brincadeira aqui é como se tivessemos um tabuleiro com 400 casas, representando nossa paisagem, com 100 delas ocupadas por peças brancas (representando nossa espécie) no início. A cada rodada, para cada casa do tabuleiro, sorteamos um número aleatório de 1 a 100 (um super dado com 100 lados!). Caso a casa esteja ocupada e o número sorteado estiver entre 5 e 1, retiramos a pedra; caso a casa esteja desocupada e o número sorteado for entre 1 e 10, colocamos uma nova peça. Caso essas condições não aconteçam mantemos a casa do tabuleiro da maneira que encontramos.
+O passo seguinte é fazer as coisas acontecerem! Mas para isso temos que ter em mente todas as possibilidades:
+  * **Manchas ocupadas**
-Agora só precisamos fazer isso  muitas vezes e observar a cada tempo quais casas foram colonizadas, quais foram extintas e quais permanecem ocupadas, anotando a proporção de manchas ocupadas (nossa variável de interesse nesse caso!).
+    * podem permanecer ocupadas (1-pe = 0,85)
+    * podem sofrer extinção local (pe = 0,15)
+  * **Manchas vazias**
+    * podem permanecer vazias (1-pi = 0,7)
+    * podem ser colonizadas (pi = 0,3)
+Calma, não se assuste. O monstrinho abaixo vai tratar apenas das manchas que estavam ocupadas no tempo inicial.
+<code>
+sum(paisag[,,1]) # numero de manchas ocupadas no tempo 1
+paisag[,,2][paisag[,,1]==1]<-sample(c(0,1),sum(paisag[,,1]),replace=T,prob=c(pe,1-pe))
+paisag[,,2]
+</code>
+Se você observar atentamente, algumas manchas que estavam ocupadas (1) ficaram vazias (0):
-===== Usando o EcoVirtual =====
+<code>
+par(mfrow=c(1,2))
+image(paisag[,,1], col=c("white", "green"))
+image(paisag[,,2], col=c("white", "green"))
+</code>
+Estes foram eventos de extinção local. Note que o número de manchas ocupadas diminuiu.
+<code>
+sum(paisag[,,1]) # total de manchas ocupadas no tempo inicial
+sum(paisag[,,2]) # total de manchas ocupadas no tempo 2
+</code>
-Como fazer isso com um tabuleiro de 400 casas tomaria muito tempo, além de ser muito tedioso, vamos usar o **//EcoVirtual//** para fazer o serviço para nós.
+Mas essa é apenas metade da história, não? Estamos esquecendo que as manchas vazias podem ser ocupadas. Vamos simular a colonização com um outro monstrinho:
+<code>
+length(paisag[,,1]) ## esse é o tamanho da nossa paisagem ou o número de manchas!
+nmanchas=length(paisag[,,1]) ## vamos guardar esse valor que não muda durante toda a simulação
+paisag[,,2][paisag[,,1]==0]<-sample(c(0,1),nmanchas-sum(paisag[,,1]),replace=T,prob=c(1-pi,pi))
+paisag[,,2]
+</code>
-Para tanto é necessário: (1) abrir o //R//, (2) baixar, instalar e carregar o pacote //Rcmdr// no R, (3) baixar e instalar os pacotes EcoVirtual e RcmdrPlugin-EcoVirtual, (4) abrir o Plugin EcoVirtual no Rcmdr.  Veja como fazer tudo isso na página de  [[ecovirt:instal|Instalação]].
+Compare a paisagem no tempo inicial e no tempo 2 e veja se no final o número de manchas ocupadas aumentou ou  diminuiu.
+<code>
+paisag[,,1]
+paisag[,,2]
-Caso tenha seguido as instruções da página de instalação, partindo da premissa que as instruções funcionam, terá na tela a janela do Rcmdr aberta e uma opção de menu //EcoVirtual//. Clique nesse menu e depois em //Metapopulation// em seguida em //Seed Rain...//. A seguinte janela de opções deve se abrir:
+sum(paisag[,,1]) # total de manchas ocupadas no tempo inicial
+sum(paisag[,,2]) # total de manchas ocupadas no tempo 2
+</code>
+Note que algumas manchas que estavam ocupadas no tempo 1 ficaram vazias no tempo 2 e vice-versa. Note também que algumas manchas continuaram vazias e outras continuaram ocupadas. O próximo passo é calcular a fração de manchas ocupadas inicialmente (f1) e no tempo 2 (f2) e depois a diferença entre elas:
+<code>
+nmanchas
+f1=sum(paisag[,,1])/nmanchas# fração das manchas ocupadas no tempo 1
+f2=sum(paisag[,,2])/nmanchas # fração das manchas ocupadas no tempo 2
-{{:ecovirt:prorain1.png?300|}}
+f1
+f2
+f2-f1
+</code>
+Eis a dinâmica de metapopulações! ;-) Mas até aqui vimos apenas a variação do f do tempo 1 para o tempo 2. Restam ainda 8 passos para chegarmos no tempo 10, só que para não termos que repetir todo o trabalho 8 vezes usaremos um comando chamado **for**, que fará o serviço por nós:
+<code>
+paisag=array(0,dim=c(10,10,10))
+nmanchas=length(paisag[,,1])
+paisag[,,1]=sample(c(rep(0,((1-fi)*nmanchas)),rep(1,fi*nmanchas)))
+resultado=numeric()
+for(t in 2:10)
+               {
+	       paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),sum(paisag[,,(t-1)]) ,replace=T, prob=c(pe,1-pe))
+	       paisag[,,t][paisag[,,(t-1)]==0]<-sample(c(0,1),nmanchas - sum(paisag[,,(t-1)]), replace=T,prob=c(1-pi,pi))
+	       resultado[t-1]=sum(paisag[,,t])/nmanchas
+	      }
+resultado
+</code>
+Vôa lá! Agora temos as frações de manchas ocupadas (f) ao longo do tempo. Só para deixar nossos resultados mais bonitos, vamos colocá-los em uma tabela:
+<code>
+resultado=data.frame(t=1:10,f=c(fi,resultado))
+resultado
+</code>
-<box 70%|Opções do Menu Chuva de Propágulos>
+Agora, para ficar mais bonito ainda, vamos criar um gráfico com esses resultados:
-  * Enter name for data set -> nome para o salvar a simulação  em um objeto no R!
+<code>
-  * Maximum time -> número de interações ($$t_max$$)
+par(mfrow=c(1,1))
-  * Coluns -> número de colunas de habitat na nossa paisagem
+plot(1:10,resultado$f,type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
-  * Rows -> número de linhas na nossa paisagem
+ylim=c(0,1),main="Dinâmica de ocupação de manchas",font.lab=2,lwd=2)
-  * Initial occupance -> fração de manchas ocupadas no inicio da simulação
+</code>
-  * Colonization probability  -> probabilidade de colonização
-  * Extintion Probability -> probabilidade de extinção
-</box>
-Produza simulações com as seguintes situações iniciais:
+Como nossa filosofia de vida é (ou pelo menos deveria ser) melhorar sempre, vamos acrescentar a esse gráfico uma informação muito importante: a fração de manchas ocupadas no equilíbrio (F). O F é calculado da seguinte forma:
-  * 1. $$ t_(max)=100, n_(col)=20, n_col=20,f = 0.2, p_i =0.5, p_e=0.2 $$
-  * 2. $$ t_(max)=100, n_(col)=20, n_col=20,f = 0.2, p_i =0.5, p_e=0.5 $$
-  * 3. $$ t_(max)=100, n_(col)=20, n_col=20,f = 0.05, p_i =0.2, p_e=0.8 $$
-  * 4. $$ t_(max)=100, n_(col)=20, n_col=20,f = 0.05, p_i =0.05, p_e=0.8 $$
-  * 5. $$ t_(max)=100, n_(col)=20, n_col=20,f = 0.10, p_i =0.05, p_e=0.9 $$
-  * 6. $$ t_(max)=100, n_(col)=10, n_col=10,f = 0.05, p_i =0.05, p_e=0.5 $$
-  * 7. $$ t_(max)=100, n_(col)=5, n_col=5,f = 0.10, p_i =0.05, p_e=0.95 $$
-  * 8. $$ t_(max)=100, n_(col)=20, n_col=20,f = 0.01, p_i =0.9, p_e=0.99 $$
-  * 9. $$ t_(max)=100, n_(col)=5, n_col=5,f = 0.01, p_i =0.9, p_e=0.99 $$
-  * 10.$$ t_(max)=100, n_(col)=20, n_col=20,f = 0.9, p_i =0.9, p_e=0.99 $$
-  * 11.$$ t_(max)=100, n_(col)=5, n_col=5,f = 0.01, p_i =0.1, p_e=0.99 $$
+$$F=\frac{p_i}{p_i+p_e}$$
+Portanto...
+<code>
+F=pi/(pi+pe)
+F
+</code>
+Agora que conhecemos o F vamos colocá-lo no gráfico na forma de uma linha horizontal:
+<code>
+plot(1:10,resultado$f,type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
+ylim=c(0,1),main="Dinâmica de ocupação de manchas",font.lab=2,lwd=2)
+abline(h=F,col=2,lwd=2,lty=2)
+</code>
-===== Outras Perguntas =====
+O serviço sujo está feito, agora é hora da diversão! Como somos bonzinhos, resolvemos criar uma função para facilitar a vida de vocês. Com essa função vocês poderão variar os parâmetros do nosso modelo à vontade, sem medo de ser feliz. Por favor, retribuam a gentileza e testem vários valores para cada parâmetro e vejam o que acontece com nossas metapopulações virtuais. Abaixo a função:
+<code>
+metapop=function(tf,c,l,fi,pi,pe){
+	paisag=array(0,dim=c(l,c,tf))
+        nmanchas=c*l
+	paisag[,,1]=matrix(sample(c(1,0),nmanchas,prob=c(fi,1-fi), replace=T),l,c)
+	resultado=numeric()
+	for(t in 2:tf){
+	       paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),sum(paisag[,,(t-1)]), replace=T, prob=c(pe,1-pe))
+	       paisag[,,t][paisag[,,(t-1)]==0]<-sample(c(0,1),c*l-sum(paisag[,,(t-1)]), replace=T, prob=c(1-pi,pi))
+	       resultado[t-1]=sum(paisag[,,t])/(c*l)
+	      }
+	F=pi/(pi+pe)
+	plot(1:tf,c(fi,resultado),type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
+	ylim=c(0,1),main=paste("Chuva de Propágulos","\n c=",c," l=",l," fi=",fi," pi=",pi," pe=",pe),font.lab=2,lwd=2)
+	abline(h=F,col=2,lwd=2,lty=2)
+      return(paisag)
+	}
+</code>
+Maravilha! A função acima tem 6 argumentos (tf, c, l, fi, pi e pe), mas para que a função funcione é preciso atribuir um valor para cada argumento, como no exemplo abaixo:
+<code>
+metapop(tf=100,c=10,l=10,fi=0.5,pi=0.3,pe=0.15)
+</code>
+<box 50%|Parâmetros da função metapop>
+  * tf - número total de passos no tempo
+  * c - número de colunas da matriz //paisag//
+  * l - número de linhas da matriz //paisag//
+  * fi - fração de manchas ocupadas inicial
+  * pi - probabilidade de colonização
+  * pe - probabilidade de extinção
+</box>
+<box 80% red |Exercício>
 Agora é com vocês! Tentem responder às seguintes perguntas:
-  * Quando você rodar mais de uma vez a função com os mesmos parâmetros o resultado é o mesmo? Por quê?((Você deve entender isto para responder várias das questões seguintes.))
+  * O que acontece você rodar mais de uma vez a função com os mesmos parâmetros ? Por quê?((Você deve entender isto para responder várias das questões seguintes.))
   * O que acontece com a variação na fração de manchas ocupadas quando:
       *  o número de manchas é muito grande? E quando é muito pequeno?
@@ Linha 77: / Linha 193: @@
   * E quando a fração inicial (fi) é muito diferente da fração em equilíbrio (F)?
   * Quais condições levam à extinção da população na paisagem neste modelo?
+</box>
-==== DICA ====
+== DICA ==
+Para responder às perguntas acima você precisa comparar resultados de diferentes simulações.  Abrir espaço para mais de um gráfico na janela gráfica do R vai ajudar bastante. Para ter 4 gráficos digite na linha de comando do R:
+<code>
+par(mfrow=c(2,2))
+</code>
+Os dois números indicam o número de linhas e colunas desejadas, na janela gráfica. No caso, pedimos duas linhas e duas colunas, ou seja, espaço para quatro gráficos. Para voltar a um gráfico apenas na janela digite:
+<code>
+par(mfrow=c(1,1))
+</code>
+Para vocês que chegaram vivos até aqui, uma recompensa. Rode os comandos abaixo e descubra:
+<code>
+anima=function(dados){
+        x11()
+	for(i in 1:dim(dados)[3])
+        {
+	image(dados[,,i], main=("Ocupação de manchas"),sub=paste("simulação no.= ", i), col=c("white","red"), bty="n",xaxt='n',yaxt='n')
+	grid(dim(dados)[1],dim(dados)[2])
+	Sys.sleep(.2)
+	}
+}
+simula1<- metapop(20,10,10,1.0, 0.4,0.2)
+anima(simula1)
-  *1. Para responder às perguntas acima você precisa comparar resultados de diferentes simulações.  Você pode salvar os gráficos do R clicando com o botão esquerdo do mouse na janela contendo o gráfico!
+</code>
-  *2. Note que no final são apresentados dois gráficos, um mostrando o cenário da paisagem em quatro tempos e outro com a variação da proporção das manchas ao longo do tempo. Um pode esconder o outro, portanto movimente as janelas para poder visualizar todos
-  *3. As janelas gráficas não são fechadas entre as simulações. De tempos em tempos, salve aquelas que interessam e feche as outras.
-  *4. Caso tenha algum conhecimento rudimentar da linguagem R e queira conhecer os códigos das funções, veja o roteiro a [[ecovirt:chuva_codigo|Código na Chuva]]
+Para fazer rodar a animação você precisa apenas salvar o resultado da função //metapop//  e chamar esse objeto na função anima. Para modificar parmâmetros do modelo é só rodar outra simulação.
+<code>
+simula2=metapop(tf=25,c=100,l=100,fi=.01,pi=0.2,pe=0.5)
+anima(simula2)
+</code>