INTRODUCTION

Structured demographic models-matrix and integral projection models (MPMs and IPMs)-are powerful tools for data-driven modeling of population and community dynamics. In contrast to MPMs for populations with discrete structure (life stage, age class, etc.), IPMs (Easterling et al., 2000) accommodate populations structured by continuous state variables, most commonly size. A related innovation of the IPM framework is its emphasis on regression-based modeling for parameter estimation, which often carries important advantages for making the most of hard-won data (Ellner et al., 2022).

A standard workflow allows ecologists to assemble an IPM from data using familiar regression tools to describe growth, survival, reproduction, and other demographic transitions as functions of size (Coulson, 2012;Ellner et al., 2016). The relative ease of regression analyses, accommodating covariates (e.g., environmental factors, experimental treatments) and complex variance structures (e.g., random effects, correlated errors), has facilitated a growing IPM literature that examines how biotic or abiotic factors affect population dynamics (e.g., Louthan et al., 2022;Ozgul et al., 2010) and explores the consequences of demographic heterogeneity associated with spatial, temporal, and individual variation (e.g., Compagnoni et al., 2016;Crone, 2016;Plard et al., 2018). The vital rate regressions (or "sub-models") are the bridge between the individual-level data and the population-level model and its predictions; it is important to get those right.

Compared with other vital rates, growth is special. The survival and reproduction sub-models only need to provide a single predicted value as functions of size (we use "size" as the name for whatever continuous variable defines the population structure). But the growth model must specify the full probability distribution of subsequent size conditional on initial size, defining the growth "kernel" G z 0 , z ð Þ that gives the probability density of future size z 0 at time t + 1 conditional on current size z at time t. Whenever survival and reproduction are size-dependent, the entire distribution of size transitions can strongly influence IPM predictions because it governs how frequently size changes are much greater or much lower than average. Easterling et al. (2000) provided the original template for modeling size transitions in IPMs. They first tried simple linear regression, assuming normally distributed size changes with constant variance. Because the residuals from this regression exhibited nonconstant variance, they used a two-step approach to estimate the size dependence in mean squared residuals (better options soon became available, such as the lme function in R). However, even after accounting for nonconstant variance, growth data may still be non-normal. Size transitions are often skewed such that large decreases are more common than large increases (Peterson et al., 2019;Salguero-G omez & Casper, 2010), or vice versa (Stubberud et al., 2019). Size transitions may also exhibit excess kurtosis ("fat tails"), where extreme growth or shrinkage is more common than predicted by the tails of the Normal distribution (Hérault et al., 2011).

The observation that the Normal (or Gaussian) distribution may poorly describe size transitions in real organisms has been made before, and several studies have emphasized that alternative distributions should be explored (Easterling et al., 2000;Peterson et al., 2019;Rees et al., 2014;Williams et al., 2012). For example, Peterson et al. (2019) showed that skewness in size transitions could be modeled through beta regression on transformed data (for reasons we describe below, this approach also has some drawbacks), or by fitting a skewed Normal distribution. They showed that incorporating skew could have important consequences for model-based inferences, and concluded that "testing of alternative distributions for growth … [should] become standard in the construction of size-structured population models." Nonetheless, default use of Gaussian growth distributions (often with nonconstant variance) remains the standard practice. The general state-of-the-art in the literature appears to remain where it was 20 or so years ago, using the default Gaussian model without examining critically whether or not it actually describes the data well. We are guilty of this, ourselves.

The persistence of Gaussian growth models is understandable. Popular packages such as lme4 (Bates et al., 2015), mgcv (Wood, 2017), and MCMCglmm (Hadfield, 2010) make it easy to fit growth models with potentially complex fixed-and random-effect structures, but the possible distributions of continuous responses are limited and default to Gaussian. Abandoning these convenient tools for the sake of more flexible growth modeling means, it may seem, sacrificing the flexibility to model diverse sources of demographic variation, some of which may be the motivation driving the study in the first place.

Our goal here is to present and illustrate a practical "recipe" that moves growth modeling past the standards set over 20 years ago. Using software tools that are now readily accessible, ecologists can escape the apparent trade-off between realistically modeling non-Gaussian size transitions and flexibly including multiple covariates and random effects. As with any recipe, users may need to make substitutions or add ingredients to suit their needs. We emphasize graphical diagnostics for developing and evaluating growth models, rather than a process centered on statistical tests or model selection. Through empirical case studies, we demonstrate how tools that were nonexistent or not readily available when IPMs first came into use now make it straightforward and relatively easy to identify when the default model is a poor fit to the data, and to then choose and fit a better growth model that is no harder to use in practice. We illustrate our approach by revisiting three published case studies (and three additional case studies in Appendix S3), including examples from our own previous work. In each case, the Gaussian assumption does not stand up to close scrutiny. We illustrate how we could have done better, and the consequences of "doing better" for our ecological inferences. All analyses were carried out in R (R Core Team, 2022) version 4.0 or higher and may be reproduced from publicly available code and data (see Data availability statement).

FLEXIBLE GROWTH MODELING

The modeling process that we suggest runs as follows (Figure 1).

Fit a "pilot" model assuming a Gaussian distribution, but allowing for nonconstant variance

