Data Augmentation para dados de Captura e Recaptura com efeito do tempo como uma variável aleatória

A algum tempo atrás, a gente falou de data augmentation aqui e depois usando data augmentation num modelo com efeito do tempo aqui.

Muito legal, mas particularmente no modelo com efeito do tempo, a gente cogitou (no livro do Marc Kery é cogitado na verdade) a possibilidade de tratar esse efeito como uma variável aleatória, então vamos tentar fazer isso e uma pequena modificação depois para ver o que acontece.

Bem os dados serão os mesmo que para o post de data augumentation com efeito do tempo. Mas vamos mudar o modelo. Na verdade vamos mudar um pouco a geração dos dados, vamos pegar o efeito de tempo de uma distribuição normal, para condizer com nossa analise.

Mas no mais é tudo igual, agora o modelo na linguagem Bugs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()

Basicamente ele continua o mesmo, só mudamos a parte do prior la em cima.

1
2
3
4
5
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo

Veja que agora temos ainda a geração dos efeitos do tempo, de cada tempo de amostragem, mas agora pegamos esses valores de uma distribuição normal, e não de uma distribuição uniforme como da outra vez.

Mas tentamos determinar agora, qual a média e o desvio padrão dessa distribuição.
Para a média partimos de uma distribuição normal com média zero e desvio 0.001, ou seja, um pico no zero, e para o desvio padrão saímos de uma distribuição uniforme de 0 a 5, lembre-se que o desvio padrão tem que ser um número maior que zero, e não pode ser negativo, senão o bugs vai reclamar.

Uma outra coisa da linguagem bugs é que a distribuição normal usa precisão, e não o desvio diretamente na função. Isso não tem nenhum impacto na estatística, tem haver com precisão numérica apenas eu acho. Até hoje não entendi esse conceito perfeitamente, mas isso não muda nada.

1
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)

No fim, tau vai ser 1 dividido pela variância, que é o desvio padrão ao quadrado, veja que o desvio só vai aumentar, e um sobre um valor entre zero e alguma coisa, vai ser cada vez mais um valor pequenininho, e bem é isso.

Fazemos umas cadeias maiores, adicionamos na função de geradores de valores iniciais nossos novos priors

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5))

E é isso, agora é so rodar o modelo.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.825 0.995 99.000 99.000 100.000 100.000 102.000 1.001 7100 p[1] 0.656 0.047 0.562 0.625 0.658 0.689 0.743 1.001 14000 p[2] 0.140 0.034 0.080 0.115 0.138 0.162 0.212 1.002 3000 p[3] 0.070 0.025 0.029 0.051 0.067 0.085 0.124 1.002 2400 p[4] 0.559 0.049 0.464 0.525 0.559 0.593 0.653 1.001 14000 p[5] 0.931 0.027 0.869 0.916 0.935 0.951 0.974 1.002 2000 med.tempo 0.470 0.306 -0.123 0.325 0.469 0.612 1.073 1.001 14000 sigma.tempo 0.579 0.381 0.236 0.359 0.471 0.663 1.608 1.001 14000 omega 0.402 0.031 0.341 0.381 0.402 0.423 0.464 1.001 5100 deviance 436.569 10.169 424.600 428.000 435.400 441.800 461.100 1.001 4200 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 51.7 and DIC = 488.3 DIC is an estimate of expected predictive error (lower deviance is better). >

E apesar da estimativa do tamanho populacional estar ok, as estimativas do tempo ficam bem ruins, confira os valores que geraram os dados com os que foram estimados.

> variacao_Tempo [1] 1.517343 -1.058690 -1.568697 1.357698 3.792000

Veja que a estimativa do desvio padrão ficou bem abaixo do real usado, que era 2. Ela ficou algo entre zero e 1.6, enquanto o desvio real era 2. A média ficou correta quanto ao usado para gerar os dados, englobando o zero, mas o desvio deixou a desejar.

