Estimando a colonização e extinção a partir de dados de censo

Ja vimos em vários posts do que se trata a teoria de metapopulações e como podemos aprender muitas coisas sobre uma espécie modelando ela como uma metapopulação. No entanto eu mesmo já me frustrei muito ao tentar fazer a ligação entre dados e teoria.
Por exemplo como estimar as taxas de colonização e extinção de uma espécie a partir de dados que podemos coletar? Que tipo de dados precisamos e principalmente onde nesses dados esta essa informação que pode nos gerar tantas conclusões.
A resposta é que se a gente coleta várias vezes em um lugar a informação esta la e é mais simples do que pelo menos o que eu esperava.
Olhando esse diagrama fica mais simples, vamos pensar em aranhas em plantinhas como estávamos fazendo nas simulações anteriores.
Então temos:

Na primeira coleta eu olho várias plantas, mas somente duas coisas podem acontecer, ou eu vejo a aranha la ou eu não vejo.
Então se eu olho 10 plantinhas e 5 tem aranhas e 5 não tem aranhas, pronto 50% de ocupação, metade das plantas estão ocupadas por aranhas.
Mas e quanto a colonização? Ue, quando eu voltar para a segunda coleta, se tinha 5 plantas não ocupadas, e agora uma dessas 5 plantas foi colonizadas, são 20% de colonização, 1 colonização de 5 plantas.
E a extinção vai ser o mesmo esquema, se eu tenho 5 plantas que tinham aranhas, e vejo uma planta que tinha aranha e passou a não ter aranha, pronto 20% de chance de extinção, de 5 plantas que tinham aranha, tem uma agora que não tem.
Mas mesmo plantas que não mudam entre uma coleta e outra são informação sobre essa colonização/extinção.
Bem se eu tenho 5 plantas colonizadas, e um aparece na próxima coleta sem aranhas, pimba, 80% de sobrevivência, e o que 1-sobrevivencia é a extinção, ou seja 20% de extinção. A outra coisa é que se eu tinha 5 plantas disponíveis para serem colonizadas e somente uma o foi, eu tenho a chance de não-colonização, ou seja no final das contas eu retirei informação de todas as plantas no intervalo entre 2 coletas. Logos essas taxas valem para esse intervalo.

Então se a gente pode fazer um censo todo ano em várias plantinhas, e ir anotando para cada planta quem esta colonizado e quem não esta colonizado, quem esta sem aranhas. Temos toda a informação que precisamos. So manter a identidade de cada plantinha.

No final nossos dados vão se parecer com uma matriz assim:

> head(yy) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 1 0 0 1 1 1 1 0 0 [2,] 0 0 0 0 0 0 0 0 0 0 [3,] 1 1 1 1 1 1 1 1 1 1 [4,] 0 0 1 1 0 0 0 1 1 1 [5,] 0 1 1 1 1 1 1 1 1 0 [6,] 1 0 0 1 1 1 1 1 1 1

Linhas são as plantas, colunas as coletas, aqui só um pedacinho dos dados que simulamos. As 6 primeiras plantas.
Nesse caso simulamos os parâmetros que sempre viemos usando, 20% de chance de colonização e 10% de chance de extinção.
O que termina com 50% de ocupação.

Mas e agora o que fazer com os dados?
Uma opção no R é usar o pacote unmarked, ele tem várias funções para trabalhar com ocupação e metapopulações.
Na verdade transformar aquilo que falamos em números exige uma modelagem um pouco mais complexa que anovas da vida. Mas usando a distribuição binomial do cara ou coroa, isso é possível, afinal os dados dizem respeito a ver ou não as aranhas nas plantas. Alias os autores desse pacote são insanos, eu chutaria que metade dos livros de estatística bayesiana voltados para ecologia ou biologia são deles, mas aqui ainda não entramos na estatística bayesiana, ajustamos os parâmetros com maior verosimilhança, ou likelihood, que já vimos a algum tempo atrás.

Mas a partir dos dados, vamos estimar então a colonização e extinção de forma simples. Sem covariaveis que podem tornar tudo bem complexo de imaginar.

Primeiro transformamos nossa planilha de observações no formato que o pacote usa…

> simUMF <- unmarkedMultFrame(y = yy,numPrimary=T)