This step is familiar to most IPM users, as it is the start and end of the standard approach. It may include model selection to identify which treatment effects or environmental drivers affect the mean and/or variance of future size. Nonconstant variance is often fitted in a two-stage process, first fitting mean growth assuming constant variance, then doing a regression relating the squared residuals to initial size or the fitted mean of subsequent size. Fitting mean and variance simultaneously as functions of initial size, as can be done with R packages mgcv and nlme, is advantageous when possible because incorrectly assuming constant variance can affect model selection for the mean. We illustrate both one-step and two-step approaches in the case studies below.

Allowing nonconstant variance removes the need for transforming the data to stabilize growth variance. Transformation may still be useful if it does not create new problems such as making some state-fate relationships highly nonlinear. In particular, log-transformation often reduces or eliminates heteroskedasticity in growth data (Ellner et al., 2016) and also helps avoid eviction at small sizes (Williams et al., 2012).

The fitted mean and variance functions should be checked before going any further. If they are perfectly correct, standardized residuals (residuals scaled by the SD) will have zero mean and unit variance overall, and will exhibit no trends in mean or variance with initial size or fitted mean value. However, estimates of the mean and variance functions are somewhat smoothed because of the inescapable bias-variance trade-off, so scaled residuals will retain some variation in location and scale. Given enough data, statistical tests will detect that variation. So instead, we take for granted the presence of F I G U R E 1 Recommended steps in growth modeling (left) and guide to common non-Gaussian distributions of size x for x ℝ that can accommodate different combinations of skewness and kurtosis (right). All of these distributions (often including multiple versions or parameterizations of each) are available in the R package gamlss.dist, except for the skewed generalized t, which is available in the package sgt (Davis, 2015).

trends and assess their importance by fitting nonparametric (NP) spline regression models for residuals (trend in mean) and absolute residuals (trend in variance) as a function of initial size or fitted value. The mean and variance functions can be accepted if the regression curves for the scaled residuals are nearly flat.

Use graphical diagnostics to identify if and how the standardized residuals deviate from Gaussian, and to choose a more appropriate distribution

If the Gaussian growth model is valid, the standardized residuals should be Gaussian with zero skewness or excess kurtosis. Growth data may deviate from this in many ways, and the nature of the deviations can guide the search for a better distribution. Tests such as the D'Agostino test of skewness (D'Agostino, 1970) and the Anscombe-Glynn test of kurtosis (Anscombe & Glynn, 1983) can be used to diagnose whether the standardized residuals, in aggregate, deviate from normality (Komsta & Novomestky, 2015). However, the aggregate distribution may be misleading if skewness or kurtosis vary with size or other covariates. Skewness changing from positive at small sizes to negative at large sizes might produce zero overall skewness, but really requires a distribution that can allow both positive and negative skew, such as the skewed Normal or Johnson S U distributions. Alternatively, growth data may exhibit leptokurtosis (in which case the t distribution may be a good choice) or may shift from platykurtosis to leptokurtosis depending on initial size (in which case the power exponential distribution may be a good choice). It is therefore essential to visualize trends in distribution properties with respect to either initial size, or expected future size for models with multiple covariates. Figure 1 includes guidance on how the skew and kurtosis properties of the standardized residuals suggest options for an appropriate growth distribution. In our case studies we exploit the many distributions in the gamlss R package (Stasinopoulos & Rigby, 2007), but other distribution families can be used.

Refit the growth model using the chosen distribution

In models with multiple covariates and/or random effects, each potentially affecting several distribution parameters, "refit the model" could entail a massive model selection process to identify the "best" non-Gaussian model. With so many options, model uncertainty may be overwhelming and over-fitting becomes a significant risk even when precautions against it are taken.

We therefore argue for adopting a more modest goal: remedy the defects evident in the standardized residuals of the Gaussian model. This recommendation is based on the finding that parameter estimation using Gaussian regression models is generally robust to deviations from normality of the residuals (Schielzeth et al., 2020). That is, the fitted mean of the Gaussian model (as a function of covariates) is probably a very good approximation for the fitted mean in the corresponding non-Gaussian model-and if it is not, the next step in the modeling process will catch that. The functional forms for skew and kurtosis of the non-Gaussian model can be guided by the qualitative features of the graphical diagnostics (e.g., that skewness switches from positive to negative with increasing size). As we demonstrate below, the mean and SD functions can often be carried over exactly from the pilot Gaussian model.

Evaluate the final growth model through graphical diagnostics comparing simulated and real growth data

A good model will generate simulated data that look like the real data. Again, it is important to inspect the properties of simulated data as a function of initial size, fitted mean, or other covariates rather than examining the aggregate distribution. We again suggest below graphical diagnostics, based mainly on quantiles, that can be used to compare simulated with real growth data. If the simulated data do not correspond well with the real data, alternative or more flexible distribution families should be considered, or more complex functions relating distribution parameters to size and other covariates.

HOW SHOULD SKEWNESS AND KURTOSIS BE MEASURED?