Se a gente olhar a figura de das estimativas do tamanho da população.

01

Veja que a gente perdeu bastante em variabilidade, já que o tempo variou menos, porque ficou “preso” dentro do nosso desvio padrão subestimado.

Mas ainda existe uma possibilidade. Lembre-se que aqui temos poucas medidas de de tempo, apenas 5. Lembra de alguma distribuição que tentamos tratar do número de amostras para controlar o desvio da média? Se pensou na distribuição t, pensou certo.

Podemos modelar a variável aleatória que gera o tempo como uma distribuição t, isso vai adicionar um parâmetro a mais no modelo, mas no entanto, teremos como ter a possibilidade de ter os valores mais longe da média como valores mais comuns para o valor do tempo. Ja que se você ler o post sobre a distribuição t, vera que isso que o novo parâmetro no t faz, o grau de liberdade, quanto mais graus de liberdade, mais parecido com uma distribuição normal, mas quanto menos, deixando a distribuição mais “uniforme”, com uma maior chance de valor mais longe da média, e como temos apenas 5 valores, isso pode ser uma propriedade interessante. É meio maluco pensar que vamos estimar esse parâmetro, já que na distribuição t ele é simplesmente o valor de amostras, mas ok, podemos estimar um valor para ele, para tentar ajustar um modelo mais realista.

De todo aquele modelo, vamos mudar apenar a parte dos prior, mais especificamente onde está a variável aleatória que é o tempo.

1
2
3
4
5
6
7
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }

Veja que usamos o dt ao invés de dnorm dentro do loop, e temos um parâmetro chamado de k.tempo agora, que vamos tirar de uma distribuição uniforme de 2 a 30, quando ele for 30, temos uma distribuição normal, quanto menor, deixamos tudo mais uniforme. Começamos no 2 porque o valor de k tem que ser no mínimo 2 para a linguagem bugs, veja a definição aqui.
E é isso, agora é só colocar ele no inicializador.

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))

E pronto, podemos rodar. Não esquecendo que temos que colocar ele entre os parâmetros que temos que acompanhar.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.866 1.010 99.000 99.000 100.000 100.000 102.000 1.001 11000 med.tempo 0.470 0.294 -0.100 0.322 0.469 0.613 1.065 1.001 14000 sigma.tempo 0.558 0.361 0.217 0.345 0.462 0.648 1.519 1.001 5900 k.tempo 16.565 7.884 2.979 9.950 16.660 23.360 29.310 1.001 5100 omega 0.402 0.031 0.343 0.380 0.402 0.423 0.465 1.001 14000 deviance 437.026 10.301 424.700 428.200 435.800 442.900 461.900 1.001 14000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 53.1 and DIC = 490.1 DIC is an estimate of expected predictive error (lower deviance is better).

E olhando para a figura temos a mesma coisa

02

E no final das contas, o valor de k ficou entre 2 e 30, ou seja, não conseguimos estimar nada para ele, ja que de 2 a 30 era as possibilidades que abrimos. As estimativas de média para a distribuição de parâmetros continua ok, mas o desvio ainda está subestimado. Então apesar dessa ideia de usar a distribuição t seja legal, ela não contribuiu em nada praticamente, aumentamos um parâmetro e não melhoramos muito o modelo. Veja que o deviance continua a mesma coisa. Bom mas valeu a tentativa.

Eu imagino que dado o pequeno número de tempos, essa ideia de usar ele como uma variável aleatória não seja uma boa ideia, ja que estamos trocando 5 parâmetros por 3, ou 2 no primeiro exemplo. Além disso, precisamos lembrar que essa é uma afirmativa forte, pensar que que o tempo é uma variável aleatória, talvez num span pequeno de tempo, pode estar tudo ok, mas num span de tempo maior, dependendo da espécie que estamos trabalhando, podemos pensar o tanto que falamos em aquecimento global, em degradação do ambiente, em tantas mudanças que falar que nada acontece no tempo é algo bem “forte”. Mas apesar de tudo, talvez para mais medidas e um span pequeno de tempo, essa seja uma possibilidade bem legal.

