Simulações

O codigo do PI parece estar funcionando bem, inclui apenas como atributo do objeto de saída os valores de argumentos de entrada na simulação. Isso facilita a vida quando estamos apresentando os resultados e fazendo gráficos. Estou pensando também em criar uma versão onde o objeto de saída trás apenas as informações que usamos. Rodei quatro simulações e o workspace está demorando para subir. Isso pode complicar nossa vida depois. Uma alternativa viável é usar os dados intermediários dos ciclos ao invés de colocar tudo no arquivo de saída.

1. arquivo com a função simula.neutra_19jul2010.r
2. entre na página para ver os resultados da simulacões: result

• $S\ = \ $ número de espécies da comunidade : inteiro positivo
• $J\ = \ $ número de indivíduos da comunidade (tamanho da comunidade): inteiro positivo
• $x_i\ = \ $ número de propágulos que o indivíduo $i$ produz por intervalo
• $c_v\ = \ $ proporção de variação na produção de propágulos em relação a sua esperança: real positivo
• $X\ = \ $ total de propágulos que um indivíduo produz (esforço reprodutivo total): inteiro positivo

• Dado que o esforço reprodutivo total é o mesmo para todos os indivíduos ($X$), o tempo médio de vida de cada indivíduo é $$E[t_i]\ = \ X/x_i $$
• A cada intervalo cada indivíduo tem uma probabilidade fixa $p_i$ de morrer. Portanto, a probabilidade de sobreviver $t$ intervalos é dada por uma distribuição geométrica¹⁾ com suporte $0,\infty$. A esperança desta distribuição é $E[t_i]=p_i^{-1}$, portanto:
  • $$p_i \ = \ x_i/X$$
• Dado $s$ igual para todas as espécies, e uma distribuição gaussiana de caracteres contínuos, o número de sementes por intervalo $x_i$ da prole de um indivíduo pode ser definido como uma variável normal, com parâmetros
  • $$\mu \ = \ x_i$$
  • $$\sigma \ = \ c_v * mu\ $$

1. Defina $J$, $S$, $X$, $c_v$
2. Estabeleça as condições iniciais de acordo com o modelo de Hubbell:
   • Comece com o mesmo número de indivíduos por espécie, que será $J/S$.

²⁾

• As probabilidades de morrer são as mesmas para todos os indivíduos, e o número esperado de mortes por intervalo é um. Portanto temos para todos os indivíduos:
  • $$p_i \ = \ 1/J $$
  • $$x_i0 \ = \ X/J $$, pois $$E[t_i] \ = \ J$$ : inteiro

1. Atribua a cada indivíduo sua espécie, e os valores de $p_i$ e $x_i$
2. Gere $x_i$ propágulos para cada indivíduo, mantendo a identidade da mãe de cada um.
3. Sorteie os indivíduos que morrem:

• Gere $u_i$, números sorteados de um distribuição uniforme entre zero e um, sendo que $i$ vai de $1$ a $J$ ;
• Morrem todos os indivíduos em que $p_i>u_i$

1. Conte o número de mortos, $D$
2. Sorteie sem reposição $D$ propágulos do total de propágulos disponíveis, para ocupar o lugar dos mortos, .
3. Para a mãe de cada propágulo sorteado, sorteie com reposição um pai, entre os indivíduos da mesma espécie.
4. Calcule para cada combinação de mãe e pai a média do número de propágulos que eles produzem por intervalo:

• $$\bar x_{ij} = \frac{x_i + x_j}{2}$$

1. Sorteie o valor de $x_i$ de cada filhote de uma normal com parâmetros:

• $$\mu = \bar x_{ij}$$
• $$\sigma = \bar x_{ij} * c_v$$

1. Armazene para cada indivíduo:

• Identidade da espécie
• $$x_i$$

1. Armazene para cada ciclo:

• $$D$$

1. Calcule para cada novo indivíduo:

• $$p_i \ = \ x_i/X$$

1. Retorne ao passo 4