Non-Gaussian growth modeling requires scrutinizing the skewness and kurtosis of standardized residuals, so measurement of these properties warrants attention. The standard measures are based on the third and fourth central moments, respectively, of the distribution: skewness = m 3 =σ 3 , excess kurtosis = m 4 =σ 4 -3 where m k ¼  X -X À Á k is the kth central moment of a random variable X and σ 2 is the variance (second central moment). A Gaussian distribution has zero skewness and zero excess kurtosis. The standard measures are simple and easy to use, but they have poor sampling properties. Because the measures involve high powers of data values, a few outliers can produce very inaccurate estimates. Figure 2 shows a simulated example, where the underlying data are samples of 200 values from a t distribution with 8 degrees of freedom, repeated 5000 times; the true skew is 0, and the true excess kurtosis is 1.5. The distance between the largest and smallest estimates (indicated by the dotted red vertical lines), relative to the distance between the 5th and 95th percentiles (dashed blue lines), shows the broad extent of extreme values that can occur even with a large sample, especially for kurtosis.

We therefore recommend nonparametric (NP) measures of skewness and kurtosis that are based on quantiles and thus are less sensitive to a few extreme values. Let q α denote the α quantile of a distribution or sample (e.g., q 0:05 is the 5th percentile). For any 0 < α < 0:5, a quantile-based measure of skewness is given by (McGillivray, 1986)

NP Skewness measures the asymmetry between the tails of the distribution above and below the median. The size of the upper tail can be measured (for any 0 < α < 0:5) by τ U ¼ q 1 -αq 0:5 ; for α ¼ 0:05 this is the difference between the 95th percentile and the median. The lower tail size is τ L ¼ q 0:5q α . The definition above is equivalent to

An NP Skewness of AE0:2 says that the difference in tail sizes is 20% of their total. The range of possible values is -1 to 1. Both α ¼ 0:25 (sometimes called "Kelly's skewness") and α ¼ 0:1 ("Bowley's skewness") are common choices. We used α ¼ 0:1. An analogous quantile-based measure of kurtosis (Jones et al., 2011) is

For α ¼ 0:05, NP Kurtosis is the difference between the 95th and 5th percentiles, relative to the interquartile range. To facilitate interpretation, we scale NP Kurtosis relative to its value for a Gaussian distribution, and subtract 1 so that the value for a Gaussian is zero. We call this "NP Excess Kurtosis." A value of AE0:2 means that the tails are on average 20% heavier than those of a Gaussian with the same interquartile range. We calculate NP Kurtosis using α ¼ 0:05, to focus on the tail edges, but again this is somewhat arbitrary. Figure 2C,D illustrate how, applied to the same simulated samples, the NP measures produce a smaller fraction of highly inaccurate estimates caused by a few extreme values. Also note that, in contrast to the moment-based measures, numerically small values of the NP measures (e.g., 0.1 or 0.2) should not be disregarded, because both measures are scaled so that a value of 1 indicates extremely large departures from a Gaussian distribution.

Using quantile-based measures carries the added value that quantile regression can be used to estimate how they vary with initial size or expected future size. In the examples below, we use the qgam package (Fasiolo et al., 2020) to fit spline quantile regression models, which accommodate nonlinear size-dependence in skewness and kurtosis. One risk of spline regression is that fitted quantiles may be excessively "wiggly" without constraints on their complexity; with realistic amounts of data, we can hope to estimate broad trends in distribution shape, but not fine-scale variation. In the examples below, we limit complexity by fitting splines with k ¼ 4 basis functions unless otherwise noted. Parametric quantile regression is also an option.

For consistency, we also use quantile-based measures of mean and SD when comparing real and simulated data, and use quantile regression to visualize their trends. Specifically, following Wan et al. (2014), NP mean ¼ q 0:25 + q 0:5 + q 0:75 3 , NP SD¼ q 0:75q 0:25 1:35 : ð4Þ

CASE STUDY: LICHEN, VULPICIDA PINASTRI

We begin with a simple example where current size is the only predictor of future size. Growth data for the epiphytic lichen Vulpicida pinastri were analyzed first by Shriver et al. (2012) and again by Peterson et al. (2019) in their study of skewed growth distributions. We therefore had an a priori expectation of deviation from normality.

The data set includes 1542 interannual transitions in thallus area (in square centimeters) observed from 2004 to 2009 in Kennicott Valley, AK. Shriver et al., 2012 used a mixture distribution that separated "normal growth or shrinkage" from "extreme shrinkage." We aimed to fit a single growth model that could realistically accommodate both types of size transition without requiring ad hoc decisions about which observations of shrinkage were "extreme" or not.

With initial size as the only predictor, a convenient way to fit a Gaussian model with nonconstant variance is the gam function in the mgcv library (Wood, 2017) using the gaulss family. Following a bit of model selection, we fit the mean and SD of future size as second-order polynomials of current size, then calculated the scaled residuals from the fitted mean and SD. (gam() is most commonly used to fit smooth splines (s()) for predictor variables, but it can also fit parametric regressions.) Here, the first argument to gam() is a two-element list that defines the linear predictors for mean and SD:

The data and fitted mean and SD are shown in Figure 3A, and the corresponding diagnostic plots are in Figure 4A,B. Our diagnostic plots are similar to plots made by R's plot.lm function, except that we use spline regression to allow data-driven choice of curve smoothness, and use absolute residuals (rather than their square roots) so that the SD of the regression curve is on the same scale as the residuals. The spline curves are not exactly flat-their SDs, given above each panel, are positive-but the trends are much too small to be worth fixing.

