Apêndice 2

Modelos populacionais matriciais aplicados à ecologia: Uma introdução à teoria e a prática - Em preparação

Gabriel S. Santos https://gabeco.wixsite.com/home (Universidade do Estado do Rio de Janeiro)
2021-05-24

O que você está vendo aqui é a versão preliminar do artigo “Modelos populacionais matriciais aplicados à ecologia: Uma introdução à teoria e a prática”. Clique em Artigo ou Apêndice 1 para acessar cada um dos respectivos documentos ou vá no meu repositório do Github aqui

Configurações iniciais

Vamos começar limpando todos os dados que porventura estejam no R. Se algo muito ruim estiver acontencendo no seu R você pode sempre voltar nessa função para limpar seja lá qual erro você tiver cometido. O problema aqui é que esse comando vai limpar tudo que você fez até o momento. Seja prudente!


rm(list=ls())

Instalar os pacotes necessáros


install.packages("diagram")
install.packages("popbio")
install.packages("cowplot")
install.packages("tidyverse")
install.packages("ggrepel")
install.packages("popdemo")

library(diagram)
library(popbio)
library(cowplot)
library(tidyverse)
library(ggrepel)
library(popdemo)

Instalando os bancos de dados COMADRE e COMPADRE

clique aqui para ir para o site do COMPADRE e COMADRE


#É mais fácil se você escolher a pasta onde o arquivo está usando a função setwd
setwd(choose.dir())

#Depois basta carregar o arquivo 
load( "COMADRE_v.3.0.0.RData" )

#OU use load(file.choose()) caso você tenha algum problema para encontrar o arquivo no seu computador

MPM

Escolhendo uma matrix

Para facilitar a nossa vida, vamos buscar no banco de dados do COMADRE o estudo referente a espécie que temos utilizado como modelo até o momento: uma população de Lontra canadensis estudada por Gorman et al. 2007.


comadre$metadata%>%
filter(SpeciesAccepted == "Lontra canadensis")%>%
  select("SpeciesAccepted","Genus","Authors","Journal","DOI.ISBN")
SpeciesAccepted Genus Authors Journal DOI.ISBN
1165 Lontra canadensis Lontra Gorman; McMillan; Erb; Deperno; Martin Am Midl Nat 10.1674/0003-0031(2008)159[98:SACMOA]2.0.CO;2

A matrix que temos interesse está na linha 1165 dos metadados do COMADRE. Então é essa a identidade da matriz que temos interesse. Vamos armazenar a identidade dessa matriz em um vetor:


ID<-1165

Agora que sabemos a identidade dessa matriz podemos chamá-la a qualquer momento. Repare que na verdade temos quatro matrizes diferente = matA,matF,matU e matC.


comadre$mat[ID]

[[1]]
[[1]]$matA
        A1    A2    A3
[1,] 0.000 0.170 0.650
[2,] 0.598 0.000 0.000
[3,] 0.000 0.598 0.765

[[1]]$matU
        U1    U2    U3
[1,] 0.000 0.000 0.000
[2,] 0.598 0.000 0.000
[3,] 0.000 0.598 0.765

[[1]]$matF
     F1   F2   F3
[1,]  0 0.17 0.65
[2,]  0 0.00 0.00
[3,]  0 0.00 0.00

[[1]]$matC
     C1 C2 C3
[1,]  0  0  0
[2,]  0  0  0
[3,]  0  0  0

A Matriz A, que é a matriz de interesse aqui, na verdade é o resultado das taxas vitais em cada classe etária (sobrevivência e fecundidade) e que podem ser representadas cada uma como uma matriz separada.É exatamente isso que temos nos bancos de dados COMADRE e COMPADRE. Lá as matrizes são diferenciadas em matriz de fecundidade (matF), matriz de sobrevivência (matU) e matriz de de clonalidade (matC).

\[A = U + F + C\] Vamos salvar nossa matriz A apenas.


A<-comadre$mat[[ID]]$matA

#E vamos também adicionar o nome de cada estágio do nosso MPM
name <- c("1 ano","2 Anos",">2 Anos")
rownames(A)<-colnames(A)<-name

A

        1 ano 2 Anos >2 Anos
1 ano   0.000  0.170   0.650
2 Anos  0.598  0.000   0.000
>2 Anos 0.000  0.598   0.765

#Tenha em mente que aqui estamos adaptando o MPM existente em Gorman et al. 2008.
#Apesar disso os nomes originais de cada classe etária pode ser adquirodo através do seguinte comando:
#colnames(A)<-do.call(rbind,strsplit(comadre$matrixClass[[ID]]$MatrixClassAuthor, "\\("))[,1]

Representação gráfica do nosso MPM

\[\begin{array}{cc} & \begin{array}{ccc} 1 ano & 2 anos & >2 anos \end{array} \\ \begin{array}{ccc} \\ A= \\ \end{array} & \left( \begin{array}{ccc} 0 & f_{2}*p_{0} & f_{3}*p_{0}\\ S_{1,2} & 0 & 0\\ 0 & S_{2,3} & S_{3,3} \end{array} \right) \end{array}\] 

Numgenerations <- dim(comadre$mat[[ID]]$matA)[1]
DiffMat <- matrix(data = 0, nrow = Numgenerations, ncol = Numgenerations)
AA <- as.data.frame(DiffMat)
#REP
 AA[[1,2]] <- "f[list(2)]*p0"
 AA[[1,3]] <- "f[list(3)]*p0"
 #
 AA[[2,1]] <- "s[list(1,2)]"
 AA[[3,2]] <- "s[list(2,3)]"
 AA[[3,3]] <- "s[list(3,3)]"
 # 
 
diagram::plotmat(A = AA, pos=3, curve = 0.4, name = name, lwd = 2,
 arr.len = 0.2, arr.width = 0.2, my = .1,
 box.size = 0.07, arr.type = "triangle", dtext = .3,cex.txt=1,
 main = NULL)

Podemos também adicionar os valores da matriz


diagram::plotmat(A = A, pos=3, curve = 0.4, name = name, lwd = 2,
 arr.len = 0.2, arr.width = 0.2, my = .1,
 box.size = 0.07, arr.type = "triangle", dtext = .3,cex.txt=1,
 main = NULL)

Usos do MPM

Projetando a população no tempo

