Comparison of marginal and mixed-effects complementary log-log regression models for predicting planted silver birch mortality Jouni Siipilehto , Daesung Lee * Natural Resources Institute Finland (Luke), Latokartanonkaari 9 00790 Helsinki, Finland A R T I C L E I N F O Keywords: Betula pendula Roth. Tree mortality Generalized linear mixed models Self-thinning Simulation A B S T R A C T Mortality is a key process in forest succession, yet modelling individual tree mortality presents significant challenges. In this study, tree mortality models for silver birch (Betula pendula Roth.) were developed to address these complexities. The modelling data comprised thinning trials for planted silver birch established between 1981 and 1991 in southern and central Finland. Thirteen experiments were established on former agricultural land and eight experiments were on forest land. The test data comprised planted silver birch stands of a spacing trial established on agricultural land in the early 1970s. The modelling options included four different types of models based on different random effect structures: a marginal model without random effects, a random site as RND_SITE, a random plot nested within the site designated as RND_PLOT(SITE), and a random year nested within the site designated as RND_YEAR(SITE) in a linear mixed-effects complementary log-log (CLL) regression. The CLL models were evaluated according to fit statistics, with the RND_YEAR(SITE) model demonstrating the best results. Furthermore, all mortality models were implemented into the MOTTI simulator to evaluate the devel opment of planted silver birch stands in terms of stem number (N, trees ha− 1) and stand basal area (G, m2 ha− 1). In the MOTTI evaluation, unthinned stands were selected, and the data were divided into density groups: initially dense (N > 2000 trees ha− 1), normal density stands (1000 trees ha− 1 ≤ N ≤ 2000 trees ha− 1), and sparse stands (N < 1000 trees ha− 1). The independent dataset demonstrated optimal performance with the RND_YEAR(SITE) model. The current MOTTI model performed generally well but underestimated N and G for the normal density stands compared to the new model options. Finally, when examining the compatibility of the RND_YEAR(SITE) model with the existing and recently introduced stand self-thinning models, the recent model demonstrated high compatibility, while the existing model showed a clear underestimation. 1. Introduction Reliable estimates of long-term forest development require a thor ough knowledge of the structure and dynamics of forests at the stand level. This includes understanding the effects of various abiotic and bi otic factors, as well as the impacts of silvicultural treatments on stand development (Fortin et al., 2008; Siipilehto et al., 2020; Sims et al., 2014). One of the main processes of forest succession is tree mortality (Jutras et al., 2003). However, modelling individual tree mortality is challenging. There are many reasons for tree mortality, with competi tion between trees being the main cause. Harvesting operations can also contribute to tree mortality. Finally, there is rare catastrophic mortality, which refers to large-scale tree death caused by natural disasters such as storms, browsing damage from herbivores, insect pests, droughts, flooding, or other unusual events (Vanclay, 1994). When all the different mortality causes are included in the modelling data, competition-induced mortality (referred to as regular mortality) shapes the survival/mortality curve, for example, in the logistic function (Laarmann et al., 2009; Sims et al., 2014). In contrast, other causes of mortality (density-independent irregular mortality) primarily affect the intercept term, which represents the overall level of survival/mortality rate. Laarmann et al. (2009) fitted simple mortality models for different causes of tree mortality in managed and semi-natural forests. Their findings showed that relatively small trees had a high probability of mortality due to competition, while larger trees were susceptible to damage from game, insects, and wind. Furthermore, they found no relationship between tree size and mortality caused by fungi and dis eases. Stand management practices, such as thinning, can affect tree mortality (Fridman and Ståhl, 2001; Sims et al., 2014; Siipilehto et al., 2021). However, stands that have been thinned between adjacent * Corresponding author. E-mail addresses: jsiipilehto@gmail.com (J. Siipilehto), daesung.lee@luke.fi (D. Lee). Contents lists available at ScienceDirect Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel https://doi.org/10.1016/j.ecolmodel.2025.111349 Received 11 January 2025; Received in revised form 20 August 2025; Accepted 16 September 2025 Ecological Modelling 510 (2025) 111349 0304-3800/© 2025 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). https://orcid.org/0000-0002-5661-8972 https://orcid.org/0000-0002-5661-8972 https://orcid.org/0000-0003-1586-9385 https://orcid.org/0000-0003-1586-9385 mailto:jsiipilehto@gmail.com mailto:daesung.lee@luke.fi www.sciencedirect.com/science/journal/03043800 https://www.elsevier.com/locate/ecolmodel https://doi.org/10.1016/j.ecolmodel.2025.111349 https://doi.org/10.1016/j.ecolmodel.2025.111349 http://creativecommons.org/licenses/by/4.0/ measurements are often excluded when modelling mortality or survival (Eid and Tuhus, 2001; Jutras et al., 2003; Boeck et al., 2014). In a growth simulator designed for forest management, mortality processes are described by empirical models. These may include stand- level models for maximum density as a self-thinning line (e.g., Hynynen, 1993; Monserud et al., 2004; Lee and Choi, 2019; Lee et al., 2025) and individual tree mortality, which is typically modelled using logistic models (Monserud, 1976; Monserud and Sterba, 1999; Eid and Tuhus, 2001; Fridman and Ståhl, 2001). In Finland, published mortality models for silver birch (Betula pendula Roth.) feature self-thinning (Hynynen, 1993) fitted for untreated control plots, which consisted of 10 stands, 17 sample plots, and 59 observations. Silver birch is a light-demanding species which benefits from quite heavy early thinnings (Niemistö, 1997; Hynynen et al., 2010). Due to its light-demanding nature, the maximum density of a birch stand is typically lower than that of Scots pine (Pinus sylvestris L.) or Norway spruce (Picea abies L.) stands, as illustrated in Fig. 9 of Hynynen et al. (2002). Existing models for individual tree mortality for birch species were based on five thinning experiments established as permanent plots on either planted or naturally regenerated stands, consisting of both silver and downy birch (Betula pubescens Ehrh.) (Hynynen et al., 2002). Much new data are now available to update mortality models, particularly for planted silver birch. Recently published self-thinning models for planted silver birch utilized new data from long-term thinning, spacing, and breeding trials in southern and central Finland (Lee et al., 2025). Nevertheless, recent studies on individual tree mortality for silver birch trees are lacking, especially concerning modelling efforts in the Nordic and Baltic countries (e.g., Maleki and Kiviste, 2016). Moreover, it is important to determine if stand-level constraints are essential compo nents when applying an individual tree mortality model. Thus, it’s necessary to cross-check and validate the stand-level complementary self-thinning models, as investigated by Hamilton (1990) and Monserud (2004). Therefore, the objective of this study was to develop a model for predicting individual tree mortality for planted silver birch, taking into account the potential impact of thinnings and site type. The fitted marginal and mixed effects models were implemented in the MOTTI simulator to assess birch stand development against an independent test dataset. Finally, the compatibility of the new individual tree mortality models with the existing self-thinning model (Hynynen, 1993) and the recent self-thinning models (Lee et al., 2025) was tested in a manner similar to that employed by Monserud et al. (2004), who used the Austrian stand simulator PROGNAUS. Specifically, the individual tree mortality models developed in this study were validated against stand-level complementary constraints. 2. Materials and methods 2.1. Study materials for tree mortality 2.1.1. Modelling dataset The data from 17 experimental areas used for planted silver birch mortality are shown in Table 1. The stand plots were located in southern and central Finland. The southernmost experiment was located in Lammi and the northernmost was in Muhos (Fig. 1). The experiments were established from 1981 to 1991 as a thinning trial for silver birch. They were followed mainly in 5-year periods until now, encompassing a total of 25 to 35 years, which corresponds to 5 to 7 growth periods of 5 years each. Experiments were established either on former agricultural land, which represents the most fertile forest site of Oxalis-Maianthemum type (OMaT, 10 sites), or on forest land, which includes herb-rich Oxalis- Myrtillus type (OMT, 2 sites) and mesic heath Myrtillus type (MT, 5 sites), according to Cajander’s (1926) forest site types. Typical plot areas were 700, 1000, and 1200 m2. The total number of tree-level observations was 91,559 in the entire dataset. Due to the thinning experiments, thinning intensity (no thinning, light, medium, and heavy), and number Table 1 Summary statistics of the tree, plot, and site-level characteristics in the mortality modelling dataset. Variable Mean SD Min Max n Tree-level characteristics ​ ​ ​ ​ ​ dbh, cm 13.9 5.2 1.5 35.5 91,558 BAL, m2 ha− 1 10.4 7.0 0.03 37.6 ​ Tree mortality 0.025 0.17 0 1 ​ plot-level characteristics ​ ​ ​ ​ ​ G, m2 ha− 1 16.2 7.2 2.6 37.6 1110 N, trees ha− 1 903 636 186 3072 ​ DQ, cm 16.5 4.4 4.1 27.7 ​ DG, cm 17.1 4.3 4.6 28.1 ​ HG, m 19.4 5.5 4.7 30.6 ​ T, year 33 9.9 13 62 ​ site-level characteristics ​ ​ ​ ​ ​ Plot area, m2 982 143 490 1250 17 Temperature sum, dd ◦C 1105 63.2 1017 1235 ​ Measurement interval, year 4.9 0.4 4 6 ​ Thinned 0.76 0.50 0 1 ​ Time since last thinning, year 6.7 6.9 0 35 ​ Note: n is the number of observations; dbh is tree diameter at breast height at a height of 1.3 m above the ground; BAL is the sum of the basal areas in larger trees; Tree mortality is the binary response variable with a code of 1 if a tree is dead or 0 if it is living. G is total stand basal area; N is number of trees; DQ is quadratic mean diameter; DG is the basal area-weighted mean diameter; HG is the basal area-weighted mean height, or Lorey’s height; T is stand age. Thinned refers to the thinning treatment as a dummy variable with a code of 1 for thinned plots and 0 for unthinned plots. Fig. 1. Locations of the silver birch experimental sites used for modelling and testing in this study. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 2 of thinnings (0–3) varied in the data. When the tree was identified as dead in a particular measurement occasion, it was also coded as a dead tree for the previous measurement to model forthcoming mortality. The number of living and dead trees was 88,797 and 2761 respectively, resulting a 3.0 % mortality rate for the 5-year growth period. A total of 1050 plot-level observations were applicable, and 76 % of the monitored plots were thinned during the survey period. Due to several (1–3) thinnings from below, 34 (3 %) study plots experienced no mortality. The plot was recorded as thinned after the first thinning. The average and the maximum time since the last thinning were 7 and 35 years, respectively. A total of 262 plot-level observations represented unthinned control plots. Since the experi ments were established for relatively young stands, the variation in stand characteristics was great; for example, stand basal area (G) varied from 2.6 to 37.6 m2 ha− 1. In the first measurement occasion, the maximum stem number (Nmax) was 3072 trees ha− 1, while the lowest density in the last period was only 186 trees ha− 1 (Table 1). According to the maximum stand age (T) of 62 years, the basal area- weighted mean diameter (DG) of 28 cm and Lorey’s height (HG) of 31 m, the oldest stand was at a stage of clearcutting at the end of the survey. The stands were aged between 46 and 67 years at the end of the last measurement period. The mean stand characteristics were very similar across site classes: OMaT had 46,919 observations, OMT had 19,518 observations, and MT had 25,122 observations. The minor discrepancy was found in stand age, which increased from OMaT (27 years) to OMT (30 years) and MT (34 years). Similarly, HG increased from 16.3 m to 17.3 m and 18.8 m, respectively. The obvious difference was found in mortality rates, which were, on average, 4.2 %, 2.8 %, and 1.0 % for OMaT, OMT and MT site, respectively. Unexpectedly, the maximum observed stand basal areas — 34.3, 35.7, and 37.6 m2 ha− 1 — decreased with increasing site fertility. 2.1.2. Test dataset The spacing trial (Niemistö, 1995), established in the beginning of 1970s on fertile, mineral-soil-based abandoned agricultural land, was used as an independent test dataset (Table 2). These trials were located in Central Finland and Ostrobothnia in the municipalities of Varkaus, Suonenjoki, and Muhos (Fig. 1). When survival models were evaluated in MOTTI, unthinned stands were selected and divided into three density groups: initially dense stands (N > 2000 trees ha− 1), normal density stands (1000 trees ha− 1 ≤ N ≤ 2000 trees ha− 1), and sparse stands (N < 1000 trees ha− 1). In the test dataset, collected from three sites, the number of dense, normal, and sparse stands was 10, 9, and 10, respec tively. The spacing trial included initially dense stands with a density of N = 5000 trees ha− 1, which exceeded the density in the modelling dataset. Trials also included initially very sparse plots (N = 500 trees ha− 1), which were not included in the modelling dataset. Even the extraordinary sparse plots exhibited some morality during the 25-year monitoring period unless the plots were thinned. Generally, the average tree and stand characteristics were similar to those in the modelling dataset (Tables 1 and 2). However, the average mortality rate was much higher at 15 % in the test dataset, compared to 3 % in the modelling dataset. The cause of mortality was unknown. A significant amount of mortality occurred between establishment and the first measurement (Niemistö, 1995), but it was not considered in this study. 2.2. Modelling approach and model structure 2.2.1. Complementary log-log regression model For a single tree, death or survival can be represented as a binary variable that takes the value of 1 if the tree dies or 0 if the tree survives over a given time interval (Fortin et al., 2008). According to Rose et al. (2006), the analysis of survival data in forestry is hindered because the data are interval-censored; the exact time of death is unknown, but the time interval during which death occurs is known (Rose et al., 2006). Methods for considering interval-censored data are described in Cox and Oakes (1984). The probability of tree mortality was expressed as an exponential function in which the period length (L) is included in the exponent, analogous to an interest rate (Fortin et al., 2008): p = 1 − exp( − exp(Xʹb))L/5 (1) where p is the probability of tree mortality, X′b is the linear model, which consists of the transposed vector of explaining variables X and the vector of estimated parameters b. The growing period length (L) was typically 5 years, although there was some variation, specifically be tween 4 and 6 years. For example, in the MOTTI simulator, the maximum and typical period length is 5 years (Hynynen et al., 2002; Salminen et al., 2005; Hynynen et al., 2014). If a 6-year period is to be simulated, this is done by first simulating 5 years (5/5) using Eq. (2), followed by simulating 1 year (1/5) separately with the same equation. The mortality model was fitted using the complementary log-log (CLL) link from Eq. (1) (e.g., Fortin et al., 2008; Penman and Johnson, 2009; Myllymäki et al., 2024) which is: CLL(p) = ln( − ln(1 − p)) = Xʹb + ln(L /5) (2) The variation in the period length (L/5) was an offset variable (i.e., the coefficient is one) in the linear model (Eq. (2)) because it functions similarly to an interest rate, as illustrated in Eq. (1). While Yang and Huang (2013) used the logit link function, the offset period length ln(L) does not follow the theory of period length as an exponent (interest rate) (Monserud and Sterba, 1999; Siipilehto and Mäkinen, 2025); rather, it functions as a multiplier, as demonstrated in Eq. (8) of Yang and Huang (2013). Note that modelling mortality (‘failures’) with the CLL link corresponds to modelling survival (‘successes’) using the log-log link function. In that case, the offset variable must be the negative logarithm of the period length, − ln(L) (Siipilehto and Mäkinen, 2025). The CLL link is particularly suitable for mortality analysis where the period length (L) varies, as shown in Eq. (2). 2.2.2. Random effects consideration Until recently, mortality models have most often been estimated as marginal models representing the average mortality of the entire dataset Table 2 Summary statistics of the tree, plot, and site-level characteristics in the test dataset. Variable Mean SD Min Max n Tree-level characteristics ​ ​ ​ ​ ​ dbh, cm 12.7 4.8 0.9 35 32,411 BAL, m2 ha− 1 11.0 7.2 0 32.2 ​ Tree mortality 0.15 0.36 0 1 ​ plot-level characteristics ​ ​ ​ ​ ​ G, m2 ha− 1 17.9 6.8 2.0 35.6 411 N, trees ha− 1 1200 749 240 4467 ​ DQ, cm 15.3 4.3 5.5 26.6 ​ DG, cm 16.0 4.4 6.2 26.7 ​ HG, m 18.0 5.3 5.1 28.4 ​ T, year 34 10.6 11 52 ​ site-level characteristics ​ ​ ​ ​ ​ Plot area, m2 1089 515 485 3520 3 Temperature sum, dd ◦C 1116 79.1 1035 1193 ​ Measurement interval, year 4.9 0.5 4 5 ​ Thinned 0.50 0.36 0 1 ​ Time since last thinning, year 6.8 5.5 0 20 ​ Note: n is the number of observations; dbh is tree diameter at breast height at a height of 1.3 m above the ground; BAL is the sum of the basal areas in larger trees; Tree mortality is the binary response variable with a code of 1 if a tree is dead or 0 if it is living. G is total stand basal area; N is number of trees; DQ is quadratic mean diameter; DG is basal area-weighted mean diameter; HG is the basal area-weighted mean height, or Lorey’s height; T is stand age. Thinned refers to the thinning treatment as a dummy variable with a code of 1 for thinned plots and 0 for unthinned plots. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 3 (e.g., Monserud, 1976; Monserud and Sterba, 1999; Yao et al., 2001; Temesgen and Mitchell, 2005). Recently, mixed-effects logistic models have become more prevalent (e.g., Fortin et al., 2008; Yang and Huang, 2013; Cailleret et al., 2016; Siipilehto et al., 2021; Myllymäki et al., 2024). The mixed-effects logistic regression for forestry applications was first presented by Jutras et al. (2003). The idea is to describe stand-level dynamics instead of the general average over the entire dataset, which is regarded as a marginal model. The term “marginal model” is preferred in this study and is also referred to as a “population-average model”, in contrast to fixed-effect and mixed-effect models. This distinction is considered important because, in practice, mixed-effect models are often applied as fixed-effect models due to the random effects being unknown and unable to be directly predicted. After Jutras et al. (2003), the observed response is represented by yijk; the status of the i-th tree in the j-th plot in the k-th experimental area at the end of the growing period is binary (0 = living, 1 = dead). The distribution for the response is binomial, yijk ~Bin(pijk,1). The model with a random effect can be written in a general form as follows in Eq. (3): pijk = f ( Xʹ ijk b+ ujk ) (3) where pijk is the expected probability of the response for the i-th tree in the j-th plot in the k-th experimental area at the end of the growth period; the function f is defined as the complementary log-log link function; the vector X′ijk is the transposed vector of independent fixed variables; b is the vector of the estimated parameters; and ujk represents the random variables for the experimental area k, plot j nested in experimental area k, or year j nested in experimental area k, depending on the tested model type, with u~N(0,su 2). The generalized linear mixed model approach was used to address the hierarchical structure of the data. Random effects were associated with the intercept of the models (Eq. (3)). In this study, a random intercept model was chosen because it provides a more parsimonious and interpretable representation of the data (Bates et al., 2018). Addi tionally, given the limited number of sites, including random slopes could lead to overfitting and convergence issues, making the model unstable and difficult to interpret. In the datasets, stand-level charac teristics represented plots within experimental areas. Within each experimental area, the site conditions are as homogenous as possible. Between experimental areas, there is variation in site conditions and locations. Therefore, the experimental area is referred to simply as site in this article. Plot-level differences were modelled using fixed effects based on the variation in tree and stand characteristics, with a random site (RND_SITE) being used rather than a random plot nested in a site (RND_PLOT(SITE)) in the model. However, models with the random effect of plots nested within sites (RND_PLOT(SITE)) and years nested within sites (RND_YEAR(SITE)) were also evaluated because of highly improved fit statistics. The plot-level random effect can be justified by the long monitoring period that includes several plot-level measurement instances. The random year effect is justified due to differences between growth periods resulting from varying climate conditions across years. Thus, the random year may explain differences in irregular mortality. Note that all the estimated models were applied as fixed-effect models because the random effects are not known in practice, especially when the model is applied to a new forest stand (Myllymäki et al., 2024). 2.2.3. Model formulation The structure of the mortality model was developed with the help of marginal stepwise logistic regression. Tree size was considered using alternatives to dbh by transforming it into various forms, such as ln(dbh), 1/dbh, ̅̅̅̅̅̅̅̅ dbh √ , and dbh2 to identify the most stable and unbiased model performance while accounting for tree allometric relationships (see Monserud and Sterba, 1999; Jutras et al., 2003). In addition to tree size, competition between individual trees was assessed using the sum of basal areas in larger trees (BAL), as defined by Wykoff (1990), and the transformed variable BAL/dbh (Kiernan et al., 2009). Stand-level competition was considered in terms of total stand basal area (G) and the number of stems per hectare (N). The developmental status was described using stand age (T) and quadratic mean diameter (DQ), with DQ typically used when modelling maximum stand density (Reineke, 1933; Eid and Øyen, 2003; Monserud et al., 2004; Condés et al., 2017; Lee and Choi, 2019; Lee et al., 2025). Because the models are intended to be used in practice, modelling data also comprised thinned stands. The effect of thinnings on silver birch mortality was tested using thinning dummy variables THIN_0_5 and THIN_6_10 for the periods of 0 to 5 years and 6 to 10 years since the last thinning, respectively. Due to the site carrying capacity and differ ences in growth rates, variation in site type was taken into consideration (e.g., Eid and Tuhus, 2001; Yao et al., 2001; Jutras et al., 2003). Because the experiments were established mainly on former agricultural land (OMaT site) and the forest site type (OMT and MT sites) was under represented in the dataset, the forest site types OMT and MT were tested in the model as dummy variables. However, it was noticed that the mortality rate in MT became unrealistically low; therefore, the forest site types were combined using a common dummy variable, denoted as Forest. The final structure of the linear part is presented as follows in Eq. (4): Xʹb = f(BAL,BAL/dbh,THIN 0 5,THIN 6 10, Forest) (4) where X′b is the linear model, which consists of the transposed vector of explaining variables X and the vector of estimated parameters b; dbh is tree diameter (cm) at breast height at a height of 1.3 m above the ground; BAL is the sum of the basal areas in larger trees (m2 ha− 1); THIN_0_5 is a dummy variable for recent thinnings from 0 to 5 years after the last thinning; THIN_6_10 is a dummy variable for the period of 6 to 10 years after the last thinning; and Forest is a dummy variable for stands on forest sites. All the tree and stand variables represented values at the beginning of the growth period. 2.3. Validation of the developed models 2.3.1. Comparison in the MOTTI simulations Comparisons were also made with the existing MOTTI models (Hynynen et al., 2002) when simulating the stand development of the unthinned planted silver birch stands of the test dataset. The Finnish MOTTI stand simulator is presented by Salminen et al. (2005). The trees of a stand are expressed as tree lists where each sample tree represents a certain number of trees per hectare. The basic tree variables are tree species, dbh, height and the number of trees per ha. Tree lists are updated with dbh and height growth models (Hynynen et al., 2014). The numbers of trees per ha are updated with the predicted tree-level sur vival. When stand mortality was simulated by applying the probability model of individual tree mortality developed in this study, the following framework was implemented using the current MOTTI simulator. The tree list in the MOTTI simulation holds a certain number of trees per hectare in each dbh class; for example, 100 trees ha− 1. If the predicted mortality rate in subsequent years is 0.1 for that dbh class, then the amount of tree mortality is 0.1 × 100 trees ha− 1 = 10 trees ha− 1. Likewise, the number of living trees is 0.9 × 100 trees ha− 1 = 90 trees ha− 1. The current tree-level survival model in MOTTI was subsequently modified from the original presented in Hynynen et al. (2002) by refitting it using the non-linear mixed-effects logistic regression (PROC NLMIXED) in SAS version 9.3 (SAS Institute Inc., 2011). The estimated linear component for silver birch survival probability is presented as follows in Eq. (5): Xʹb = 8.5158 − 4.8386RDFL − 34.9224RDFL/dsh (5) J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 4 where RDFL is the relative density factor in larger trees, and dsh refers to stump height diameter (dsh = 2.0 + 1.25 × dbh) according to Laasase naho (1975). The relative density factor (RDF) is based on the estimated self-thinning line (Hynynen, 1993). For a more detailed description of the models and modelling dataset, see Hynynen et al. (2002). Both models exhibited a similar structure including tree-level competition (either BAL or RDFL) and tree-level competition divided by tree diameter (BAL/dbh). The alternative marginal and random- effects CLL mortality models (Eq. (2)) were fitted using PROC GLIM MIX in SAS version 9.4 (SAS Institute Inc., 2017). Using PROC GLIMMIX, a marginal (or population-average) model is fitted with the syntax “Random residual”, which provides the residual variance as an output term. Mixed-effect models, on the other hand, include terms such as random site or nested terms like random plot nested within site, which estimate the variance associated with the respective random effects. For clarity and further reference on using PROC GLIMMIX, the SAS code syntax is provided in the Appendix. The complementary log-log (CLL) link function for tree mortality was employed in the PROC GLIMMIX in SAS version 9.4 due to the varying period lengths L (see Eq. (2)). The PROC GLIMMIX offers flexibility in incorporating multiple random effects, which is advantageous for our analysis (SAS Institute Inc., 2017). Alternatively, a nonlinear logistic survival model could be used (e.g., Siipilehto et al., 2021); however, the PROC NLMIXED (SAS Institute Inc., 2011) allows only a single random effect, limiting its applicability in this context. All the models were evaluated using the fit statistics of − 2 log- likelihood (− 2 logLik) and Akaike Information Criterion (AIC), along with the significance of parameter estimates. When the predicted values from the MOTTI simulation were compared, additional evaluation sta tistics were calculated: Bias, Bias %, SD, RMSE, and RMSE %. Here, Bias is the mean of the observed values minus the predicted (MOTTI) values; Bias % is the Bias divided by the mean observed value and multiplied by 100 to express it as a percentage; SD is the standard deviation; RMSE is the root mean square error; and RMSE % is calculated by dividing the RMSE by the mean observed value and then multiplying by 100 to ex press the error as a percentage. 2.3.2. Self-thinning models as the stand-level constraints Self-thinning models of silver birch stands, used here as the stand- level constraints, can be calculated using the following equations: the self-thinning rule (STR) by Hynynen (1993) (Eq. (6)), Reineke’s STR applying the linear mixed-effects quantile regression models fitted at the 99 % quantile (Eq. (7)), and Nilson’s stand sparsity index (SSI) fitted with 1 % quantile regression (Eq. (8)), both based on Lee et al. (2025), as shown below: ln(Nmax ) = 14.198 − 2.218ln(DGMs) (6) ln(Nmax) = 11.9013 − 1.5848ln(DQ) (7) Nmax = (100/(0.1457DQ + 0.1112))2 (8) where DGMs is the basal area-median diameter (cm) at stump height, given by the formula DGMs = 2.0 + 1.25 × DGM; DQ is quadratic mean diameter (cm), given by the formula DQ = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ G/(qN) √ , where G is the stand basal area (m2 ha− 1) and q is the conversion factor q = π 4(10000) ≈ 7.854 × 10− 5. When calculating the maximum basal area (Gmax), the equation is Gmax = DQ2 × Nmax × q. To calculate Gmax for the STR model of Hynynen (1993), the average relationship between DGM and DQ is described by the equation DGM=0.855 + 0.987 × DQ. 3. Results 3.1. Estimated models for five-year tree mortality The estimated CLL regression models for mortality are shown in Table 3. High correlations could be avoided even when BAL and BAL/ dbh were included in the model. In particular, when the random plot effect was incorporated into the model, the two highest correlations between predictor variables (or random effects) were only 0.44 and 0.43. The random effects were associated with the intercept terms, affecting only the level of the mortality model. The estimated parame ters for the mixed models changed notably compared with those of the marginal model (Table 3). The absolute value of the estimated param eters for tree-level competition (BAL) was found to be higher in the mixed models than in the marginal model. An opposite pattern was observed for the two thinning dummy variables. Finally, the model with RND_YEAR(SITE) rendered the thinning dummy variable THIN_6_10 insignificant, leading to its exclusion from that model option. The model fit statistics for − 2 log-likelihood and AIC improved significantly from the marginal model to the RND_SITE model (Table 3). Further improvement was found using the RND_PLOT(SITE) effect and, finally, the RND_YEAR(SITE) effect (Table 3). Fig. 2 shows an example of the RND_YEAR(SITE) model predictions for a dense stand plot on an agricultural field (the initial stage of the plot no. 133 in Varkaus), assuming recent thinning that decreased mortality by about 10 %. If the stand is on forest land, the mortality is much lower (Fig. 2). The smallest tree had a dbh of 4.2 cm, and the predicted probability of mortality was 0.83. If the stand was recently thinned, the probability of mortality was 0.73. If, instead of agricultural land, the site were on forest land, the probability of mortality would be much lower at Table 3 Parameter estimates and fit statistics for the 5-year planted silver birch mortality models using linear marginal and mixed-effects complementary log-log regression approaches. The marginal model represents the model fitting without a random effect; the RND_SITE model includes a random effect of sites; the RND_PLOT(SITE) model includes a random effect of plots within sites; and the RND_YEAR(SITE) model includes a random effect of years within sites. Marginal RND_SITE RND_PLOT(SITE) RND_YEAR(SITE) Category Variable Estimate S.E. t-value Estimate S.E. t-value Estimate S.E. t-value Estimate S.E. t-value Fixed effects Intercept − 5.1036 0.0580 − 87.97 − 5.4573 0.1800 − 30.33 − 5.8872 0.1365 − 43.14 − 5.8032 0.1373 − 42.26 ​ BAL 0.0940 0.0031 30.80 0.1125 0.0036 31.05 0.1328 0.0042 31.87 0.1119 0.0050 22.12 ​ BAL/dbh 0.5654 0.0142 39.79 0.5710 0.0155 36.83 0.5759 0.0164 35.09 0.6172 0.0184 33.42 ​ THIN_0_5 − 0.6861 0.0712 − 9.64 − 0.6269 0.0722 − 8.68 − 0.7510 0.0944 − 7.96 − 0.2938 0.0783 − 3.75 ​ THIN_6_10 − 0.6636 0.1191 − 5.57 − 0.5954 0.1208 − 4.93 − 0.6950 0.1326 − 5.24 ​ ​ ​ ​ Forest − 1.2042 0.0431 − 27.97 − 1.3779 0.2649 − 5.20 − 1.2394 0.1773 − 6.99 − 1.3368 0.1769 − 7.56 Random parts Variance 0.9764 ​ ​ 0.2774 0.1000 ​ 0.9753 0.1478 ​ 0.6914 0.1146 ​ Fit statistics − 2logLik 18,285 ​ ​ 17,804 ​ ​ 17,373 ​ ​ 17,059 ​ ​ ​ AIC 18,297 ​ ​ 17,818 ​ ​ 17,387 ​ ​ 17,071 ​ ​ Note: BAL represents the sum of the basal areas in larger trees (m2 ha− 1); dbh refers to tree diameter (cm) at breast height, measured at a height of 1.3 m above the ground; THIN_0_5 and THIN_6_10 are the dummy variables indicating the periods of 0 to 5 years and 6 to 10 years since the last thinning, respectively; Forest is a dummy variable denoting a forest site as opposed to former agricultural land; AIC refers to the Akaike Information Criterion; − 2logLik denotes the − 2 × log-likelihood value. Note that the marginal model included random residual for the population-average model, and the PROC GLIMMIX syntax in SAS reported a residual variance of 0.9764 for it. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 5 0.37. For the largest tree, with a dbh of 11.8 cm, the predicted mortality was practically zero (0.0009). 3.2. Model evaluation The advantage of the marginal model is that it is unbiased over the whole modelling dataset. If representative data (e.g., National Forest Inventory data) are available for modelling, the marginal model is a highly suitable option that can be utilised in practice. The mixed effects model using the RND_SITE effect underestimated mortality by 0.12 %, while the model using the RND_PLOT(SITE) effect underestimated mortality only marginally, at 0.05 %. The model using the RND_YEAR (SITE) effect resulted in the highest underestimation for mortality at 0.56 % in the modelling dataset. In the spacing trial test dataset, all the models underestimated mortality by about 9–10 %, with the least un derestimation occurring in the RND_SITE model. All the model evalua tions were performed using only the fixed effects. All the models underestimated mortality in the smallest 2.5 cm dbh class in the modelling data (Fig. 3a). The random effect models under estimated mortality in this class more than the marginal model, and the model with the RND_YEAR(SITE) effect underestimated it the most (Fig. 3a). In the next 7.5 cm class, the RND_PLOT(SITE) effect provided almost correct mortality (0.056 vs 0.058) while the other models underestimated mortality, with the RND_YEAR(SITE) model under estimating it the most (0.045) (Fig. 3a). For the 12.5 cm class, mortality was generally slightly overestimated; however, it was accurate with the RND_YEAR(SITE) model (Fig. 3a). For the 17.5 cm class, all the models provided quite accurate prediction for mortality while the RND_PLOT (SITE) model was the best. For the largest two dbh classes, all the pre dictions were underestimations: least by the marginal model and the RND_YEAR(SITE) model (Fig. 3a). In general, the mortality rate in test dataset was much higher (0.15) than in the modelling dataset (0.03). Thus, each model underestimated the mortality in the test dataset in all but the smallest 2.5 cm dbh class (Fig. 3b). The observed mortality rate was unexpectedly lower in the smallest dbh class (0.10) than in the subsequent 7.5 and 12.5 classes, which had rates of 0.169 and 0.195, respectively (Fig. 3b). In general, the marginal model provided the best predictions for the test dataset (Fig. 3b). Two exceptions were the 2.5 and 7.5 cm dbh classes, where the RND_YEAR(SITE) model predicted the lowest mortality and the RND_PLOT(SITE) model predicted the highest mortality, respectively; therefore, these models were closest to the observed mortality (Fig. 3b). Generally, the models performed well against BAL. All models except the RND_PLOT(SITE) effect provided similar results until the last two BAL classes, 27.5 and 32.5 m2 ha− 1 (Fig. 4a). The RND_PLOT(SITE) model slightly overestimated the mortality rate (0.23 vs 0.22) in the BAL class of 27.5 m2 ha− 1, whereas the other models slightly underestimated it (0.16–0.20) (Fig. 4a). For the last BAL class of 32.5 m2 ha− 1, all the models overestimated mortality in both datasets (Fig. 4, plots a and b). Clearly, the highest overestimations were detected in the RND_PLOT (SITE) model, specifically 0.421 and 0.800 in the modelling and test datasets, respectively. The RND_YEAR(SITE) model resulted in pre dictions that were close to those of marginal model for high BAL classes; these were the best predictions in both datasets (Fig. 4). 3.3. Evaluation of the models in the MOTTI simulator To evaluate the models in practice, a total of 29 untinned stands were simulated from the test dataset using the MOTTI simulator to compare the existing survival model with the new mortality models. The output variables of N and G were followed in order to characterise and examine stand density. The differences in simulated DG were so marginal Fig. 2. Probability of silver birch tree mortality by tree diameter at breast height (dbh, cm) according to the RND_YEAR(SITE) model. An example of a dense stand from Plot no. 133 in Varkaus was chosen for demonstration pur poses. The selected sample plot was classified as a dense stand with a stem density (N) of 4380 trees ha− 1, a stand basal area (G) of 24.7 m2 ha− 1, and a basal area-weighted mean diameter (DG) of 9.0 cm. Applications of the dummy variables THIN_0_5 and Forest are depicted by line type. THIN_0_5 denotes 1 if the period is between 0 and 5 years since the last thinning, or 0 if not; Forest denotes 1 if the forest site type is classified as OMT or MT according to Cajander (1926), or 0 if not. Fig. 3. Probability of mortality against tree diameter at breast height (dbh, cm) in the modelling dataset (a) and the test dataset (b), categorized by observed probability and model type: marginal, RND_SITE, RND_PLOT(SITE), and RND_YEAR(SITE). The marginal model represents a model fit without the random effect, while the mixed models represent the model fits with a random site effect, a random plot nested within site, and a random year nested within site. Only fixed-effect terms from each model, as shown in Table 3, were applied in this evaluation. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 6 between model options that it was excluded from the detailed evalua tion. In general, DG was overestimated by an average of 11.1 % in initially dense stands (15.8 % with the current survival model in the MOTTI simulator as described in Eq. (5)), underestimated by 0.4 % in initially normal density stands (1.7 % with the current MOTTI model), and overestimated by 2.3 % in initially sparse stands (2.7 % with the current MOTTI model). The initial state is always correct in MOTTI; thus, it was excluded from the calculations of evaluation characteristics. The simulation results for dense stands indicated that the greater the underestimation in stem number, the smaller the overestimation in stand basal area (Table 4). Notably, a 6.5 % underestimation in N was connected to the highest 7.4 % overestimation in G, while a 9.3 % un derestimation in N resulted in the smallest 4.2 % overestimation in G. The marginal model provided the lowest standard deviation (SD) of 324 trees ha− 1, and the RND_YEAR(SITE) model provided the best RMSE of 17.1 % in N (Table 4). The smallest prediction error for G was provided by the model with the RND_SITE effect, but a smaller standard deviation (SD) and RMSE were found using the model with the RND_PLOT(SITE) effect (Table 4). These two models provided less than 5 % over estimation of G for the overstocked stands in the MOTTI simulation. The existing MOTTI model (Eq. (5)) underestimated G by 4.5 % and N by 23 % (Table 4). Despite the similar absolute bias in G, the SD of 2.35 m2 ha− 1 and RMSE of 11.2 % were even smaller with the current MOTTI model compared with the new models. The evaluation of stands with normal initial density showed an un biased prediction of N using the RND_PLOT(SITE) model (Table 5). The marginal and the RND_SITE models resulted in slight underestimations of N (0.6 %–0.8 %), while the RND_YEAR(SITE) model overestimated it by 0.4 % (Table 5). Using the RND_YEAR(SITE) model provided the least biased estimate of G (1.6 %), while the other models underestimated G by 1.7 % to 2.6 %. The smallest SD and RMSE for G were achieved by the RND_PLOT(SITE) model (Table 5). Here, the RND_YEAR(SITE) model provided results that were remarkably close to those of the RND_PLOT (SITE) model. The existing RDFL-based MOTTI model underestimated N and G by 5.4 % and 5.1 %, respectively. The SD and RMSE in G were similar to those of the new models, but they were approximately twice as much in N. The RND_YEAR(SITE) model had considerably lower SD and RMSE for N in the initially dense and normal density stands (Tables 4 and 5). This was due to the fact that the RND_YEAR(SITE) model provided the most accurate predictions for the last simulation steps, while the other models showed that N dropped substantially, resulting in considerable underestimations. This is shown in Fig. 5, which is an example of the observed and predicted N in the MOTTI simulation for Plot no. 6 in Suonenjoki. The predicted N was almost the same for the first two 5-year growth steps. However, after the third growth step, the RND_YEAR (SITE) model provided a higher N, which was closest to the observed N after the fourth growth step (observed N vs. predicted N: 1125 vs. 1126 trees ha− 1) and the fifth growth step (observed N vs. predicted N: 1058 vs. 989). The lowest predicted N was 905 from the RND_PLOT(SITE) model in the fifth growth step (Fig. 5). The current MOTTI model pro vided accurate N for the first growth step (observed N vs. predicted N: 1475 vs. 1439 trees ha− 1), but thereafter, the underestimation was obvious, and the predicted N was only 766 trees ha− 1 by the last growth step. The evaluation of the new models on the initially sparse stands Fig. 4. Probability of mortality against the sum of the basal areas in larger trees in a plot (BAL, m2 ha− 1) in the modelling dataset (a) and the test dataset (b), categorized by observed probability and model type: marginal, RND_SITE, RND_PLOT(SITE), and RND_YEAR(SITE). The marginal model represents the model fit without a random effect, while the mixed models represent the model fits with a random site effect, a random plot nested within site, and a random year nested within site. Only fixed-effect terms from each model, as shown in Table 3, were applied in this evaluation. Table 4 Evaluation statistics for dense stands (N > 2000 trees ha− 1) in the test dataset after applying existing and new MOTTI simulations equipped with alternative models: marginal, RND_SITE, RND_PLOT(SITE), and RND_YEAR(SITE). The best evaluation characteristics among the new models and the current MOTTI model are high lighted in bold. The marginal model represents the model fit without a random effect, while the mixed models represent the model fits with a random site effect, a random plot nested within site, and a random year nested within site. MOTTI represents the current version of a MOTTI simulation (Salminen et al., 2005). N is the number of stems per hectare (trees ha− 1) and G is the stand basal area (m2 ha− 1). Model Marginal RND_SITE RND_PLOT(SITE) RND_YEAR(SITE) MOTTI N G N G N G N G N G Bias 186.44 − 1.40 213.92 ¡1.07 231.35 − 1.14 148.92 − 1.90 533.06 1.15 Bias % 8.09 − 5.45 9.28 ¡4.16 10.04 − 4.44 6.46 − 7.40 23.13 4.46 SD 323.80 2.72 345.72 2.88 364.25 2.62 324.61 2.80 490.41 2.35 RMSE 412.37 3.37 448.80 3.39 476.40 3.15 393.90 3.73 801.45 2.88 RMSE % 17.89 13.12 19.47 13.17 20.67 12.27 17.09 14.53 34.78 11.21 Note: Bias is the mean of the observed values minus the predicted (MOTTI) values; Bias % is the Bias divided by the mean observed value and multiplied by 100; SD is the standard deviation; RMSE is the root mean square error; and RMSE % is calculated by dividing the RMSE by the mean observed value and then multiplying by 100 to express the error as a percentage. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 7 showed the best performance for G using the marginal model (Table 6). While the marginal model overestimated G by only 1.7 %, the prediction errors in G ranged from 2.0 % to 2.7 % overestimation by the mixed- effects models. Simultaneously, the bias in N was only one tree, which corresponds to − 0.14 % using the RND_YEAR(SITE) model, and ranged from 0.2 % to 0.8 % with the other models. The existing RDFL-based MOTTI model resulted in prediction biases of 2.6 % for N and − 0.9 % for G, respectively (Table 6). Thus, for the sparse stands, the existing MOTTI model showed a smaller bias in G than any of the new BAL-based mortality models. However, the SD and RMSE in G were slightly lower using the new models. All in all, the differences in the simulation results were mostly slight between the marginal and mixed-effect model options. Among the new models, the RND_YEAR(SITE) model demonstrated the best perfor mance in dense stands, as measured by RMSE for the development of N; however, it simultaneously resulted in the highest overestimation of 7.4 % in G (Table 4). In normal density stands, it showed the lowest bias in G and the smallest SD and RMSE in N (Table 5). Finally, in sparce stands, it exhibited the lowest bias and RMSE for the development of N in the MOTTI simulation (Table 6). In the initially sparse stand, the differences between the new models in their development of G were only slight; however, the marginal model provided the best performance. The existing RDFL-based MOTTI model for silver birch survival performed well in the development of G, particularly in dense and sparse stands. However, in normal density stands, all the new models provided smaller bias than the current MOTTI model. Overall, the underestimation of N was considerable when using the current MOTTI model. The average and maximum G observations for the last measurements of the dense stands in the test dataset were 31.6 m2 ha− 1 and 35.6 m2 ha− 1 (Table 7). The maximum G during long-term MOTTI simulations was searched for each plot and for the alternative mortality model in initially dense and normal stands. The marginal model resulted in a slightly lower mean for maximum G of 31.0 m2 ha− 1, while the RND_SITE and RND_PLOT(SITE) models provided even lower maximum G values of 30.6 m2 ha− 1 and 30.4 m2 ha− 1, respectively. In contrast, the RND_YEAR(SITE) model resulted in the highest G of 32.0 m2 ha− 1 (Table 7). None of the maximum G values were regarded as unrealisti cally high because they didn’t reach the maximum observed value of 35.6 m2 ha− 1. Note that the maximum G of 28.1 m2 ha− 1 using the Table 5 Evaluation statistics in normal density stands (1000 ≤ N ≤ 2000 trees ha− 1) of the test dataset after applying existing and new MOTTI simulations equipped with alternative models: marginal, RND_SITE, RND_PLOT(SITE), and RND_YEAR(SITE). The best evaluation characteristics among the new models and the current MOTTI model are highlighted in bold. The marginal model represents a model fit without a random effect, while the mixed models represent model fits with a random site effect, a random plot nested within site, and a random year nested within site. MOTTI represents the current version of a MOTTI simulation (Salminen et al., 2005). N is the number of stems per hectare (trees ha− 1) and G is the stand basal area (m2 ha− 1). Model Marginal RND_SITE RND_PLOT(SITE) RND_YEAR(SITE) MOTTI N G N G N G N G N G Bias 8.90 0.52 6.43 0.47 ¡0.79 0.34 − 5.04 0.32 62.32 1.02 Bias % 0.77 2.59 0.56 2.38 ¡0.07 1.73 − 0.44 1.62 5.40 5.14 SD 53.02 2.45 55.92 2.43 59.76 2.43 48.91 2.43 98.21 2.24 RMSE 57.65 2.69 60.34 2.66 64.06 2.63 52.71 2.63 125.00 2.64 RMSE % 4.99 13.50 5.22 13.34 5.55 13.19 4.56 13.20 10.82 13.27 Note: Bias is the mean of the observed values minus the predicted (MOTTI) values; Bias % is the Bias divided by the mean observed value and multiplied by 100; SD is the standard deviation; RMSE is the root mean square error; and RMSE % is calculated by dividing the RMSE by the mean observed value and then multiplying by 100 to express the error as a percentage. Fig. 5. Development of observed and predicted stem number (N, trees ha− 1) for Plot no. 6 in Suonenjoki, with stem number of 1558 trees ha− 1, stand basal area (G) of 15.6 m2 ha− 1, and basal area-weighted mean diameter (DG) of 12.7 cm, in the MOTTI simulation using alternative mortality models over five 5-year growth steps. The marginal model represents the model fitting without a random effect; the RND_SITE model includes a random effect of sites; the RND_PLOT(SITE) model includes a random effect of plots within sites; and the RND_YEAR(SITE) model includes a random effect of years within sites. Table 6 Evaluation statistics in sparse stands (N < 1000 trees ha− 1) of the test dataset after applying existing and new MOTTI simulations equipped with alternative models: marginal, RND_SITE, RND_PLOT(SITE), and RND_YEAR(SITE). The best evaluation characteristics among the new models and the current MOTTI model are high lighted in bold. The marginal model represents model fitting without a random effect, while the mixed models represent model fits with a random site effect, a random plot nested within site, and a random year nested within site. MOTTI represents the current version of a MOTTI simulation (Salminen et al., 2005). N is the number of stems per hectare (trees ha− 1) and G is the stand basal area (m2 ha− 1). Model Marginal RND_SITE RND_PLOT(SITE) RND_YEAR(SITE) MOTTI N G N G N G N G N G Bias 5.98 ¡0.30 3.32 − 0.36 − 1.77 − 0.47 ¡1.04 − 0.44 19.47 ¡0.15 Bias % 0.79 ¡1.71 0.44 − 2.03 − 0.23 − 2.68 ¡0.14 − 2.51 2.59 ¡0.86 SD 26.60 2.27 27.29 2.27 27.22 2.28 26.65 2.27 46.30 2.34 RMSE 29.88 2.50 30.11 2.52 29.89 2.55 29.21 2.53 55.09 2.57 RMSE % 3.97 14.24 4.00 14.31 3.97 14.53 3.88 14.41 7.32 14.61 Note: Bias is the mean of the observed values minus the predicted (MOTTI) values; Bias % is the Bias divided by the mean observed value and multiplied by 100; SD is the standard deviation; RMSE is the root mean square error; and RMSE % is calculated by dividing the RMSE by the mean observed value and then multiplying by 100 to express the error as a percentage. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 8 existing MOTTI model did not reach the average basal area from the last measurements, even though those observations might not represent the absolute maximum. 3.4. Evaluation of tree-level mortality models against self-thinning lines Presented models were compared to the existing self-thinning line for silver birch by Hynynen (1993) and to the new self-thinning models for planted silver birch (Lee et al., 2025). This examination was con ducted using the MOTTI stand simulator to simulate the development of unthinned stands of our test dataset as an example (Fig. 6). In addition to the observations and the MOTTI simulations, the applied stand-level complementary constraints applied were based on the self-thinning models calculated by Eqs. (6–8). The RND_YEAR(SITE) model was selected due to its generally good performance and the highest G values achieved in the MOTTI simulations. The idea was to test if the individual tree mortality models can be used without the self-thinning constraint (see Monserud et al., 2004). The ideal situation would be for the pre dicted Reineke’s STR and Nilson’s SSI, as presented by Lee et al. (2025), to coincide with the development of stand density using the new tree-level mortality models presented by this current study. Lee et al. (2025) concluded that the most relevant maximum density for planted silver birch was a combination of Reineke’s STR for lower DQ and Nil son’s SSI for higher DQ. Here, the models for maximum stem number (Nmax) are expressed as maximum basal area (Gmax) using the formula Gmax = DQ2 × Nmax × q, where q = p/2002. In terms of maximum density, the differences in Gmax are easier to observe than in highly variable Nmax. Results in Fig. 6 showed compatibility between stand- and tree-level mortality models. The RND_YEAR(SITE) model represented the tree- level mortality in the MOTTI simulation. Whether the initial density was dense (4000 trees ha− 1), normal (1600 trees ha− 1), or sparse (900 trees ha− 1), the simulated G did not reach the maximum density rep resented by Reineke’s STR fitted using 99 % or Nilson’s SSI fitted using 1 % quantile regression (Lee et al., 2025). The simulation of initially high density resulted in the maximum density that was closest to Reineke’s STR when DQ was 13 cm, and thereafter, it was slightly below Nilson’s SSI at higher DQ values (Fig. 6). Overall, the new tree-level mortality models can be used without a self-thinning constraint because the simulated G did not reach the modelled stand-level Gmax. The existing STR by Hynynen (1993) showed too low Gmax (< 30 m2 ha− 1) for the planted silver birch. According to example plots, the observed G reached about 30 m2 ha− 1, but it is expected to continue increasing in the future (Fig. 6). 4. Discussion 4.1. Overall performance of the developed models When constructing a mortality model, several alternative model structures were evaluated. The most relevant variables representing tree vigor are the crown ratio (Monserud and Sterba, 1999) and diameter growth (Monserud, 1976; Yao et al., 2001). Those are preferred not to be used because when the model is applied in a forest decision support system for growth and yield simulation, such as a MOTTI simulator, it would rely entirely on model-based data with no measured observa tions, even in the initial stage. Instead, when using tree dbh and the stand characteristics based on it, the initial state can be accurately described using a tree list. The final model of this study was rather simple, comprising only the key variables of competition (BAL) and tree dbh, along with the thinning effects (THIN_0_5 and THIN_6_10) and the Forest dummy variable. Thus, there was no concern of over-parameterisation or high correlations between estimated parameters, and the resulting model was robust and stable (Monserud and Sterba, 1999; Jutras et al., 2003). In addition, the model structure used behaved such that during the long simulation run, the Gmax reached by the model in this study remained for an extended period, staying slightly below the modelled Gmax as reported by Lee et al. (2025). If additional stand characteristics were included (e.g., DQ, G or N), this behavior was lost; instead, the density of a stand in terms of G collapsed more rapidly after reaching its maximum value. In general, the data did not support this kind of behaviour because there were only a few such observations. Otherwise, the observed G was continuously increasing, similar to the values shown in the example (Fig. 6). Mixed effects in the generalized linear models improved the fit sta tistics considerably. The improvement in the fit statistics did not demonstrate similar behaviour when tested in the MOTTI simulator. Still, the RND_YEAR(SITE) model provided the best fit statistics and generally performed well in MOTTI simulations. The differences be tween the models were rather negligible in practice because all the mixed-effects models had to be treated as fixed-effect models, at least when used with the independent test data (Myllymäki et al., 2024). Table 7 The maximum stand basal area (Gmax, m2 ha− 1) at the last measurement observation of dense stands and the absolute maximum stand basal area reached during the long-term MOTTI simulation. This includes results from new mar ginal and mixed-effects models incorporating random effects of RND_SITE, RND_PLOT(SITE), or RND_YEAR(SITE), as well as from the current RDFL-based survival model in MOTTI. The marginal model represents a model fit without a random effect, while the mixed models represent model fits with a random site effect, a random plot nested within site, and a random year nested within site. MOTTI represents the current version of a MOTTI simulation (Salminen et al., 2005). RDFL is the relative density factor in larger trees. Category Mean Min Max Last observation 31.6 28.5 35.6 Marginal 31.0 28.2 34.1 RND_SITE 30.6 28.0 33.4 RND_PLOT(SITE) 30.4 28.0 33.0 RND_YEAR(SITE) 32.0 29.5 34.6 MOTTI 28.1 26.9 29.1 Fig. 6. Comparison of stand basal area between preceding self-thinning models from other studies, the tree-level mortality model from this study, and obser vations. The preceding models for maximum basal area (Gmax) were based on self-thinning lines according to Reineke’s self-thinning rule (STR) and Nilson’s stand sparsity index (SSI), both as described by Lee et al. (2025), along with the existing STR by Hynynen (1993). The stand basal area from the MOTTI simu lation (—●—), using the RND_YEAR(SITE) model for the individual tree mor tality, in initially dense, normal, and sparse stands of the test dataset, is presented alongside the observations (-▴-). The RND_YEAR(SITE) model is the tree-level mortality model fitted in this study, incorporating random effects of years nested within sites. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 9 4.2. Mixed-effects modelling and thinning effects on tree mortality All the fitted mortality models (marginal and mixed effects models) performed well. The independent test dataset showed the overall best performance (in terms of N and G) using the RND_YEAR(SITE) model. In addition, the RND_YEAR(SITE) model performed the best in terms of G for initially normal density stands. The accuracy in G is much more meaningful compared to N when estimating stand volume characteris tics (total volume and assortment volume). When the previous marginal model for silver birch in MOTTI (Hynynen et al., 2002) was refitted using mixed-effects logistic modelling, the performance improved significantly. This was shown, for example, by making the stand-level constraint (self-thinning) unnecessary. In this study, the model fit sta tistics improved significantly from the marginal model to the alternative mixed-effects models. However, only the fixed part of the models can be used when predicting mortality for trees in a new stand. Lastly, despite a much better fit (− 2logLik and AIC), some of the mixed-effects models performed worse than the marginal model. As the RND_YEAR(SITE) model showed the best performance, this suggests that natural mortality varies across both time and space, and is not solely driven by competition or stand dynamics. It could therefore be argued that, for broader applicability, proxies for year (and potentially site) should be included as fixed effects rather than random effects, since year may reflect climatic variation during the observation period. Given that climate data (e.g., precipitation, temperature) are increasingly available and expected to become more important under climate change, omitting them as fixed effects could reduce model performance in future predictions. However, treating year as a fixed effect is not appropriate in the context of our study, as it would not facilitate reliable prediction of mortality beyond the observed period. The variation be tween years is largely stochastic and cannot be adequately captured by simply including year as a categorical fixed effect. To account for tem poral variation more meaningfully, it would be necessary to incorporate relevant environmental covariates—such as maximum wind speed, precipitation, or temperature—that reflect climatic influences. With larger datasets and more comprehensive climate data, integrating these factors as fixed effects could improve model performance and predictive capacity, particularly under changing climate conditions. Unfortu nately, due to the relatively small size of our dataset, including such variables is not feasible at this stage. Planted silver birch responded positively to thinnings, as demon strated by increased survival rates over a period of 10 years following the last thinning from below (see Niemistö, 1997; Simard et al., 2004; Hynynen et al., 2010). Unlike in the model for Norway spruce (Siipilehto et al., 2021), the dummy for the most recent thinning (THIN_0_5) was significant for silver birch. The effect of thinnings (THIN_0_5 and THIN_6_10) for silver birch were greatest in the marginal model. In contrast, when applying the RND_YEAR(SITE) model, only THIN_0_5 remained significant, while THIN_6_10 became insignificant (Table 3). 4.3. Recommendations and limitations for field application The new mortality models presented here showed reasonable compatibility with the new self-thinning models by Lee et al. (2025). The recommended initial density of 1600 trees ha− 1 showed good per formance in the MOTTI simulation. The survival rates in the test dataset revealed some surprising results. For example, the highest mortality rate of 20 % occurred in the dbh class of 12.5 cm, which was clearly higher than the mortality rates in the two lowest dbh classes (Fig. 3b). Addi tionally, the initial spacing appeared to affect tree mortality. According to Niemistö (1995), using a row-to-row distance of 5 m, mortality was about twice that found for even spacing at the same density. Thus, there were other types or causes of mortality besides competition-induced regular mortality in the test dataset. Niemistö (1995) further noted that mortality in sparse stands occurred mainly at the seedling stage, and subsequent observations indicated that mortality continued despite sparse and uneven tree locations. The developed individual tree mortality models are applicable to silver birch plantations, taking into account the potential effects of thinnings and site type. This limitation arises because the behaviour of any Finnish mortality model for birch, including those developed in our study, is questionable in natural forest conditions. Regarding this, Sii pilehto et al. (2021) noted that in near-natural forests in Finland, the proportion of birch was small—approximately 4 % silver birch and 3 % downy birch—due to their light-demanding nature. The models were subsequently rigorously tested on agricultural land. However, due to the limited data availability on forest land and the absence of an indepen dent test dataset, direct application of these models to forest sites is not recommended. The significant reduction in mortality observed with the Forest dummy variable raises concerns about potential underestimation of mortality when applying the models to forest sites. The models were tested on forest land over only one 5-year period across six breeding study plots. In this case, the current MOTTI model provided the best estimate of N, with an overestimation of only 4 trees ha− 1, while the new models overestimated N by 16 to 19 trees ha− 1. Therefore, it is recom mended that the models be further tested once suitable data from forest sites become available. In the previous study, Hamilton (1990) suggested that incorporating stand-level complementary constraints such as Reineke’s STR with in dividual tree mortality models improved performance and extended the range of applicability. Monserud et al. (2004) concluded that stand-level constraints are not always necessary if the modelling dataset is suffi cient; however, they suggested the constraints should be imposed if the modelling dataset is inadequate. In the present study, it is concluded that stand-level self-thinning complementary constraints such as Rein eke’s STR and Nilson’s SSI are not mandatory based on the examination conducted (Fig. 6). 5. Conclusion In this study, individual tree mortality models were developed using a complementary log-log regression approach with mixed effects. Due to growing period length as an offset variable, these models are available for any period length, not just the subsequent five years that were typical in the modelling dataset. The finally selected key variables for fixed effects were dbh and BAL, which effectively represented tree size and competition. BAL explained the higher mortality probability associated with a severely competitive state characterised by high BAL. Small dbh increased the probability in the BAL/dbh term, as indicated by positive parameter estimates in all model types. THIN_0_5 and THIN_6_10 were selected, indicating that the thinning treatment had a converse impact on the probability of tree mortality; specifically, the probability decreased when thinning dummy value was 1. The dummy variable for site type, Forest, significantly indicated a negative parameter, suggesting that the probability of tree mortality decreased in forest sites compared to the probability in abandoned agricultural land. Among all the types of models developed in this study using different classes of random effects, the RND_YEAR(SITE) model demonstrated the best performance. Although there was some variation when the RND_YEAR(SITE) model was tested and validated, it is considered applicable in practice, given that the best fit statistics and the reliably predicted maximum G. However, considering the relatively limited test dataset examined, it is prudent to use the model cautiously in forest sites; therefore, formal agricultural land would be more suitable and accurate for applying the developed model. Lastly, the individual tree mortality model demonstrated its compatibility with the recently developed self-thinning at the stand level, based on the concepts of Reineke’s STR and Nilson’s SSI. In contrast, the old self-thinning model of MOTTI was found to be incompatible with the model developed in this study. Consequently, it is concluded that the models developed in this study are more adequate for predicting the mortality probability of individual trees in the near subsequent years, especially in silver birch J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 10 plantations in abandoned fields in southern and central Finland. Funding information This research was funded by the Academy of Finland (Suomen Akatemia) under the UNITE (Forest-Human-Machine Interplay) flagship ecosystem (Grant No. 337655) and the ForClimate (Managing Forests for Climate Change Mitigation) project (Grant No. 347793); and by the Natural Resources Institute Finland (LUKE) through the strategically funded project ‘ProBirch - Possibilities and challenges of increasing birch with special attention to productivity and risk management’. CRediT authorship contribution statement Jouni Siipilehto: Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Formal analysis, Conceptualization. Daesung Lee: Writing – review & editing, Writing – original draft, Visualization, Validation, Data curation, Conceptualization. Declaration of competing interest The authors declare that they have no conflict of interest. Acknowledgements This study was carried out based on the empirical data provided by Natural Resources Institute Finland (LUKE) in Finland. The authors wish to thank all associated members for their invaluable efforts in fieldwork and data maintenance support. Appendix For the modelling of Eq. (2), the following PROC GLIMMIX syntax fits the marginal model in SAS version 9.4 (SAS Institute Inc., 2017). The mixed models are shown as alternatives and are commented out using *, as follows: PROC GLIMMIX Method=Laplace; CLASS site plot year; MODEL dead = BAL BAL_per_dbh THIN_0_5 THIN_6_10 Forest / dist=BINARY LINK=CLL OFFSET=lnL SOLUTION; RANDOM residual; *RANDOM site; *RANDOM plot(site); *RANDOM year(site); RUN; Data availability The data are stored at Natural Resources Institute Finland (LUKE) and are not publicly available due to their ownership. The data in this study are available on reasonable request by contacting the authors. References Bates, D., Kliegl, R., Vasishth, S., Baayen, H., 2018. Parsimonious mixed models. htt ps://doi.org/10.48550/arXiv.1506.04967. Boeck, A., Dieler, J., Biber, P., Pretzsch, H., Ankerst, D.P., 2014. Predicting tree mortality for European beech in southern Germany using spatially explicit competition indices. For. Sci. 60, 613–622. https://doi.org/10.5849/forsci.12-133. Cailleret, M., Bigler, C., Bugmann, H., Camarero, J.J., Cufar, K., Davi, H., Mészáros, I., Minunno, F., Peltoniemi, M., Robert, E.M.R., Suarez, M.L., Tognetti, R., Martínez- Vilalta, J., 2016. Towards a common methodology for developing logistic tree mortality models based on ring-width data. Ecol. Appl. 26, 1827–1841. https://doi. org/10.1890/15-1402.1. Cajander, A.K., 1926. The theory of forest types. Acta For. Fenn. 29, 108. Condés, S., Vallet, P., Bielak, K., Bravo-Oviedo, A., Coll, L., Ducey, M.J., Pach, M., Pretzsch, H., Sterba, H., Vayreda, J., Del Río, M., 2017. Climate influences on the maximum size-density relationship in Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) stands. For. Ecol. Manage. 385, 295–307. https://doi.org/ 10.1016/j.foreco.2016.10.059. Cox, D.R., Oakes, D., 1984. Analysis of Survival Data. Chapman and Hall/CRC, New York. https://doi.org/10.1201/9781315137438. Eid, T., Øyen, B., 2003. Models for prediction of mortality in even-aged forest. Scand. J. For. Res. 18, 64–77. https://doi.org/10.1080/02827581.2003.10383139. Eid, T., Tuhus, E., 2001. Models for individual tree mortality in Norway. For. Ecol. Manag. 154, 69–84. https://doi.org/10.1016/S0378-1127(00)00634-4. Fortin, M., Bédard, S., DeBlois, J., Meunier, S., 2008. Predicting individual tree mortality in northern hardwood stands under uneven-aged management in southern Québec. Canada. Ann. For. Sci. 65, 205. https://doi.org/10.1051/forest:2007088. Fridman, J., Ståhl, G., 2001. A three-step approach for modelling tree mortality in Swedish forests. Scand. J. For. Res. 16, 455–466. https://doi.org/10.1080/ 02827580152632856. Hamilton Jr., D.A., 1990. Extending the range of applicability of an individual tree mortality model. Can. J. For. Res. 20, 1212–1218. https://doi.org/10.1139/x90-160. Hynynen, J., 1993. Self-thinning models for even-aged stands of Pinus sylvestris, Picea abies and Betula pendula. Scand. J. For. Res. 8, 326–336. https://doi.org/10.1080/ 02827589309382781. Hynynen, J., Niemisto, P., Vihera-Aarnio, A., Brunner, A., Hein, S., Velling, P., 2010. Silviculture of birch (Betula pendula Roth and Betula pubescens Ehrh.) in northern Europe. Forestry 83, 103–119. https://doi.org/10.1093/forestry/cpp035. Hynynen, J., Ojansuu, R., Hökkä, H., Siipilehto, J., Salminen, H., Haapala, P., 2002. Models for predicting stand development in MELA System. In. Research Papers, 835. Finnish Forest Research Institute, p. 116. Hynynen, J., Salminen, H., Ahtikoski, A., Huuskonen, S., Ojansuu, R., Siipilehto, J., Lehtonen, M., Rummukainen, A., Kojola, S., Eerikäinen, K., 2014. Scenario analysis for the biomass supply potential and the future development of Finnish forest resources. Finnish Forest Research Institute. Res. Pap. 302, 106. Jutras, S., Hökkä, H., Alenius, V., Salminen, H., 2003. Modeling mortality of individual trees in drained peatland sites in Finland. Silva Fenn. 37, 504. https://doi.org/ 10.14214/sf.504. Kiernan, D., Bevilacqua, E., Nyland, R., Zhang, L., 2009. Modeling tree mortality in low- to medium-density uneven-aged hardwood stands under a selection system using generalized estimating equations. For. Sci. 55, 343–351. https://doi.org/10.1093/ forestscience/55.4.343. Laarmann, D., Korjus, H., Sims, A., Stanturf, J.A., Kiviste, A., Köster, K., 2009. Analysis of forest naturalness and tree mortality patterns in Estonia. In: For. Ecol. Manag., IUFRO Conference on Biodiversity in Forest Ecosystems and Landscapes, 258, pp. S187–S195. https://doi.org/10.1016/j.foreco.2009.07.014. Laasasenaho, J., 1975. Runkopuun saannon riippuvuus kannon korkeudesta ja latvan katkaisulpimitasta [Dependence of the amount of harvestable timber upon the stump height and the top-logging diameter] (In Finnish with English summary). Folia Forestalia 223, 1–20. Lee, D., Choi, J., 2019. Evaluating maximum stand density and size–density relationships based on the Competition Density Rule in Korean pines and Japanese larch. For. Ecol. Manag. 446, 204–213. https://doi.org/10.1016/j.foreco.2019.05.017. Lee, D., Siipilehto, J., Hynynen, J., 2025. Comparison and analysis of self-thinning models based on diameter-based maximum size-density relationships. For. Ecol. Manag. 575, 122374. https://doi.org/10.1016/j.foreco.2024.122374. Maleki, K., Kiviste, A., 2016. Individual tree mortality of silver birch (Betula pendula Roth) in Estonia. iForest 9, 643–651. https://doi.org/10.3832/ifor1672-008. Monserud, R.A., 1976. Simulation of forest tree mortality. For. Sci. 22, 438–444. https:// doi.org/10.1093/forestscience/22.4.438. Monserud, R.A., Ledermann, T., Sterba, H., 2004. Are self-thinning constraints needed in a tree-specific mortality model? For. Sci. 50, 848–858. https://doi.org/10.1093/ forestscience/50.6.848. Monserud, R.A., Sterba, H., 1999. Modeling individual tree mortality for Austrian forest species. For. Ecol. Manag. 113, 109–123. https://doi.org/10.1016/S0378-1127(98) 00419-8. Myllymäki, M., Kuronen, M., Bianchi, S., Pommerening, A., Mehtätalo, L., 2024. A Bayesian approach to projecting forest dynamics and related uncertainty: an application to continuous cover forests. Ecol. Model. 491, 110669. https://doi.org/ 10.1016/j.ecolmodel.2024.110669. Niemistö, P., 1997. Ensiharvennuksen ajankohdan ja voimakkuuden vaikutus istutetun rauduskoivikon kasvuun ja tuotokseen [Effect of the timing and intensity of the first thinning on growth and yield of a Betula pendula stand (in Finnish only)]. Metsätieteen Aikakauskirja 1997, 439–454. https://doi.org/10.14214/ma.6232. Niemistö, P., 1995. Influence of initial spacing and row-to-row distance on the growth and yield of silver birch (Betula pendula). Scand. J. For. Res. 10, 245–255. https:// doi.org/10.1080/02827589509382890. Penman, A.D., Johnson, W.D., 2009. Complementary log–log regression for the estimation of covariate-adjusted prevalence ratios in the analysis of data from cross- sectional studies. Biom. J. 51, 433–442. https://doi.org/10.1002/bimj.200800236. Reineke, L.H., 1933. Perfecting a stand-density index for even-aged forests. J. Agric. Res. 46, 627–638. Rose, C.E., Hall, D.B., Shiver, B.D., Clutter, M.L., Borders, B., 2006. A multilevel approach to individual tree survival prediction. For. Sci. 52, 31–43. https://doi.org/ 10.1093/forestscience/52.1.31. Salminen, H., Lehtonen, M., Hynynen, J., 2005. Reusing legacy FORTRAN in the MOTTI growth and yield simulator. Comput. Electron. Agric. 49, 103–113. https://doi.org/ 10.1016/j.compag.2005.02.005. SAS Institute Inc., 2017. SAS/STAT® 14.3 User’s Guide: The GLIMMIX Procedure. SAS Institute Inc., Cary, NC, USA. SAS Institute Inc., 2011. SAS/STAT® 9.3 User’s Guide: The NLMIXED Procedure. SAS Institute Inc., Cary, NC, USA. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 11 https://doi.org/10.48550/arXiv.1506.04967 https://doi.org/10.48550/arXiv.1506.04967 https://doi.org/10.5849/forsci.12-133 https://doi.org/10.1890/15-1402.1 https://doi.org/10.1890/15-1402.1 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0004 https://doi.org/10.1016/j.foreco.2016.10.059 https://doi.org/10.1016/j.foreco.2016.10.059 https://doi.org/10.1201/9781315137438 https://doi.org/10.1080/02827581.2003.10383139 https://doi.org/10.1016/S0378-1127(00)00634-4 https://doi.org/10.1051/forest:2007088 https://doi.org/10.1080/02827580152632856 https://doi.org/10.1080/02827580152632856 https://doi.org/10.1139/x90-160 https://doi.org/10.1080/02827589309382781 https://doi.org/10.1080/02827589309382781 https://doi.org/10.1093/forestry/cpp035 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0014 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0014 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0014 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0015 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0015 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0015 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0015 https://doi.org/10.14214/sf.504 https://doi.org/10.14214/sf.504 https://doi.org/10.1093/forestscience/55.4.343 https://doi.org/10.1093/forestscience/55.4.343 https://doi.org/10.1016/j.foreco.2009.07.014 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0019 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0019 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0019 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0019 https://doi.org/10.1016/j.foreco.2019.05.017 https://doi.org/10.1016/j.foreco.2024.122374 https://doi.org/10.3832/ifor1672-008 https://doi.org/10.1093/forestscience/22.4.438 https://doi.org/10.1093/forestscience/22.4.438 https://doi.org/10.1093/forestscience/50.6.848 https://doi.org/10.1093/forestscience/50.6.848 https://doi.org/10.1016/S0378-1127(98)00419-8 https://doi.org/10.1016/S0378-1127(98)00419-8 https://doi.org/10.1016/j.ecolmodel.2024.110669 https://doi.org/10.1016/j.ecolmodel.2024.110669 https://doi.org/10.14214/ma.6232 https://doi.org/10.1080/02827589509382890 https://doi.org/10.1080/02827589509382890 https://doi.org/10.1002/bimj.200800236 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0030 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0030 https://doi.org/10.1093/forestscience/52.1.31 https://doi.org/10.1093/forestscience/52.1.31 https://doi.org/10.1016/j.compag.2005.02.005 https://doi.org/10.1016/j.compag.2005.02.005 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0033 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0033 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0034 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0034 Siipilehto, J., Allen, M., Nilsson, U., Brunner, A., Huuskonen, S., Haikarainen, S., Subramanian, N., Antón-Fernández, C., Holmström, E., Andreassen, K., Hynynen, J., 2020. Stand-level mortality models for Nordic boreal forests. Silva Fenn. 54, 10414. https://doi.org/10.14214/sf.10414. Siipilehto, J., Mäkinen, H., 2025. A generalized marginal and mixed-effect models for predicting tree-level mortality with unequal measurement intervals for Scots pine in Finland. For. Sci. 71. https://doi.org/10.1007/s44391-025-00028-6. Siipilehto, J., Mäkinen, H., Andreassen, K., Peltoniemi, M., 2021. Models for integrating and identifying the effect of senescence on individual tree survival probability for Norway spruce. Silva Fenn. 55, 10496. https://doi.org/10.14214/sf.10496. Simard, S.W., Blenner-Hassett, T., Cameron, I.R., 2004. Pre-commercial thinning effects on growth, yield and mortality in even-aged paper birch stands in British Columbia. For. Ecol. Manag. 190, 163–178. https://doi.org/10.1016/j.foreco.2003.09.010. Sims, A., Mändma, R., Laarmann, D., Korjus, H., 2014. Assessment of tree mortality on the Estonian Network of Forest Research Plots/Puude väljalangevuse hindamine metsa kasvukäigu püsiproovitükkidel. Forestry Stud. 60, 57–68. https://doi.org/ 10.2478/fsmu-2014-0005. Temesgen, H., Mitchell, S.J., 2005. An individual-tree mortality model for complex stands of southeastern British Columbia. West. J. Appl. For. 20, 101–109. https:// doi.org/10.1093/wjaf/20.2.101. Vanclay, J.K., 1994. Modelling Forest Growth and Yield: Applications to Mixed Tropical Forests. CAB International, Wallingford, U.K. Wykoff, W.R., 1990. A basal area increment model for individual conifers in the northern Rocky Mountains. For. Sci. 36, 1077–1104. https://doi.org/10.1093/forestscience/ 36.4.1077. Yang, Y., Huang, S., 2013. A generalized mixed logistic model for predicting individual tree survival probability with unequal measurement intervals. For. Sci. 59, 177–187. https://doi.org/10.5849/forsci.10-092. Yao, X., Titus, S.J., MacDonald, S.E., 2001. A generalized logistic model of individual tree mortality for aspen, white spruce, and lodgepole pine in Alberta mixedwood forests. Can. J. For. Res. 31, 283–291. https://doi.org/10.1139/x00-162. J. Siipilehto and D. Lee Ecological Modelling 510 (2025) 111349 12 https://doi.org/10.14214/sf.10414 https://doi.org/10.1007/s44391-025-00028-6 https://doi.org/10.14214/sf.10496 https://doi.org/10.1016/j.foreco.2003.09.010 https://doi.org/10.2478/fsmu-2014-0005 https://doi.org/10.2478/fsmu-2014-0005 https://doi.org/10.1093/wjaf/20.2.101 https://doi.org/10.1093/wjaf/20.2.101 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0041 http://refhub.elsevier.com/S0304-3800(25)00335-7/sbref0041 https://doi.org/10.1093/forestscience/36.4.1077 https://doi.org/10.1093/forestscience/36.4.1077 https://doi.org/10.5849/forsci.10-092 https://doi.org/10.1139/x00-162 Comparison of marginal and mixed-effects complementary log-log regression models for predicting planted silver birch mortality 1 Introduction 2 Materials and methods 2.1 Study materials for tree mortality 2.1.1 Modelling dataset 2.1.2 Test dataset 2.2 Modelling approach and model structure 2.2.1 Complementary log-log regression model 2.2.2 Random effects consideration 2.2.3 Model formulation 2.3 Validation of the developed models 2.3.1 Comparison in the MOTTI simulations 2.3.2 Self-thinning models as the stand-level constraints 3 Results 3.1 Estimated models for five-year tree mortality 3.2 Model evaluation 3.3 Evaluation of the models in the MOTTI simulator 3.4 Evaluation of tree-level mortality models against self-thinning lines 4 Discussion 4.1 Overall performance of the developed models 4.2 Mixed-effects modelling and thinning effects on tree mortality 4.3 Recommendations and limitations for field application 5 Conclusion Funding information CRediT authorship contribution statement Declaration of competing interest Acknowledgements Appendix Data availability References