Quantile regression on the scaled residuals generates the skewness and kurtosis diagnostics shown in Figure 3B. As expected based on previous analyses, the graphical analysis of the standardized residuals indicates negative skew, especially at larger sizes (Figure 3B). We also find positive excess kurtosis for all sizes.

We turned to the Johnson's S-U (JSU) distribution for improvement. The JSU is a four-parameter leptokurtic distribution allowing positive or negative skew, with the convenient property that its location and scale parameters mu and sigma are the mean and SD, respectively, which greatly facilitates the transition from a pilot Gaussian model. JSU is not available in any standard linear or additive modeling packages, to our knowledge. But that is not a barrier because we can write a likelihood function using the dJSU() function in the gamlss.dist package. Following the best fit Gaussian model, we defined mu and sigma of the JSU as quadratic polynomials of initial size and, based on Figure 3B we define the skewness parameter nu as a linear function of size and kurtosis parameter tau as a positive constant. The likelihood function therefore has nine parameters to estimate. We fit the model using the maxLik package with starting coefficient values for mu and sigma based on the pilot Gaussian model. We chose maxLik because it offers the BHHH optimization method, which works well for non-Gaussian likelihoods in our experience.

Simulating data from the fitted JSU model indicates a compelling improvement over the best Gaussian model, not only in skewness and kurtosis (Figure 5C,D) but also the NP SD (Figure 5B). Note, in these data simulation figures Gaussian and non-Gaussian data are offset by an arbitrary amount to more easily visualize their correspondence to the real data (black lines in Figure 5).

To understand the practical consequences of improved growth modeling, we assembled the remainder of the lichen IPM following Shriver et al. (2012). The asymptotic population growth rate λ based on Gaussian growth differs from the JSU growth model by about 1% annual population growth (Table 1), in line with results of Peterson et al. (2019). However, even this modest difference can lead to biased estimates of extinction risk from the Gaussian model, particularly over longer time horizons (Figure 6). We also explored differences in other life history metrics (Table 1) using functions from Hern andez et al. (2024). For example, the JSU growth model predicts values for mean lifespan, mean lifetime reproductive success, and generation time that are 15%-25% lower than the Gaussian growth model. In this case study, properly modeling non-normal size transitions-which was easy to do with a few extra lines of code-can influence ecological inferences, at least based on point estimates. However, Table 1 also provides bias-corrected, bootstrapped CIs (Diciccio & Efron, 1996), and these are heavily overlapping between the Gaussian and JSU models for all life history traits, suggesting that effects of "improved" growth modeling are small relative to our uncertainty in model parameters.

One could argue that this example was a convenient "straw man" to disqualify Gaussian growth, because it was recognized by the original and subsequent analysts that size transitions are strongly skewed (Peterson et al., 2019;Shriver et al., 2012). In all remaining case studies, including those in Appendix S3, we reexamine growth data that were modeled as Gaussian in the original published analysis.

CASE STUDY: TREE CHOLLA CACTUS, CYLINDROPUNTIA IMBRICATA

The next case study, focused on the tree cholla cactus Cylindropuntia imbricata at the Sevilleta Long-Term Ecological Research site in central New Mexico, adds a new feature to the simple size-dependent regressions in the previous study: random effects associated with temporal (year) and spatial (plot) environmental heterogeneity. This long-term study was initiated in 2004, and different subsets of the data have been analyzed in various IPM studies, all using Gaussian growth kernels (Compagnoni et al., 2016;Czachura & Miller, 2020;Elderd & Miller, 2016;Miller et al., 2009;Ohm & Miller, 2014). In fact, Elderd and Miller (2016) presented a Gaussian growth model as an example of a well-fit growth function, based on an overall distribution of residuals that appeared Gaussian and posterior predictive checks (PPCs) of a Bayesian model that suggested consistency between the real data and data simulated from the fitted model (fig. 4 in Elderd & Miller, 2016). While PPCs and the associated "Bayesian p-value" are popular diagnostic tools, they are often too conservative (Conn et al., 2018;Zhang, 2014), failing to reject marginally bad Miller (2020). Following previous studies, we quantified size as the natural logarithm of plant volume (in cubic centimeters), derived from height and width measurements.

We begin growth modeling, as above, with a generalized additive model with the mean and SD of size in year t + 1 modeled as smooth function of size in year t, with random intercepts for year and plot and assuming normally distributed residuals:

Note that here we fitted the SD function with k ¼ 6 basis functions rather than our default of k ¼ 4 because, in a preliminary analysis, we found a moderate variance trend in the standardized residuals using k ¼ 4, suggesting a need for greater flexibility. With k ¼ 6, spline regression detected essentially no trend in the mean of the resulting standardized residuals (Figure 4C,D).