Vamos primeiro criar uma condição inicial para a nossa população configurando a quantidade de organismos em cada classe etária


#N0=round(stable.stage(A)*(100+round(rnorm(3,70,50),0)),0)
N0=round(10*abs(rnorm(3,1,10)),0)
Classes N0
1 ano 17
2 anos 106
>2 anos 4

Projeta a população do tempo t para o tempo t+1


N1= A%*%N0

relembre que isso é uma simplificação então tudo bem ter 0,5 individuos! Para deixarmos o processo mais realista cabe a nós arrendondarmos os valores quando for pertinente! Para isso utilizandos a função "round


round(N1,0)
Classes N0 N1
1 ano 17 21
2 anos 106 10
>2 anos 4 66

E podemos projetar mais uma vez


N2= A%*%N1
round(N2)
Classes N0 N1 N2
1 ano 17 21 45
2 anos 106 10 12
>2 anos 4 66 57

Também podemos projetar o tamanho populacional quantas vezes quizermos utilizando o pacote popbio


futuro<-pop.projection(A,N0, 16)
plot(futuro$pop.sizes,type="l",
     xlab="Anos projetados",
     ylab="Tamanho populacional")

Classes N0 N1 N2 N3 N4 N5 N6 N7 N8 N9 N10 N11 N12 N13 N14 N15
1 ano 1 ano 17 21 45 39 38 40 40 41 41 42 43 43 44 45 45 46
2 Anos 2 anos 106 10 12 27 23 23 24 24 24 25 25 25 26 26 27 27
>2 Anos >2 anos 4 66 57 51 55 56 56 57 58 59 60 61 62 63 64 65

O gráfico e a tabela acima são excelentes exemplos para discutirmos a idéia de estrutura etária estável. Isso porque tanto o crescimento do gráfico começa a estabilizar a partir do 5-6 ano assim como a proporção de indivíduos em cada classe etária passa a sofrer quase nenhuma mudança.Podemos utilizar a função stable.stage() se quisermos verificar a proporção de individíduos em cada classe etária e plotar essa proporção ao longo do tempo usando stage.vector.plot()


stable.stage(A)

    1 ano    2 Anos   >2 Anos 
0.3336842 0.1965536 0.4697621
Classes N0 N1 N3 N4 N5 N6 N13 N14 N15 Estrutura estável
1 ano 1 ano 0.13 0.18 0.31 0.29 0.29 0.35 0.37 0.35 0.34 0.33
2 Anos 2 anos 1.09 0.09 0.21 0.17 0.18 0.21 0.21 0.21 0.20 0.20
>2 Anos >2 anos 0.04 0.55 0.40 0.40 0.58 0.48 0.50 0.48 0.47 0.47

Também podemos avaliar o crescimento populacional de forma bastante objetiva aqui.


lambda(A)

\[\lambda = 1.01\]

Uma vez que o \(\lambda >1\) podemos concluir que a população está crescendo apesar de que numa velocidade bastante baixa.

As mesmas operações feitas por meio de uma única linha de comando podem ser feitas de forma detalhada passo a passo caso haja interesse de se aprofundar nos processos matemáticos explicados no apêndice 1 desse artigo.


#ou utilizando as fórmulas matriciais
ev <- eigen(A)                      #Calcula os autovalores e autovetores associados a cada autovalor
lmax <- which.max(Re(ev$values))    #Identifica o autovalor dominante (maior valor)


#Retorna o valor de lambda
print("Lambda")

[1] "Lambda"

round(ev$values[lmax],3)

[1] 1.015+0i

W <- ev$vectors                     #separa os autovetores dos autovalores
w <- abs(Re(W[, lmax]))             #extrai o autovetor associado ao autovalor dominante
print("Estrutura etária estável")

[1] "Estrutura etária estável"

round(w/sum(w),2)

[1] 0.33 0.20 0.47

                            #Uma vez extraido o autovetor precisamos reescalona-lo para que ele represente de fato a proporção de indivíduos em cada classe etária

par(mrow=c(1,2))
stage.vector.plot(futuro$stage.vectors,proportions = FALSE, col=2:4,main="Valores absolutos",xlab="Anos projetados",ylab="Número de indivíduos")
stage.vector.plot(futuro$stage.vectors,proportions = TRUE, col=2:4,main="Proporção",xlab="Anos projetados",ylab="Número de indivíduos")

Vamos sobrepor a projeção populacional utilizando a estrutura etária estável contra a projeção do crescimento populacional observada. A idéia aqui é relembrar como chegamos no conceito de autovalor e autovetor apresentado na equação 3 e 4 do apêndice 1.

\(A*w = \lambda * w\) ou \((A-\lambda I)*w = 0\)


#Vamos utilizar a quantidade de indivíduos projetadas até o 15º ano. 
#Aqui temos certeza que a população já alcançaou a estrutura etária estável!
futuro$stage.vectors[,15]

   1 ano   2 Anos  >2 Anos 
45.27621 26.66961 63.74022 

#Vamos criar um objeto para armazenar a projeção dessa população para os próximos 100 anos
trajetoria.observada<-pop.projection(A,futuro$stage.vectors[,15],100)$stage.vectors


#E por fim vamos criar um objeto para armazenar a projeção dessa populaçao utilizando apenas o conceito matemático do autovalor dominante e autovetor.
#Aqui é importante ter em mente que a estrutura etária soma 1 então precisamos colocar essa estrutura etária numa escala semelhante.
#Isso foi feito multiplicando a estrutura etária pelo número de indivíduos totais projetados até o ano 15 com o comando abaixo:
#"stable.stage(A)*sum(futuro$stage.vectors[,15]),100)$stage.vectors"

trajetoria.autovalor.estruturaestavel<-pop.projection(lambda(A),
      stable.stage(A)*sum(futuro$stage.vectors[,15]),100)$stage.vectors


par(mfrow=c(1,2))
stage.vector.plot(trajetoria.observada,proportions = FALSE, col=2:4,
      xlab="Anos projetados",ylab="Número de indivíduos")
title(expression(""~A~"* x"))

stage.vector.plot(trajetoria.observada,proportions = FALSE, col=2:4,
      xlab="Anos projetados",ylab="Número de indivíduos")
title(expression(""~lambda~"* w"))


