Elasticity and Sensitivity of the Vital rates

Script

#Code to extract vital rates, vital rate sensitivities and vital rate elasticities from projection matrix models.
#Code developed by R. Salguero-Gomez (Max Planck Institute) modified after  E. Jongejans (Nijmegen University)
#Last modified: Mar 29th, 2012
#Modified by Maria Paula Balcazar-Vargas

library(fBasics)
library(Matrix)
library(MASS)
library(popbio)

DataSP<-read.table(matrix_file,header=T,sep=";", na.strings="NA", dec=".", strip.white=TRUE)

for (i in 1:dim(DataSP)[2]) {  
 if (i<4) 
DataSP[,i]=as.character(DataSP[,i]) 
 else 
DataSP[,i]=as.numeric(DataSP[,i]) 
}  
sapply(DataSP,mode)

startRows=1 
matDim=dimensions


SurSens=matrix(NA,matDim,matDim)
  GRSCSens=matrix(NA,matDim,matDim)
  GRSens=matrix(NA,matDim,matDim)
  SexSens=matrix(NA,matDim,matDim)
  Clo1Sens=matrix(NA,matDim,matDim)
  
  
  SurElas=matrix(NA,matDim,matDim)
  GRSCElas=matrix(NA,matDim,matDim)
  GRElas=matrix(NA,matDim,matDim)
  SexElas=matrix(NA,matDim,matDim)
  Clo1Elas=matrix(NA,matDim,matDim)
    
  ASurIndependent=matrix(NA,matDim,matDim)

sp=1   
#This reads the Sur (survival and change in stage, including falling probability), Sex (sexual reprod), and the three Clo (clonal reprod types: Clo1, Clo2, Clo3) matrices, following new nomenclature by Caswell and adapted to Maria Paula's species' crazy life:
 Sur=as.matrix(DataSP[,5:(4+matDim[sp])]) 
Sex=as.matrix(DataSP[,(5+matDim):((4+matDim)+matDim[sp])]) 
Clo1=as.matrix(DataSP[,(5+matDim+matDim):((4+matDim+matDim)+matDim[sp])])

 #Total matrix A (Note that in this life cycle one has to add all matrices... more than Sur, Rep and Clo!)
 A=Sur+Sex+Clo1

 #Proportion of each cell that corresponds to each vital rate
SurProp=(Sur/A) 
SexProp=(Sex/A) 
Clo1Prop=(Clo1/A)  

#Calculates the proportion of each matrix element that corresponds to each vital rate
SurProp=as.character(SurProp) 
SexProp=as.character(SexProp) 
Clo1Prop=as.character(Clo1Prop) 
 
for (i in 1:(matDim^2)) {
    if (SurProp[i]=="NaN") (SurProp[i]="0") else (SurProp[i]=SurProp[i])
    if (SexProp[i]=="NaN") (SexProp[i]="0") else (SexProp[i]=SexProp[i])
    if (Clo1Prop[i]=="NaN") (Clo1Prop[i]="0") else (Clo1Prop[i]=Clo1Prop[i])
}


SurProp=matrix(as.numeric(SurProp),matDim,matDim) 
SexProp=matrix(as.numeric(SexProp),matDim,matDim) 
Clo1Prop=matrix(as.numeric(Clo1Prop),matDim,matDim) 

#For convenience, I'm also separating things that are strictly reproduction (either sexual or asexual):
Reprod=Sex+Clo1

#Retrieves info about the broad life stage categories
lifeStages=DataSP[,"classOrganized"]

 #Add life stage descriptors to all matrices: 
  colnames(Sur)=colnames(Sex)=colnames(Clo1)=colnames(A)=colnames(Reprod)=lifeStages
  rownames(Sur)=rownames(Sex)=rownames(Clo1)=rownames(A)=rownames(Reprod)=c(1:matDim)
  
 #Extracting eigenstructure of matrix A
lambda=eigen.analysis(A)$lambda1 
Sens=eigen.analysis(A)$sensitivities 
Elas=eigen.analysis(A)$elasticities 

