Crescimento denso-independente com estocasticidade ambiental - Roteiro no R

No modelo de crescimento discreto independente de densidade projetamos o tamanho população multiplicando o seu tamanho inicial $N_0$ pela mesma taxa de crescimento $\lambda$ em todas as $t$ gerações:

\begin{equation} N_t=N_0 \lambda^t \label{eq1} \end{equation}

Com isso, estamos assumindo que a taxa de crescimento é a mesma a cada geração, o que é pouco realista. Condições e recursos variam ao longo do tempo, o que deve fazer $\lambda$ flutuar. Se aceitamos isso, e chamando de $\lambda_i$ os valores da taxa de crescimento a cada geração¹⁾, nosso modelo passa a ser:

\begin{equation} N_t=N_0 \lambda_1 \lambda_2 \lambda_3 \ldots \lambda_t \ = \ N_0 \prod_{i=1}^t \lambda_i \label{eq2} \end{equation}

Ainda assim, teremos uma taxa média de crescimento, que poderíamos estimar com a média das taxas observadas ao longo das gerações. Se o ambiente é pouco variável, as taxas observadas deveriam estar próximas dessa média a maior parte do tempo, ou seja, os valores $\lambda_i$ teriam pouca variância²⁾. Se há muita variação ambiental, isso deve se refletir em taxas de crescimento mais variáveis. Chamamos de estocasticidade ambiental a incerteza no valor de taxas demográficas, devido às flutuações nas condições e recursos experienciadas por todos os indivíduos da população.

De que maneira a estocasticidade ambiental afeta as projeções de tamanho populacional? Vamos responder isso simulando o modelo de crescimento discreto (equação $\ref{eq2}$), com valores de $\lambda_i$ que mudam a cada geração. Escolhemos este modelo por razões didáticas e de implementação computacional, mas as conclusões são as mesmas para modelos em tempo contínuo.

Abaixo os parâmetros do nosso modelo:

Cada linha colorida é a projeção no tempo do tamanho de uma população de acordo com o modelo ($\ref{eq2}$). Todas as populações têm o mesmo valor médio de taxa de crescimento e mesmo tamanho inicial. As projeções diferem porque a cada intervalo de tempo é sorteado um valor da taxa de crescimento independente para cada população.

Os sorteios são feitos de uma distribuição de probabilidades chamada lognormal, que garante que as taxas serão sempre positivas, como exige o modelo. Você pode mudar a média e variância da lognormal usada no sorteio (veja a seguir). Quanto maior a estocasticidade ambiental, maior a variância da taxas de crescimento.

A linha preta indica a projeção de uma população sob uma taxa constante igual à da média da distribuição de taxas. Apesar do modelo ser discreto, representamos as projeções com linhas para facilitar seguir a trajetória das populações no tempo.

Nosso objetivo é entender os efeitos da estocasticidade demográfica sobre o crescimento discreto independente da densidade. Vamos começar avaliando o comportamento de uma população por vez:

1. Altere o valor de Maximum time (tmax) para 10;
2. Altere o valor da opção Number of stochastic simulations (npop) para 1;
3. Altere o valor de Population growth rate (lambda) para 2;

1. Repita algumas vezes a operação e observe o gráfico.
2. Repita com valores mais altos de estocasticidade ambiental: para isso, aumente o valor em lambda variance(varr) e observe o gráfico.

A cada vez que você roda a função, o EcoVirtual sorteia uma valor de $\lambda$ para cada intervalo de tempo e os utiliza para projetar o tamanho da população ao longo do tempo, conforme a equação $\ref{eq2}$.

Os sorteios são tomados da mesma distribuição de probabilidade, com a média e variância especificadas. Os sorteios são independentes, isto é, o valor em um momento não afeta o valor sorteado para os intervalos seguintes. Por isso dizemos que as taxas sorteadas são variáveis aleatórias independentes mas de igual distribuição de probabilidades.

1. O que é preciso para simular ausência de estocasticidade ambiental?
2. Qual o comportamento da projeção na ausência de estocasticidade ambiental? Verifique sua previsão com o EcoVirtual