Bem é isso ai, a gente tem que tentar fazer as coisas para aprender, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

John K. Kruschke 2014 Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan.2nd Edition. Academic Press / Elsevier.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = rnorm(tempo,0,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                  sigma.tempo = runif(1,0, 5))
 
# Parametros para saida
params <- c("N", "p","med.tempo","sigma.tempo","omega")
 
# Configurações do MCMC
ni <- 5000
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
jpeg("01.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()
#####################################
## Efeito aleatorio com distribuiçao t
#####################################
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
 
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }
 
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                         sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))
inits()
# Parametros para saida
params <- c("N","med.tempo","sigma.tempo","k.tempo","omega")
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura2
jpeg("02.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()

Fraude científica?

Hoje eu vi um artigo muito legal na revista Science, revisando a controvérsia Mendel-Fisher.
Eu particularmente nem sabia disso, mas foi uma discussão interessante apresentada no artigo, que está na referência ali em baixo.

A controvérsia foi que um cara chamado W.F.R. Weldon achou que os dados do Mendel eram bom demais para serem de verdade.

So lembrando que de acordo “Law of Segregation”, os gametas são formados por um alelo apenas de cada progenitor, e precisamos de dois alelos dominantes para exibirmos o fenótipo dominante.
Ou seja, para uma característica “simples”, se cruzarmos dois indivíduos heterozigotos, esperamos que 3:1 seja a razão entre dominantes e recessivos, 75% de dominantes e 25% de recessivos.

O trabalho do Mendel foi a base para a genética, e mesmo antes da era do DNA e máquinas de sequenciar DNA, ele já tinha a teve a sacada que seriamos diploides e tal. Acho que não há muita controvérsia sobre as leis de Mendel hoje.

Mas a critica do Weldon foi a seguinte, ele olhou os dados do Mendel:

figura1

De um dos experimentos dele, para uma característica das ervilhas, mas olhando todos os dados, e pensou, “Esses dados são bons demais para ser de verdade”. Claro que estamos hoje em outra época, e temos um conhecimento muito maior dos mecanismos que fazem os dados saírem desse jeito, no experimento do Mendel. Quando eu estava na graduação, eu fiz cruzamentos com moscas Drosophilas melanogaster, e via o resultado dele certinho, é batata, mas na época, esse Weldon achou que tudo era bom demais.

Realmente, hoje, vira e mexe a gente vê artigos sendo retratados por má conduta científica. Mas é muito difícil provar que dados são gerados. Mesmo a gente tentando gerar dados para os resultados do Mendel, esperamos mais ou menos o que ele encontrou.

figura2

Mas depois o Ronald Aylmer Fisher deu um chega na discussão. Ainda que podemos reclamar dele, já que ele usou extensivamente os resultados do Mendel, combinando com a Teoria da Evolução para formar a Teoria NeoDarwinista, que junto com mais uma galera, deu uma descrição matemática para a teoria do Darwin usando a ideia do Mendel como um dos pontos de apoio, mas ele observou que Mendel conseguiu abstrair sobre sua teoria, pensar o que seria esperado e so então fez seus experimentos, muito mais que talvez ter “inventado dados”. Acho que essa é uma boa observação, porque de tantas características ou espécies que ele podia avaliar, ele escolheu o modelo certo junto com as características certas, ja que ele não pegou características co-dominantes ou genes deletérios, que alteram as razões com que esperamos os descendentes. Talvez ele tenha sido seletivo no que observar, mas é bastante razoável, essa precisão, dado o número de sementes que ele contou. Não posso falar muito mais também porque não li o trabalho do Mendel, não sei nem onde achar o original, mas acho que como o Fisher disse, foi no máximo uma pseudo-replicação do Experimento, segundo o Hulbert.

Hoje não existe muito bem como dar uma resposta definitiva se uma sequência é aleatória ou não, podemos procurar padrões nelas, mas é outros quinhentos afirmar a ausência deles. No gerador de números aleatórios que usamos toda hora no R mesmo, o Mersenne Twister, a gente não consegue ali no R, nos nossos notebooks demostrar a presença de padrões, porque eles são muito longos, mas eles existem.

Mas ainda sim, a gente tem muitas historias de quem faz fraudes cientificas e é pego, mas o que não deve ser o caso do Mendel, mas fica a historia de o que é legal na ciência, que o que importa é a ideia e não quem a propos, e todo mundo está apto a ser criticado.

Referência

Beyond the “Mendel-Fisher controversy”Gregory RadickScience 9 October 2015: 159-160

Hulbert 1984 – PSEUDOREPLICATION AND THE DESIGN OF ECOLOGICALFIELD EXPERIMENTS Ecological Monographs 54(2),187-192

Evolução – Fisherian Optimality Models – Parte 5

Descendo a quantidade de um post por mês quase, mas ainda tenho esperança de reverter esse quadro.

Agora, voltando a evolução básica de Fisher, lembrando que é básica para ele, pro Fisher, não que eu ache isso básico, de forma nenhuma, gostaria, mas não dou conta de achar isso básico ainda. Mas ok, até agora nos assumimos que a medida apropriada para o fitness é o sucesso reprodutivo ao longo da vida, R0. Apesar que essa medida pode ser apropriada para uma população estável, uma medida de fitness mais geral é o parâmetro maltusiano, r, que é igual a taxa de crescimento da população na distribuição de idades estável.

 \int_0^\infty e^{-rt} l(t)m(t)dt =1

 \sum\limits_{t=1}^\infty  e^{-rt}l_t m_t = 1

onde l(t), l_t são as probabilidades de sobrevivência na idade t e m(t), m_t são os nascimentos específicos da idade de fêmeas (= fecundidade/2, assumindo uma razão sexual de 1:1). A diferença da notação usada nas duas equações são geralmente de pouco ou nenhuma consequência ( a equação da diferença pode ser igualmente escrita como a equação da integral) e usada simplesmente para ilustrar como a diferença na notação não deve implicar em diferença na interpretação. Note que a equação da diferença (do somatório) começa no fim do primeiro período de tempo, desde que nos assumimos que a fecundidade é zero no nascimento (claro que poderíamos iniciar do zero se nos simplesmente colocarmos m_0=0). Apesar que essas equações não englobam diretamente a contribuição dos machos, nos poderíamos escrever de forma similar equações para os machos relacionando seu sucesso em encontrar fêmeas a taxa de crescimento. Esse pressuposto referente ao r é que qualquer mutação que aumenta o r devera aumentar em frequência na população. A algum tempo atras, escrevemos sobre isso aqui, esse post foi muito legal, para ver como o r é determinante ao longo do tempo, para ver como um fenótipo ou característica muda de proporção nas populações, mesmo quando por exemplo ta todo mundo decaindo. Outra coisa para ter claro na cabeça olhando esse modelo, que sempre estamos levando em conta o tamanho, é lembrar da seleção sexual. Da para dar uma olhada no artigo do Ronald Aylmer Fisher aqui para ter uma ideia de como é um cara que escreveu seu nome na historia rediz um artigo, alias esses artigos meio que escritos a maquina de datilografar, sem figuras, porque acho que o povo desenhava a mão, mas uma leitura uma vez na vida de uns paper desse não mata ninguém.

Mas ok, vamos voltar ao que interessa, nos vamos considerar os mesmos pressupostos, com exceção que agora a medida de fitness é o r.

Pressupostos gerais.

1. O organismo é iteroparo.
2. Fecundidade, F, aumenta com o tamanho do corpo, x, o qual não muda depois da maturidade
3. Sobrevivência, S, diminui com o tamanho do corpo, x.
4. Fitness, W, é uma função da fecundidade e da sobrevivência.

Pressupostos matemáticos

1. Maturidade ocorre na idade 1 depois da qual, não há mais crescimento.
2. Fecundidade aumenta linearmente com o tamanho na maturidade, resultando na fecundidade sendo uma função uniforme da idade

 F_t = a_F + b_F x

3. A taxa instantânea de mortalidade aumenta linearmente com o tamanho do corpo obtido na idade 1 e é constante por unidade de tempo. Sobre esse pressuposto, a sobrevivência na idade t é dada por

 S_t = e^{-(a_S+b_S x) t}

Note que para fazer a sobrevivência declinar em função do tamanho do corpo, dado que temos uma função exponencial, nos trocamos o a_S-b_Sx por a_S+b_Sx, outra coisa que pode ser confuso aqui, é falar que é linear e a equação ser uma exponencial, veja que mesmo que tenhamos uma chance p de morrer no tempo 1, no tempo 2, vai ser o p da primeira e o p da segunda, no tempo 3, vai ser o p do tempo 1, o p do tempo 2 e o p do tempo 3, por isso aquele e elevado a equação ali.

4. O Fitness, W, é o parâmetro maltusiano r, pegando o r para ser a medida de fitness, a função de fitness é dada pela solução da equação característica:

 \int_1^\infty e^{-rt} (a_F+b_Fx)e^{-(a_S+b_Sx)t} dt =1

onde o valor inicial da integral é 1, ja que essa é a primeira idade reprodutiva.
Os dois expoentes podem ser absorvidos em um único termo, ja que os dois são “e” elevado a algo, dando

 \int_1^\infty (a_F+b_Fx)e^{-(a_S+b_Sx+r)t}dt=1

Agora, a equação acima tem a mesma forma geral da equação da da parte 4, e a integração dela vai ser mais ou menos a mesma coisa, a gente so esta adicionando a parte do crescimento populacional exponencial.

 1= (a_F+b_Fx) [\frac{-e^{-(a_S+b_Sx+r)t}}{a_S+b_Sx+r}]_1^\infty = 0 + \frac{a_F+b_Fx}{a_S+b_Sx+r} e^{-(a_S+b_Sx+r)}

Plotando a função do fitness

Para plotar r, nossa medida de fitness, como uma função de x, nos precisamos resolver a integral anterior ou a função acima para cada valor de x, ao gosto do cliente. Podemos primeiro examinar os métodos de estimar o r sem nos mesmos resolvermos a integral, e depois usando a integral resolvida.

Usando o R para fazer a integração.

Nos primeiro consideramos a integração numérica para estimar r. A estratégia é de iterar pelo espaço de x e para cada valor e calcular r.

1. integral produz a função a ser integrada, o  (a_F+b_Fx)e^{-(a_s+b_Sx+r)t}

1
2
3
4
5
6
7
8
9
10
11
 
# Função para retornar a função a ser integrada
integral <- function(idade,x,r) {
    #Parametros
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    #Função
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}

2. cal_integral chama a função do R integrate para obter a integral, passando ela para o integral que acabamos de definir

1
2
3
4
cal_integral <- function(r,x) {
    #retornamos somente o valor
    return(1-integrate(integral,1,Inf,x,r)$value)
}

3. A função encontra_raiz usa a função do R uniroot para achar o valor de r que satisfaz a equação característica para um dado x.

1
2
3
encontra_raiz <- function(x) {
    uniroot(cal_integral, interval=c(1e-7,10),x)$root
}

4. Agora é só criar uma matriz de com uma única coluna valores de x e usa a função apply do R para usar encontra_raiz em cada valor de x (linha da matrix x) para pegar o valor requisito de r, guardando no vetor r.

1
2
3
4
5
6
#Figura
#Valor de x a explorar
x <- matrix(seq(0.5,3, length=100)) 
length=100
# Calcular r dado x
r <- apply(x,1,encontra_raiz)

5.O vetor x e r são então usados para plotar a relação entre r(fitness) e x(tamanho do corpo). Rm sequencia, para um dado x, RCALC chama uniroot que chama INTEGRAL que então chama INTEGRAND

1
2
#Figura
plot(x,r,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, r",las=1,lwd=4)

Usando nossa solução para a integral

Nos podemos usar a nossa solução de integral também, é mais simples no quesito de código, entretanto, como muitas funções não podem ser integradas analiticamente, ela é menos geral, além do que estamos assumindo que você consegue integrar bem as funções, eu sou péssimo com integrais, saindo daquilo que são os exemplos “triviais”.

#Fazendo a figura com nossa equação integrada

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#Nossa função de integral
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S)
}
 