• Código inicial do Alexandre, com as seguintes modificações feitas pelo Paulo Inácio:
  • A fertilidade dos filhos é tomada de uma normal discretizada truncada, com parâmetro $$mu$$ = média das fertilidades dos pais, e coeficiente de variação definido pelo usuário (ver discussão). Os limites do truncamento são 1 e X.
  • Pai sorteado entre os co-específicos (código original soretava entre todos os indivíduos).
  • Ao invés do usuário definir o total de indivíduos (J), agora há um argumento para definir o número de indivíduos por espécie (j). O total de indivíduos é calculado pela função como o produto $$J = jS$$, garantindo que J seja um múltiplo de S.
  • Opção de gravar os resultados a intervalos regulares de tempo, que é útil para não estourar memória nas simulações longas. O argumento ciclo controla o tempo total da simulação, e o argumento step o intervalo para gravação.
• Argumentos da função:
  • S : Número de espécies inicial
  • j : número de individuos por espécie
  • X : esforço reprodutivo total³⁾
  • cv : coeficiente de variação fenotípica da fertilidade individual
  • ciclo : número total de ciclos que devem ocorrer. Em cada ciclo os adultos estão sujeitos a morrer, de acordo com a probabilidade que é fertilidades/esforço reprodutivo.
  • step : intervalos para gravação dos resulatdos. O dados da comunidade são gravados a cada step ciclos. Para gravar em cada ciclo, faça step=1

• Saída da Função: uma lista com so seguintes objetos:
  • tempo: vetor com rótulos do tempo passado em ciclos a cada instante gravado
  • sp.list: matriz de J linhas por ciclos/step colunas, com a identidade de cada indivíduso a cada tempo gravado.
  • sementes : matriz de J linhas por ciclos/step colunas, com a fertilidade (n de propágulos produzidos/ciclo) de cada indivíduso a cada tempo gravado.
  • prob.morte: matriz de J linhas por ciclos/step colunas, com a probabilidade de morte de cada indivíduso a cada tempo gravado.
  • n.mortes: vetor com o número total de mortes ao longo de cada intervalo de step ciclos que é gravado.