Uma maneira de lidar com processos variáveis é descrever seu comportamento médio. Por exemplo, aqui está um gráfico em tempo real do congestionamento na cidade de São Paulo ao longo das horas do dia de hoje, em comparação com as médias históricas dos menores e maiores valores de cada horário. O percentual de vias congestionadas varia a cada dia, mas isso não o torna completamente imprevisível. Se acumulamos dados de muitos dias podemos ter uma ideia do que esperar, e mesmo avaliar se estamos em um dia atípico. Este valor esperado, que designamos $E[X]$, é a média teórica de uma variável aleatória $X$ ³⁾, no caso o percentual da extensão de vias que estão congestionadas numa dada hora.

Em nossas simulações, os tamanhos populacionais são resultado da multiplicação de taxas de crescimento sorteadas sempre de uma mesma distribuição de probabilidades. Se esses produtos tiverem um valor esperado definido poderíamos prever o comportamento das populações, mesmo com a incerteza da estocasticidade ambiental. Será possível?

Vamos avaliar isso projetando várias populações independentes, para ter muitos valores de tamanhos populacionais a cada tempo. Para isso entre os seguintes valores na função:

# grave o resultado da função com o nome sim1

tmax = 51
npop = 1000
N0 = 10
lambda = 1.1
varr = 0.03
ext = FALSE

Rode a função e verá um gráfico com a evolução independente de 1000 populações ao longo do tempo, segundo o modelo ($\ref{eq2}$), para um $\lambda$ médio de 1,1 com variância de 0,03.

Ao colocar um nome para o resultado da simulação, você gravou as projeções das populações em cada tempo na memória do ambiente R. É uma tabela de dados em que as colunas são os tamanhos populacionais e as linhas os intervalos de tempo. É o que precisamos para calcular o tamanho médio das populações a cada tempo!

## Elimina as duas primeiras colunas da tabela (tempo e valor deterministico)
sim1b <- sim1[,-(1:2)]

## Calcula a media para cada tempo
medias=apply(sim1b,1,mean)

## Vetor de tempo de zero ao maximo
tempo <- length(medias)-1

## Grafico com medias
plot(0:tempo,medias, xlab="Tempo", ylab="N")

e agora rode estes comandos no R. Se tudo deu certo a janela gráfica vai mostrar uma figura como esta:

As médias parecem ter um crescimento geométrico, similar ao modelo de crescimento discreto sem estocasticidade. Vamos verificar isso para os cinco primeiros valores da média das projeções. Para vê-los, copie o comando abaixo na janela Script e rode-o:

# Visualizando as médias das projecoes de t=0 a t=5. 

medias[1:6]

Agora calcule as projeções para os mesmo tempos usando o $\lambda$ médio. O comando abaixo executa essa conta em único comando. O codigo 0:5 cria o vetor de valores 0, 1, 2, 3, 4, 5 …

## lembre-se que a populacional que simulou tinha 10 indivíduos inicialmente e o lambda de ''1.1''
 
10*1.1^(0:5)

Para ajudar na comparação entre estas duas projeções, podemos adicionar a que usa o $\lambda$ médio ao gráfico com os comandos

points(0:tempo, 10*1.1^(0:tempo), pch=19, col="blue", cex=0.5)

legend("topleft", c("medias observadas", "projecao com lambda medio"),
        pch=c(1,19), col=c("black", "blue"))

Dados os resultados que você obteve e o modelo de crescimento discreto com estocasticidade ambiental da equação $\ref{eq2}$:

$$N_t \ = \ N_0 \prod_{i=1}^t \lambda_i$$

e se definimos a média dos valores de $\lambda_i$ como

$$\bar \lambda \ = \ E[\lambda_i]$$

proponha a equação para o valor esperado do tamanho populacional $E[N_t]$ em função de $\bar \lambda$.

A seção anterior mostra que mesmo com estocasticidade ambiental as projeções em média são as mesmas do modelo sem estocasticidade.

É um resultado reconfortante mas, como toda síntese, a média pode ser traiçoeira, como diria o estatístico que se afogou no lago com profundidade média de 90 cm. A figura desse pobre homem da ciência, ao lado, é um link para um ótimo post com mais exemplos de médias enganadoras.

