mastSim

mastSim generates inputs needed for model fitting with mast. These objects are simulated data, including the four data.frames ( treeData, seedData, xytree, xytrap ) and formulas ( formulaFec, formulaRep ). Other objects are “true” parameter values and latent states in the list truevalues used to simulate the data ( fec, repr, betaFec, betaRep, upar, R ). I want to see if mast can recover these parameter values for the type of data I simulated.

Here is a mapping of objects created by mastSim to the model in Clark et al. ( 2019 ):

Table 4. Some of the objects from the list mastSim.

mastSim assumes that predictors of maturation and fecundity are limited to intercept and diameter. Thus, the formulaFec and formulaRep are identical. Here is fecundity:

formulaFec

Before fitting the model I say a bit more about input data.

treeData and xytree

The information on individual trees is held in the tree-years by variables data.frame treeData. Here are a few lines of treeData generated by mastSim,

head( treeData )

treeData has a row for each tree-year. Several of these columns are required:

data.frame treeData columns

There can be additional columns, but they are not required for model fitting. The variable repr is reproduction status, often NA.

There is a corresponding data.frame xytree that holds tree locations.

Table 5. data.frame xytree columns

xytree has fewer rows than treeData, because it is not repeated each year–it assumes that tree locations are fixed. They are cross-referenced by plot and tree:

head( xytree, 5 )

##seedData and xytrap

Seed counts are held in the data.frame seedData as trap-years by seed types. Here are a few lines of seedData,

head( seedData )

Required columns are in Table 6.

Table 6. data.frame seedData columns

Seed trap location data are held in the data.frame xytrap.

Table 7. data.frame xytrap columns

There is no year in xytrap, because locations are assumed to be fixed. If a seed trap location is moved during a study, simply assign it a new trap name. Then seed counts for years when it is not present are assigned NA ( missing data ) will be imputed.

data summary and maps

Because they are stochastic, not all simulations from mastSim generate data with sufficient pattern to allow model fitting. I can look at the relationship between tree diameters and seeds on a map for plot mapPlot and year mapYear:

dataTab <- table( treeData$plot, treeData$year )

w <- which( dataTab > 0, arr.ind=T ) # a plot-year with observations
w <- w[sample( nrow( w ), 1 ), ]

mapYears <- as.numeric( colnames( dataTab )[w[2]] )
mapPlot  <- rownames( dataTab )[w[1]]
inputs$mapPlot    <- mapPlot
inputs$mapYears   <- mapYears
inputs$treeSymbol <- treeData$diam
inputs$SCALEBAR   <- T

mastMap( inputs )

Depending on the numbers and sizes of maps I adust treeScale and trapScale to see how seed accumulation compares with tree size or fecundity. In the above map, trees are shown as green circles and traps as gray squares. The size of the circle comes from the input variable treeSymbol, which is set to tree diameter. I could instead set it to the ‘true’ number of seeds produced by each tree, an output from mastSim.

inputs$treeSymbol <- trueValues$fec
inputs$treeScale  <- 1.5
inputs$trapScale  <- 1
mastMap( inputs )

Which are reproductive (mature)? Here are true values, showing just trees that are mature:

inputs$treeSymbol <- trueValues$repr
inputs$treeScale  <- .5
mastMap( inputs )

To fit the data, there must be sufficient seed traps, and sources must not be so dense such that overlapping seed shadows make the individual contributions undetectable.

Here is the frequency distribution of seeds ( observed ) and of fecundities ( unknown ) from the simulated data:

par( mfrow=c( 1, 2 ), bty='n', mar=c( 4, 4, 1, 1 ) )
seedData  <- inputs$seedData
seedNames <- inputs$seedNames

hist( as.matrix( seedData[, seedNames] ) , nclass=100, 
      xlab = 'seed count', ylab='per trap', main='' )
hist( trueValues$fec, nclass=100, xlab = 'seeds produced', ylab = 'per tree', main = '' )

In data of this type, most seed counts and most tree fecundities are zero.

output summary, lists, and plots

Sample size, parameter estimates, goodness of fit are all provides as tables by summary.

summary( output )

By default, this summary is sent to the console. It can also be saved to a list, e.g., outputSummary <- summary( output ).