# Função para achar r dado x
acha_raiz <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Explorando x
x <- matrix(seq(0.5,3, length=100))
r <- apply(x,1,acha_raiz)
#Figura
plot(x, r, type="l", xlab="Tamanho, x", ylab="Fitness, r",las=1,lwd=4)

E então temos nossa figura:

01

E temos ai, como o fitness varia em relação ao tamanho do corpo, podemos ver que temos aquele bom e velho ponto máximo, o melhor tamanho de corpo para se ter, agora é so iniciar a busca para saber qual valor é esse.

Encontrando o máximo usando cálculo.

Nossa medida de fitness não está mais em um lado da equação, e ela não pode ser simplificada para ser assim. Mas ainda podemos usar a diferenciação implícita ja que o r depende do x. Por “conveniência” (porque senão não anda), nos usamos logaritmos para converter a função em uma que é aditiva.

 0 = ln(a_F+b_Fx)-(a_S+b_Sx +r) - ln(a_S+b_Sx+r)

Veja que temos 3 termos separados por um sinal de menos

 0 = T1 - T2 - T3

Agora, trabalhando com cada termo separadamente temos:

T1 :  \frac{dr}{dx} = \frac{b_F}{a_F+b_Fx}
Esse é fácil, não tem r na parada

T2:  \frac{dr}{dx} = b_S + \frac{dr}{dx}
Aqui veja que temos o r e x

