Ternary plots
Rob Salguero-Gomez and Owen Jones
2021-05-05
Source:vignettes/a05_TernaryPlots.Rmd
a05_TernaryPlots.RmdHere we produce a ternary plot a la Silvertown & Franco (1993) with population growth rate as the “fourth” dimension. We will use functions from the popdemo, Rage and ggtern packages, so first we need to load those packages:
Next, we will load the dataset directly from the COMPADRE online repository. Then, we will subset making sure that we only select matrix population models (MPMs, hereafter) fulfilling two conditions. The first one is that the components of the MPM A (e.g. A = U + F + C; U: survival-dependent matrix; F: sexual reproduction matrix; C: clonal reproduction matrix) have already been split so we can automatically calculate vital-rate dependent processes (e.g. elasticities of population growth rate (*$). The second condition is that the MPM A must be non-ergodic (See Caswell 2001 for further details), so there is a unique solution for the dominant eigenvalue of A.
data(Compadre)
Compadre <- cdb_flag(Compadre)
Compadre <- subset(Compadre, MatrixSplit == "Divided" & check_ergodic == TRUE)We can calculate elasticities of population growth rate to changes in the matrix elements using the elas function from the popdemo package. For example, the elasticity matrix for the first A matrix in this example database looks like this:
popdemo::elas(matA(Compadre)[[1]])
#> A1 A2 A3 A4
#> A1 0.1645631 0.06524073 0.02430743 0.03131899
#> A2 0.1208672 0.16721803 0.04892799 0.01682975
#> A3 0.0000000 0.07802604 0.11046976 0.02794895
#> A4 0.0000000 0.04335812 0.03273957 0.06818440Because the subset of MPMs here has split A into the matrices U, F, and C, we can now also automatise the calculation of elasticities, regardless of the varying degree of life cycle complexities represented in these MPMs. Therefore, we can classify the elements of the matrix to changes in stages (which we can further identify as stasis, progression, and retrogression), sexual reproduction, and clonal reproduction.
We can use the matrixElementPerturbation function from the Rage package to conduct this element-by-element elasticity analysis. The function outputs both sensitivities (prefixed with S) and elasticities (prefixed with E) as follows:
mat_a <- matA(Compadre)[[1]]
perturb_matrix(mat_a, type = "elasticity")
#> A1 A2 A3 A4
#> A1 0.1645633 0.06524073 0.02430742 0.03131899
#> A2 0.1208670 0.16721815 0.04892799 0.01682975
#> A3 0.0000000 0.07802598 0.11046990 0.02794895
#> A4 0.0000000 0.04335804 0.03273956 0.06818449
#Rage::matrixElementPerturbation(matU = matU(Compadre)[[1]], This function used to be able to accommodate partioning
# matF = matF(Compadre)[[1]], into U F and C so we coudl partition this... but it's now gone
# matC = matC(Compadre)[[1]]) in perturb_matrix.... why?We can use a for loop to run through each set of matrices in turn to calculate the summed elasticity for survival (S), growth (G) and reproduction (R) like this:
Amats <- matA(Compadre)
Umats <- matU(Compadre)
Fmats <- matF(Compadre)
Cmats <- matC(Compadre)
output <- data.frame(S=rep(NA,length(Umats)),G=NA,R=NA,lam=NA)
for(i in 1:length(Umats)){
temp <- perturb_vr(Umats[[i]], #The same issue as above
Fmats[[i]],
Cmats[[i]],
type = "elasticity")
output$S[i] <- temp$survival
output$G[i] <- temp$growth + temp$shrinkage
output$R[i] <- temp$fecundity + temp$clonality
#Calculate growth rate
output$lam[i] <- popdemo::eigs(Amats[[i]], "lambda")
}
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zeroLet’s take a look at this output:
head(output)
#> S G R lam
#> 1 1.0000003 0.077203676 0.137696890 1.0580286
#> 2 0.9999160 0.192840443 0.154347265 1.4118127
#> 3 1.0000008 0.001738428 0.007431082 0.9962044
#> 4 0.9999998 0.194812694 0.284476663 1.1241625
#> 5 1.0000001 -0.000474047 0.003098288 0.9659322
#> 6 0.9999999 -0.351429957 0.000000000 0.9154951Now we have elasticities for our three demographic processes we can place these onto a ternary plot. However, we should first scale the elasticities so that they sum to 1 - this is necessary because of possible rounding issues:
Now for the plot:
B<-ggtern::ggtern(data = output,
aes(x = R,
y = G,
z = S,
colour = lam)) +
geom_point() +
scale_color_viridis_c()+
theme_showarrows() +
theme_clockwise()
B
Now you can try this, try using another variable, such as reactivity, or life expectancy as the “fourth” dimension instead of lambda.
IN THE TEXT BELOW I RE-DO THE FUNCTIONALITY ABOVE AT THE VITAL RATE LEVEL. WE SHOULD FIRST FIX THE FUNCTION ABOVE SO IT DOES THE PARTITIONING OF MATRICES INTO SUBCOMPONENTS U F AND C
We can use the perturb_vr function from the Rage package to conduct the perturbations at the level of underlying vital rates rather than matrix elements. It is worth noting that in a MPM, most matrix elements are composites of vital rates. For instance, the element a[2,1] (second row, first column in the matrix A) is composed by the survival probability of the first stage times the probability that individuals in the first stage transition to the second one in the life cycle of the species. To proceed with the calculation of the vital rate perturbations:
data(mpm1)
perturb_vr(mpm1$matU,
mpm1$matF,
mpm1$matF,
type = "elasticity")
#> $survival
#> [1] 0.9999983
#>
#> $growth
#> [1] 0.3843816
#>
#> $shrinkage
#> [1] -0.02546099
#>
#> $fecundity
#> [1] 0.2227753
#>
#> $clonality
#> [1] 0.2227753We can use a for loop to run through each set of matrices in turn to calculate the summed elasticity for (i) survival (\(\sigma), (ii) growth (\)) and shrinkage (\(\rho), and (iii) sexual (\)) and clonal reproduction ($) like this:
Amats <- matA(Compadre)
Umats <- matU(Compadre)
Fmats <- matF(Compadre)
Cmats <- matC(Compadre)
output <- data.frame(survival=rep(NA,length(Umats)),growthShrinkage=NA,reproduction=NA,lam=NA)
for(i in 1:length(Umats)){
temp <- perturb_vr(Umats[[i]],
Fmats[[i]],
Cmats[[i]],
type = "elasticity")
output$survival[i] <- temp$survival
output$growthShrinkage[i] <- temp$growth + temp$shrinkage
output$reproduction[i] <- temp$fecundity + temp$clonality
#Calculate growth rate
output$lam[i] <- popdemo::eigs(Amats[[i]], "lambda")
}
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zero
#> Warning: All elements of matF are zeroLet’s take a look at this output:
head(output)
#> survival growthShrinkage reproduction lam
#> 1 1.0000003 0.077203676 0.137696890 1.0580286
#> 2 0.9999160 0.192840443 0.154347265 1.4118127
#> 3 1.0000008 0.001738428 0.007431082 0.9962044
#> 4 0.9999998 0.194812694 0.284476663 1.1241625
#> 5 1.0000001 -0.000474047 0.003098288 0.9659322
#> 6 0.9999999 -0.351429957 0.000000000 0.9154951Now we have elasticities for our three collapsed vital rate processes we can place these onto a ternary plot. However, we should first scale the elasticities so that they sum to 1 - this is necessary because of possible rounding issues:
Now for the plot:
B<-ggtern::ggtern(data = output,
aes(x = reproduction,
y = growthShrinkage,
z = survival,
colour = lam)) +
geom_point() +
scale_color_viridis_c()+
theme_showarrows() +
theme_clockwise()
B