Nesse momento, nos também declararíamos covariaveis, como isolamento, tamanho das plantas e qualquer variável que pode afetar a colonização e extinção, já que esperamos que por exemplo o isolamento faça o acesso se tornar difícil, logo diminua a chance de colonização. Mas não trataremos disso aqui, imaginemos que todas as plantas são perfeitamente iguais.
Outra coisa que achei até maluco a primeira vez, foi como é possível lidar com a detecção imperfeita na coleta, por exemplo num arbusto, ele pode não estar com aranha porque elas não estão la ou ainda porque falhamos em ver elas. Estimar esse erro é possível aqui, o que traz estimativas de colonização e extinção muito melhores, mas ainda assim não simulamos isso.
A nossa simulação foi simplesmente colonização de 20%, extinção de 10% com detectabildiade das aranhas de 100%, vamos ver se realmente conseguimos recuperar esses parâmetros dos dados, usando a função colext() do pacote unmarked do jeito que falamos la em cima.

> summary(simUMF) unmarkedFrame Object 250 sites Maximum number of observations per site: 10 Mean number of observations per site: 10 Number of primary survey periods: 10 Number of secondary survey periods: 1 Sites with at least one detection: 234 Tabulation of y observations: 0 1 988 1512 0

Antes ainda podemos dar um summary pra ver o que criamos, e vemos com o que estamos lidando, são 250 plantas que coletamos durante 10 anos, so visitamos essas plantas 1 vez por coleta, e 234 plantas vimos aranhas pelo menos uma vez nelas.

> #estimando colonização e extinção > m1 <- colext(psiformula = ~1, # Ocupação do primeiro Ano + gammaformula = ~ 1, # Colonização + epsilonformula = ~ 1, # Extinção + pformula = ~ 1, # Deteção + data = simUMF)

Usamos o comando colext(), e nela poderíamos declarar tudo que achamos que pode influencias na colonização, extinção e detectabilidade, e na ocupação inicial, deixamos tudo 1, ou seja sem influencia nenhuma de covariaveis.
Mas assim obtemos o resultado.

> m1 Call: colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, pformula = ~1, data = simUMF) Initial: Estimate SE z P(>|z|) 0.0164 0.127 0.13 0.897 Colonization: Estimate SE z P(>|z|) -1.43 0.0855 -16.8 5.56e-63 Extinction: Estimate SE z P(>|z|) -2.18 0.0936 -23.3 2.36e-120 Detection: Estimate SE z P(>|z|) 7.8 6.21 1.26 0.209 AIC: 2125.146

E ai estão as estimativas de Colonização e Extinção. Usando a função logística podemos obter os resultados em porcentagem, que é o que sempre falamos e que é mais fácil de pensar, então destransformando temos:

> plogis(m1@opt$par) [1] 0.5041035 0.1928461 0.1012372 0.9995916

Ohhh que lindo, ocupação inicial de 50%, praticamente 20% de colonização, levando em conta o erro para estabelecer o intervalo de confiança ali em cima, mas olhando também a extinção temos 10%, o que usamos na simulação, ohhhhhhhh muito bom, e como não simulamos falhas na detectabilidade, ai esta ela 100% a mão.
Realmente é muito legal ver dados saindo da forma bruta para esses parâmetros, que podem nos gerar tantos insights e conclusões de como essa espécie esta agora e qual o seu futuro. E é um exercício também para se aprender que nem toda analise de dados te da um valor p, aqui não tem nenhum.

Bem segue o código usado, essa simulação de dados esta no vignnete do pacote, apenas mudei os parâmetros e tirei as covariaveis.
Mas os vignnetes desse pacote são especialmente bem escritos diga-se de passagem, pelo menos bem mais simples de serem entendidos por mentes não inclinadas a matemática.

Bem e tem também material muito legal aqui nesse site da aula do Richard Chandler. Alias das apresentações dele que eu roubei aquele diagrama do inicio do post.

#######################################################
### Estimando Colonização e Extinção a partir de dados#
#######################################################
set.seed(666)
 
M <- 250 # Número de sítios
J <- 1 # Coletas dentro de um periodo
T <- 10 # Periodos de Coleta
 
psi <- rep(NA, T) # Probabilidade de Ocupação
muZ <- z <- array(dim = c(M, T)) # Ocorrencia observada e verdadeira
y <- array(NA, dim = c(M, J, T)) # Historico de Detecção
 