The growth variance is estimated to peak at small to medium sizes (Figure 3C). The standardized residuals show clear signals of negative skew and positive excess kurtosis across most of the size distribution, but strongest in the middle (Figure 3D). We therefore need a distribution family allowing negative skew and positive excess kurtosis, both of which may be negligible at some sizes. We first tried Johnson's S U and then the skewed t distributions, which provided some improvements but there were still visible discrepancies between simulated and real data. We next turned to the SHASH distribution, which allows a greater range of kurtosis for a given amount of skew, and vice versa (Jones and Pewsey 2009; T A B L E 1 Life history attributes derived from integral projection model (IPM) kernels that included Gaussian or "improved" growth sub-models for six case studies.

Species

Growth model λ Lifespan Lifetime reproductive output Age at reproduction Generation time Lichen (V. pinastri) Gaussian 1.01 (0.99, 1.04) 6.4 (3.6, 11.1) 1.4 (0.5, 3.1) 6.5 (5.7, 7.3) 40.8 (30.5,57.4) Improved 1.00 (0.98, 1.03) 5.4 (3.1, 9.7) 1 (0.4, 2.4) 6.4 (5.4, 7.3) 36.6 (27.5, 48.6) Cactus (C. imbricata) Gaussian 0.994 (0.99, 0.996) 6.11 (3.66, 8.63) 21.8 (8.27, 49.4) 17.6 (1.75, 22.7) 189 (131, 266) Improved 0.993 (0.991, 0.998) 5.38 (3.34, 16.3) 13.4 (5.72, 251) 20.3 (1.21, 22.2) 179 (133, 298) Orchid (O. purpurea) Gaussian 1.09 (1.08, 1.1) 1.08 (1.06, 1.11) 20.0 (12.6, 31.0) 5.07 (4.78, 5.31) 104 (73.1, 150) Improved 1.09 (1.08, 1.1) 1.08 (1.06, 1.1) 19.3 (12.0, 29.9) 5.03 (4.75, 5.3) 100.7 (71.0, 145.0) Pike (E. Lucius) Gaussian 1.62 (1.35, 1.89) 1.2 (1.09, 1.35) 5.75 (2.9, 9.7) 1.09 (1.03, 1.18) 4.96 (4.26, 5.84) Improved 1.62 (1.35, 1.88) 1.2 (1.09, 1.35) 5.76 (2.91, 9.73) 1.09 (1.03, 1.18) 4.94 (4.30, 5.84) Creosote (L. tridentata) Gaussian 1.033 (1.029, 1.04) 4.52 × 10 6 (2.14 × 10 5 , 1.82 × 10 8 ) 3.19 × 10 5 (1.27 × 10 4 , 1.24 × 10 7 ) 32.7 (29.2, 36.0) 5.27 × 10 6 (2.50 × 10 5 , 1.95 × 10 8 ) Improved 1.034 (1.03, 1.04) 3.26 × 10 5 (1.98 × 10 3 , 1.66 × 10 7 ) 2.31 × 10 4 (5.83 × 10 2 , 1.27 × 10 6 ) 32.8 (29.3,36.0) 3.7 × 10 5 (2.63 × 10 3 , 1.93 × 10 7 ) Coral (G. ventalina) Gaussian … 17.3 (11.9, 24.3) … 10.5 (9.3, 11.8) 31.6 (28.3, 36.7) Improved … 17.5 (12.1, 24.3) … 10.7 (9.4, 12.2) 30.9 (27.4, 35.3) Note: The improved distributions were Johnson's S-U (JSU) (lichen, creosote), SHASH (cactus, pike, coral), and skewed t (orchid). Pike, creosote, and coral case studies are presented in Appendix S3. The original coral case study assumed an open population with constant recruitment from a large source region, so some life history attributes cannot be computed from the published model. Values in parenthesis are 95% bootstrap CIs, specifically the bias-corrected (BC) bootstrap CIs (Diciccio & Efron, 1996). Table can be reproduced from scripts crossspp_growth.R, Vulpicida_boot.R, Akumal_corals_boot.R (Ellner & Miller, 2025). Appendix S1). This flexibility proved necessary to generate simulated data that compared favorably with the real data, so we proceeded with the SHASH. Conveniently, SHASH is available as an mgcv family, allowing for flexible size-dependence in skewness and kurtosis without having to select specific size-dependent functions.

Here, the first argument to gam() is now a four-element list specifying the linear predictors for the four parameters of the SHASH distribution.

Data simulated from the SHASH model compared favorably with the real data (Appendix S4: Figure S1). Similar to the lichen case study, we see that correctly modeling skewness and kurtosis improved estimation of the NP mean and SD (Appendix S4: Figure S1A,B), yielding a growth model that is truer to the data.