#Você pode simplesmente remover o termo #"*sum(futuro$stage.vectors[,15]),100)$stage.vectors".
#Você vai ver que a trajetória continua exatamente a mesma 
#apesar dos eixos y estarem em escalas diferentes.

Outras aplicações do MPM

Por fim, vamos demonstar alguns comandos representados na tabela 1 do apêndice 1.

você pode buscar muito mais aplicações utilizando a função help do pacote popbio.Lá é possivel consultar várias aplicações demonstradas nos livros Caswell (2001) e Morris e Doak (2002).


help(popbio)

R0


net.reproductive.rate(A)

[1] 1.090777

Generation time


generation.time(A)

[1] 5.756156

Damping ratio


damping.ratio(A)

[1] 2.600902

Sensitivitidade e elasticidade

Opção 1: por derivada parcial
Sensitividade

sens(A)
\[\begin{array}{cc} & \begin{array}{ccc} 1 ano & 2 anos & >2 anos \end{array} \\ \begin{array}{ccc} \\ S= \\ \end{array} & \left( \begin{array}{ccc} 0 & 0.10 & 0.25\\ 0.3 & 0 & 0\\ 0 & 0.27 & 0.65 \end{array} \right) \end{array}\]
Elasticidade

elasticity(A)
\[\begin{array}{cc} & \begin{array}{ccc} 1 ano & 2 anos & >2 anos \end{array} \\ \begin{array}{ccc} \\ E= \\ \end{array} & \left( \begin{array}{ccc} 0 & 0.02 & 0.16\\ 0.18 & 0 & 0\\ 0 & 0.16 & 0.49 \end{array} \right) \end{array}\]
Opção 2: Por iteração

#Criamos um vetor com as perturbações que iremos aplicar nos elementos da matriz
perturb<-seq(0.01,.99,0.01)

#Criamos também algumas matrizes onde iremos armazenar as mudanças feitas
new.mx<-SensOfLambda<-matrix(NA,nrow=dim(A)[1]^2,ncol=length(perturb))

for (i in 1:dim(A)[1]^2){
    for(j in 1:length(perturb)){
A2<-A
if (A[[i]]==0) {
NULL
} else {
A2[[i]]<-perturb[j]
 new.mx[[i]]<-A2[[i]]
 SensOfLambda[i,j]<-lambda(A2)}
#print(j)
}
#print(i)
}

colnames(SensOfLambda)<-perturb
rownames(SensOfLambda)<-paste(c("s11","s21","s31","F12","s22","s32", "F13", "F23", "s33"))

Pronto, agora podemos demonstrar gráficamente essas mudanças


SensOfLambda

#Primeiro precisamos construir uma tabela para armazenar todas as informações em um único ambiente
as.data.frame(              #Transforma o contjunto de dados em uma tabela
  cbind(                    #Concatena (junta) as seguintes tabelas: 
    SensOfLambda,           #SensOfLambda = aqui foi armazenada o resultado de lambda calculado para cada observação
    Observado=as.vector(A)))%>%   #as.vector(A) serve para extrair os valores observados em cada elemento da matriz que será importante para marcar os valores observados em cada um dos gráficos que serão gerados abaixo
tibble::rownames_to_column(var = "Elementos")%>%
pivot_longer(!c(Elementos,Observado),names_to = "Perturbação", values_to = "Lambda")%>%
filter(Lambda != "NA")%>%
as.data.frame()%>%
mutate(Perturbação=as.numeric(Perturbação))%>%
mutate(Lambda.observado=lambda(A))%>%
ggplot(., aes(x=Perturbação,y=Lambda))+
geom_line()+
  geom_hline(yintercept = lambda(A),linetype="longdash",col="gray70")+
  geom_vline(aes(xintercept = Observado),linetype="longdash",col="gray70")+
ylab(expression(""~lambda~""))+
  facet_grid(.~Elementos,scales="free_y")+
theme_classic()+
 labs(subtitle = "Linha tracejada horizontal indica o valor de lambda observado \n e a linha tracejada vertical indica o valor observado do elemento da matriz")

Olhando a inclinação de cada uma das retas fica claro que a sobrevivência dos indivíduos com >2 anos é a que mais gera mudanças no crescimento populacional.

Note também que a relação entre a mudança de um parámetro frequentemente é não linear. A magnitude e forma dessa curvatura (côncava - \(\cap\) ;Convexa - \(\cup\) ) pode ser acessada a partir da derivada parcial de segunda ordem. Scripts para isso podem ser acessado no pacote demogR e algumas aplicações ecológicas e para manejo podem ser encontradas em Carslake et al. 2009 e Santos et al.2021 2021 para aplicações evolutivas.

NOTA: O uso da derivada parcial para calcular a sensitividade e elasticidade assume uma relação linear entre os valores dos elementos da matriz e a sua contribuição para o crescimento populacional. Como vimos no gráfico acima essa relação raramente é linear e por isso é importante ter em mente que as análises de perturbação são realistas apenas quando a perturbação é muito pequena (geralmente infinitesimal).

Figuras

Abaixo fornecemos os códigos para a elaboração de cada figura. As legendas para cada figura podem ser encontradas no texto principal e por isso não serão repetidas aqui

Figura 2

#=========================================================================
#       FIGURA 2
#=========================================================================
estoc<-round(rnorm(3,50,20),0)
#estoc<-c(28, 64, 59, -2)   #Gerado com a função acima!

lambda(A)

[1] 1.01521

N0<-stable.stage(A)*(100+estoc)
A%*%floor(N0)

          [,1]
1 ano   48.480
2 Anos  29.302
>2 Anos 68.597

#plot(pop.projection(A,N0, 15)$pop.changes,type="l")
#abline(h=lambda(A))
projecao<-pop.projection(A,N0, 15)
proj.vec<-t(projecao$stage.vectors)

colnames(proj.vec)<-c("1 ano","2 Anos",">2 Anos")

proj.vec<-reshape2::melt(proj.vec)
    colnames(proj.vec)<-c("t","Estágio","N")

#-------------------------------------------------------------
#   TOP
#-------------------------------------------------------------

TOP<-ggplot(proj.vec,aes(y=N,x=t,group=Estágio))+
geom_line(aes(color=Estágio),size=.5,type=2,linetype="longdash")+
geom_point(aes(color=Estágio,shape=Estágio),size=2)+
geom_segment(aes(x = 4.5, y = 110, xend = 0.2, yend = 110),
                  arrow = arrow(length = unit(0.5, "cm")))+