psi[1] <- 0.5 # Ocupação Inicial
p <- c(1,1,1,1,1,1,1,1,1,1) # Detectabilidade
phi <- runif(n=T-1, min=0.9, max=0.9) # Sobrevivencia (1-extinção)
gamma <- runif(n=T-1, min=0.2, max=0.2) # Colonização
 
z[,1] <- rbinom(M, 1, psi[1])
for(i in 1:M){
for(k in 2:T){
muZ[k] <- z[i, k-1]*phi[k-1] + (1-z[i, k-1])*gamma[k-1]
z[i,k] <- rbinom(1, 1, muZ[k])
}
}
for(i in 1:M){
for(k in 1:T){
prob <- z[i,k] * p[k]
for(j in 1:J){
y[i,j,k] <- rbinom(1, 1, prob)
}
}
}
for (k in 2:T){
psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
}
yy <- matrix(y, M, J*T)
 
head(yy)
#Transformando para a classe da dos usado pelo Unmarked
simUMF <- unmarkedMultFrame(y = yy,numPrimary=T)
summary(simUMF)
 
#estimando colonização e extinção
m1 <- colext(psiformula = ~1, # Ocupação do primeiro Ano
gammaformula = ~ 1, # Colonização
epsilonformula = ~ 1, # Extinção
pformula = ~ 1, # Deteção
data = simUMF)
 
#estimativas
m1
 
#transformando para probabildiade
plogis(m1@opt$par)

Analise de sobrevivência

Ja vimos como algumas vezes mesmo um valor p significativo não significa muita coisa.
Um tipo de modelo muito interessante e o que não tive muitas oportunidades de explorar são os de sobrevivência.

Para fazer um pequeno teste onde ele é melhor que uma análise de variância, vamos fazer um pequeno exemplo.

Uma aranha, vivendo sua vidinha numa flor qualquer corre o risco todo dia de ser predada. Podemos pensar nesse risco como uma chance.
Então podemos simular 100 dias da vida dessa aranha e ver se ela morre algum dia dado que ela tem 10% de chance de morrer por dia.

> rnos which(rnos<= 0.1) [1] 22 24 47 62 70 84 91

Então vemos aqui quais dias ela não passou na prova da sobrevivência, e quais dias ela morreu. Mas como todo mundo só morre uma vez, podemos dizer na verdade que essa aranha morreu no dia 22 de sua vida.

> which(rnos<= 0.1)[1] [1] 22

Certo e qual a relação do dia da morte com a chance de morte da aranha?
Bem se observamos varias aranhas, simulando 30 aranhas com o mesmo procedimento acima vemos que a chance de morte vai ser 1 dividido pela média em dias de vida

> death1 [1] 7 10 1 1 12 9 1 24 3 1 21 5 1 12 3 4 11 19 16 7 15 2 2 6 6 36 6 1 2 4 > 1/mean(death1) [1] 0.1209677

Usamos 10% na simulação, para 30 aranhas obtivemos 12% como estimativa grosseira de sobrevivência. Estimativa razoavelmente boa eu diria.
Agora poderíamos pensar em outra espécie com uma chance maior de morrer, ou ser predada, uma espécie com uma camuflagem ruim, o que aconteceria é que a média em dias seria menor correto? Os indivíduos tenderiam a morrer mais jovens. Vamos simular isso e ver o que acontece, então agora com 20% de chance de morrer.

> death2 [1] 1 1 1 2 8 5 2 3 20 1 4 1 5 9 6 3 1 3 2 3 2 5 11 7 3 8 2 5 2 2 > 1/mean(death2) [1] 0.234375

Certo vemos uma estimativa próxima do que foi usado na simulação, obtivemos 23%, usamos 20%, e vemos aranhas morrendo mais cedo.
Podemos perguntar então se em média aranhas da segunda espécie simulada aqui morrem mais cedo que da primeira e fazendo um modelo linear (teste t) temos.

> summary(model0) Call: lm(formula = death ~ factory) Residuals: Min 1Q Median 3Q Max -7.267 -3.267 -1.767 1.983 27.733 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.267 1.188 6.956 3.48e-09 *** factory2 -4.000 1.681 -2.380 0.0206 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.51 on 58 degrees of freedom Multiple R-squared: 0.08896, Adjusted R-squared: 0.07326 F-statistic: 5.664 on 1 and 58 DF, p-value: 0.02063

