data(lathyrus)dim(lathyrus)#> [1] 1119   38#summary(lathyrus)

Step 1. Life history model development

To begin, we will create astageframe. A stageframe isa data frame that describes all stages in the life history of theorganism, in a way usable by the functions in this package and usingstage names and classifications that completely match those used in thedataset. It must include complete descriptions of all stages that occurin the dataset, with each stage defined uniquely. Since this object canbe used for automated classification of individuals, all sizes,reproductive states, and other characteristics defining each stage inthe dataset need to be accounted for explicitly. This can be difficultif a few data points do not fit neatly into any stage description, sogreat care must be taken to include all relevant size values and valuesof other descriptor variables occurring within the dataset. The finaldescription of each stage occurring in the dataset must not completelyoverlap with any other stage, although partial overlap is allowed. Wewill base our stageframe on the life history model provided inEhrlén (2000), but use a different sizeclassification based on leaf volume to allow IPM construction and makeall mature stages other than vegetative dormancy reproductive, as shownin Figure 6.1.

Figure 6.1. Life history model ofLathyrus vernus. Not alladult classes are shown. Survival transitions are indicated with solidarrows, while fecundity transitions are indicated with dashedarrows.

In the stageframe code below, we show that we want an IPM by choosingtwo stages that serve as the size limits for IPM size classification.These two size classes should have exactly the same characteristics inthe stageframeother than size. By choosing these twosize limits, we can skip adding and describing the many size classesthat will fall between these limits - functionsf_create()will create all of these for us. We mark these limits in the vector thatwe load into thestagenames option using the lower-casestringipm. Packagelefko3 will then createand name all IPM size classes according to its own conventions. Thedefault number of size classes is 100 bins, which can be changed usingtheipmbins option.

sizevector<-c(0,100,0,1,7100)stagevector<-c("Sd","Sdl","Dorm","ipm","ipm")repvector<-c(0,0,0,1,1)obsvector<-c(0,1,0,1,1)matvector<-c(0,0,1,1,1)immvector<-c(1,1,0,0,0)propvector<-c(1,0,0,0,0)indataset<-c(0,1,1,1,1)binvec<-c(0,100,0.5,1,1)comments<-c("Dormant seed","Seedling","Dormant","ipm adult stage","ipm adult stage")lathframeipm<-sf_create(sizes = sizevector,stagenames = stagevector,repstatus = repvector,obsstatus = obsvector,propstatus = propvector,immstatus = immvector,matstatus = matvector,comments = comments,indataset = indataset,binhalfwidth = binvec,ipmbins =100,roundsize =3)#> Warning: Values supplied in vectors sizemin, sizemax, and binhalfwidth will be overwritten for ipm bins.#lathframeipm

This stageframe has 103 stages. The IPM portion technically starts withthe fourth stage and keeps going through the 103rd stage. Stage nameswithin this range are concatenations of the size centroid (designatedwithsz), and, given a limited string length, thereproductive status, maturity status, and observation status. The firstthree stages, which fall outside of the IPM classification, are leftunaltered. Note that we are using the midpoint approach to determiningthe size bins here, using the default bin halfwidths of 0.5 (since allsize bins are using this default, we do not include abinhalfwidth vector option in thesf_create()input). However, we could have used thesizemin andsizemax options to more deliberately set the size binminima and maxima instead, although that would have been a more tediousapproach in a high dimensonal object like an IPM.

Step 2a. Dataset standardization

To work with this dataset, we first need to format the data intovertical format, in which each row corresponds to the state ofa single individual in two (if ahistorical) or three (if historical)consecutive time intervals. In the input toverticalize3()below, we utilize a repeating pattern of variable names arranged in thesame order for each monitoring occasion. This arrangement allows us toenter only the first variable in each set, as long asnoyears andblocksize are set properly and thecolumns have no gaps or shuffles. The data management functions inlefko3 do not require such repeating patterns, but they domake the required input in the function shorter and more succinct.Because this is an IPM, we will need to estimate linear models of vitalrates. This will require us to avoid NAs in size and fecundity(fortunately, in this dataset, NA is equivalent to 0), so we will setNAas0 = TRUE to inform R to treat NAs as zeroes. We willalso setNRasRep = TRUE because we will assume that alladult stages other than dormancy are reproductive, and there are matureindividuals in the dataset that do not reproduce but need to be includedin reproductive stages.