We next explored how improved growth modeling influenced IPM results. The λ values predicted by Gaussian and SHASH growth functions, corresponding to the average plot and year, were nearly identical (Table 1) but we could also leverage structure of the study design to quantify demographic variance associated with temporal and spatial heterogeneity. We used the fitted random effects from the vital rate models to estimate the asymptotic growth rate for each year (λ t ), centered on the average plot, and for each plot (λ p ), centered on the average year. Estimates of λ t from the Gaussian growth model were often greater than estimates from the SHASH growth model, particularly in some of the harshest years (Figure 7A), and therefore the Gaussian model predicted lower temporal variance in fitness (SD

Plot-to-plot variation was more similar between the two models (SD

, although spatial variation in fitness was much lower than temporal variation (Figure 7B). The difference in temporal variance would suggest that Gaussian growth modeling would predict a higher stochastic growth rate λ S , because temporal variance has a negative effect on λ S . However, the stochastic growth rate from the Gaussian growth model (λ S ¼ 0:992) was nearly identical to that of the SHASH growth model (λ S ¼ 0:991). This is likely because temporal fluctuations in vital rates, which is where the SHASH growth model would make a difference, have a weaker influence on λ S than the temporal fluctuations in size structure that they generate (Compagnoni et al., 2016;Ellis & Crone, 2013). The SHASH and Gaussian growth models predicted differences in other life history traits but, as in the lichen case study, these differences were small relative to the uncertainty captured by bootstrapped CIs (Table 1). Interestingly, the SHASH model had wider uncertainty intervals, particularly for lifespan and lifetime reproductive output, presumably because the additional parameters it requires introduce additional sources of uncertainty in these estimates. Thus, in this case study, modeling non-Gaussian size transitions with a Gaussian growth model may or may not influence IPM results depending on the target of the analysis, and whether the emphasis is on point estimates or uncertainty intervals.

CASE STUDY: LADY ORCHID, ORCHIS PURPUREA

Our final case study examines selection on life history strategies in the lady orchid Orchis purpurea. In a prior study, Miller et al. (2012) analyzed how costs of reproduction (flowering or not in year t) affected growth from year t to t + 1. The two growth kernels for flowering and nonflowering were then used in an IPM to quantify the optimal flowering size that balances the benefits of waiting to flower at larger sizes against the greater risk of death before flowering. The original study assumed Gaussian size transitions with nonconstant variance depending on initial size. Here we revisit that analysis to derive improved growth kernels. We use this case study to illustrate several new elements and challenges, including modeling skewness and kurtosis as functions of expected future size.

The data, originated by Dr. Hans Jacquemyn, come from 368 plants in a Belgian population censused annually from 2003 through 2011. Here we use data only from the "light" habitat in the original study. We used the natural logarithm of total leaf area as the size variable in the IPM.

As a variation on software, we fitted the pilot Gaussian model using the lmer function in the lme4 package, as in the original study. We fit three candidate linear models that included fixed effects of size in year t (model 1), additive effects of size and flowering status in year t (model 2), or an interaction between size and flowering (model 3), all including random intercepts for year. The interaction model was strongly favored (difference in AIC between the best and next-best models [ΔAIC] = 10.5). Unlike our previous case studies, here we have multiple fixed effects (initial size and flowering status) that may influence the variance of future size. In cases such as this it is convenient to model variance as a function of expected future size, rather than initial size as we did with the lichens and cacti. The expected (or "fitted") values reflect the combined influence of all fixed and random effects, and therefore implicitly account for multiple sources of variation in the variance.

Models where error variance is a function of fitted values cannot be fitted directly with lme4 (nor in the mgcv functions for generalized additive models). But it can still be done with lmer through an iterative re-weighting approach, as follows. In lmer, weights w i can be used to indicate that the observations y i have error variance proportional to 1=w 2 i . The iterative steps are as follows, and code that executes these steps is in orchid_growth_modeling.R.

1. Fit the expected value assuming Gaussian-distributed residuals with constant variance. 2. Fit the SD of the residuals as a function of the corresponding fitted value. 3. Refit the model, with weights equal to the inverse of the SD estimated in step 2.

We iterated steps 2 and 3 until the root mean square change in weights was below 10 -6 . This is not elegant, but it works and converges quickly. In step 2, we modeled the log of the SD (because SDs cannot be negative) as a quadratic polynomial in the fitted mean. In exploratory analyses we found that the quadratic term was necessary to fit the SD. We did this for all candidate models and, for a fair Akaike information criterion (AIC) comparison, we then refit all candidate models with the weights estimated from the top model.

The updated model selection continued to favor the size × flowering interaction model (3), but now with a weaker improvement over the next best model (ΔAIC = 6.7). The fitted mean (a function of initial size and flowering status) and fitted SD (a function of the fitted mean) are shown in Figure 3E. Spline regression found no trend in the mean of the resulting standardized residuals, and only small variation in the variance (Figure 4E,F).

The best Gaussian model indicated a growth cost associated with flowering at the start of the census interval and a decline in growth variance with increasing expected values (Figure 3E). The standardized residuals indicated negative skewness (10%-20% difference in tail weight) and excess kurtosis (10%-40% fatter than Gaussian) across much of the size distribution but both were negligible at large expected sizes (Figure 3F).

As possible improvements, we explored the skewed t and JSU distributions, both leptokurtic distributions with flexible skewness. Based on comparisons between real and simulated data we were happier with the skewed t, which we fit with a custom likelihood function similar to the JSU growth model for the lichen data. However, rather than refitting all parameters of the skewed t model, as we did with the lichen JSU, we built a "hybrid" likelihood function that uses the fitted mean and SD from the best Gaussian model, and estimates parameters Ellner & Miller, 2025).

that control skewness and kurtosis as linear functions of expected future size. This is easy because the gamlss.dist package provides a parameterization of the skewed t in which the location parameter μ is the mean and scale parameter σ is the SD (Rigby et al., 2019). The hybrid likelihood looks like this:

Based on diagnostics of the standardized residuals, parameters that control skewness and kurtosis are defined as linear functions of the mean (note that the tau parameter uses a log x -2 ð Þlink function). This approach relies on the robustness of fitted Gaussian models to deviations from normality, which implies that the fitted mean and variance from a Gaussian model are good approximations for the mean and variance of the corresponding non-Gaussian model. If one is skeptical of this approach, it is possible to simultaneously refit all parameters of the skewed t. However, recall that unlike the lichen case study, the pilot Gaussian model here includes random year effects, and the expected values getting passed into dSST account for this source of variation. Estimating random effects "from scratch" with a custom likelihood model is possible (we provide guidance on doing this with a "shrinkage" approach, in Appendix S2), but generally should not be necessary. Instead, a key advantage of the hybrid approach is retention of the fitted random effects and associated variance components, which get shuttled from the Gaussian model into the non-Gaussian model without any fuss (though it was critical to use a parameterization of the skewed t for which mu is the mean and sigma is the SD). And, if this approach does not "work" (i.e., deviations from normality biased the fitted values of the Gaussian model) one would quickly find out when comparing simulated with real data. In this case, size transition data simulated from this model corresponded favorably to the real data, much better than the pilot Gaussian model, including improvements in the SD, skewness, and kurtosis of future size (Appendix S4: Figure S2).

Finally, we used the improved growth model to revisit key results of the original study. Miller et al. (2012) used the orchid IPM to estimate the evolutionarily stable strategy (ESS) as the mean size at flowering that maximizes lifetime reproductive success (R 0 ), given the constraint that flowering when small reduces growth and thus elevates mortality risk. Repeating that analysis here, we found that improved growth modeling has virtually no influence on predictions for optimal life history strategies (Figure 8). ESS flowering sizes were nearly identical between IPMs with Gaussian versus skewed t growth models, and both aligned well with the observed mean flowering size (dashed vertical line in Figure 8). Similarly, there were very small differences between growth functions in other metrics of orchid life history and, again, these differences were overwhelmed by uncertainty associated with parameter estimation (Table 1).

DISCUSSION

Much of the appeal of IPMs has stemmed from their embrace of continuous size structure through regression-based approaches, and the potentially complex fixed-and random-effect structures that those approaches allow. Using familiar statistical tools and with relatively few parameters to estimate, IPM users can incorporate important sources of variation in demography and interrogate their influence on ecological and evolutionary dynamics. With this opportunity comes the burden of getting it right: an IPM is only as good as the statistical sub-models for the underlying data. The growth sub-model is the trickiest part because it defines a distribution of future size conditional on current size. Distributions have many properties-"moments"-and a good growth model should recapitulate the properties of real size transitions. The default assumption of Normally distributed size transitions, employed overwhelmingly across 20+ years of IPM studies, is an arbitrary historical precedent. In our case studies and, we suspect, more broadly, skewness and excess kurtosis were common features of size transitions. Our most important message is that the assumption of normally distributed size transitions can easily be abandoned, and a more inquisitive process of growth modeling should take its place.

We have attempted to lay out what that process should look like, emphasizing visual diagnostics to characterize how data deviate from Gaussian. One implication of relying on visual diagnostics is that goodness of fit is in the eye of the beholder. This empowers IPM users to make informed choices, but it is not very prescriptive; we have not suggested any hard rules for choosing among distributions, only that a good growth model should generate data that look like the real thing. Alternatively, model selection could be used to identify best fitting growth distributions and best fitting functions for higher moments. However, model selection among growth distributions with 3-5 parameters, each of which may be functions of multiple state variables or fitted values, can quickly explode in complexity, and we are not convinced it is worth the trouble.

Our work follows the important contribution of Peterson et al. (2019), who were similarly motivated by inadequacy of the Gaussian model but arrived at different recommendations. These authors developed a creative approach in which size data are transformed onto a 0, 1 ½ scale and size transitions on that scale are modeled using beta regression. The beta distribution can accommodate positive, negative, or zero skew, potentially varying with size, so the Peterson et al. approach is a flexible option for skewed growth data. However, beta regression also has some limitations: common beta regression packages do not fit random effects (e.g., betareg; Cribari-Neto & Zeileis, 2010) or do not do so reliably (in our experience gamlss regressions with random effects are numerically unstable); and the two-parameter beta distribution does not allow skewness and kurtosis to be fitted independently. Additionally, the initial transformation onto 0, 1 ½ scale requires estimating extreme quantiles of the growth distribution (e.g., 0.01 and 0.99) as a function of initial size. In our experience those quantile estimates can be very sensitive to how size dependence is modeled, and model selection is challenging for extreme quantiles where data are (by definition) very sparse. Rather than picking one distribution as a new default, users can leverage the vast arsenal of continuous probability distributions-all at one's fingertips with a few lines of code-so that the data and their particular deviations from normality can guide the choice of a better distribution. It is also possible to use mixtures of multiple growth distributions, as done by Shriver et al., 2012 to model "normal" and "extreme" types of lichen shrinkage. In reanalyzing that data set, we found that a single, flexible distribution (Johnson's S-U) could recapitulate the observed size transitions, but there may be other cases where mixture distributions are preferable or necessary.