#Stage specific survival
  surStage=colSums(Sur)
  
  #Survival-independent A matrix
  for (i in 1:matDim) {
    for (j in 1:matDim) {
      ASurIndependent[j,i]=A[j,i]/surStage[i]
    }
  }

#Sensitivity matrix of survival
  SurSens=ASurIndependent*Sens
  SurSensStage=colSums(SurSens)
  SurElas=SurSens*surStage/lambda
  SurElasStage=SurSensStage*surStage/lambda
  
  #Sensitivity of growth, retrogression and reproduction:
  for (i in 1:matDim) {
    for (j in 1:matDim) {
      if (A[j,i]==0) (GRSCSens[j,i]=0)
      if (Sur[j,i]>0) (GRSCSens[j,i]=Sens[i,i]*(-surStage[i])+Sens[j,i]*(surStage[i]))
      if (Reprod[j,i]>0) (GRSCSens[j,i]=Sens[j,i]*(surStage[i]))
    }
  }


GRSCElas=GRSCSens*ASurIndependent/lambda
  
  GRSens=SurProp*GRSCSens
  SexSens=SexProp*GRSCSens
  Clo1Sens=Clo1Prop*GRSCSens
    
  GRElas=SurProp*GRSCElas
  SexElas=SexProp*GRSCElas
  Clo1Elas=Clo1Prop*GRSCElas
  
#Getting rid of NaNs and Infs
  GRSens=as.character(GRSens)
  SexSens=as.character(SexSens)
  Clo1Sens=as.character(Clo1Sens)
  
  GRElas=as.character(GRElas)
  SexElas=as.character(SexElas)
  Clo1Elas=as.character(Clo1Elas)
  
 
for (i in 1:(matDim^2)) {
    if (GRSens[i]=="NaN" | GRSens[i]=="Inf") (GRSens[i]="0") else (GRSens[i]=GRSens[i])
    if (SexSens[i]=="NaN" | SexSens[i]=="Inf") (SexSens[i]="0") else (SexSens[i]=SexSens[i])
    if (Clo1Sens[i]=="NaN" | Clo1Sens[i]=="Inf") (Clo1Sens[i]="0") else (Clo1Sens[i]=Clo1Sens[i])
       
    if (GRElas[i]=="NaN" | GRElas[i]=="Inf") (GRElas[i]="0") else (GRElas[i]=GRElas[i])
    if (SexElas[i]=="NaN" | SexElas[i]=="Inf") (SexElas[i]="0") else (SexElas[i]=SexElas[i])
    if (Clo1Elas[i]=="NaN" | Clo1Elas[i]=="Inf") (Clo1Elas[i]="0") else (Clo1Elas[i]=Clo1Elas[i])
}
  


GRSens=matrix(as.numeric(GRSens),matDim,matDim)
  diag(GRSens)=0
  SexSens=matrix(as.numeric(SexSens),matDim,matDim)
  Clo1Sens=matrix(as.numeric(Clo1Sens),matDim,matDim)
    
  VRSens=cbind(SurSens,GRSens,SexSens,Clo1Sens)
  rownames(VRSens)=lifeStages
  colnames(VRSens)=c(rep("Sur",matDim),rep("GrowShrink",matDim),rep("Sex",matDim),rep("Clo1",matDim))
  
  GRElas=matrix(as.numeric(GRElas),matDim,matDim)
  diag(GRElas)=0
  SexElas=matrix(as.numeric(SexElas),matDim,matDim)
  Clo1Elas=matrix(as.numeric(Clo1Elas),matDim,matDim)
    
  VRElas=cbind(SurElas,GRElas,SexElas,Clo1Elas)
  rownames(VRElas)=lifeStages
  colnames(VRElas)=c(rep("Sur",matDim),rep("GrowShrink",matDim),rep("Sex",matDim),rep("Clo1",matDim))
  

#dput(write.csv(VRSens, paste(species_name, "- VitalRateSens.csv")), sens)
#dput(write.csv(VRElas, paste(species_name, "- VitalRateElas.csv")), elas) 
write.csv(VRSens, sens)
write.csv(VRElas, elas)

R Server