geom_segment(aes(x = 5.5, y = 110, xend = 14.5, yend = 110),
                  arrow = arrow(length = unit(0.5, "cm")))+
annotate("text", x = 10, y = 124, label = "Proporção estável \nentre estágios \n (crescimento assintótico)")+
annotate("text", x = 2, y = 124, label = "Proporção não estável \nentre estágios \n ")+
xlab(NULL)+
ylab("Nº de Indivíduos")+
ylim(0,130)+
 scale_x_continuous(limits = c(0, 14),breaks=c(0:14))+
geom_vline(xintercept=5,linetype="dotdash",color="gray70")+
theme_classic()+
theme(text=element_text(size=20))

#TOP
#-------------------------------------------------------------
# BOTTOM
#-------------------------------------------------------------
df.changes<-data.frame(Changes=c(pop.projection(A,N0, 15)$pop.changes),Time=c(1:14))

BOTTOM<-
ggplot(df.changes,aes(y=Changes,x=Time-.5))+
geom_line(size=.1,type=2,linetype="longdash")+
geom_point(size=2)+
 scale_x_continuous(limits = c(0, 14),breaks=c(0:14))+
xlab("Tempo projetado")+
ylab(expression(""~lambda~""))+
#ylim(0,120)+
geom_vline(xintercept=5,linetype="dotdash",color="gray70")+
geom_hline(yintercept=lambda(A),linetype="longdash",color="red")+
annotate("text", size = 7,col="red",,x = 0.01, y = lambda(A)+0.002, parse = TRUE, label =as.character(expression(paste("",lambda[1],""))))+
theme_classic()+
theme(text=element_text(size=20))

#BOTTOM
#-------------------------------------------------------------

legend <- get_legend(
  # create some space to the left of the legend
  TOP + theme(legend.box.margin = margin(0, 0, 0, 12))
)

plots <- align_plots(TOP+theme(legend.position="none"), BOTTOM, align = 'v', axis = 'l')
bottom_row <- plot_grid(plots[[2]],  nrow = 1)

plot_grid(plots[[1]], legend,bottom_row,NULL,rel_widths = c(1, .3),rel_heights = c(1, .4),ncol = 2)

cowplot::plot_grid(align_plots(TOP, BOTTOM, align = 'v', axis = 'l'))

Figura 6

#=========================================================================
#       Figura 6 - LTRE x Pfister plot
#=========================================================================
df<-data.frame(as.numeric(as.vector(sensitivity(A))),
as.numeric(as.vector(elasticity(A))),
(c("s11","s21","s31","F12","s22","s32", "F13", "F23", "s33")))
colnames(df)<-c("sensibilidade", "elasticidade","classes")
df$sensibilidade[df$elasticidade==0]<-0
df<-subset(df,elasticidade>0)
rownames(df)<-df$classes


nonzero<-length(A[A>0])
labels<-(c("s11","s21","s31","F12","s22","s32", "F13", "F23", "s33"))
nonzero.labels<-labels[AA>0]

projecoes.menos<-projecoes.mais<-matrix(0,ncol=5,nrow=50)


for (i in 1: nonzero){ 
  A2<-A
if(A2[i]>0){
A2[i]<-A2[i]+(A2[i]/10)
    }else{(NULL)}
projecoes.mais[,i]<-pop.projection(A2,N0, 50)$pop.sizes
}   
colnames(projecoes.mais)<-paste(nonzero.labels,"+10%")


#------------------------------------------------------------
#   SIMULANDO REDUÇÃO DE -10%
#------------------------------------------------------------
for (i in 1: nonzero){ 
  A2<-A
if(A2[i]>0){
A2<-A2[i]-(A2[i]/10)
    }else{(NULL)}
projecoes.menos[,i]<-pop.projection(A2,N0, 50)$pop.sizes
}   

colnames(projecoes.menos)<-paste(nonzero.labels,"-10%")

df.proj<-data.frame(
    pop.projection(A,N0, 50)$pop.sizes,
    cbind(projecoes.menos,projecoes.mais))

colnames(df.proj)<-c("projecao",
            colnames(projecoes.menos),
            colnames(projecoes.mais))

df.proj<-reshape2::melt(df.proj)
    colnames(df.proj)<-c("Cenario","N")

df.proj$time<-as.vector(replicate(11,0:49))


df.proj2<-df.proj
for(i in 1:dim(df.proj)[1]){
if(df.proj$time[i] > 0){
    df.proj2$N[i]<-df.proj$N[i]+rnorm(1,10,5)
    }
}

barplot(df[,2],names.arg=df$classes,ylab=expression("Contribução para o cresc. pop. ( "~lambda~")"))

ggplot(subset(df.proj2[-c(1:21),],time<21),aes(y=N,x=time,group=Cenario))+
geom_line(aes(color=Cenario),size=1,linetype="longdash",
    position=position_jitter(w=.1, h=.1))+
  xlab("Tempo projetado")+
  ylab("Tamanho populacional (N)")+
geom_line(data=df.proj2[1:21,],color="black",size=1,linetype="longdash")+
  theme_classic()

Figura 7

#=========================================================================
#       Figura 7 - LTRE
#=========================================================================
data(calathea)

#-----------------------------------
#Contribuição dos sítios amostrais
#-----------------------------------
plots.lamb<-lapply(calathea[-17],lambda)

df<-data.frame(lambda=unlist(
lapply(lapply(split(plots.lamb, rep(1:4,each=4)),unlist),mean)),
            sd=
unlist(lapply(lapply(split(plots.lamb, rep(1:4,each=4)),unlist),sd)))

lambs<-ggplot(df,aes(y=lambda,x=rownames(df)))+
geom_point()+
geom_linerange(aes(ymin=lambda-sd,ymax=lambda+sd))+
theme_classic()+xlab("Sítio amostral")+ylab("Taxa de crescimento ("~lambda~")")

#-----------------------------------
#Contribuição da variação anual
#-----------------------------------

lambs.anual<-data.frame(
lambda=unlist(
    lapply(lapply(
        split(plots.lamb, rep(1:4,4)),unlist),mean)),
sd=unlist(
    lapply(lapply(
        split(plots.lamb, rep(1:4,4)),unlist),sd)))%>%
