Allow collapseMatrix to collapse non-adjacent stages & remove rearrangeMatrix

[patrickbarks]

collapseMatrix is currently written in a way that prevents a user from collapsing non-adjacent stages. Specifically, the collapse argument is a character vector whose elements must reflect single stages or sequences of adjacent stages, e.g.

collapse <- c('1', '2-4', '5')

I suggest we change collapse to a list, to allow the user to collapse non-adjacent stages, e.g.

collapse <- list(1, 2:4, 5)           # same as above
collapse <- list(1, c(2, 4), c(3, 5)) # collapse non-adjacent stages (e.g. stage 3 inter-reprod.)

I think this would make rearrangeMatrix unnecessary, because it's only function is to move any inter-reproductive stages to the right side of a matrix, to segregate reproductive and post/non-reproductive stages ahead of collapsing.

If we remove rearrangeMatrix, reprodStages can be updated to identify inter-reproductive stages (it already takes the required arguments, i.e. reproStages).

[RobSalGo]

Hi @patrickbarks,
I agree with keeping the functionality within collapseMatrix to let the author collapse non-contiguous classes, but I also note that letting the authors re-arrange the matrices is a valuable thing per se - some times authors do not follow an ever-increasing level of comlplexity in stages from left to right, and that leads to the probabilities of retrogression being incorrectly calculated unless re-arranged

[patrickbarks]

That makes sense, I'll leave rearrangeMatrix.

I think there are a few MPMs where inter-reproductive stages are legitimate, and appropriately ordered by the author (e.g. doi.org/10.2307/3061069), in which case moving them to the right inappropriately swaps growth and retrogression. Perhaps we can just add a warning to the documentation to this effect.

(Side note: I think there are other scenarios where the distinction between growth/retrogression is tricky, such as when an MPM includes a dormant stage, or when the MPM incorporates additional states such as sex, management status, disease status, etc. Still thinking about how to deal with some these — particularly dormancy.)

[RobSalGo]

That's a good point, but luckily a rather easy one to take care of. See what I've done here, where dorm is:

dorm <- which(compadre$matrixClass[[i]]$MatrixClassOrganized=="dorm")

one could do the same thing about seedbank dynamics as I do with going into, staying dormant and awakening from dormancy (below), but going into, staying in seedbank, and germinating from seedbank w something like:

propdorm <- which(compadre$matrixClass[[i]]$MatrixClassOrganized=="prop")

vitalRate <- function(matU, matF, matC, dorm){

matU[is.na(matU)]=0
matF[is.na(matF)]=0
matC[is.na(matC)]=0

matA=matU+matF+matC
matDim=dim(matA)[1]
dorm=dorm

out = data.frame("Survival"=NA,"Progression"=NA,"Retrogression"=NA,"Reproduction"=NA,"Clonality"=NA,"GoDorm"=NA,"StayDorm"=NA,"Awake"=NA,
                 "SurvivalSSD"=NA,"ProgressionSSD"=NA,"RetrogressionSSD"=NA,"ReproductionSSD"=NA,"ClonalitySSD"=NA,"GoDormSSD"=NA,"StayDormSSD"=NA,"AwakeSSD"=NA)

#Extracting SSD corrected vital rate values
SSD=eigen.analysis(matA)$stable.stage
f=colSums(matF)
out$Reproduction=mean(f)
out$ReproductionSSD=mean(f*SSD)
c=colSums(matC)
out$Clonality=mean(c)
out$ClonalitySSD=mean(c*SSD)

#Preparing survival-independent matrix to calculate survival and dormancy
uDistrib=matrix(NA,ncol=matDim,nrow=matDim)
u=colSums(matU)
out$Survival=mean(u[which(!(1:matDim)%in%dorm)])
out$SurvivalSSD=mean((u[which(!(1:matDim)%in%dorm)])*SSD[which(!(1:matDim)%in%dorm)])
if (length(dorm)>0) {
  out$StayDorm=mean(u[dorm])
  out$StayDormSSD=mean((u[dorm])*(SSD[dorm]))
}

#Making matrix for transitions conditional on survival
for (j in which(u>0)) uDistrib[,j]=matU[,j]/u[j]
UPrime=uDistrib
UPrime[is.na(UPrime)]=0

#Extracting proxy to individual progressive growth rate
UPrimeGrowth=UPrime
UPrimeGrowth[upper.tri(UPrime, diag = T)]=NA
UPrimeGrowth[matDim,matDim]=UPrime[matDim,matDim]  #Putting back the last element of stasis bc there is likely growth on the top of class
out$Progression=mean(colSums(UPrimeGrowth,na.rm=T)[which(!(1:matDim)%in%dorm)])
out$ProgressionSSD=mean(colSums(UPrimeGrowth,na.rm=T)[which(!(1:matDim)%in%dorm)]*(SSD[which(!(1:matDim)%in%dorm)]))

#Extracting proxy to individual retrogressive growth rate
UPrimeShrinkage=UPrime
UPrimeShrinkage[lower.tri(UPrime, diag = T)]=NA
out$Retrogression=mean(colSums(UPrimeShrinkage,na.rm=T)[which(!(1:matDim)%in%dorm)])
out$RetrogressionSSD=mean((colSums(UPrimeShrinkage,na.rm=T)[which(!(1:matDim)%in%dorm)])*(SSD[which(!(1:matDim)%in%dorm)]))

if (length(dorm)>0) {
  
  #Extracting proxy to going into dormancy
  UPrimeGoDorm=UPrime
  UPrimeGoDorm[upper.tri(UPrime, diag = T)]=NA
  UPrimeGoDorm[matDim,matDim]=UPrime[matDim,matDim]  #Putting back the last element of stasis bc there is likely growth on the top of class
  out$GoDorm=mean(colSums(UPrimeGrowth,na.rm=T)[dorm])
  out$GoDormSSD=mean(colSums(UPrimeGrowth,na.rm=T)[dorm]*(SSD[dorm]))
  
  #Extracting proxy to awaking from dormancy
  UPrimeAwake=UPrime
  UPrimeAwake[lower.tri(UPrime, diag = T)]=NA
  out$Awake=mean(colSums(UPrimeAwake,na.rm=T)[dorm])
  out$AwakeSSD=mean((colSums(UPrimeAwake,na.rm=T)[dorm])*(SSD[dorm]))
}  

return(out)  

}

[patrickbarks]

Well alright then! Now the challenge is to incorporate this idea into our functions (potentially including vitalRates, matrixElementPerturbation, vitalRatePerturbation, reprodStages, etc.).

[patrickbarks]

In addition to converting the 'collapse' argument of collapseMatrix to a list, I suggest we correspondingly convert the output of reprodStages to a list, e.g.:

list(propStages = propStage,
     preRepStages = preRep,
     repStages = Rep,
     postRepStages = postRep)

This would allow the output of reprodStages to be directly input into collapseMatrix.

(This idea is already suggested within reprodStages as a #FIXME note attributed to Rob)