T3:  \frac{dr}{dx} = ( \frac{1}{a_S+b_Sx+r} )(\frac{dr}{dx}) + \frac{b_S}{a_S+b_S+r}

Juntando todos os termos, ficamos com

 \frac{dr}{dx} (1+ \frac{1}{a_S+B_Sx +r}) = \frac{b_F}{a_F+b_Fx} -b_S - \frac{b_S}{a_S+b_Sx+r}

Agora quando  \frac{dr}{dx}=0 , toda a parte da direita se iguala a zero, dado que o termo no parenteses do lado esquerdo não se iguala a zero também, o que geralmente não sera o caso. Assim podemos rearranjar a parte da direita para fazer o r uma função de x:

 r^*= \frac{b_S(a_F+b_Fx^*)}{b_F-b_Sa_F-b_Sb_Fx^*} -b_Sx^* - a_S

onde r^* é o máximo valor de r, obtido no x^*. Para achar x^* nos substituímos essa equação la na primeira com os logaritmos.

 0=ln(a_F+b_Fx^*)-(a_S+b_Sx^*+r^*)-ln(a_S+b_Sx^*+r^*)

 0=ln(a_F+b_Fx^*)-[a_S+b_Sx^*+\frac{b_S(a_F+b_Fx^*)}{b_F-b^Sa_F-b_Sb_Fx^*} -b_Sx^*-a_S]  -ln[a_S+b_Sx^*+\frac{b_S(a_F+b_Fx^*)}{b_F-b_Sa_F-b_Sb_Fx^*}-b_Sx^*-a_S]