simula.neutra.step=function(S= 100, j=10, X=1000, cv=0.1, ciclo=1e6, step=100){ 
  ## Tamanho da comunidade
  J <- S*j
  ##Verifica se X e multiplo de J
  if(abs(X/J - round(X/J)) > .Machine$double.eps^0.5)
    {stop("\n\tO potencial reprodutivo (X) precisa ser multiplo do tamanho da comunidade (J). Tente novamente!\n\n")} 
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step) 
  ##Matriz de propagulos produzidos por individuo em cada ciclo
  prop.mat=matrix(nrow=J,ncol=1+ciclo/step)
  ## Matriz de probabilidade de morte de cada individuo, por ciclo
  dead.mat=matrix(nrow=J,ncol=1+ciclo/step)
  ##Vetor com numero de mortos por ciclo
  n.dead <- c()
  n.dead[1] <- 0
  ##CONDICOES INICIAIS##
  ##Deduzidas de acordo com o modelo de Hubbell:
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ind.mat[,1] <- rep(1:S,each=j)
  cod.sp <- ind.mat[,1]
  ##A probabilidade inicial de morte, tomada de uma geométrica,
  ##dado que o numero de mortes esperado por ciclo eh a mesma para todos (p=1/J)
  dead.mat[,1] <- 1/J
  p.death <- dead.mat[,1]
  ##O esforco reprodutivo inicial e deduzido da esperanca de tempo de vida da geometrica E[t]=J
  ## o que portanto resulta em X/J propagulos produzidos a cada ciclo, por todos
  prop.mat[,1] <- X/J
  n.propag <- prop.mat[,1]
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    n.mortes <- 0
    for(j in 1:step){
      ## Sorteio dos que morrerao
      morte=rbinom(J, 1, prob=p.death)
      ##Total de mortos, que e armazenado em n.dead
      D=sum(morte)
      n.mortes <- n.mortes+D
      ## Se ha mortos comeca aqui a rotina de substituicao dos valores
      if(D>0) 
        {
           ##vetor de propagulos: cada propagulo tem o codigo numerico do individuo
          seed.bank <- rep(1:J,n.propag)
          ##Indices dos individuos que morreram
          nascer= which(morte==1)
          ##Sorteio dos propagulos que irao repor os mortos
          mami=sample(seed.bank,D)
          ##Um vetor para armazenar o fenotipo do pai
          papi <- c()
          ##Um loop para sortear o pai entre os individuos da especie de cada
          ## propagulo-mae sorteado
          for(w in 1:D){
            papi[w] <- sample(n.propag[ cod.sp==cod.sp[mami[w]] ],1)
          }
          ##O valor esperado de propagulos dos filhotes eh
          ## a media do numero medio de propagulos produzidos pelo pai e pela mae
          medias.prop=(n.propag[mami]+papi)/2
          ##Codigos das especies dos mortos sao substituidos pelos nascidos
          ##dos propagulos sorteados
          cod.sp[nascer]<-cod.sp[mami]
          ## Numero medio de propagulos dos novos individuos nascidos e sorteada
          ## de uma normal discretizada e truncada entre 1 e J,
          ## com o coeficiente de variacao estabelecido
          n.propag[nascer] <- sapply(medias.prop,rnormt,cv=cv,min=1,max=X)
          ##A matriz de probabilidades de morrer eh atualizada para os novos individuos
          p.death[nascer] <- n.propag[nascer]/X
        }
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
    dead.mat[,i] <- p.death
    prop.mat[,i] <- n.propag
    n.dead[i] <- n.mortes
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  colnames(dead.mat) <- tempo
  colnames(prop.mat) <- tempo
  names(n.dead) <- tempo
  resulta=list(tempo=tempo,sp.list=ind.mat,sementes=prop.mat,prob.morte=dead.mat,n.mortes=n.dead)
  return(resulta)
}

Esta é a função simula.neutra.step editada para retirara o tradeoff reprodutivo e a variação neste caracter. O resulatdo é o modelo neutro do Hubbell, com a diferença que o número esperado de mortes a cada ciclo é um, mas há variação estocástica nisto. no modelo original, há sempre uma morte por ciclo.

simula.neutra.hub=function(S= 100, j=10, ciclo=1e4, step=100){ 
  ## Tamanho da comunidade
  J <- S*j
  ##Matrizes para guardar os resultados
  ## matriz da especie de cada individuo por ciclo
  ind.mat=matrix(nrow=J,ncol=1+ciclo/step) 
  ##Vetor com numero de mortos por ciclo
  n.dead <- c()
  n.dead[1] <- 0
  ##CONDICOES INICIAIS##
  ##Deduzidas de acordo com o modelo de Hubbell:
  ## Todas as especies comecam com o mesmo numero de individuos (j=J/S)
  ind.mat[,1] <- rep(1:S,each=j)
  cod.sp <- ind.mat[,1]
  ##A probabilidade inicial de morte, tomada de uma geométrica,
  ##dado que o numero de mortes esperado por ciclo eh a mesma para todos (p=1/J)
  p.death <- 1/J
  ##Aqui comecam as simulacoes
  for(i in 2:(1+ciclo/step)){
    n.mortes <- 0
    for(j in 1:step){
      ## Sorteio dos que morrerao
      morte=rbinom(J, 1, prob=p.death)
      ##Total de mortos, que e armazenado em n.dead
      D=sum(morte)
      n.mortes <- n.mortes+D
      ## Se ha mortos comeca aqui a rotina de substituicao dos valores
      if(D>0) 
        {
          ##Indices dos individuos que morreram
          nascer= which(morte==1)
          ##Quem substitui os mortos
          novos <- sample(1:J,D,replace=T)
          ##Substituindo
          cod.sp[nascer]<-cod.sp[novos]
        }
    }
    ## A cada step ciclos os resultados sao gravados
    ind.mat[,i] <- cod.sp
    n.dead[i] <- n.mortes
  }
  tempo <- seq(0,ciclo,by=step)
  colnames(ind.mat) <- tempo
  names(n.dead) <- tempo
  resulta=list(tempo=tempo,sp.list=ind.mat,n.mortes=n.dead)
  return(resulta)
}

