Nesta aula, iremos simular a dinâmica de uma metapopulação onde as taxas de extinção e de colonização são constantes (modelo de chuva de propágulos). Neste caso há sempre uma fonte constante de imigrantes que podem colonizar qualquer mancha vazia.
Primeiro vamos conhecer um novo comando do R que cria uma matriz de dados tridimensional. O comando é array. Teste a seguinte linha para ver um array de 3 dimensões.
array(NA,dim=c(4,4,3))
Veja que o R criou 3 matrizes de 4 linhas por 4 colunas. Essa matriz assemelha-se a um bolo com 3 camadas, porém quando o objeto é chamado o R mostra as camadas separadamente. No exercício abaixo a primeira dessas 3 matrizes representará o tempo t = 1, e cada célula da matriz
representará uma mancha no espaço, que poderá ou não estar ocupada.
## crie esses objetos para ir testando as linhas da sua função enquanto cria tmax = 20; nc = 10; nl = 10; psi0 = 0.5; cln = 0.3; ext = 0.15 ## Para criar a função. metapop<-function(nl,nc,tmax,psi0,cln,ext){ ## Isso vai criar tmax matrizes com nl linhas e nc colunas. resulta<-array(0,dim=c(nl,nc,tmax)) ## como vamos usar muito o nl*nc (número de manchas), vamos criá-lo
Agora, para cada mancha você deverá jogar uma moeda para decidir se a mancha está ou não ocupada no tempo t =0. A probabilidade de estar ocupada deve ser psi0 (a probabilidade de ocorrência inicial). Para isso use a função sample()
.
Usando sample(c(0,1),1)
obtemos uma amostra de um valor que pode ser 0 ou 1, com igual probabilidade. Mas nós precisamos amostrar nl*nc valores de 0 ou 1, pois temos esse número de manchas.
Então vamos usar:
sample(c(0,1),nl*nc,replace=T )
onde replace=T
significa amostragem com reposição.
Porém, ainda falta informarmos a probabilidade de ocorrência. Se amostrarmos com igual probabilidade de obter 0 (mancha desocupada) ou 1 (mancha ocupada) a probabilidade de ocorrência será sempre 0.5, mas não é isso que pretendemos. Queremos ter a opção de alterar a probabilidade de ocorrência. Para isso usamos o argumento “prob=c(1-psi0,psi0)” que especifica a probabilidade de obter os valores indicados dentro da função sample
. Ou seja, 1-psi0
será a probabilidade de um 0 ser amostrado e psi0
será a probabilidade de um 1 ser amostrado.
Então vamos simular a ocupação no tempo t=0 de acordo com psi0:
resulta[ , ,1]<-sample(c(0,1),nl*nc,replace=T,prob=c(1-psi0,psi0))
A probabilidade de ocorrência psi0 irá variar nos próximos passos, então crie um objeto chamado de resultados “res.psi” e salve o psi0 na posição 1 (psi0 é o valor inicial de fs). No Gotelli f é a fração de manchas ocupadas (número de manchas ocupadas / número de manchas).
Use um for()
para simular o que acontece com a metapopulação em cada tempo subseqüente, até que tmax seja atingido. Repare que já preencheu o array de resultado na primeira posição relativa ao tempo inicial
Agora precisamos criar uma condição para saber se uma mancha que estava ocupada continuará ocupada e se uma mancha que estava desocupada continuará desocupada, ou se será ocupada.
Veja o que apareçe no objeto resulta[,,1]: apenas valores 0 e 1, significando ausência ou ocupação de mancha no tempo t=1. No tempo t+1 ou no caso t=2 as possibilidades são:
Para o primeiro caso a probabilidade de 0 passar a 1 é cln e a probabilidade de 1 passar a 0 é ext. Uma forma de fazer isso é:
num1=sum(resulta[,,t-1]) num0=(nl*nc)-num1 resulta[,,2][resulta[,,1]==1]<-sample(c(0,1),num1,replace=T,prob=c(ext,1-ext)) ...
Esse comando pode ser lido da seguinte forma: nas posições do resulta t2 onde resulta t1 são iguais a 1, coloque 0 ou 1 com probabilidade ext e 1-ext respectivamente.
Faça o mesmo para as posições em t1 iguais a 0, pensando na colonização. Pronto, agora e só mudar o valores de tempo pelo seu indexador para que em cada ciclo um novo tempo seja calculado. Para finalizar o ciclo só falta calcular o res.psi de cada tempo.
res.psi[t]=sum(resulta[,,t])/(nl*nc)
<box 100% |Dica avançada> Sugerimos que use a solução anterior por ser mais intuitiva. Caso queira se aprofundar nos comandos do R, veja a solução abaixo para o mesmo problema: substituir os dados de ocupação de manchas de um tempo para outro, utilizando as probabilidades de extinção e colonização:
É possível fazer o mesmo com uma única linha de comando utilizando a função
ifelse()
. Ela funciona assim:
ifelse(condição, “faça isso”, “se não faça aquilo”)
Então tente algo assim:
resulta[,,t]<-ifelse(Passo 1, Passo 2, Passo3)
Veja abaixo as instruções para os Passos 1, 2 e 3.
ifelse
Instruções para os passos 1, 2 e 3 acima.
Na função acima você deverá, no Passo 1, criar uma condição para saber se uma determinada mancha estava ou não ocupada no tempo t-1.
Caso a mancha estivesse desocupada você deverá, no Passo 2, dizer o que acontecerá com ela no tempo t. Para isso gere um número aleatório de uma distribuição uniforme e veja se esse valor é menor que a taxa de colonização, se o valor for menor a mancha deverá ser ocupada, se for maior ela permanece desocupada. Vamos ver como fazer isso:
Atenção, as duas linhas abaixo são para você testar no R e entender os comandos, estas linhas não devem estar no seu script final.
runif(1)<cln
Veja que o resultado será TRUE ou FALSE, mas nós queremos que seja 0 ou 1
Precisamos transformar o TRUE e o FALSE em valores numéricos, usando a função as.numeric()
as.numeric(runif(1)<cln)
veja que o resultado agora será 0 ou 1.
runif(1'') você não deverá ter apenas 1 valor, e sim nc*nl valores, ou seja, um valor pra cada mancha.
Agora você já sabe como fazer o Passo 2, mas lembre que no
Passo 3: Aprenda com o Passo 2.
Fim das instruções para os passos 1, 2 e 3.
</box>
Agora faça um gráfico, dentro da sua função, que mostre a trajetória da proporção de manchas ocupadas ao longo do tempo. Acrescente uma reta mostrando o valor da fração de manchas ocupadas no equilíbrio f.equilibrio.
Equação 4.4 do Gotelli: f=Pi/(Pi+Pe)
Primeiro calcule o f.equilibrio
f.equilibrio<- ??
plot(1:tmax),fs ,type="l",xlab="tempo",ylab="fs",ylim=c(0,1),main="Metapopulação ") abline(h=f.equilibrio, col=2)
Agora vamos criar uma animação para olhar nossos resultados ao longo do tempo. Este mostra uma animação da dinâmica da metapopulação ao longo do tempo.
Primeiro salve o resultado de uma simulação da função metapop e crie um gráfico do tempo 1.
image(dados[,,i], main=("Ocupação de manchas"),col=c("white","red")) grid(nl,nc)
Agora é só criar o ciclo com for para que os gráficos de cada tempo sejam sobrepostos no mesmo dispositivo gráfico. Para isso vamos criar uma função anima.meta:
Sys.sleep(.5)
Use sua função e altere os parâmetros iniciais para ver se a metapopulação persiste, dependendo das condições iniciais. Se tiver tempo, procure alterar a sua função para implementar os outros modelos descritos na aula téorica.
<box 90% blue| Uma animação mais informativa> Uma forma de representar a dinâmica espacial da metapopulação é incluir no mapa não só os espaços ocupados e vagos mas entre esses aqueles que foram recentemente colonizados ou extintos. Ou seja, teremos quatro tipos de informações em cada momento: ocupado antigo, ocupado recente, vago, extinção recente.
meta.anima2=function(dados) { nsim=dim(dados)[3] image(dados[,,1], main=("Ocupacao de manchas"),col=c("white","green"),sub="tempo=1") for(i in 2:nsim) { conta12=dados[,,(i-1)]+ (2*dados[,,i]) image(conta12, main=("Ocupacao de manchas"),col=c("white", "red","green", "darkgreen"),sub=paste("tempo= ",i)) Sys.sleep(.5) } }
</box>
Podemos eliminar do modelo anterior o pressuposto de uma chuva de propágulos constante e fazer com que a colonização seja uma função do número de lugares ocupados. Em uma formulação simples desse modelo, a fonte de propágulos é unicamente interna (sistema fechado) e a probabilidade de colonização varia de forma linear à proporção de lugares ocupados. Nossa função começa então diferente com relação aos argumentos necessários:
meta.inter<-function(nl,nc,tmax,psi0,ext, ef) { ...
... for(t in 2: tmax) { cln=ef*res.psi[t-1] ... ... }
... f.equilibrio=1-(ext/ef) ...