que pode ser finalmente resolvido numericamente, mas deus me livre ter que fazer isso na ponta do lápis, acho que desde a parte 2, já não rola mais heheh. Mas o R ta ai pra isso.

Usando uniroot

A primeira abordagem é o uniroot. Note que os limites são setados ligeiramente perto dos valores requeridos. Se os limites estiverem muito longes um do outro () a função pode falhar, retornando uma mensagem de erro.

> uniroot(f=funcao,interval=c(1.2,3))$root Error in uniroot(f = funcao, interval = c(1.2, 3)) : f.upper = f(upper) é NA Além disso: Warning message: In log(As + Bs * x + r) : NaNs produzidos

Esse erro mostra a importância do plot original, no “olhômetro”, da para ver quais seriam valores razoáveis e mais próximos do “ótimo”.

Então definimos a função que chegamos usando o calculo

1
2
3
4
5
6
7
8
funcao <- function(x) {
    Af <- 0;
    Bf<-4*4;
    As<-1;
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As # r from eqn (2.40)
    return(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r))
}

e usamos o uniroot

> uniroot(f=funcao,interval=c(1.2,1.8))$root [1] 1.389974

Usando o nlm

Uma forma alternativa é usar o nlm, usando o valor absoluto na função, que no caso o mínimo tem de ser zero.