Diferença significativa entre as espécies. Mas temos um problema, estamos quebrando fortemente uma premissa dos modelos lineares, a da homogeneidade das variâncias.
Olhem as variâncias das espécies:

> var(death1) [1] 68.75402 > var(death2) [1] 15.9954

Olhando o gráfico, vejam como a espécie um, que tinha 10% de chance de morrer tem uma variância muito maior que a aranha que tinha 20% de morrer. Isso não esta sendo considerado no modelo, o que o deixa muito ruim.

Uma solução é fazer um modelo geral linear com distribuição gamma.
Mas porque essa complicação? O que vamos ver sempre, é que quanto maior a chance de morrer, menor vai ser a média e a variação.
Não é isso que estamos vendo no gráfico ali em cima, dobramos a chance de morte, de 10% para 20% e a média e variância caíram vertiginosamente, podemos fazer mais simulações e ver se isso se confirma, caso não consigamos deduzir isso matematicamente, eu particularmente não consigo, mas deixando isso para depois, vamos ajustar o modelo geral linear e ver o que vamos obter.

> model1 summary(model1) Call: glm(formula = death ~ factory, family = Gamma) Deviance Residuals: Min 1Q Median 3Q Max -1.5705 -0.9383 -0.3327 0.3153 2.0701 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.12097 0.02144 5.642 5.26e-07 *** factory2 0.11341 0.04675 2.426 0.0184 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for Gamma family taken to be 0.9423881) Null deviance: 56.867 on 59 degrees of freedom Residual deviance: 50.422 on 58 degrees of freedom AIC: 337.43 Number of Fisher Scoring iterations: 6

As estimativas são nossas probabilidades, a da primeira espécie (tomado como base, intercepto) é 0.12, próximo do usado na simulação, e para a espécies dois temos uma diferença significativa (p de 0.0184), sendo essa maior que o intercepto 0.1231.
0.12 mais 0.11 da 0.23, próximo dos 20% que usamos na simulação, se pensarmos nos erros dessas medidas (Std. Error) temos que o valor da simulação esta ai dentro.

Se compararmos os resíduos, o que a parte determinística do modelo não consegue descrever, podemos ter uma noção de qual abordagem é mais apropriada.

Vemos como usando o modelo linear temos resíduos gigantes, além de uma tendência a sempre sobrar mais resíduos para baixo, enquanto o que queremos é mais o que vemos no resíduos gerados pelo glm gamma, pequenos e sem nenhuma tendência, esperamos errar igual tanto para mais quanto para menos.

Mas existe ainda uma terceira abordagem para esse problema, que é finalmente a analise de sobrevivência do titulo do post.

Então usamos uma função de sobrevivência para ajustar uma chance de sobrevivência para nossas arahas. Isso e muitas possibilidades estão disponíveis no pacote survival no R.

Então vamos la e primeiro ajustamos as curvas de sobrevivências para as duas espécies de aranhas que estamos trabalhando:

E temos uma escadinha de como a sobrevivência vai diminuindo ao longo do tempo. Lembre-se a cada dia você acumula mais chance de morrer, então o inverso, diminui sua chance de estar vivo a cada dia.

Bem, pelo gráfico já podemos ver como a diferença entre as aranhas, as do exemplo 2 (cor azul) tem uma mortalidade maior, logo sobrevivem menos, logo a curva cai bem mais rápido.
Podemos ajustar parâmetros, no caso essa chance de sobrevivência usando a função survreg() do pacote survival e temos:

> summary(model) Call: survreg(formula = Surv(death, event = rep(1, length(death))) ~ factory, dist = "exponential") Value Std. Error z p (Intercept) 2.112 0.183 11.57 5.91e-31 factory2 -0.661 0.258 -2.56 1.04e-02 Scale fixed at 1 Exponential distribution Loglik(model)= -166.9 Loglik(intercept only)= -170.1 Chisq= 6.45 on 1 degrees of freedom, p= 0.011 Number of Newton-Raphson Iterations: 4 n= 60

Vemos as estimativas, 2.112 para aranhas do exemplo 1 e as aranhas do exemplo 2 diferentes por -0.66. Bem esses valores precisam ser convertidos pela função logística de volta para probabilidade (chance) para podermos comparar com os valores da simulação, então convertendo…