ggplot(aes(y=lambda,x=rownames(df)))+
geom_point()+
geom_linerange(aes(ymin=lambda-sd,ymax=lambda+sd))+
theme_classic()+xlab("Ano")+ylab("Taxa de crescimento ("~lambda~")")

esquerda<-cowplot::plot_grid(lambs,lambs.anual,ncol=1,labels=c("A","C"))

#-----------------------------------
#LTRE
#-----------------------------------
#script retirado de help(LTRE) do pacote popbio

data(calathea)
calathea_pool<-calathea[['pooled']]

plots<- split(calathea[-17], rep(1:4,each=4))
plots<- lapply(plots, mean)
Cm<-LTRE(plots, calathea_pool)
pe<-sapply(Cm, sum)

##YEARS -- split recycles vector
yrs<-split(calathea[-17], 1:4)
yrs <- lapply(yrs, mean)
names(yrs)<-1982:1985
Cm<-LTRE(yrs, calathea_pool)
ye<-sapply(Cm, sum)

#=========================================================================
#       INTERVALO DE CONFIANÇA PARA AS ANÁLISES DE LTRE
#=========================================================================
##YEARS -- split recycles vector

yrs<-split(calathea[-17], 1:4)

years.LTRE<-matrix(nrow=100,ncol=length(yrs),dimnames=list(NULL,paste ("Ano",1982:1985)))
 for(i in 1:length(yrs)){
years.LTRE[,i]<-replicate(100,sum(LTRE(mean(sample(yrs[[i]])[-1]), calathea_pool)))
}


##PLOTS -- split recycles vector

plots<- split(calathea[-17], rep(1:4,each=4))

plots.LTRE<-matrix(nrow=100,ncol=length(plots),dimnames=list(NULL,paste ("Sítio",1:4)))
 for(i in 1:length(plots)){
plots.LTRE[,i]<-replicate(100,sum(LTRE(mean(sample(plots[[i]])[-1]), calathea_pool)))
}


plot.LTRE.sitios<-gather(data.frame(plots.LTRE))%>%
ggplot(aes(x=key, y=value)) + 
#    geom_col()+
 stat_summary(fun.data = mean_cl_normal, geom = "errorbar", mult = 1,size=.1)+
 stat_summary(fun.y = mean, geom = "bar") +
xlab("Plots")+
ylab("Contribuição para " ~lambda~ "")+
theme_classic()


plot.LTRE.anos<-gather(data.frame(years.LTRE))%>%
ggplot(aes(x=key, y=value)) + 
#    geom_col()+
 stat_summary(fun.data = mean_cl_normal, geom = "errorbar", mult = 1,size=.1)+
 stat_summary(fun.y = mean, geom = "bar") +
xlab("Anos")+
ylab("Contribuição para " ~lambda~ "")+
theme_classic()


cowplot::plot_grid(esquerda,
cowplot::plot_grid(plot.LTRE.sitios,plot.LTRE.anos,ncol=1,labels=c("B","D")))

Figura 8

#=========================================================================
#       Figura 8 - LTRE
#=========================================================================
x <- dim(calathea[[1]])[1]      #dimensão da matriz
data.frame(
    Parametros = as.vector(calathea_pool),
    Efeito=as.vector(mean(LTRE(lapply(yrs, mean),calathea_pool))),
    Variacao=as.vector(mean(lapply(yrs, var2))),
    elementos.matriz=paste("a", 1:x, rep(1:x, each = x), sep = ""),
    #variable="Efeito",
    Sensitivity=as.vector(popdemo::sens(calathea_pool)),
    Elasticity=as.vector(elasticity(calathea_pool)))%>%
filter(Parametros != 0)%>%
ggplot(aes(x=Variacao,y=abs(Efeito),label=elementos.matriz))+
geom_point(aes(size=Elasticity),alpha=.2)+
#geom_text(size=3)+
geom_text_repel(size=3)+
scale_x_log10(breaks=c(0.0001,0.01,0.1,1,10,100),labels=c(0.001,0.01,0.1,1,10,100))+
#scale_y_log10()+
xlab("Variação temporal")+
labs(size="Elasticidade")+
geom_hline(yintercept=0,linetype="dashed")+     #geom_hline(yintercept, linetype, color, size)
ylab("| Efeito da variação para o " ~lambda~ "|")+theme_classic()

Figura 9

#=========================================================================
#       Figura 9 - Demographic buffering hypothesis
#=========================================================================
data.frame(
    Parametros = as.vector(calathea_pool),
    Efeito=as.vector(mean(LTRE(lapply(yrs, mean),calathea_pool))),
    Variacao=as.vector(mean(lapply(yrs, var2))),
    elementos.matriz=paste("a", 1:x, rep(1:x, each = x), sep = ""),
    #variable="Efeito",
    Sensitivity=as.vector(popdemo::sens(calathea_pool)),
    Elasticity=as.vector(elasticity(calathea_pool)))%>%
filter(Parametros != 0)%>%
ggplot(aes(x=Sensitivity,y=Variacao,label=elementos.matriz))+
geom_point(aes(size=Elasticity),alpha=.2)+
#geom_text(size=3)+
geom_text_repel(size=3)+
scale_x_log10(breaks=c(0.0001,0.01,0.1,1,10,100),labels=c(0.001,0.01,0.1,1,10,100))+
scale_y_log10()+

  xlab("Sensitividade")+
labs(size="Elasticidade")+
#geom_hline(yintercept=0,linetype="dashed")+        #geom_hline(yintercept, linetype, color, size)
ylab("Variação temporal")+theme_classic()

Metadados COMADRE & COMPADRE