Sorteia n números inteiros de uma aproximação discreta da normal truncada em seus limites inferiores e superiores dados sua média e coeficiente de variação.

• Pseudocódigo:
  • Dados os limites de truncamento (no caso um mínimo de um filhote e máximo de X filhotes), divida este intervalo contínuo em trechos de tamanho um, cujo ponto médio são os valores inteiros (e.g. 0.5-1.5, 1.5-2.5 …).
  • A probabilidade associada a cada um destes intervalos é calculada para uma normal com os parâmetros especificados.
  • Use estas probabilidades para sortear com reposição os números inteiros que são os pontos médios de cada intervalo.

• Argumentos:
  • mean : parâmetro $$mu$$ da normal, que corresponde à média
  • n : número de sorteios
  • cv : coeficiente de variação (desvio padrão / média)
  • min , max : limites mínimo e máximo de truncamneto

##Funcao de sorteio de uma normal discretizada e truncada
rnormt <- function(mean,n=1,cv,min,max){
  vals <- (min-0.5):(max+0.5)
  p <- pnorm(vals,mean=mean,sd=mean*cv)
  p2 <- diff(p)
  sample(min:max,n,prob=p2,replace=T)
}

Função simples que retorna o número de elementos únicos de um vetor alfanumérico. Aplicado ao vetor de coódigo de espécie de cada indivíduo retorna a riqueza.

rich <- function(x)length(unique(x))

Dado um vetor de abundãncia de espécies, faz um gráfico simples de Whitaker (log abundâncias x ranque).

w.plot <- function(x,...){
  plot(1:length(x),sort(x,decreasing=T),
       xlab="Rank order",ylab="Abundance",log="y",...)
}

Função feita para explorar graficamente a distribuição da fertilidade. Plota box-plots da fertlidade dos indivíduso de cada espécie, sobre um gr[áfico de rank-abundância da fertilidade de cada espécie, num dado instante gravado da simulação (argumento tempo).

box.w <- function(tempo,cod.sp,n.propag){
  sp <- cod.sp[,tempo]
  prop <- n.propag[,tempo]
  ab <- tapply(sp,sp,length)
  par(mfrow=c(2,1))
  boxplot(prop~factor(sp,levels=names(ab)[order(ab,decreasing=T)]))
  w.plot(ab)
  par(mfrow=c(1,1))
  }

A simulação grava uma matriz com as fertilidades de cada indivíduo a cada step ciclos. Esta função resume estes dados, calculando uma estatística descritiva (o default é a média) da fertilidade de cada espécie, em cada instante de tempo.

• Argumentos:
  • cod.sp: matriz retronada pela simulação de identidade de cada indivíduo x tempo
  • n.propag: matriz retronada pela simulação de fertilidade de cada indivíduo x tempo
  • fun: função a aplicar a cada espécie por tempo.

• Saída
  • Uma matriz das médias de fertilidade de cada espécie em cada tempo. Atenção: os valores nesta matriz são não existentes (NA) a patir do tempo em que a espécie de sextingue.

fert.t <- function(cod.sp,n.propag,fun=mean){
  especie <- unique(cod.sp[,1])
  medias <- matrix(ncol=ncol(n.propag),nrow=length(especie))
  colnames(medias) <- colnames(n.propag)
  rownames(medias) <- especie
  for(i in 1:ncol(n.propag)){
    medias[,i] <- tapply(n.propag[,i],factor(cod.sp[,i],levels=especie),fun)
  }
  medias
}

Novos tópicos com “new blog entry”. Para discutir um tópico, use o link comments.

blog:neutral1