> plogis(coefficients(model)[1]) (Intercept) 0.8920863 > plogis(coefficients(model)[1]+coefficients(model)[2]) (Intercept) 0.8101266

Hummmm 89% e 81%, ahhhh, mas aqui são estimativas da sobrevivência, e simulamos a chance de morrer, logo 1-sobrevivencia será a chance de mortalidade.

> 1-plogis(coefficients(model)[1]) (Intercept) 0.1079137 > 1-plogis(coefficients(model)[1]+coefficients(model)[2]) (Intercept) 0.1898734

Ae, 11% e 19%, valores muito próximos do usados na simulação.

Bem parece que analise de sobrevivência funciona muito bem recuperando os parâmetros que usamos para simular os dados.
Mas existem dois motivos que acho muito muito interessantes para usar analise de sobrevivências.

O primeiro, que está relacionado até com a motivação de criação desta técnica, é a de que para traçar a curva de sobrevivência, não precisamos acompanhar todos os indivíduos até a morte, pode parecer pouca coisa isso, mas se pensarmos por exemplo em pessoas e remédios, muitas vezes acompanhamos pessoas utilizando ou não um determinado tratamento, mas só podemos acompanhar essas pessoas durante o tratamento, então mesmo que não saibamos quando elas morrem, sabemos que ela sobreviveram o período de tempo que as acompanhamos, então temos ai informação a cerca da sobrevivência. Ou ainda temos que realizar um experimento num tempo determinado, então aquilo que não seria informação usando glm, aqui será informação, porque la esperamos todos os indivíduos morrerem.

Outra coisa é que em todos os casos, consideramos uma chance constante de mortalidade (ou sobrevivência se preferir), mas facilmente sabemos que isso pode ser diferente. O exemplo clássico é comparar pólipos com humanos com tartarugas, pólipos sim tem uma chance constante de morrer no mar, já que morrem quando são predados, nos humanos, muito diferente disso, quando nascemos, temos o cuidado de nossos pais, o que diminui muito a chance de morrer quando jovens, mas conforme envelhecemos, a partir de certa idade, doenças e muitos problemas vão aparecendo, diminuição da capacidade física logo podemos esperar um aumento gradativo da chance de mortalidade.
E o ultimo extremo seriam tartarugas, que quando nascem sozinhas na areia tem uma chance altíssima de mortalidade, mas esta vai diminuindo com o tempo, já que ela vai ficando grande e dai consegue se defender cada vez melhor. Isso seria uma complicação enorme nos modelos lineares e glms, mas podemos usando a regressão de sobrevivencia por exemplo, ajustar uma distribuição e modelo para essa alteração de sobrevivência, vejam que aqui usamos dist=”exponential”, ou seja a mortalidade é constante. Mas podemos ajustar vários modelos e comparar o ajuste por exemplo para ter uma visão mais clara de como a sobrevivência muda ao longo do tempo, se um tratamento demora a começar a fazer efeito, etc.

Bem é isso, uma primeira visita a possibilidade do uso de modelos de sobrevivência dentro da biologia. Segue o script, a maior parte dele esta no capitulo de analise de sobrevivência do livro do Crawley “The R Book” para quem quiser ler mais sobre o assunto.

##################################################
set.seed(3)
rnos<-runif(100)
which(rnos<= 0.1)
which(rnos<= 0.1)[1]
 
#exemplo 1
death1<-numeric(30)
for (i in 1:30) {
rnos<-runif(100)
death1[i]<- which(rnos<= 0.1)[1]
}
death1
 
1/mean(death1)
 
#exemplo 2
death2<-numeric(30)
for (i in 1:30) {
rnos<-runif(100)
death2[i]<- which(rnos<= 0.2)[1]
}
death2
1/mean(death2)
 
#juntando os dados
death<-c(death1,death2)
factory<-factor(c(rep(1,30),rep(2,30)))
 
plot(factory,death,ylab="Dia da Morte",frame.plot=F)
 
#Modelo linear
model0<-lm(death~factory)
summary(model0)
 
var(death1)
var(death2)
 
#Modelo gamma
model1<-glm(death~factory,Gamma)
summary(model1)
 