plotPars <- list( trueValues = trueValues )
mastPlot( output, plotPars )

output   <- mastif( inputs = output, ng = 5000, burnin = 3000 )

output$predList <- list( mapMeters = 3, plots = mapPlot, years = mapYears )
output   <- mastif( inputs = output, ng = 2000, burnin = 1000 )

output$mapPlot    <- mapPlot
output$mapYears   <- mapYears
output$treeScale  <- 1.2
output$trapScale  <- .7
output$PREDICT    <- T
output$scaleValue <- 10
output$plotScale  <- 1
output$COLORSCALE <- T
output$LEGEND     <- T
  
mastMap( output )

plotPars$RMD <- 'pdf'
mastPlot( output, plotPars )

slow convergence?

Seed shadow models confront convoluted likelihood surfaces, in the sense that we expect local optima. These surfaces are hard to traverse with MCMC ( and impossible with HMC ), because maturation status is a binary state that must be proposed and accepted together with latent fecundity. Posterior simulation can get bogged down when fecundity estimates converge for a combination of mature and immature trees that is locally but not globally optimal. Especially when there are more trees of a species than there are seed traps, we expect many iterations before the algorithm can ‘find’ the specific combination of trees that together best describe the specific combination of seed counts in many seed traps. Compounding the challenge is the fact that, because both maturation and fecundity are latent variables for each individual, the redistribution kernel must be constructed anew at each MCMC step. Yet, analysis of simulated data shows that convergence to the correct posterior distribution is common. It just may depend on finding reasonable prior parameter values ( see below ) and long chains.

Examples in this vignette assign enough interations to show this progress toward convergence may be happening, but sufficiently few to avoid long waiting times. To evaluate convergence, consider the plots for chains of sigma ( the residual variance \(\sigma\) on log fecundity ), the rspse ( the seed count residual ), and the deviance. Finally, the plot of seed observed vs predicted gives a sense of progress.

As demonstrated above, restart mastif with the fitted object.

multiple seed types per species

Often a seed type could have come from trees of more than one species. Seeds that are only identified to genus level include in seedNames the character string UNKN. In this simulation the three species pinuTaeda, pinuEchi, and pinuVirg contribute most seeds to the type pinuUNKN. Important: there can be only one element in seedNames containing the string UNKN.

Here are some inputs for the simulation:

specNames <- c( 'pinuTaeda', 'pinuEchi', 'pinuVirg' )
seedNames <- c( 'pinuTaeda', 'pinuEchi', 'pinuVirg', 'pinuUNKN' )
sim    <- list( nyr=4, ntree=25, nplot=10, ntrap=50, specNames = specNames, 
                  seedNames = seedNames )

There is a specNames-by-seedNames matrix M that is estimated as part of the posterior distribution. For this example seedNames containing the string UNKN refers to a seed type that cannot be differentiated beyond the level of the genus pinu.

Here is the simulation with output objects with 2/3 of all seeds identified only to genus level:

inputs <- mastSim( sim )        # simulate dispersal data
R      <- inputs$trueValues$R   # species to seedNames probability matrix
round( R, 2 )