Note that prior to standardizing the dataset, we will create a newvariable to code individual identity, since different plants indifferent subpopulations use the same identifiers. To see a moredetailed summary of the standardized dataset, please addfull = TRUE to thesummary_hfv() call.

lathyrus$indiv_id<-paste(lathyrus$SUBPLOT, lathyrus$GENET)lathvertipm<-verticalize3(lathyrus,noyears =4,firstyear =1988,individcol ="indiv_id",blocksize =9,juvcol ="Seedling1988",sizeacol ="Volume88",repstracol ="FCODE88",fecacol ="Intactseed88",deadacol ="Dead1988",nonobsacol ="Dormant1988",stageassign = lathframeipm,stagesize ="sizea",censorcol ="Missing1988",censorkeep =NA,censorRepeat =TRUE,censor =TRUE,NAas0 =TRUE,NRasRep =TRUE)summary_hfv(lathvertipm)#>#> This hfv dataset contains 2527 rows, 54 variables, 1 population,#> 1 patch, 1053 individuals, and 3 time steps.

Before we move on to the next steps in analysis, let’s take a closerlook at fecundity. In this dataset, fecundity is mostly a count ofintact seeds, and only differs in six cases where the seed output wasestimated based on other models. To see this, try the following code,which focuses on fecundity in timet.

# Length of fecundity variable in t:length(lathvertipm$feca2)#> [1] 2527# Number of non-integer entries:length(which(lathvertipm$feca2!=round(lathvertipm$feca2)))#> [1] 6

We see that we have quite a bit of fecundity data, and that it isoverwhelmingly but not exclusively integer. So, we can either treatfecundity as a continuous variable, or round the values and treatfecundity as a count variable. We will choose the latter approach inthis analysis.

lathvertipm$feca3<-round(lathvertipm$feca3)lathvertipm$feca2<-round(lathvertipm$feca2)lathvertipm$feca1<-round(lathvertipm$feca1)

Although we wish to treat fecundity as a count, it is still not clearwhat underlying distribution we should use. The Poisson distributionassumes that the mean and variance are equal, and so we can test thisassumption using a chi-squared test. If it is not significantlydifferent, then we may use some variant of the Poisson distribution. Ifthe data are significantly over-dispersed, then we should use thenegative binomial distribution. If fecundity of 0 is possible inreproductive stages, as in cases where reproductive status is defined byflowering rather than by offspring production, then we should also testwhether the number of zeroes is significantly greater than expectedunder these distributions, and use a zero-inflated distribution if so.

Let’s formally test our assumptions. We will use thehfv_qc() function, which will allow us to assess thequality of the data from the viewpoint of the linear models to be builtlater. We will use most of the input from themodelsearch()call that we will conduct later.