#grafico dos residuos
par(mfrow=c(2,1))
plot(resid(model0),ylim=c(-10,30),main="Modelo Linear (Distribuição normal)"
,ylab="Residuos",xlab="Aranhas")
abline(h=0,lwd=2,lty=2)
plot(resid(model1),main="Modelo Geral Linear (Distribuição gamma)",
ylim=c(-10,30),ylab="Residuos",xlab="Aranhas")
abline(h=0,lwd=2,lty=2)
 
#curva de sobrevivencia
library(survival)
plot(survfit(Surv(death,event=rep(1,length(death)))~factory),lwd=2,
col=c("red","blue"),frame.plot=F,xlab="Dias",ylab="Sobrevivencia")
legend("topright",legend=c("Exemplo 1","Exemplo 2"),lwd=2,col=c("red","blue"))
 
#Regressão de sobrevivencia
model<-survreg(Surv(death,event=rep(1,length(death)))~factory,dist="exponential")
summary(model)
 
#Probabilidades de sobrevivencias
plogis(coefficients(model)[1])
plogis(coefficients(model)[1]+coefficients(model)[2])
 
#Probabilidades de mortalidade
1-plogis(coefficients(model)[1])
1-plogis(coefficients(model)[1]+coefficients(model)[2])

Rosalind – Problema 2 “Transcribing DNA into RNA”

O segundo problema ainda é bastante trivial, converter uma sequência de DNA em RNA.
Na verdade o problema só pede que dentro da sequência de DNA, todos os “T” sejam trocados por U.
Pelo R isso é bem simples de fazer usando Regular Expressions.
Olhando a documentação da função gsub (digite ?gsub) da para ver a família de funções do R mais básicas para tratar de expressões regulares. Apesar de o uso aqui ser simples, Expressões regulares podem ser uma ferramenta muito legal para organizar dados ou extrair informações de um texto gigante.
Dentro da documentação da função também tem um link para uma explicação básica de como montar padrões para usar Expressões Regulares.

Mas nosso intuito aqui é trocar “T” por “U”, então é gsub(achar padrão, substituir por esse, dados para trabalhar)

subs.t.u<-function(sequencia) {
gsub("T","U",sequencia)
}

Com isso o funcionamento fica assim:

> sequencia subs.t.u(sequencia) [1] "GAUGGAACUUGACUACGUAAAUU"

Simples, mas cumpri o papel. E no fórum pessoas inclusive falaram de usar um editor de texto e mandar localizar e substituir uma letra pela outra, solução bem genérica.

No python também, modificar o string é bem simples com a função replace.
É só falar o que é para ser modificado, e o que colocar no lugar.
Só fica diferente a sintaxe, que aqui você tem que colocar string.comando(antigo,novo).
Diferente de usar um argumento dentro da função somente, mas nada muito difícil de digerir com o tempo.

Então minha função ficou assim:

def substseq(sequencia):
    return sequencia.replace("T","U")

Como eu ainda não entrei em abrir arquivos pelo python, colocamos dentro do programa mesmo a entrada de dados, vai ser uma linha, nada muito grande, então tudo ficou assim:

def substseq(sequencia):
    return sequencia.replace("T", "U")
 
sequencia = "GATGGAACTTGACTACGTAAATT"
 
print substseq(sequencia)

E assim temos nossa resposta:

>>> GAUGGAACUUGACUACGUAAAUU

Até o quarto ou quinto desafio da até para se enganar que tudo é simples assim 🙂

Explorando mais um pouco a teoria de metapopulações.

Após a primeira exploração ao modelo de Levins, podemos ainda tentar uma outra abordagem mais “visual” que é resolver o modelo dele para vários valores de colonização e extinção e ver o que acontece. Quando eu digo visual, é porque tudo pode ser deduzido diretamente da equação. Equações de primeiro grau são retas, de segundo grau são tortas e assim vai. Equações logísticas sempre tem essas forma além de outras características. Mas eu particularmente sou péssimo para entender isso, então eu tento resolver a equação para vários valores diferentes dos parâmetros dele e vejo o que acontece, como primeira exploração acho muito legal fazer isso.

Aqui começaremos sempre com a população com 50% de ocupação e usaremos os seguintes valores, 10%, 20%, 40% e 80% para ocupação e extinção, esses valores eu escolhi, eles não tem uma motivação especial. Bem isso da 16 possibilidades, com colonização 10% e extinção 10%, colonização 20% e extinção de 10% e assim até completar todas as 16 possibilidades, isso vai dar 16 gráficos que montaremos numa tabelinha de gráficos.