A informação que falta é quanta variação há nos dados. Uma das maneiras de medir isso é a variância, que expressa o desvio esperado dos dados em relação à média. Para uma variável aleatória $X$ a variância teórica é:

$$VAR[X] \ = \ E[(X - E[X])^2]$$

Note que a variância é a média do desvio em relação à média. Este desvio está elevado ao quadrado para ter algumas propriedades convenientes⁵⁾. Por isso a variância também é chamada desvio quadrático médio ou desvio quadrático esperado.

Mas há uma inconveniência de expressar a variação com desvios quadráticos: as unidades da variância são o quadrado das unidades da medida original. Se a média estiver em centímetros, a variância desta média estará em centímetros quadrados . Isso é resolvido tomando-se a raiz quadrada da variância, que chamamos desvio-padrão.

A expressão da variância acima é uma definição teórica. Na prática não sabemos valores esperados teóricos, mas podemos estimá-los. Fazemos isso de um conjunto de medidas, supondo que estas medidas são uma amostra representativa do processo cujo valor esperado queremos estimar. Assim, em uma amostra de $n$ medidas estimamos a variância com a expressão:

$$s^2 \ = \ \sum_{i=1}^n \frac{(x_i - \bar x)^2}{n-1}$$

onde $\bar x$ é a estimativa da média

$$\bar x \ = \ \frac{1}{n}\sum_{i=1}^n x_i $$

e $n$ o tamanho da amostra (seu número de elementos medidos). A estimativa do desvio-padrão é obtida com a raiz quadrada de $s^2$:

$$s \ = \ \sqrt{\sum_{i=1}^n\frac{(x_i - \bar x)^2}{n-1}} $$

Com isso podemos avaliar como as projeções de nosso modelo com estocasticidade variaram ao longo do tempo. Vamos usar as projeções da última simulação, que já guardados em uma tabela na memória do ambiente R. Para criar um gráfico da média e desvio-padrão em função do tempo copie e rode os comandos abaixo:

## calcula os desvios-padrão para cada tempo
desvios <- apply(sim1b,1,sd)

## Medias em funcao do tempo
plot(0:tempo,medias, xlab="Tempo", ylab="N")

#3 Acrescenta desvios-padrao ao grafico
points(0:tempo,desvios, xlab="Tempo", ylab="N", col="blue")

## Adiciona legenda ao grafico
legend("topleft", c("medias", "desvios-padrão"),
        pch=c(1,1), col=c("black", "blue"))

O desvio-padrão também cresce exponecialmente com o tempo! Este crescimento pode ser até mais rápido que o crescimento do tamanho populacional esperado, como para a simulação que usamos aqui.

Tomando o valor do desvio-padrão como uma medida de incerteza de nossa expectativa, podemos dizer que projeções a longo prazo são extremamente imprecisas, apesar de terem um valor esperado conhecido. Agora que você já sabe isso, rode de novo a simulação e veja o gráfico, que deve ser parecido com este

Note como a amplitude das projeções se abre como um funil à medida que o tempo passa. Isso acontece porque as trajetórias divergem com a multiplicação de taxas variáveis, mesmo que as populações partam do mesmo valor inicial.

A variância e o desvio-padrão são médias (dos desvios) e, portanto, também podem não caracterizar completamente a estrutura de variação dos dados. Considere este dois conjuntos de valores⁶⁾:

para os dois conjuntos a média é exatamente 9 e a variância exatamente 11. O que escapa aos dois índices é a distribuição dos valores em relação à média, que é bem mais assimétrica no conjunto B, em que apenas um valor é maior que a média.

Agora olhe de novo o gráfico das projeções da seção anterior e avalie a simetria da distribuição dos valores projetados em relação à média, indicada pela linha preta. Para ajudar, execute o código do R abaixo para fazer histogramas das projeções nos tempos 5, 10, 20 e 30:

## valor maximo para definir a escala do histograma
sim1b.m <- (max(sim1b[c(6,11,21,31),]))

##define uma janela para 4 gráficos
par(mfrow=c(2,2))