1
2
3
4
5
6
7
8
funcao <- function(x) {
    Af<-0
    Bf<-4*4
    As<-1
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(abs(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r)))
}

Veja que o que mudou para a função anterior é somente o abs de valor absoluto ali na saída.
agora usamos a função nlm

> nlm(funcao, p=1.2)$estimate [1] 1.389943 Warning messages: 1: In log(As + Bs * x + r) : NaNs produzidos 2: In nlm(funcao, p = 1.2) : NA/Inf substituido pelo máximo valor positivo

Como x não é restrito a uma região, esse método não é satisfatório e apesar de termos a resposta correta, uma mensagem de aviso é gerada.

Usando optimize

Usar optimize é bem simples no entanto, basta usarmos a mesma função do nlm. Sendo que aqui definimos o mínimo e máximo, que é fácil estimar a partir da figura do fitness.

> optimize(f=funcao, interval = c(1.2,1.8),maximum= FALSE)$minimum [1] 1.389952

O intervalo de nlm de busca é muito grande, veja também que apesar de muito próximos, encontramos valores diferentes, diferentes a partir da quinta casa depois do zero, mas valores diferentes.

Usando métodos numéricos

Nos vamos ver duas abordagens aqui, a primeira, usando a equação integral resolvida e a segunda usando métodos númericos para a integração do modelo original. Esse é o caso mais geral usando métodos númericos mas também o que consome mais tempo. Mesmo que a integral possa ser resolvida, essa é uma boa prática para checar se a solução da integral está realmente correta.

Usando a função integrada.

Aqui nos usamos a função optimize, dando a ele a função que calcula a raiz (raiz_funcao) usando uniroot para achar a raiz da função especificada em “funcao”.

1
2
3
4
5
6
7
8
9
10
11
12
13
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S) #lembrando que aqui ja é a integral
}
 
# Achando r em função de x
raiz_funcao <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}

E assim nos temos

> optimize(f = raiz_funcao, interval = c(.5,3),maximum = TRUE)$maximum [1] 1.389946

Usando integração númerica para a função.