Recapitulando o que vimos no post anterior, vimos que quando a colonização é maior que a extinção não teremos uma extinção, já quando a colonização é igual a extinção a ocupação ira diminuir ao longo do tempo, e essa diminuição será mais intensa quando a extinção é maior que a colonização.

Então agora resolvendo tudo veremos:

Primeira vamos pensar como os gráficos foram organizados.
Da esquerda para direita, sempre a colonização esta aumentando, e de cima para baixo a extinção esta aumentando.
Logo, na diagonal são sempre valores de colonização e extinção iguais. E nessa diagonal vemos como a tendência ao longo do tempo é a extinção, mas conforme a extinção é maior, ela fica mais intensa, comparem os gráficos mais extremos onde a colonização (ci) e extinção (e) são 10% cada e 80% cada.

Outra coisa que vemos devido a essa organização dos gráficos é que fora dessa diagonal onde “ci” e “e” são sempre iguais, temos na parte superior casos onde a “ci” é maior que a “e”, na parte inferior casos onde “ci” é menor que “e”.

Na diagonal inferior vemos a confirmação da outra afirmação do post anterior que é sempre a tendência é extinção caso a “ci” seja menor que “e”. E vemos como essa tendência é mais forte, o declínio da ocupação ocorre mais rápida conforme a extinção é muito maior que a colonização.

E por último, vamos olhar a diagonal superior, onde sempre o valor de colonização é maior que a extinção.Primeira coisa legal é olhar como a ocupação se manteve em 50%, igual o valor inicial quando a colonização era o dobro da extinção, ou seja com colonização e extinção respectivamente com os valores de 20% e 10%, 40% e 20%, 80% e 40%, sempre a ocupação se manteve em 50%, ou seja para a gente ver como a ocupação depende de um “trade off” da Colonização versus extinção. Nos outros três casos vimos a ocupação aumentar rapidamente.

Bem essas são mais algumas características das metapopulações de Levins visto de um jeito, podemos dizer um pouco mais heurístico, de sair resolvendo a equação para tudo quanto é valor e ver o que acontece.

E no próximo post, chega de teoria e veremos como transformar observações e censos em nossas parâmetros vistos aqui, colonização e extinção.

###############################################################
#Explorando os possiveis resultados da metapopulação de Levins#
###############################################################
#Função de metapopulações de Levins
levins <- function(t, y, parms) {
p <- y[1]
with(as.list(parms), {
dp <- ci * p * (1 – p) – e * p
return(list(dp))
})
}
 
library(deSolve)
 
#Parametros usados
valores<-c(0.10,0.20,0.40,0.80)
#criando uma matrix com todas as combinações
prms<-expand.grid(ci=valores,e=valores)
 
#Ocupação Inicial (50%)
Initial.p <- 0.5
 
#Criando uma matriz para armazenar os resultados
out<-matrix(NA,nrow=100,ncol=16)
 
#Resolvendo para 100 eventos
for(i in 1:nrow(prms)) {
out[,i]<-ode(y = Initial.p, times = 1:100,func = levins,parms = prms[i,])[,2]
}
 
#gráfico
par(mfrow=c(4,4),mar=c(4,4,2,2))
for(i in 1:16) {
plot(out[,i]~c(1:100),type="l",ylim=c(0, 1),xlab="Tempo",ylab="Ocupação",lwd=3)
text(60,0.95,labels=c(paste("Ci=",prms[i,1]," e=",prms[i,2],sep="")),cex=1)
}

G. H. Hardy reclamando da vida…

Todo biólogo, de modo geral já ouviu falar do Equilíbrio de Hardy-Weinberg.

G. H. Hardy  , um dos autores, escreveu em 1940 um ensaio chamado A Mathematician’s Apology. Escrito quando ele tinha 62 anos, lamentando o declínio de sua capacidade matemática. Neste ensaio ele diz…

I have never done anything “useful”. No discovery of mine has made, or is likely to make, directly or indirectly, for good or ill, the least difference to the amenity of the world.

Traduzindo.

Eu nunca fiz nada de “útil”. Nenhuma descoberta minha fez, ou é provável que fará, direta ou indiretamente, para o bem ou para o mal, a menor diferença para o mundo.

Parece que ele errou nessa.

Vi isso na minha aulinha de genética, achei interessante.