The matrix M stacks the length-\(M\) vectors \(\mathbf{m}'_s\) discussed in Clark et al. ( 2019 ) as a species-by-seed type matrix. There is a matrix for each plot, here stacked as a single matrix. This is the matrix of values used in simulation. mastif will estimate this matrix as part of the posterior distribution.

Here is a model fit:

output <- mastif( inputs = inputs, ng = 2000, burnin = 1000 ) 

The model summary now includes estimates for the unknown M as the “species to seed type matrix”:

summary( output )

Output plots will include chains for estimates in M. There will also be estimates for each species included in specNames:

plotPars <- list( trueValues = inputs$trueValues, RMAT = TRUE ) 
mastPlot( output, plotPars )

The inverse of \(\mathbf{M}\) is the probability that an unknown seed of type \(m\) was produced by a tree of species \(s\). These estimates are included in the plot Species to undiff seed.

Again, the .Rmd file can be knitted to a pdf file.

Here is a restart, now with plots and years specified for prediction in predList:

tab   <- with( inputs$seedData, table( plot, year ) )
mapPlot <- 'p1'
mapYears <- as.numeric( colnames( tab )[tab[1, ] > 0] )   # years for 1st plot
output$predList <- list( plots = mapPlot, years = mapYears ) 
output <- mastif( inputs = output, ng = 3000, burnin = 1500 )
mastPlot( output, plotPars = list( trueValues = inputs$trueValues, MAPS = TRUE )  )

Chains for the matrix M will converge for plots where there are sufficient trees and seeds to obtain an estimate. They will not be identified on plots where one or the other are rare or absent.

To map just the seed rain prediction maps, do this:

output$PREDICT  <- T
output$LEGEND   <- T
output$mapPlot  <- mapPlot
output$mapYears <- mapYears
mastMap( output )

specNames and seedNames

In the many data sets where seeds of multiple species are not differentiated the undifferentiated seed type occupies a column in seedData having a column name that includes the string UNKN. For a concrete example, if specNames on a plot include trees in treeData$species with the specNames pinuTaed and pinuEchi, then many seeds might be indentified as seedNames pinuUNKN. In the foregoing example, estimates of the matrix M reflect the contributions of each species to the UNKN type. This could be specified in several ways:

specNames <- c( 'pinuTaed', 'pinuEchi' )

#seeds not differentiated:
seedNames <- c( 'pinuUNKN' ) 

#one species sometimes differentiated:
seedNames <- c( 'pinuTaed', 'pinuUNKN' )    

#both species sometimes differentiated:
seedNames <- c( 'pinuTaed', 'pinuEchi', 'pinuUNKN' )

diameter effect

Here I fit the model.

formulaFec <- as.formula( ~ diam )    # fecundity model
formulaRep <- as.formula( ~ diam )    # maturation model
inputs   <- list( specNames = specNames, seedNames = seedNames, 
                 treeData = treeData, seedData = seedData, xytree = xytree, 
                 xytrap = xytrap )
output <- mastif( inputs, formulaFec, formulaRep, ng = 3000, burnin = 1000 )

Following this short sequence, I fit a longer one and predict one of the plots:

output$predList <- list( mapMeters = 10, plots = 'DUKE_EW', years = 2010:2015 ) 
output <- mastif( inputs = output, ng = 4000, burnin = 1000 )

Here are some plots followed by comments on display panels. Notice that when it progresses to the maps for plot DUKE_HW, they will include the predicted seed shadows:

mastPlot( output )

From fecundity chains, still more iterations are needed for convergence.

From seed shadow, seed rain beneath a 30-cm diameter tree averages near 0.2 seeds per m\(^2\).

Although a combination of maturation/fecundity can be found to predict the seed data in prediction ( part a ), the process model does not well predict the maturation/fecundity combination ( part b ). In other words, the maturation/fecundity/dispersal aspect of the model is effective, whereas the design does not yet include variables that help explain maturation/fecundity.

From the many maps in predicted fecundity, seed data, note that the DUKE_EW site includes prediction surfaces, as specfied in predList.

Here is the summary:

summary( output )

At this point in the the Gibbs sampler, the DIC and root mean square prediction error are:

output$fit

As mentioned above, prediction scores for seed-trap observations are based on the estimated fecundities \([\mathbf{y} | \phi, \rho]\) as scoreStates and on the estimated parameters \([\mathbf{y} | \phi, \rho][\phi, \rho| \boldsymbol{\beta}^x, \boldsymbol{\beta}^w, \dots]\) as predictScore. The TGc and TGd parameters are the affine-transformation parameters ( intercept, slope ) relating the predictive variance to the error variance ( Thorarinsdottir and Gneiting, 2010 ).

year and random effects

Individual ( tree differences ) can be random and year-to-year differences can be fixed; the latter are not ‘exchangeble’. In cases where there are no predictors that explain variation among individuals the formula can be limited to an intercept, using ~ 1. In this case, it can be valuable to estimate fecundity even when there are no good predictors. For this intercept-only model, random effects and/or year effects can allow for variation. Prior parameter values are supplied as discussed below.

#group plots in regions for year effects
region <- rep( 'sApps', nrow( treeData ) )
region[ as.character( treeData$plot ) == 'DUKE_EW' ] <- 'piedmont'

treeData$region <- region

formulaFec   <- as.formula( ~ diam )
formulaRep   <- as.formula( ~ diam ) 
yearEffect   <- list( groups = 'region' )
randomEffect <- list( randGroups = 'treeID', formulaRan = as.formula( ~ 1 ) )
inputs <- list( specNames = specNames, seedNames = seedNames, 
               treeData = treeData, seedData = seedData, 
               xytree = xytree, xytrap = xytrap, 
               priorDist = 28, priorVDist = 15, maxDist = 50, minDist = 15, 
               minDiam = 25, maxF = 1e+6, 
               randomEffect = randomEffect, yearEffect = yearEffect )
output <- mastif( inputs, formulaFec, formulaRep, ng = 2000, burnin = 1000 )
mastPlot( output )

Without predictors, the fitted variation is coming from random effects and year effects.

To substitute year effects for AR( \(p\) ) effects, simply remove the argument specification of p in the yearEffect list:

yearEffect   <- list( groups = c( 'species', 'region' ) )   # year effects

prior parameter values

In previous examples, the prior parameter values were past in the list inputs. I can also specify prior values as inputs through the inputs$priorList or inputs$priorTable holding parameters named in Table 8.

Table 8. Components of inputs with default values.

The mean and variance of the dispersal kernel, priorDist and priorVDist, refer to the dispersal parameter \(u\), transformed to the mean distance for the dispersal kernel. The parameters minDist and maxDist place minimum and maximum bounds on \(d\).
Fecundity \(\phi\) has a prior maximum value of maxF.

minDiam and maxDiam bound diameter ranges where a tree of unknown maturation status can be mature or not. This prior range is overridden by values in the treeData$repr column that establish an observed maturation state for a tree-year ( 0 - immature, 1 - mature ).

In general, large-seeded species, especially those that can be dispersed by vertebrates, generate noisy seed-trap data, despite the fact that the bulk of the counts still occur close to the parent, and the mean dispersal distance is relatively low. Large-seeded species produce few fruits/seeds. Long-distance dispersal cannot be estimated from inventory plots, regardless of plot size, because the fit is dominated by locally-derived seed. Consider values for maximum fecundity as low as maxF = 10000, minimum dispersal minDist = 2, and maximum dispersal maxDist = 12. Recall that the latter values refer to the dispersal kernel parameter value, not the maximum distance a seed can travel, which is un-constrained.

Prior parameter values can be passed as a list, e.g.,

inputs$priorList <- list( minDiam = 15, maxDiam = 60 )

Alternatively, prior parameter values can be passed as a table, with parameters for column names for specNames for rownames. Here is a file holding parameter values and the function mastPriors to obtain a table for several species in the genus \(Pinus\),

d <- "https://github.com/jimclarkatduke/mast/blob/master/priorParametersPinus.csv?raw=True"
download.file( d, destfile="priorParametersPinus.csv" )

specNames <- c( "pinuEchi", "pinuRigi", "pinuStro", "pinuTaed", "pinuVirg" )
seedNames <- c( "pinuEchi", "pinuRigi", "pinuStro", "pinuTaed", "pinuVirg", "pinuUNKN" )
priorTable <- mastPriors( "priorParametersPinus.csv", specNames, 
                         code = 'code4', genus = 'pinus' )

The argument code = 'code8' specifies the column in priorParametersPinus.csv holding the codes corresponding to specNames. In mastif the format code4 means that names consist of the first 4 letters from the genus and species, like this: pinuTaed. The genus is provided to allow for genus-level parameters in cases where priors for individual species have not been specified in priorParametersPinus.csv.

For coefficients in fecundity, betaPrior specifies predictor effects by sign. The use of prior distributions that are flat, but truncated, is two-fold ( Clark et al. 2013 ). First, where prior information exists, it often is limited to a range of values. For regressions coefficients ( e.g., \(\boldsymbol{\beta}\) ) we typically have a prior belief about the sign of the effect ( positive or negative ), but not its magnitude or appropriate weight. Second, within bounds placed by the prior, the posterior has the shape of the likelihood, unaffected by a prior weight. Prior weight is hard to specify sensibly for this hierarchical, non-linear model.

Here is an example using betaPrior to specify a quadratic diameter effect. First, here is a data set for Pinus:

d <- "https://github.com/jimclarkatduke/mast/blob/master/pinusExample.rdata?raw=True"
repmis::source_data( d )

Here I specify formulas, prior bounds for diam ( quadratic ), and fit the model:

formulaRep <- as.formula( ~ diam )
formulaFec <- as.formula( ~ diam + I( diam^2 ) )   
betaPrior  <- list( pos = 'diam', neg = 'I( diam^2 )' )
                  
inputs <- list( treeData = treeData, seedData = seedData, xytree = xytree, 
                xytrap = xytrap, specNames = specNames, seedNames = seedNames, 
                betaPrior = betaPrior, priorTable = priorTable, 
                seedTraits = seedTraits )
output <- mastif( inputs = inputs, formulaFec, formulaRep, 
                 ng = 500, burnin = 200 )

# restart
output <- mastif( output, formulaFec, formulaRep, ng = 5000, burnin = 1000 )
mastPlot( output )

Here again, many of the seed ID maxtrix elements are well identified ( trace plots labeled “M matrix plot”, others are uncertain due to small numbers of the different species that could contribute to a seed type. These are presented as bar and whisker plots for elements of matrix \(\mathbf{M}\) ( fraction of seed from a species that is counted as each seed type ) in the plot labeled species -> seed type. Bar plots show the fraction of unknown seeds that derive from each species, labeled species -> seed type. The distribution of diameters in the data will limit the ability to estimate its effect on maturation and fecundity.

Note that the estimates for fecundity coefficients, \(\boldsymbol{\beta}^x\), are positive for diam and negative for diam^2. If there is no evidence in the data for a decrease \(\frac{\partial{\psi}}{\partial{diam}}\), then the quadratic term is near zero. In this case, the predicted fecundity in the plot labelled maturation, fecundity by diameter will increase exponentially.

maturation/fecundity data, if available

Most individuals produce no seed, due to resource limitation, especially light, in crowded stands. As trees increase in size and resource access they mature, and become capable of seed production. Maturation is hard to observe, so it must be modelled. In our data sets, maturation status is observed, but most values are NA; it is only assigned if certain. A value of 1 means that reproduction is observed. A value of 0 means that the entire crown is visible during the fruiting season, and it is clear that no reproduction is present. Needless to say, most observations from beneath closed canopies are NA. The observations are entered in the column repr in the data.frame treeData. Above I showed an example where minDiam is set as a prior that trees of smaller diameters are immature and maxDiam is the prior that trees above this diameter are mature.

mastif further admits estimates of cone production. The column treeData$cropCount refers to the cones that will open and contribute to trap counts in the year treeData$year. There is another column treeData$cropFraction, which refers to the fraction of the tree canopy that is represented by the cone count. When crop counts are used, a matrix seedTraits must be supplied with one row per species and a column with seeds per fruit, cone, or pod. The rownames of seedTraits are the specNames. The colname of seedTraits with seeds per fruit is labeled seedsPerFruit.

when trees are sampled less frequently than seeds

Seed studies often include seed collections in years when trees are not censused. For example, tree data might be collected every 2 to 5 years, whereas seed data are available as annual counts. In this case, there there are years in seedData$year that are missing from treeData$year. mastif needs a tree year for each trap year. If tree data are missing from some trap years, I need to constuct a complete treeData data.frame that includes these missing years. This amended version of treeData must be accessible to the user, to allow for the addition of covariates such as weather variables, soil type, and so forth.

The function mastFillCensus allows the user to access the filled-in version of treeData that will be fitted by mastif. mastFillCensus accepts the same list of inputs that is passed to the function mastif. The missing years are inserted for each tree with interpolated diameters. The list inputs is returned with objects updated to include the missing census years and modified slightly for analysis by mastif. inputs$treeData can now be annotated with the covariates that will be included in the model. Here is the example from the help file. First I read a file that has complete years, so I randomly remove years:

# randomly remove years for this example:
years <- sort( unique( treeData$year ) )
sy    <- sample( years, 5 )
treeData <- treeData[treeData$year %in% sy, ]
treeData[1:10, ]

Note the missing years, as is typical of mast data sets. Here is the file after filling missing years:

inputs   <- list( specNames = specNames, seedNames = seedNames, 
                  treeData = treeData, seedData = seedData, 
                  xytree = xytree, xytrap = xytrap, priorTable = priorTable, 
                  seedTraits = seedTraits )
inputs <- mastFillCensus( inputs, beforeFirst=10, afterLast=10 )
inputs$treeData[1:10, ]

The missing plotYr combinations in treeData have been filled to match thos in seedData. Now columns can be added to inputs$treeData, as needed in formulaFec or formulaRep.

In cases where tree censuses start after seed trapping begins or tree censuses end before seed trapping ends, it may be reasonable to assume that the same trees are producing seed in those years pre-trapping or post-trapping years. In these cases, mastFillCensus can accommodate these early and/or late seed trap data by extrapolating treeData beforeFirst years before seed trapping begins or afterLast years after seed tapping ends.

adding covariates

Continuing with the \(Pinus\) example, here I want to add covariates for missing years as needed for an AR( \(p\) ) model, in this example \(p\) = 4. First, consider the tree-years included in this sample, before and after in-filling:

# original data
table( treeData$year )

# filled census
table( inputs$treeData$year )
table( seedData$year )

Note that the infilled version of tree data in the list inputs has the missing years, corresponding to those included in seedData. Here I set up the year effects model and infill to allow for p = 3,

p <- 3
inputs   <- mastFillCensus( inputs, p = p )
treeData <- inputs$treeData

Here are regions for random effects in the AR( 3 ) model:

region <- c( 'CWT', 'DUKE', 'HARV' )
treeData$region <- 'SCBI'
for( j in 1:length( region ) ){
  wj <- which( startsWith( treeData$plot, region[j] ) )
  treeData$region[wj] <- region[j]
}
inputs$treeData <- treeData
yearEffect <- list( groups = c( 'species', 'region' ), p = p )

I add the climatic deficit ( monthly precipitation minus PET ) as a covariant. First, here is a data file,

d <- "https://github.com/jimclarkatduke/mast/blob/master/def.csv?raw=True"
download.file( d, destfile="def.csv" )

The format for the file dev.csv is plot by year_month. I have saved it in my local directory. The function mastClimate returns a list holding three vectors, each having length equal to nrow( treeData ). Here is an example for the cumulative moisture deficit for the previous summer. I provide the file name, the vector of plot names, the vector of previous years, and the months of the year. To get the cumulative deficit, I use FUN = 'sum':

treeData <- inputs$treeData
deficit  <- mastClimate( file = 'def.csv', plots = treeData$plot, 
                         years = treeData$year - 1, months = 6:8, 
                         FUN = 'sum', vname='def' )
treeData <- cbind( treeData, deficit$x )
summary( deficit )

The first column in deficit is the variable itself for each tree-year. The second column holds the site mean value for the variable. The third column is the difference between the first two columns All three can be useful covariates, each capturing different effects ( Clark et al. 2014 ). I could append any or all of them as columns to treeData.

Here is an example using minimum temperature of the preceeding winter. I obtain the minimum with two calls to mastClimate, first for Dec of the previous year ( years = treeData$year - 1, months = 12 ), then from Jan, Feb, Mar for the current year ( years = treeData$year, months = 1:3 ). I then take the minimum of the two values:

# include min winter temperature
d <- "https://github.com/jimclarkatduke/mast/blob/master/tmin.csv?raw=True"
download.file( d, destfile="tmin.csv" )

# minimum winter temperature December through March of previous winter

t1 <- mastClimate( file = 'tmin.csv', plots = treeData$plot, 
                   years = treeData$year - 1, months = 12, FUN = 'min', 
                   vname = 'tmin' )
t2 <- mastClimate( file = 'tmin.csv', plots = treeData$plot, 
                   years = treeData$year, months = 1:3, FUN = 'min', 
                   vname = 'tmin' )
tmin <- apply( cbind( t1$x[, 1], t2$x[, 1] ), 1, min )
treeData$tminDecJanFebMar <- tmin
inputs$treeData <- treeData

Here is a model using several of these variables:

formulaRep <- as.formula( ~ diam )
formulaFec <- as.formula( ~ diam + defJunJulAugAnom + tminDecJanFebMar )   
inputs$yearEffect <- yearEffect
output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2500, burnin = 1000 )

Data for future years can be a scenario, based on assumptions of status quo in the mean and variance, an assumed rate of climate change, and so on.

random individual effects

Random individual effects currently can include random intercepts for the fecundities of each individual tree that is imputed to be in the mature state for at least 3 years. [There are no random effects on maturation, because they would be hard to identify from seed trap data for this binary response.]

I include in the inputs list the list randomEffect, which includes the column name for the random group. Typically this would be a unique identifier for a tree within a plot, e.g., randomEffect$randGroups = 'tree'. However, randGroups could be the plot name. This column is interpreted as a factor, each level being a group. Random effects will not be fitted on individuals that are in the mature state less than 3 years.

The formulaRan is the random effects model. Because individual time series tend to short, it is currently implemented only for intercepts. Here is a fit with random effects.

formulaFec <- as.formula( ~ diam )    # fecundity model
formulaRep <- as.formula( ~ diam )    # maturation model 
inputs$randomEffect <- list( randGroups = 'tree', formulaRan = as.formula( ~ 1 ) )
output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2000, burnin = 1000 )