In all of our case studies, non-Gaussian growth models always yielded more satisfying fits to size transition data than the Gaussian models published in those papers. However, to our relief, none of these reanalyses yielded a "gotcha" result that overturned results of the original study. In fact, in this small sampling of case studies, improved growth modeling had weak to modest effects on IPM results, similar in magnitude to the results of Peterson et al. (2019), and for most species and life history metrics, these effects were overwhelmed by the uncertainty associated with parameter estimation (Table 1). For some case studies, one might argue that non-Gaussian modeling was not worth the trouble-only it was almost no trouble at all, and we could not have known whether or not a non-Gaussian model would have made a difference before fitting it. Even where Gaussian and improved growth models differed in IPM results, we do not know "true" values to compare them against, only that one describes the data better than another, and a closer match to the data does not necessarily translate to better predictive ability due to the bias-variance trade-off.

We caution against taking too much comfort in weak effects of "improved" growth modeling; in other scenarios, the choice of growth distribution could be more consequential. We focused on life history metrics such as mean lifespan and mean lifetime reproductive success (Table 1). It is possible that higher moments of those traits (e.g., variance or skewness of lifespan and lifetime reproduction) are more sensitive to the tails of the growth distribution. It is also worth noting that most of our case studies focused on perennial life histories (perennial plants and lichens) characterized by relatively slow growth, heavy losses during recruitment, and high survival once established, and these species all had mean lifespans between one and six years and generation times on the order of decades. Life histories such as these may be relatively robust to subtle features of the growth kernel. In Appendix S3 we present three additional case studies that broaden our life history coverage, including pike (Esox lucius), a fish with a generation time of four to five years, and creosote bush (Larrea tridentata), a desert shrub that is virtually immortal once established. Life history metrics from these "fast" and "slow" populations were no more sensitive to improved growth modeling than those of the perennial plants and lichens (Table 1). More systematic comparative analyses may provide insight into which types of species and life histories are more likely to exhibit strong skewness and kurtosis, and which demographic quantities are more or less sensitive to these features of size transition.