comadre$metadata%>%
as_tibble()%>%
filter(Continent == "S America")%>%
distinct(Authors,.keep_all = TRUE)%>%
select(-c( Genus:AngioGymno, AdditionalSource:MatrixPopulation,Altitude,Ecoregion:SurvivalIssue))
Table 1: Metadados COMADRE filtrado para América do Sul
SpeciesAuthor SpeciesAccepted CommonName Authors Journal YearPublication DOI.ISBN Lat Lon Country Continent
Astroblepus_ubidiai Astroblepus ubidiai Andean catfish Velez-Espino Ecol Freshwater Fish 2005 10.1111/j.1600-0633.2005.00084.x 0.1005833 -78.21724 ECU S America
Brachyrhaphis_rhabdophora Brachyrhaphis rhabdophora Live-bearing fish Johnson; Zuniga-Vega Ecology 2009 10.1890/07-1672.1 10.4015833 -85.08469 CRI S America
Genypterus_blacodes Genypterus blacodes Pink cusk-eel Gonzales-Olivares; Aranguiz-Acuna; Ramos-Jiliberto; Rojas-Palma Fish Res 2009 10.1016/j.fishres.2008.11.006 -49.0000000 -77.00000 CHL S America
Agaricia_agaricites Agaricia agaricites Tan lettuce-leaf coral Hughes; Tanner Ecology 2000 10.1890/0012-9658(2000)081[2250:RFLHAL]2.0.CO;2 NA NA JAM S America
Orbicella_annularis Montastraea annularis Caribbean star coral Edmunds Limnol Oceanogr 2015 10.1002/lno.10075 18.3166667 -64.71722 VIR S America
Plexaura_A Plexaura A Gorgonian coral Lasker Oecologia 1991 10.1007/BF00318316 9.5666667 -78.81667 PAN S America
Diomedea_melanophris Thalassarche melanophris Black-browed albatross Arnold; Brault; Croxall Ecol Appl 2006 10.1890/03-5340 -54.0006111 -38.03356 SGS S America
Forpus_passerinus Forpus passerinus Green-rumped parrotlets Sandercock; Beissinger J Appl Stats 2002 10.1080/02664760120108818 8.5666667 -67.58333 VEN S America
Homo_sapiens_subsp._sapiens Homo sapiens sapiens Human Keyfitz; Flieger Book 1990 0-226-43237-8 10.6500000 -61.51667 TTO S America
Panstrongylus_geniculatus Panstrongylus geniculatus NA Rabinovich; Feliciangeli J Med Entomol 2015 10.1093/jme/tjv112 NA NA VEN S America
Brachyteles_hypoxanthus Brachyteles hypoxanthus Northern muriqui Morris; Altmann; Brockman; Cords; Fedigan; Pusey; Stoinski; Bronikowski; Alberts; Strier Am Nat 2011 10.1086/657443 -19.7181111 -41.41689 BRA S America
Didelphis_aurita Didelphis aurita Common marsupial Ferreira; Kajin; Vieira; Zangrandi; Cerqueira; Gentile Mammal Biol 2013 10.1016/j.mambio.2013.03.002 -22.0346111 -42.68392 BRA S America
Hippocamelus_bisulcus Hippocamelus bisulcus Huemul deer Corti; Wittmer; Festa-Bianchet J Mamm 2010 10.1644/09-MAMM-A-047.1 47.2000000 -72.50000 CHL S America
Lagothrix_lagothricha Lagothrix lagotricha Humboldt’s woolly monkey Defler Book 2014 978-1-4939-0697-0 2.6666667 -74.16667 COL S America
Physeter_macrocephalus Physeter macrocephalus Sperm whale Whitehead; Gero Endangered Spp Res 2015 10.3354/esr00657 NA NA NODC-27 S America
Pseudalopex_culpaeus Lycalopex culpaeus Andean fox Novaro; Funes; Walker J Appl Ecol 2005 10.1111/j.1365-2664.2005.01067.x -40.0000000 -71.00000 ARG S America
Saguinus_fuscicollis Saguinus fuscicollis Saddlebacked tamarin Watsa PhD Thesis 2013 10.7936/K7DB7ZTD -12.5668611 -70.08492 PER S America
Oithona_hebes Oithona hebes Copepod Torres-Sorando; Zacarias; Zoppi de Roa; Rodriguez Ecol Model 2003 10.1016/S0304-3800(02)00355-1 10.2508333 -65.80083 VEN S America
Stratiodrilus_aeglaphilus Stratiodrilus aeglaphilus Freshwater crayfish Pardo; Vila; Bustamante Hydrobio 2008 10.1007/s10750-007-9136-8 -33.7340833 -70.90000 CHL S America
Podocnemis_lewyana Podocnemis lewyana Magdalena River Turtle Páez; Bock; Espinal-García; Rendón-Valencia; Alzate-Estrada; Cartagena-Otálvaro; Heppell Copeia 2015 10.1643/CE-14-191 NA NA COL S America
Podocnemis_expansa Podocnemis expansa Arrau turtle Mogollones; Rodríguez; Hernández; Barreto Chel Cons Biol 2010 10.2744/CCB-0778.1 6.4180278 -67.15125 VEN S America
Alouatta_seniculus_2 Alouatta_seniculus Red howler Wiederholt; Fernandez-Duque; Diefenbach; Rudran Ecol Model 2010 10.1016/j.ecolmodel.2010.06.026 8.5666667 -67.58333 VEN S America
Ara_glaucogularis_2 Ara_glaucogularis Blue-throated Macaw Maestri; Ferrati; Berkunsky Ecol Model 2017 10.1016/j.ecolmodel.2017.07.023 -14.3678333 -65.08458 BOL S America
Arctocephalus_australis Arctocephalus_australis South american fur seal Lima; Páez J Mamm 1997 10.2307/1382951 -35.0177222 -54.86819 URY S America 