hfv_qc(data = lathvertipm,vitalrates =c("surv","obs","size","fec"),juvestimate ="Sdl",indiv ="individ",year ="year2",age ="obsage")#> Survival:#>#>   Data subset has 54 variables and 2246 transitions.#>#>   Variable alive3 has 0 missing values.#>   Variable alive3 is a binomial variable.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 931   (singleton categories: 181)#>   year2: 3   (singleton categories: 0)#>#> Observation status:#>#>   Data subset has 54 variables and 2121 transitions.#>#>   Variable obsstatus3 has 0 missing values.#>   Variable obsstatus3 is a binomial variable.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 858   (singleton categories: 143)#>   year2: 3   (singleton categories: 0)#>#> Primary size:#>#>   Data subset has 54 variables and 1916 transitions.#>#>   Variable sizea3 has 0 missing values.#>   Variable sizea3 appears to be a floating point variable.#>   1256 elements are not integers.#>   The minimum value of sizea3 is 3.4 and the maximum is 6646.#>   The mean value of sizea3 is 512.8 and the variance is 507200.#>   The value of the Shapiro-Wilk test of normality is 0.7134 with P = 3.014e-49.#>   Variable sizea3 differs significantly from a Gaussian distribution.#>#>   Variable sizea3 is fully positive, lacking even 0s.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 845   (singleton categories: 185)#>   year2: 3   (singleton categories: 0)#>#> Fecundity:#>#>   Data subset has 54 variables and 2246 transitions.#>#>   Variable feca2 has 0 missing values.#>   Variable feca2 appears to be an integer variable.#>#>   Variable feca2 is fully non-negative.#>#>   Overdispersion test:#>     Mean feca2 is 1.282#>     The variance in feca2 is 23.21#>     The probability of this dispersion level by chance assuming that#>     the true mean feca2 = variance in feca2,#>     and an alternative hypothesis of overdispersion, is 0#>     Variable feca2 is significantly overdispersed.#>#>   Zero-inflation and truncation tests:#>     Mean lambda in feca2 is 0.2774#>     The actual number of 0s in feca2 is 1980#>     The expected number of 0s in feca2 under the null hypothesis is 623#>     The probability of this deviation in 0s from expectation by chance is 0#>     Variable feca2 is significantly zero-inflated.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 931   (singleton categories: 181)#>   year2: 3   (singleton categories: 0)#>#> Juvenile survival:#>#>   Data subset has 54 variables and 281 transitions.#>#>   Variable alive3 has 0 missing values.#>   Variable alive3 is a binomial variable.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 281   (singleton categories: 281)#>   year2: 3   (singleton categories: 0)#>#> Juvenile observation status:#>#>   Data subset has 54 variables and 210 transitions.#>#>   Variable obsstatus3 has 0 missing values.#>   Variable obsstatus3 is a binomial variable.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 210   (singleton categories: 210)#>   year2: 3   (singleton categories: 0)#>#> Juvenile primary size:#>#>   Data subset has 54 variables and 193 transitions.#>#>   Variable sizea3 has 0 missing values.#>   Variable sizea3 appears to be a floating point variable.#>   127 elements are not integers.#>   The minimum value of sizea3 is 2.1 and the maximum is 61.#>   The mean value of sizea3 is 11.23 and the variance is 50.81.#>   The value of the Shapiro-Wilk test of normality is 0.5997 with P = 5.72e-21.#>   Variable sizea3 differs significantly from a Gaussian distribution.#>#>   Variable sizea3 is fully positive, lacking even 0s.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 193   (singleton categories: 193)#>   year2: 3   (singleton categories: 0)#>#> Juvenile maturity status:#>#>   Data subset has 54 variables and 210 transitions.#>#>   Variable matstatus3 has 0 missing values.#>   Variable matstatus3 is a binomial variable.#>#>   Numbers of categories in data subset in possible random variables:#>   indiv id: 210   (singleton categories: 210)#>   year2: 3   (singleton categories: 0)

All of the probability variables look like binomials, so we have noproblem there. Size and juvenile size appear to be continuous variablesthat differ significantly from the assumptions of the Gaussiandistribution. However, we will still use a Gaussian distribution in thiscase, for instructional purposes. Fecundity is a count variable withsignificant overdispersion and significantly more zeros than expected,so we will use a zero-inflated negative binomial distribution. Note thatthere might be some problems in modeling in the juvenile cases, where wesee that individual identity has as many levels as data points.

Step 2b: Develop supplemental information for matrix estimation

Now we will createsupplement tables, which provideextra data for matrix estimation that is not included in the maindemographic dataset. Specifically, we will provide the seed dormancyprobability and germination rate, which are given as transitions fromthe dormant seed stage to another year of seed dormancy or to thegerminated seedling stage, respectively. We assume that the germinationrate is the same regardless of whether seed was produced in the previousyear or has been in the seedbank for longer. We will incorporate theseterms both as fixed constants for specific transitions within theresulting matrices, and as multipliers for fecundity, since ultimatelyfecundity will be estimated as the production of seed multiplied by theseed germination rate or the seed dormancy rate. The fecunditymultipliers will also serve to tell R which transitions are thefecundity transitions. Because some individuals stay in the seedlingstage for only 1 year, and the seed stage itself cannot be observed andso does not exist in the dataset, we will also set a proxy set oftransitions assuming that transitions from seed in occasiont-1to seedling in occasiont to all mature stages in occasiont+1 are equal to the equivalent transitions from seedling inboth occasionst-1 andt. We will start with theahistorical case, and then move on to the historical case, where we alsoneed to input the corresponding stages in occasiont-1 andtransition types from occasiont-1 tot for eachtransition.