Here is a restart:

output <- mastif( inputs = output, ng = 4000, burnin = 1000 )
mastPlot( output )

There is a new panel for fixed plus random effects, showing the individual combinations of intercepts.

• the prediction panel ( b ) shows some improvement, indicating that even with random effects, diameter struggles to predict maturation/fecundity.

output$fit

random groups for years and lags

As discussed previously, the model admits year effects and lag effects, the latter as an AR( \(p\) ) model. Year effects assign a coefficient to each year \(t = [1, \dots, T_j]\). Lag effects assign a coefficient to each of \(j \in \{1, \dots, p\}\) plags, where the maximum lag \(p\) should be substantially less than the number of years in the study. Year effects can be organized in random groups. Specification of random groups is done in the same way for year effects and for lag effects.

I define random groups for year and lag effects by species, by plots, or both. When there are multiple species that contribute to the modeled seed types, I expect the year effects to depend on which species is actually producing the seed. When there are multiple plots sufficiently distant from one another, I might allow for the fact that year effects or lag effects differ by group; yearEffect$groups allows that they need not mast in the same years. In the example below, I’ll use the term region for plots in the same plotGroup.

Here is a breakdown for this data set by region:

with( treeData, colSums( table( plot, region ) ) )

Here are year effects structured by random groups of plots, given by the column region in treeData:

yearEffect <- list( groups = c('species', 'region' ) )

This option will fit a year effect for both provinces and years having sufficient individuals estimated to be in the mature state. Here is the model with random year effects,

\[\log \psi_{ij, t} \sim N \left( \mathbf{x}'_{ij, t} \mathbf{\beta}^x + \gamma_t + \gamma_{g[i], t}, \sigma^2 \right )\] where the year effect \(\gamma_{g[i], t}\) is shared by trees in all plots defined by yearEffect$province, \(ij \in g\). Year effects are sampled directly from conditional posteriors.

inputs$treeData      <- treeData 
inputs$randomEffect  <- randomEffect
inputs$yearEffect    <- yearEffect
output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2500, burnin = 500 )

Here is a restart, with predictions for one of the plots:

predList <- list( mapMeters = 10, plots = 'DUKE_BW', years = 1998:2014 ) 
output$predList <- predList
output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
mastPlot( output )

Note that still more iterations are needed for convergence. Here are some comments on mastPlot:

• In the dispersal parameter u panel there are now year effects plotted for the two random groups mtn and piedmont.

• The subsequent panel dispersal mean and variance shows the mean and variance of random effects

• There is a dispersal by group panel showing posterior estimate for the two random groups, with scales for the parameter \(u\) on the left ( m\(^2\) ) and mean parameter \(d\) on the right ( m ).

• There has been some improvement in the prediction, panel ( b ).

• The predicted fecundity, seed data maps for the plot DUKE_BW show seed prediction surfaces.

• The year effect groups shows year effects for random groups in treeData$region.

• partial ACF shows partial autocorrelation by species and plot.

random effects when there are multiple species

Most data sets have multiple seed types that complicate estimation of mast production by each species. This example considers Pinus spp, seeds of which cannot be confidently assigned to species. Here I load the data and generate a sample of maps from several years, including all species and seed types.

d <- "https://github.com/jimclarkatduke/mast/blob/master/pinusExample.rdata?raw=True"
repmis::source_data( d )

mapList <- list( treeData = treeData, seedData = seedData, 
                 specNames = specNames, seedNames = seedNames, 
                 xytree = xytree, xytrap = xytrap, mapPlot = 'DUKE_EW', 
                 mapYears = c( 2007:2010 ), treeSymbol = treeData$diam, 
                 treeScale = .6, trapScale=1.4, 
                 plotScale = 1.2, LEGEND=T )
mastMap( mapList )

Note the tendency for high seed accumulation ( large green squares ) near dense, large trees ( large brown circles ).

In this example, I again model seed production as a function of log diameter, diam, now for multiple species and seed types. This is an AR( p ) model, because I include the number of lag terms yearEffect$p = 5,

\[\log \psi_{ij, t} \sim N \left( \mathbf{x}'_{ij, t} \mathbf{\beta}^x + \sum^p_{l=1} ( \alpha_l + \alpha_{g[i], l} ) \psi_{ij, t-l}, \sigma^2 \right )\] Only years \(p < t \le T_i\) are used for fitting. Samples are drawn directly from the conditional posterior distribution.

In the table printed at the outset are trees by plot and year, i.e., the groups assigned in inputs$yearEffect. The zeros indicate either absence of trees or that no plots were sampled in those years. ( These are not the same thing, mastif knows the difference ).

In the code below I specify formulas, AR( \(p\) ), and random effects, and some prior values. Due to the large number of trees, convergence is slow. Because I do not assume that trees of different species necessarily mast in the same years, I allow them to differ through random groups on the AR( \(p\) ) terms.

formulaFec <- as.formula( ~ diam )   # fecundity model
formulaRep <- as.formula( ~ diam )   # maturation model

yearEffect   <- list( groups = 'species', p = 4 )   # AR( 4 )
randomEffect <- list( randGroups = 'tree', 
                     formulaRan = as.formula( ~ 1 ) )

inputs   <- list( specNames = specNames, seedNames = seedNames, 
                  treeData = treeData, seedData = seedData, 
                  yearEffect = yearEffect, 
                  randomEffect = randomEffect,
                  xytree = xytree, xytrap = xytrap, priorDist = 20, 
                  priorVDist = 5, minDist = 15, maxDist = 30, 
                  minDiam = 12, maxDiam = 40, 
                  maxF = 1e+6, seedTraits = seedTraits )
output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 500, burnin = 100 )

Here is a restart:

output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
mastPlot( output, plotPars = plotPars )

Again, convergence will require more iterations. The AR( \(p\) ) coefficients in the lag effect group panel shows the coefficients by random group. They are also shown in a separate panel, with each group plotted separately. In the ACF eigenvalues panel are shown the eigenvalues for AR lag coefficients on the real ( horizontal ) and imaginary ( vertical ) scales with the unit circle, within which oscillations are damped. The imaginary axis describes oscillations.

Here’s a restart with predictions:

plots <- c( 'DUKE_EW', 'CWT_118' )
years <- 1980:2025
output$predList <- list( mapMeters = 10, plots = plots, years = years ) 
output <- mastif( inputs = output, ng = 3000, burnin = 1000 )

and updated plots:

mastPlot( output, )

Note that convergence requires additional iterations ( larger ng ). The predictions of seed production will progressively improve with convergence.

mapList <- output
mapList$mapPlot <- 'DUKE_EW'
mapList$mapYears <- c( 2011:2012 )
mapList$PREDICT <- T
mapList$treeScale <- .5
mapList$trapScale <- .8
mapList$LEGEND <- T
mapList$scaleValue <- 50
mapList$plotScale  <- 2
mapList$COLORSCALE <- T
mapList$mfrow <- c( 2, 1 )

mastMap( mapList )

Or a larger view of a single map:

mapList$mapPlot <- 'CWT_118'
mapList$mapYears <- 2015
mapList$PREDICT <- T
mapList$treeScale <- 1.5
mapList$trapScale <- .8
mapList$LEGEND <- T
mapList$scaleValue <- 50
mapList$plotScale <- 2
mapList$COLORSCALE <- T
mapList$mfrow <- c( 1, 1 )
mastMap( mapList )

Here is a summary of parameter estimates:

summary( output )