Lembrando que a gente ja teve uma noção do que é integração númerica aqui e aqui, nos podemos pular a parte da integral feita a mão (caso lidemos com um modelo difícil) e integrar e depois achar o valor de r em função de x e depois maximizar o fitness. Como dito anteriormente, mesmo com uma função integrável, essa pode ser uma garantia que tudo está correto.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
funcao_integrar <- function(idade,x,r) {
    Af <- 0
    Bf <- 4*4
    As <- 1
    Bs <- 0.5
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função que calcular a integral
integral <- function(r,x){
    1-integrate(funcao_integrar,1,Inf,x,r)$value
}
# Função para achar r dado x
raiz_integral <- function(x) {
    uniroot( integral, interval=c(1e-07,10),x)$root
}

E agora usamos o optimize para achar o máximo da função de fitness.

> optimize(f = raiz_integral, interval = c(1.2,1.8),maximum = TRUE)$maximum [1] 1.389934

Veja que aqui, estamos começando a esbarrar com problemas de precisão númerica, que não tem muito efeito por enquanto, mas pode vir a ser um problema.

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail. Esse livro do Derek Roff é realmente muito legal, estamos ainda no segundo capítulo para quem quiser acompanhar, e queria eu ter a oportunidade de desde a graduação, estar exposto a um conteúdo desses, mas nunca é tarde para estudar e aprender coisas novas. Nas próxima partes, os modelos vão começar a ficar realmente legais com a entrada de variação e mais variáveis a maximizar. Até a próxima espero eu que a próxima chegue mais rápido.

Referência:
Derek A. Roff 2010 – Modeling Evolution – an introduction to numerical methods. Oxford University Press 352pp.

R. A. Fisher 1915 The evolution of sexual preference Eugenics Reviews 7(3): 184–192.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#Plotando a função de fitness
 
 
# Função para retornar a função a ser integrada
integral <- function(idade,x,r) {
    #Parametros
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    #Função
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função para integrar a equação e retornar os valores
cal_integral <- function(r,x) {
    #retornamos somente o valor
    return(1-integrate(integral,1,Inf,x,r)$value)
}
 
 
encontra_raiz <- function(x) {
    uniroot(cal_integral, interval=c(1e-7,10),x)$root
}
 
#Figura
#Valor de x a explorar
x <- matrix(seq(0.5,3, length=100)) 
length=100
# Calcular r dado x
r <- apply(x,1,encontra_raiz)
#Figura
jpeg("01.jpg")
plot(x,r,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, r",las=1,lwd=4)
dev.off()
 
 
#Fazendo a figura com nossa equação integrada
 
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S)
}
 
# Função para achar r dado x
acha_raiz <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Explorando x
x <- matrix(seq(0.5,3, length=100))
r <- apply(x,1,acha_raiz)
#Figura
plot(x, r, type="l", xlab="Tamanho, x", ylab="Fitness, r",las=1,lwd=4)
 
#####################################
## Achando o maximo com calculo
#####################################
 
##uniroot
funcao <- function(x) {
    Af <- 0;
    Bf<-4*4;
    As<-1;
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r))
}
#
uniroot(f=funcao,interval=c(1.2,1.8))$root
 
## Usando nlm
funcao <- function(x) {
    Af<-0
    Bf<-4*4
    As<-1
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(abs(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r)))
}
#
nlm(funcao, p=1.2)$estimate
 
## Optimize
optimize(f=funcao, interval = c(1.2,1.8),maximum= FALSE)$minimum
 
#####################################
## Usando métodos númericos
#####################################
 
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S) #lembrando que aqui ja é a integral
}
 
# Achando r em função de x
raiz_funcao <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Achando o maximo
optimize(f = raiz_funcao, interval = c(.5,3),maximum = TRUE)$maximum
 
#####################################
## Usando integração númerica
#####################################
 
funcao_integrar <- function(idade,x,r) {
    Af <- 0
    Bf <- 4*4
    As <- 1
    Bs <- 0.5
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função que calcular a integral
integral <- function(r,x){
    1-integrate(funcao_integrar,1,Inf,x,r)$value
}
# Função para achar r dado x
raiz_integral <- function(x) {
    uniroot( integral, interval=c(1e-07,10),x)$root
}
 
#Achando o maximo
optimize(f = raiz_integral, interval = c(1.2,1.8),maximum = TRUE)$maximum