compadre$metadata%>%
as_tibble()%>%
filter(Continent == "S America")%>%
distinct(Authors,.keep_all = TRUE)%>%
select(-c( Genus:AngioGymno, AdditionalSource:MatrixPopulation,Altitude,Ecoregion:SurvivalIssue))
Table 2: Metadados COMPADRE filtrado para América do Sul
SpeciesAuthor SpeciesAccepted CommonName Authors Journal YearPublication DOI.ISBN Lat Lon Country Continent
Andropogon_brevifolius Schizachyrium brevifolium NA Canales; Trevisan; Silva; Caswell Acta Oeco 1994 NA 9.4166667 -68.25000 VEN S America
Aspasia_principissa Aspasia principissa NA Zotz; Schmidt Biol Cons 2006 10.1016/j.biocon.2005.07.022 9.1666667 -79.85000 PAN S America
Heteropsis_flexuosa Heteropsis flexuosa NA Balcazar PhD thesis 2013 978-90-9027402-7 -4.0000000 -69.88333 COL S America
Tillandsia_flexuosa Tillandsia flexuosa NA Wester, Zotz J Trop Ecol 2010 10.1017/S0266467409990459 NA NA PAN S America
Werauhia_sanguinolenta Vriesea sanguinolenta NA Zotz Acta Oeco 2005 10.1016/j.actao.2005.05.009 9.1666667 -79.85000 PAN S America
Aechmea_nudicaulis Aechmea nudicaulis Nurse bromeliad Sampaio; Pico; Scarano Am J Bot 2005 10.3732/ajb.92.4.674 -22.2004444 -41.50003 BRA S America
Andropogon_semiberbis Schizachyrium sanguineum NA Silva; Raventos; Caswell; Trevisan J Ecol 1991 10.2307/2260717 8.6333333 -70.20000 VEN S America
Calathea_marantifolia Calathea marantifolia NA Matlaga PhD thesis 2008 NA 8.4680278 -83.58394 CRI S America
Calathea_micans Calathea micans NA Le Corff; Horvitz Ecol Model 2005 10.1016/j.ecolmodel.2005.05.009 10.4333333 -84.00000 CRI S America
Festuca_gracillima Festuca gracillima Tussock grass Oliva; Collantes; Humano Rang Ecol Manag 2005 10.2111/1551-5028(2005)58[466:DOGTGP]2.0.CO;2 -51.4003778 -69.61725 ARG S America
Heliconia_acuminata Heliconia acuminata NA Bruna Ecology 2003 10.1890/0012-9658(2003)084[0932:APPIFH]2.0.CO;2 -2.5000000 -60.00000 BRA S America
Heliconia_metallica Heliconia metallica NA Schleuning; Huamán; Matthies J Ecol 2008 10.1111/j.1365-2745.2008.01416.x -12.3505556 -70.70122 PER S America
Leptocoryphium_lanatum Leptocoryphium lanatum NA Raventós; Segarra; Acevedo Ecol Model 2004 10.1016/j.ecolmodel.2003.12.044 8.6333333 -70.20000 VEN S America
Periandra_mediterranea Periandra mediterranea NA Hoffmann; Solbrig Forest Ecol Manag 2003 10.1016/S0378-1127(02)00566-2 -15.9333333 -47.88333 BRA S America
Phaseolus_lunatus Phaseolus lunatus Lima bean Degreef; Baudoin; Rocha Gen Res Crop Evol 1997 10.1023/A:1008623521755 NA NA CRI S America
Syngonanthus_nitens Syngonanthus nitens Golden-grass Schmidt; Ticktin Biol Cons 2012 10.1016/j.biocon.2012.03.018 -1.0500000 -46.96667 BRA S America
Machaerium_cuspidatum Machaerium cuspidatum NA Nabe-Nielsen J Trop Ecol 2004 10.1017/S0266467404001609 0.6666667 -76.38333 ECU S America
Astrocaryum_aculeatissimum Astrocaryum aculeatissimum NA Portela; Bruna; dos Santos Biodivers Conserv 2010 10.1007/s10531-010-9846-5 NA NA BRA S America
Attalea_humilis Attalea humilis NA Souza; Martins Biodivers Conserv 2004 10.1023/B:BIOC.0000029326.44647.7f -22.5166667 -42.28333 BRA S America
Euterpe_edulis Euterpe edulis NA Silva-Matos; Freckleton; Watkinson Ecology 1999 10.1890/0012-9658(1999)080[2635:TRODDI]2.0.CO;2 -22.8179167 -47.10092 BRA S America
Euterpe_edulis_2 Euterpe edulis NA Freckleton; Matos; Bovi; Watkinson J Appl Ecol 2003 10.1046/j.1365-2664.2003.00842.x -22.8179167 -47.10092 BRA S America
Euterpe_oleracea Euterpe oleracea NA Arango; Duque; Muñoz Int J Trop Biol 2010 10.15517/rbt.v58i1.5222 NA NA COL S America
Euterpe_precatoria Euterpe precatoria NA Zuidema Book 2000 NA -10.9833333 -65.71667 BOL S America
Euterpe_precatoria_2 Euterpe precatoria NA Otárola; Avalos Am J Bot 2014 10.3732/ajb.1400089 10.4166667 -84.00000 CRI S America
Geonoma_deversa Geonoma deversa NA Zuidema; de Kroon; Werger Ecol Appl 2007 10.1890/1051-0761(2007)017[0118:TSBPAR]2.0.CO;2 -11.7500000 -67.33333 BOL S America
Geonoma_macrostachys Geonoma macrostachys NA Svenning Plant Ecol 2002 10.1023/A:1015520116260 0.6666667 -76.38333 ECU S America
Geonoma_orbignyana Geonoma orbignyana NA Rodríguez-Buriticá; Orjuela; Galeano Forest Ecol Manag 2005 10.1016/j.foreco.2005.02.052 4.3678650 -74.05129 COL S America
Geonoma_schottiana Geonoma schottiana NA Sampaio; Scariot J Trop Ecol 2010 10.1017/S0266467409990599 -15.6666667 -47.98333 BRA S America
Iriartea_deltoidea Iriartea deltoidea Paxiuba barriguda Pinard Biotrop 1993 10.2307/2388974 NA NA BRA S America
Lepidocaryum_tenue Lepidocaryum tenue NA Navarro; Galeano; Bernal Trop Cons Sci 2011 10.1177/194008291100400104 -4.0005556 -69.88486 COL S America
Mauritia_flexuosa Mauritia flexuosa Canangucho; Morete Holm; Miller; Cropper Biotrop 2008 10.1111/j.1744-7429.2008.00412.x 0.1166667 -76.18333 ECU S America
Phytelephas_seemannii Phytelephas seemannii NA Bernal J Appl Ecol 1998 10.1046/j.1365-2664.1998.00280.x NA NA COL S America
Adesmia_volckmannii Adesmia volckmannii NA Cipriotti; Aguiar Appl Veg Sci 2012 10.1111/j.1654-109X.2011.01138.x -45.6833333 -70.26667 ARG S America
Fabiana_imbricata Fabiana imbricata NA Ruete BSc thesis 2006 NA -41.0500000 -71.03333 ARG S America
Miconia_albicans Miconia albicans NA Hoffmann Ecology 1999 10.2307/177080 -15.9333333 -47.88333 BRA S America
Espeletia_spicata Coespeletia spicata NA Silva; Trevisan; Estrada; Monosterio Global Ecol Biogeogr 2000 10.1046/j.1365-2699.2000.00187.x 8.8666667 -70.91667 VEN S America
Acacia_pennatula Acacia pennatula NA Somarriba Agrofor Syst 2011 10.1007/s10457-011-9447-7 13.1833333 -86.26667 NIC S America
Actinostemon_concolor Actinostemon concolor NA Bianchini; de Araújo; Green; Pimenta Braz Arch Biol Technol 2013 10.1590/S1516-89132013000100009 -23.4500000 -51.25000 BRA S America
Araucaria_araucana Araucaria araucana NA Bekessy; Newton; Fox; Lara et al. Book 2004 NA NA NA CHL S America
Avicennia_germinans Avicennia germinans NA López-Hoffman; Ackerly; Anten; Denoyer;Ramos J Ecol 2007 10.1111/j.1365-2745.2007.01298.x 10.9678333 -71.73372 VEN S America
Bertholletia_excelsa Bertholletia excelsa Brazil nut tree Zuidema; Boot J Trop Ecol 2002 10.1017/S0266467402002018 NA NA BOL S America
Carapa_guianensis_2 Carapa guianensis NA Klimas; Cropper Jr.; Kainer; Wadt Ecol Model 2012 10.1016/j.ecolmodel.2012.07.022 -10.0174444 -67.70053 BRA S America
Cedrela_odorata Cedrela odorata Spanish cedar Zuidema; Brienen; During; Guneralp Am Nat 2009 10.1086/605981 -11.4000000 -68.71667 BOL S America
Chlorocardium_rodiei Chlorocardium rodiei Greenheart ter Steege; Boot; Brouwer; Hammond; Vanderhout; Jetten; Khan; Polak; Raaimakers; Zagt Ecol Appl 1995 10.2307/2269341 5.2166667 -58.80000 GUY S America
Dicorynia_guianensis Dicorynia guianensis NA Picard; Mortier; Chagneau Ecol Model 2010 10.1016/j.ecolmodel.2010.06.010 5.2666667 -52.91667 GUF S America
Dicymbe_altsonii Dicymbe altsonii NA Zagt; Boot PhD thesis 1997 NA 5.2175000 -58.80083 GUY S America
Eperua_falcata Eperua falcata NA Chagneau; Mortier; Picard J R Stat Soc C 2009 10.1111/j.1467-9876.2008.00657.x 5.2666667 -52.91667 GUF S America
Grias_peruviana Grias peruviana NA Peters Book 1991 NA -4.2166667 -74.21667 PER S America
Guettarda_viburnoides Guettarda viburnoides NA Loayza; Knight Ecology 2010 10.1890/09-0480.1 -14.6666667 -66.41667 BOL S America
Himatanthus_drasticus Himatanthus drasticus Plumel Baldauf; Correa; Ferreira; Santos Forest Ecol Manag 2015 10.1016/j.foreco.2015.06.022 -7.1833333 -39.21667 BRA S America
Pentaclethra_macroloba Pentaclethra macroloba NA Hartshorn PhD thesis 1972 NA 10.4333333 -83.98333 CRI S America
Prioria_copaifera Prioria copaifera El cativo Condit Forest Ecol Manag 1993 10.1016/0378-1127(93)90045-O 9.1502103 -79.83465 PAN S America
Prosopis_<U+FB02>exuosa Prosopis flexuosa NA Aschero; Morris; Vázquez; Alvarez; Villagra For Ecol Manag 2016 10.1016/j.foreco.2016.03.028 -34.3333333 -67.96667 ARG S America
Rhizophora_mangle Rhizophora mangle NA López-Hoffman; Ackerly; Anten; DeNoyer;Ramos J Ecol 2007 10.1111/j.1365-2745.2007.01298.x 10.9678333 -71.73372 VEN S America
Swietenia_macrophylla Swietenia macrophylla Big leaf mahogany Verwer; Peña-Claros; van der Staak; Ohlson-Kiehn; Sterck J Appl Ecol 2008 10.1111/j.1365-2664.2008.01564.x -15.7833333 -62.91667 BOL S America
Tachigali_vasquezii Tachigali vasquezii NA Poorter; Zuidema; Peno-Claros; Boot J Ecol 2005 10.1111/j.1365-2745.2005.00958.x -10.9833333 -65.71667 BOL S America
Vochysia_ferruginea Vochysia ferruginea NA Boucher; Mallono Forest Ecol Manag 1997 10.1016/S0378-1127(96)03890-X NA NA NIC S America
Broughtonia_cubensis Broughtonia cubensis NA Raventos; Gonzalez; Mujica; Bonet Biotropica 2015 10.1111/btp.12231 21.9336556 -84.50153 CUB S America
Euterpe_edulis_4 Euterpe edulis NA de Souza; Portela; de Mattos Ecol Evol 2018 10.1002/ece3.4686 -22.4174732 -43.00011 BRA S America
Lepanthes_acuminata Lepanthes acuminata NA Raventós; García-González; Riverón-Giró; Damon Plant Ecol Divers 2018 10.1080/17550874.2018.1444110 15.0842083 -92.15158 MEX S America
Lepanthes_caritensis_2 Lepanthes caritensis NA Tremblay Selbyana 1997 NA 18.0844306 -66.01777 PRI S America
Lepanthes_rubripetala_2 Lepanthes rubripetala NA Schödelbauerová; Tremblay; Kindlmann Biodivers Conserv 2010 10.1007/s10531-009-9724-1 18.2844889 -65.80001 PRI S America
Melocactus_bahiensis Melocactus bahiensis NA Hughes; Figueira; Jacobi; Borba Braz J Bot 2018 10.1007/s40415-018-0483-7 -16.5671325 -41.46823 BRA S America
Betula_pendula Betula pendula Silver Birch Maillette J Appl Ecol 1982 10.2307/2403005 NA NA BRA S America
Spathoglottis_plicata Spathoglottis plicata Philippine ground orchid; Philippine orchid; Large purple orchid Falcón; Ackerman; Tremblay Biol Invas 2017 10.1007/s10530-016-1318-8 18.3349444 -66.68344 PRI S America
Vellozia_aff._sincorana Vellozia sinconara Candomba Souza; Schmidt; Conceicao Flora 2017 10.1016/j.<U+FB02>ora.2017.02.007 -12.6179028 -41.50046 BRA S America