lathsupp2<-supplemental(stage3 =c("Sd","Sdl","Sd","Sdl"),stage2 =c("Sd","Sd","rep","rep"),givenrate =c(0.345,0.054,NA,NA),multiplier =c(NA,NA,0.345,0.054),type =c(1,1,3,3),stageframe = lathframeipm,historical =FALSE)#> Warning: NA values in argument multiplier will be treated as 1 values.lathsupp3<-supplemental(stage3 =c("Sd","Sd","Sdl","Sdl","npr","Sd","Sdl"),stage2 =c("Sd","Sd","Sd","Sd","Sdl","rep","rep"),stage1 =c("Sd","rep","Sd","rep","Sd","mat","mat"),eststage3 =c(NA,NA,NA,NA,"npr",NA,NA),eststage2 =c(NA,NA,NA,NA,"Sdl",NA,NA),eststage1 =c(NA,NA,NA,NA,"Sdl",NA,NA),givenrate =c(0.345,0.345,0.054,0.054,NA,NA,NA),multiplier =c(NA,NA,NA,NA,NA,0.345,0.054),type =c(1,1,1,1,1,3,3),type_t12 =c(1,2,1,2,1,1,1),stageframe = lathframeipm,historical =TRUE)#> Warning: NA values in argument multiplier will be treated as 1 values.#lathsupp2#lathsupp3

Step 3. Tests of history, and vital rate modeling

Integral projection models (IPMs) require functions of vital rates topopulate them. Here, we will develop these functions as linear modelsusingmodelsearch(). This function automates severalcrucial and complex tasks in MPM analysis. Specifically, it automates 1)the building of global models for each vital rate requested, 2) theexhaustive construction of all reduced models, and 3) the selection ofthe best-fit models. This function also allows us to test whetherindividual history affects demography. Settinghistorical = TRUE fits size and/or reproductive status inoccasionst andt-1 into global models, while settinghistorical = FALSE limits testing to the impacts of statein occasiont only. The model building and selection protocolscan then be used to see if history has a significant impact on a vitalrate, by assessing whether a historical term is retained in the best-fitmodel.

First, we will create the historical models to assess whether history isa significant influence on vital rates.

lathmodels3ipm<-modelsearch(lathvertipm,historical =TRUE,approach ="mixed",suite ="main",vitalrates =c("surv","obs","size","fec"),juvestimate ="Sdl",bestfit ="AICc&k",sizedist ="gaussian",fecdist ="negbin",fec.zero =TRUE,indiv ="individ",year ="year2",year.as.random =TRUE,juvsize =TRUE,quiet ="partial")#>#> Developing global model of survival probability...#>#> Global model of survival probability developed. Proceeding with model dredge...#>#> Developing global model of observation probability...#>#> Global model of observation probability developed. Proceeding with model dredge...#>#> Developing global model of primary size...#>#> Global model of primary size developed. Proceeding with model dredge...#>#> Developing global model of fecundity...#>#> Global model of fecundity developed. Proceeding with model dredge...#>#> Developing global model of juvenile survival probability...#>#> Global model of juvenile survival probability developed. Proceeding with model dredge...#> Warning: Juvenile maturity status in time t+1 appears to be constant (1). Setting to constant.#>#> Developing global model of juvenile observation probability...#>#> Global model of juvenile observation probability developed. Proceeding with model dredge...#>#> Developing global model of juvenile primary size...#>#> Global model estimation failed. Dropping individual identity term.#>#> Redeveloping global model of juvenile primary size...#>#> Global model of juvenile primary size developed. Proceeding with model dredge...#>#> Finished selecting best-fit models.#summary(lathmodels3ipm)

We see here, as before, that status in occasiont-1 exerts aninfluence on some vital rates, particularly survival to occasiont+1, size in occasiont+1, and the conditional portionof fecundity. So, the historical IPM is the correct choice here.Accuracy is also quite high for adult survival and juvenile and adultobservation, but quite poor for primary size and fecundity.

We will also create an ahistorical IPM for comparison. For that purpose,we will create the ahistorical linear model set.

lathmodels2ipm<-modelsearch(lathvertipm,historical =FALSE,approach ="mixed",suite ="main",vitalrates =c("surv","obs","size","fec"),juvestimate ="Sdl",bestfit ="AICc&k",sizedist ="gaussian",fecdist ="negbin",fec.zero =TRUE,indiv ="individ",year ="year2",year.as.random =TRUE,juvsize =TRUE,quiet ="partial")#>#> Developing global model of survival probability...#>#> Global model of survival probability developed. Proceeding with model dredge...#>#> Developing global model of observation probability...#>#> Global model of observation probability developed. Proceeding with model dredge...#>#> Developing global model of primary size...#>#> Global model of primary size developed. Proceeding with model dredge...#>#> Developing global model of fecundity...#>#> Global model of fecundity developed. Proceeding with model dredge...#>#> Developing global model of juvenile survival probability...#>#> Global model of juvenile survival probability developed. Proceeding with model dredge...#> Warning: Juvenile maturity status in time t+1 appears to be constant (1). Setting to constant.#>#> Developing global model of juvenile observation probability...#>#> Global model of juvenile observation probability developed. Proceeding with model dredge...#>#> Developing global model of juvenile primary size...#>#> Global model estimation failed. Dropping individual identity term.#>#> Redeveloping global model of juvenile primary size...#>#> Global model of juvenile primary size developed. Proceeding with model dredge...#>#> Finished selecting best-fit models.#summary(lathmodels2ipm)