Our case studies illustrate a diversity of software packages and computational approaches to reflect the diversity of preferences and habits that the community of IPM analysts bring to their own problems. We like spline generalized additive models (gams) for their flexibility and for mgcv's numerous options for distribution families and overall speed and reliability. However, there are some applications for which classical parametric regression would be preferable because the coefficients carry biological meaning. For example, regression coefficients may be targets of natural selection (Rees & Ellner, 2016) and may combine to influence traits of interest such as the expected size at flowering (e.g., in Figure 8A), a function of the intercept and slope of the size-dependent flowering function (Metcalf et al., 2003). Some potentially useful distributions are not available in linear modeling software packages, but that should not be a barrier to their use: as in several of our case studies, custom likelihood functions allow non-Gaussian models without sacrificing the complex, multilevel features that one might be accustomed to fitting in lme4, for example. Bayesian analysis may further broaden the options for non-Gaussian candidate distributions and may help estimate hard-to-fit parameters through the brute force of sampling algorithms. Bayesian analysis also provides a natural way to propagate uncertainty from vital rate sub-models to full model predictions (Elderd & Miller, 2016), as an alternative to our bootstrapping approach. However, as of this writing, most of the non-Gaussian distributions that we have discussed are not available in popular Bayesian software packages such as Stan or JAGS. While user-defined distributions can always be coded from scratch, this may be a significant technical barrier for many IPM analysts.

From the outset there have been concerns about "how well these methods [IPM growth kernels] can deal with different patterns of growth, stasis, and shrinkage" (Morris & Doak, 2002, p. 200), compared with "binning" methods that use observed transition frequencies between user-defined size classes as the transition probabilities in a (possibly large) matrix model (Doak et al., 2021). The non-Gaussian models that we have considered here are not a panacea. For example, none of them allow bimodal growth, such as might occur if herbivore-or pathogen-attached individuals experience rapid tissue loss. When the shape of the growth distribution is nearly the same for all initial sizes, a nonparametric IPM growth kernel can be defined from a kernel density estimate for scaled residuals (Ellner et al., 2016, p. 288). Outside that special situation, nonparametric approaches require choosing multiple smoothing parameters, which is very challenging. We are currently exploring whether "targeted learning" approaches developed for causal inference (van der Laan & Rose, 2011) can be used to circumvent smoothing parameter selection. Targeted learning starts with a pilot model and updates it iteratively to achieve unbiased estimates and valid CIs for a particular "target" quantity, such as λ or mean lifespan. Preliminary results suggest that targeted learning with a deliberately under-smoothed pilot model works well for complex growth patterns (Zhou & Hooker, 2024). But nonparametric methods are data-hungry, so when departures from Gaussian are quantitative rather than qualitative, parametric modeling as developed here will make more efficient use of limited data.

CONCLUSION

Gaussian-distributed size transitions are probably the exception in nature, not the rule, yet two decades of IPM studies have relied overwhelmingly on Gaussian growth models. Using tools not available when IPMs were first developed, it should often be possible now to make major improvements over a Gaussian model, without worrying about finding the "best" alternative. By generating predicted size transitions that are truer to the data, IPM analysts can narrow the gap between model and nature.