##Cria os 4 histogramas
hist(sim1b[6,], xlab="N", ylab="Frequencia", main="T=5",
     breaks=seq(0,max(sim1b[6,]+50),by=50), xlim=c(0,sim1b.m))
hist(sim1b[11,], xlab="", ylab="", main="T=10",
     breaks=seq(0,max(sim1b[11,]+50),by=50), xlim=c(0,sim1b.m))
hist(sim1b[21,], xlab="", ylab="", main="T=20",
     breaks=seq(0,max(sim1b[21,]+50),by=50), xlim=c(0,sim1b.m))
hist(sim1b[31,], xlab="", ylab="", main="T=30",
     breaks=seq(0,max(sim1b[31,]+50),by=50), xlim=c(0,sim1b.m))
     
## Volta a janela para 1 gráfico     
par(mfrow=c(1,1))

As distribuições das projeções tornam-se cada vez mais assimétricas com o passar do tempo. Isso acontece porque poucas populações crescem muito, enquanto a maioria permanece em tamanhos pequenos. Dependendo da média e variância das taxas de crescimento, muitas populações podem até diminuir de tamanho.

Assim, mesmo com taxa média de crescimento maior que um ($\bar \lambda>1$) a estocasticidade ambiental pode fazer com que a maioria das populações não cresça, ou mesmo diminua! Vamos verificar isso contando quantas populações passam do tamanho inicial em função do tempo em uma simulação em que $\lambda$ tem uma variância alta.

Execute as simulações abaixo na função de estocasticidade ambiental:

# salve os resultados da simulação no objeto 'sim2'
tmax = 51
npop = 1000
N0 = 10
lambda = 1.05
varr = 0.2
ext = FALSE

Os resultados estão guardados em uma nova tabela na memória do R, chamada sim2. Com ela podemos calcular a proporção das 1000 populações simuladas que tiveram tamanhos maiores do que $N_0=10$ a cada tempo. Para isso, copie e execute os comandos abaixo:

maiorN0 <- apply(sim2[,-(1:2)],1,
                 function(x)sum(x>10)/length(x))
                 
plot(sim2[-1,1],maiorN0[-1], xlab="Tempo", 
     ylab="Proporcao projecoes > N0")

A proporção de projeções maiores do que $N_0$ cai com o tempo. Isso indica que a probabilidade de uma população efetivamente crescer é cada vez menor, mesmo que em média haja crescimento .

Este aparente paradoxo se explica pela forte assimetria das distribuições das projeções populacionais: a maioria das populações está abaixo da média, o que se acentua com o tempo ⁷⁾. E mais: se a variância ambiental ⁸⁾ é grande, cada vez menos populações estarão acima de um valor abritrário qualquer, como o tamanho inicial da população.

Na função de estocasticidade ambiental, o parâmetro $ext$ se marcado (ext=TRUE) executas as simulações com $N_{min}=1$, ou seja, faz com que toda população que caia abaixo de um extinga-se⁹⁾. O número de populações que se extinguiram ao final da simulação é indicado abaixo do eixo X. Use esta opção para descobrir:

1. Utilize o tamanho inicial da população como $N_0 = 100$ e $N_0 = 1000$ para invetigar se o tamanho inicial das populações modificam o resultado das questões abaixo.
2. Qual o efeito do aumento de estocasticidade ambiental sobre a probabilidade de extinção?
3. Dado um certo nível de estocasticidade, qual o efeito de mais tempo de simulação sobre a probabilidade de extinção?
4. Deduza uma consequência desses resultados para a conservação de populações ameaçadas.

• On population growth in a randomly varying environment: neste artigo clássico, Richard Lewontin e Dan Cohen chamaram a atenção dos biólogos para as propriedades surpreendentes dos modelos que exploramos aqui.
• tutorial em R de multiplicação de variáveis aleatórias: veja como o modelo que usamos demonstra que o resultado desta multiplicação tende a uma distribuição lognormal. Da disciplina de de modelagem probabilística da pós-graduação em ecologia da USP.
• Life is log-normal : argumentação provocadora a favor do uso da log-normal ao invés da normal, com link para o artigo e uma demonstração animada.