Step 4. IPM estimation

We will now create the historical suite of matrices covering the yearsof study. These matrices will be extremely large, but they are sparseand we will build them in sparse format, so they should not take toomuch memory.

lathmat3ipm<-flefko3(stageframe = lathframeipm,modelsuite = lathmodels3ipm,supplement = lathsupp3,data = lathvertipm,reduce =FALSE,sparse_output =TRUE)#summary(lathmat3ipm)

These are giant matrices. With 10,609 rows and columns, there are atotal of 112,550,881 elements per matrix. But they are also amazinglysparse - with approximately 960,209 elements estimated per matrix, only0.9% of elements per matrix are non-zero. The survival probability sumsall look good, so we appear to have no problems with overly large givenand proxy survival transitions provided through our supplemental tables.Let’s now build the ahistorical IPMs.

lathmat2ipm<-flefko2(stageframe = lathframeipm,modelsuite = lathmodels2ipm,supplement = lathsupp2,data = lathvertipm,reduce =FALSE)#summary(lathmat2ipm)

The ahistorical IPMs are certainly smaller than the historical IPMs, butare nonetheless quiet large. Although huge, these matrices are notsparse - an average of 9,180.3 elements out of 10,609 per matrix areestimated (86.5%; note that the summary() function only counts elementsas estimated if they equal a value other than 0, leading the exactnumber of estimated elements to vary among matrices when elements areestimated at near 0).

Let’s do a further comparison - let’s view a matrix plot of each kind ofMPM, ahistorical and then historical. First, the ahistorical plot.

image3(lathmat2ipm,used =1)#> [[1]]

Figure 6.2. Image of ahistorical projectionmatrix

Now the historical plot.

image3(lathmat3ipm,used =1)#> [[1]]

Figure 6.3. Image of historical projectionmatrix

The plots above show that the ahistorical matrix is large but dense,mostly full of non-zero entries (the colored elements correspond tonon-zero elements). In contrast, the historical matrix is huge andsparse, mostly full of zeroes with a general pattern to the distributionof non-zero elements.

Now let’s estimate the mean IPM matrices. We will estimate one meanmatrix each, because we did not separate patches in the datareorganization and vital rate modeling. As a check, let’s also look overthe summaries.

lath2ipmmean<-lmean(lathmat2ipm)lath3ipmmean<-lmean(lathmat3ipm)#summary(lath2ipmmean)#summary(lath3ipmmean)

All looks fine! Let’s also take a look at a portion of one of theconditional historical matrices, particularly the matrix conditional onvegetative dormancy in occasiont-1. This matrix can becompared to the ahistorical mean to assess the impacts of history on thematrix elements themselves. Please remember to remove the hashtag.

l3mcond<-cond_hmpm(lath3ipmmean)#l3mcond$Mcond[[1]]$Dorm

Step 5. MPM analysis

Now let’s estimate and plot the deterministic population growth rate,\(\lambda\), for each year. We findthat the ahistorical and historical estimates of\(\lambda\) are roughly in line with eachother. Please note we may also estimate standard errors for lambda, andmethods to do so are illustrated in other vignettes (see particularlythe vignette on Cypripedium candidum function-based MPMs). We will notdo so here because the historical IPM will be too large in bootstrappedform.

ipm2lambda<-lambda3(lathmat2ipm)ipm3lambda<-lambda3(lathmat3ipm)plot(lambda~ year2,data = ipm2lambda,xlab ="Year",ylab ="Lambda",ylim =c(0.65,1.00),type ="l",lwd =2,bty ="n")lines(lambda~ year2,data = ipm3lambda,lwd =2,lty =2,col ="red")legend("bottomleft",c("ahistorical","historical"),lty =c(1,2),col =c("black","red"),lwd =2,bty ="n")

Figure 6.4. Ahistorical vs. historicallambda

Let’s now look at\(\lambda\) for themean matrices and compare with the stochastic growth rate,\(a = \text{log} \lambda _{S}\). We will setthe number of simulations low in the historical case in order to keepthe amount of memory used and computation time low, because the size ofthe historical matrices will use up plenty of both (the default is10,000). We will also set the seed for the random number generator tomake our results reproducible.

lambda3(lath2ipmmean)#>   pop patch    lambda#> 1   1  <NA> 0.8858591lambda3(lath3ipmmean)#>   pop patch    lambda#> 1   1  <NA> 0.9529018set.seed(42)slambda3(lathmat2ipm)#>   pop patch          a        var        sd          se#> 1   1  <NA> -0.1252668 0.05120795 0.2262917 0.002262917set.seed(42)slambda3(lathmat3ipm,times =1000)#>   pop patch           a        var       sd         se#> 1   1  <NA> -0.05142838 0.04279148 0.206861 0.00654152

The historical growth rate is larger than the ahistorical in bothdeterministic and stochastic analyses, although all four numbers suggesta declining population. Now let’s compare the stable stage distributionfrom both the ahistorical and historical mean MPMs.

ipm2ss<-stablestage3(lath2ipmmean)ipm3ss<-stablestage3(lath3ipmmean)ipm2ss_s<-stablestage3(lathmat2ipm,stochastic =TRUE,seed =42)ipm3ss_s<-stablestage3(lathmat3ipm,stochastic =TRUE,times =1000,seed =42)ss_put_together<-cbind.data.frame(ipm2ss$ss_prop, ipm3ss$ahist$ss_prop,  ipm2ss_s$ss_prop, ipm3ss_s$ahist$ss_prop)names(ss_put_together)<-c("det ahist","det hist","sto ahist","sto hist")rownames(ss_put_together)<- ipm2ss$stage_idlty.o<-par("lty")par(lty =0)barplot(t(ss_put_together),beside=T,ylab ="Proportion",xlab ="Stage",col =c("black","orangered","grey","darkred"),bty ="n")legend("topright",c("det ahist","det hist","sto ahist","sto hist"),col =c("black","orangered","grey","darkred"),pch =15,bty ="n")

Figure 6.5. Ahistoricalvs. historically-corrected stable stage distribution

par(lty = lty.o)

Both ahistorical and historical approaches show the stable stagedistribution dominated by the dormant seed stage. The next highest isthe seedling stage, followed by the adult dormant stage. Stochasticanalysis shows similar patterns.

We will now move on to elasticity analysis. We will reduce the number ofoccasion steps simulated in the stochastic historical analysis to reducethe amount of computer processing time (this might take several hourseven with the number of steps reduced).

lath2ipmelas<-elasticity3(lath2ipmmean)#> Running deterministic analysis...lath3ipmelas<-elasticity3(lath3ipmmean)#> Running deterministic analysis...lath2ipmelas_s<-elasticity3(lathmat2ipm,stochastic =TRUE,seed =42)#> Running stochastic analysis...lath3ipmelas_s<-elasticity3(lathmat3ipm,stochastic =TRUE,times =200,sparse =TRUE,seed =42)#> Running stochastic analysis...elas_put_together<-cbind.data.frame(colSums(lath2ipmelas$ah_elasmats[[1]]),  Matrix::colSums(lath3ipmelas$ah_elasmats[[1]]),colSums(lath2ipmelas_s$ah_elasmats[[1]]),  Matrix::colSums(lath3ipmelas_s$ah_elasmats[[1]]))names(elas_put_together)<-c("det ahist","det hist","sto ahist","sto hist")rownames(elas_put_together)<- lath2ipmelas$ah_stages$stage_idlty.o<-par("lty")par(lty =0)barplot(t(elas_put_together),beside=T,ylab ="Elasticity",xlab ="Stage",col =c("black","orangered","grey","darkred"),bty ="n")legend("topright",c("det ahist","det hist","sto ahist","sto hist"),col =c("black","orangered","grey","darkred"),pch =15,bty ="n")

Figure 6.6. Ahistoricalvs. historically-corrected elasticity of lambda to stage

par(lty = lty.o)

The plot of these distributions a strong influence of history, and alsoshows the strong importance of vegetative dormancy and small adults.Both environmental stochasticity and individual history increase theelasticity associated with larger, though still small, adults.

Packagelefko3 also includes functions to conduct manyother analyses, including deterministic and stochastic life tableresponse experiments, and general projection including quasi-extinctionanalysis and density dependent analysis. Users wishing to conduct theseanalyses should see our free e-manual calledlefko3: agentle introduction, as well as through the resources found onthe{r}evolutionarydemography website.