model summary

The basic model is detailed in Clark et al. (2017) and summarized at the end of this vignette (see reference notes). An observation consists of two vectors \((\mathbf{x}_i, \mathbf{y}_i)^n_{i = 1}\), where \(\mathbf{x}_i\) is a length-\(Q\) design vector (intercept and predictors) and \(\mathbf{y}_i\) is a length-\(S\) vector of response variables, each potentially measured in different ways. \(\mathbf{y}_i\) may include both continuous and discrete variables. There is a latent vector \(\mathbf{w}_i\) that represents response variables all on a continuous scale. The vector \(\mathbf{w}_i\) has the joint distribution \(\mathbf{w}_i \sim MVN(\boldsymbol{\mu}_i, \boldsymbol{\Sigma})\), where \(\boldsymbol{\mu}_i\) is the length-\(S\) mean vector, and \(\boldsymbol{\Sigma}\), is an \(S \times S\) covariance matrix.

As a data-generating mechanism the model can be thought of like this: There is a vector of continuous responses \(\mathbf{w}_{i}\) generated from mean vector \(\boldsymbol{\mu}_{i}\) and covariance \(\boldsymbol{\Sigma}\) (Fig. 1a). A partition \(\mathbf{p}_{is} = (-\infty, \dots, \infty)\) segments the real line into intervals, some of which are censored and others not. Each interval is defined by two values, \((p_{is,k}, p_{is,k+1}]\). For a value of \(w_{is}\) that falls within a censored interval \(k\) the observed \(y_{is}\) is assigned to discrete interval \(z_{is} = k\). For a value of \(w_{is}\) that falls in an uncensored interval \(y_{is}\) is assigned \(w_{is}\).

Of course, data present us with the inverse problem: the observed \(y_{is}\) are continuous or discrete, with known or unknown partition \((p_{is,k}, p_{is,k+1}]\) (Fig. 1b). Depending on how the data are observed, we must impute the elements of \(n \times S\) matrix \(\mathbf{W}\) that lie within censored intervals. Unknown elements of \(\mathcal{P}\) will also be imputed in order to estimate \(\mathbf{B}\) and \(\boldsymbol{\Sigma}\).

Censoring in gjam. As a data-generating model (a), a realization \(w_{is}\) that lies within a censored interval is translated by the partition \(\mathbf{p}_{is}\) to discrete \(y_{is}\). The distribution of data (bars at left) is induced by the latent scale and the partition. For inference (b), observed discrete \(y_{is}\) takes values on the latent scale from a truncated distribution.

data types

The different types of data that can be included in the model are summarized in Table 1, assigned to the character variable typeNames that is included in the modelList passed to gjam:

Table 1. Partition for each data type

\(^1\) For 'DA' and 'CC' data the second element of the partition is not zero, but rather depends on effort. There is thus zero-inflation. The default partition for each data type can be changed with the function gjamCensorY (see specifying censored intervals).
\(^2\) For 'CAT' data species \(s\) has \(K_s - 1\) non-reference categories. The category with the largest \(w_{is,k}\) is the ‘1’, all others are zeros.

effort and weight of discrete data

The partition for a discrete interval \(k\) depends on effort for sample \(i\)

\[(p_{i,k}, p_{i,k+1}] = \left(\frac{k - 1/2}{E_{i}}, \frac{k + 1/2}{E_{i}}\right]\]

Effort affects the partition and, thus, the weight of each observation; wide intervals allow large variance, and vice versa. For discrete abundance ('DA') data on plots of a given area, large plots contribute more weight than small plots. Because plots have different areas one might choose to model \(w_{is}\) on a ‘per-area’ scale (density) rather than a ‘per-plot’ scale. If so, plot area becomes the ‘effort’. Here is a table of variables for the case where counts represent the same density of trees per area, but have different effort due to different plot areas:

The wide partition on the 0.1-ha plot admits large variance around the observation of 10 trees per 0.1 ha plot. Wide variance on an observation decreases its contribution to the fit. Conversely, the narrow partition on the 1.0-ha plot constrains density to a narrow interval around the observed value.

For composition count ('CC') data effort is represented by the total count. For \(0 < y_{is} < E_i\) the variable \(0 < w_{is} < 1\), i.e., the composition scale. Using the same partition as previously the table for two observations that represent the fraction 0.10 with different effort (e.g., total reads in PCR data) looks like this:

Again, on the composition scale \([0, 1]\), weight of the observation is determined by the partition width and, in turn, effort.

simulated data

Simulated data are used to check that the algorithm can recover true parameter values and predict data, including underlying latent variables. To illustrate I simulate a sample of size \(n = 500\) for \(S = 10\) species and \(Q = 4\) predictors. To indicate that all species are continuous abundance data I specify typeNames as 'CA'. CA data are continuous above zero, with point mass at zero.

library(gjam)
f <- gjamSimData(n = 500, S = 10, Q = 4, typeNames = 'CA')
summary(f)

The object f includes elements needed to analyze the simulated data set. f$typeNames is a length-\(S\) character vector. The formula follows standard R syntax. It does not start with y ~, because gjam is multivariate. The multivariate response is supplied as a \(n \times S\) matrix or data.frame ydata. Here is the formula for this example:

f$formula

The simulated parameter values are returned from gjamSimData in the list f$trueValues, shown in Table 2, with the corresponding names of estimates from function gjam, which I discuss below:

Table 2. Variable names and dimensions in simulation and fitting

\(^1\) simulated object from gjamSimData.

\(^2\) fitted object from gjam.

\(^3\) coefficients are \(\mathbf{B}_u\): unstandardized; \(\mathbf{B}\): standardized for \(\mathbf{X}\); \(\tilde{\mathbf{B}}\): standardized for \(\mathbf{X}\), correlation scale for \(\mathbf{W}\) and factor levels relative to the mean for each factor (see section factors in \(\mathbf{X}\)**).

\(^4\) sensitivities based on unstandardized \(\mathbf{B}\) and covariance scales \(\Sigma\).

\(^5\) Only when ydata includes ordinal types "OC".

\(^6\) Only with dimension reduction, reductList is included in modelList.

\(^7\) Only for trait analysis, traitList is included in modelList (Trait vignette).

The dimension for response \(Y\) is \(W \times E\), where \(W\) is the latent variable scale, and \(E\) is sample effort. When effort \(E\) = 1, then \(Y\) and \(W\) have the same dimension.

The matrix \(\mathbf{F}\) contains the covariance between predictors in \(\mathbf{X}\) in terms of the responses they elicit from \(\mathbf{Y}\). The diagonal vector \(\mathbf{f} = diag( \mathbf{F} )\) is the sensitivity of the entire response matrix to each predictor in \(\mathbf{X}\).

The matrix \(\mathbf{E}\) is the correlation among species in terms of their responses to \(\mathbf{X}\). Relationships to outputs are discussed in the Reference notes.

Simulated data are typical of real data in that there is a large fraction of zeros.

par(bty = 'n', mfrow = c(1,2), family='')
h <- hist(c(-1,f$y), nclass = 50, plot = F)
plot( h$counts, h$mids,type = 's' )
plot( f$w,f$y,cex = .2 )

Here is a short Gibbs sampler to estimate parameters and fit the data. The function gjam needs the formula for the model, the data.frame xdata, which includes the predictors, the response matrix or data.frame ydata, and a modelList specifying number of Gibbs steps ng, the burnin, and typeNames.

ml  <- list( ng = 1000, burnin = 100, typeNames = f$typeNames )
out <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )
summary(out)

The print function acts like summary

print(out)

Among the objects to consider initially are the design matrix out$inputs$xUnstand, response matrix out$inputs$y, and the MCMC out$chains with these names and sizes:

summary(out$chains)

$chains is a list of matrices, each with ng rows and as many columns as needed to hold parameter estimates. For example, each row of $chains$bgibbs is a length-\(QS\) vector of values for the \(Q \times S\) matrix \(\mathbf{B}\), standardized for \(\mathbf{X}\). In other words, a coefficient \(\mathbf{B}_{qs}\) has the units of \(w_s\). $chains$sgibbs holds either the \(S(S + 1)/2\) unique values of \(\boldsymbol{\Sigma}\) or the \(N \times r\) unique values of the dimension reduced covariance matrix. A summary of the chains is given in Table 2.

Additional summaries are available in the list inputs:

summary(out$inputs)

The matrix classBySpec shows the number of observations in each interval. For this example of continuous data censored at zero, the two labels are \(k \in \{0, 1\}\) corresponding to the intervals \((p_{s,0}, p_{s,1}] = (-\infty,0]\) and \((p_{s,1}, p_{s,2}) = (0, \infty)\). The length-\((K + 1)\) partition vector is the same for all species, \(\mathbf{p} = (-\infty, 0, \infty)\). Here is classBySpec for this example:

out$inputs$classBySpec

The first interval is censored (all values of \(y_{is}\) = 0). The second interval is not censored (\(y_{is} = w_{is}\)).

The fitted coefficients in $parameters, as summarized in Table 2. For example, here is posterior mean estimate of unstandardized coefficients \(\mathbf{B}_u\),

out$parameters$betaMu

Here are posterior summaries,

out$parameters$betaMu           # S by M coefficient matrix unstandardized
out$parameters$betaSe           # S by M coefficient SE
out$parameters$betaStandXmu     # S by M standardized for X
out$parameters$betaStandXWmu    # (S-F) by M standardized for W/X, centered factors

out$parameters$betaTable        # SM by stats posterior summary
out$parameters$betaStandXTable  # SM by stats posterior summary
out$parameters$betaStandXWTable # (S-F)M by stats posterior summary

out$parameters$sensBeta         # sensitivity to response variables
out$parameters$sensTable        # sensitivity to predictor variables
 
out$parameters$sigMu            # S by S covariance matrix omega
out$parameters$sigSe            # S by S covariance std errors

Again, check Table 2 for names of all fitted coefficients.

The data are also predicted in gjam, summarized by predictive means and standard errors. These are contained in \(n \times Q\) matrices $prediction$xpredMu and $prediction$xpredSd and \(n \times S\) matrices $prediction$ypredMu and $prediction$ypredSd. These are in-sample predictions.

The output can be viewed with the function gjamPlot:

f   <- gjamSimData( n = 500, S = 10, typeNames = 'CA' )
ml  <- list( ng = 1000, burnin = 200, typeNames = f$typeNames )
out <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )
pl  <- list( trueValues = f$trueValues, GRIDPLOTS = T )
gjamPlot( output = out, plotPars = pl )

gjamPlot creates a number of plots comparing true and estimated parameters (for simulated data).

par( bty = 'n', mfrow = c(1,3) )
plot( f$trueValues$beta, out$parameters$betaMu, cex = .2 )
plot( f$trueValues$corSpec, out$parameters$corMu, cex = .2 )
plot( f$y,out$prediction$ypredMu, cex = .2 )

To process the output beyond what is provided in gjamPlot I can work directly with the chains.

my data

gjam uses the standard R syntax in the formula that I would use with functions like lm and glm. Because gjam uses inverse prediction to summarize large multivariate output, it is important to abide by this syntax. For example, to analyze a model with quadratic and interaction terms, I might simply construct my own design matrix with these columns included, i.e., side-stepping the standard syntax for these effects that can be specified in formula. This would be fine for model fitting. However, without specifying this in the formula there is no way for gjam to know that these columns are in fact non-linear transformations of other columns. Without this knowledge there is no way to properly predict them. The prediction that gjam would return would include silly variable combinations.

To illustrate proper model specification I use a few lines from the data.frame of predictors in the forestTraits data set:

library(repmis)
d <- "https://github.com/jimclarkatduke/gjam/blob/master/forestTraits.RData?raw=True"
source_data(d)
xdata <- forestTraits$xdata[,c(1,2,7,8)]

Here are a few rows:

xdata[1:5,]

Here is a simple model specification with as.formula() that includes only main effects:

formula <- as.formula( ~ temp + deficit + soil )

The design matrix x that is generated in gjam has an intercept, two covariates, and four columns for the multilevel factor soil:

##   (Intercept)  temp deficit soilreference
## 1           1  1.22    0.04             1
## 2           1  0.18    0.21             1
## 3           1 -0.94    0.20             0
## 4           1  0.64    0.82             1
## 5           1  0.82   -0.18             1

To include interactions between temp and soil I use the symbol ‘*’:

formula <- as.formula( ~ temp*soil )

Here is the design matrix that results from this formula with interaction terms indicated by the symbol ':'

##   (Intercept)  temp soilreference temp:soilreference
## 1           1  1.22             1               1.22
## 2           1  0.18             1               0.18
## 3           1 -0.94             0               0.00
## 4           1  0.64             1               0.64
## 5           1  0.82             1               0.82

For a quadratic term I use the R function I():

formula <- as.formula( ~ temp + I(temp^2) + deficit )

Here is the design matrix with linear and quadratic terms:

##   (Intercept)  temp I(temp^2) deficit
## 1           1  1.22    1.4884    0.04
## 2           1  0.18    0.0324    0.21
## 3           1 -0.94    0.8836    0.20
## 4           1  0.64    0.4096    0.82
## 5           1  0.82    0.6724   -0.18

Here is a quadratic response surface for temp and deficit:

formula <- as.formula( ~ temp*deficit + I(temp^2) + I(deficit^2) )

Here is the design matrix with all combinations:

##   (Intercept)  temp deficit I(temp^2) I(deficit^2) temp:deficit
## 1           1  1.22    0.04    1.4884       0.0016       0.0488
## 2           1  0.18    0.21    0.0324       0.0441       0.0378
## 3           1 -0.94    0.20    0.8836       0.0400      -0.1880
## 4           1  0.64    0.82    0.4096       0.6724       0.5248
## 5           1  0.82   -0.18    0.6724       0.0324      -0.1476

These are examples of the formula options available in gjam. Using them will allow for proper inverse prediction of x. To optimize MCMC, gjam does not predict x for higher order polynomials–they are rarely used, being both hard to interpret and a cause of unstable predictions. To accelerate MCMC I can set PREDICTX = F in the modelList.

I can use this model to analyze a tree data set. For my data set I use the tree data contained in forestTraits. It is stored in de-zeroed (sparse matrix) format, so I extract it with the function gjamReZero. Here are dimensions and the upper left corner of the response matrix \(\mathbf{Y}\),

y  <- gjamReZero(forestTraits$treesDeZero)  # extract y
treeYdata  <- gjamTrimY(y,10)$y             # at least 10 plots
dim(treeYdata)
treeYdata[1:5,1:6]

In code that follows I treat responses as discrete counts, typeNames = 'DA'. Because of the large number of columns (98) I speed things up by calling for dimension reduction, passed as \(N \times r = 20 \times 8\):

ml   <- list( ng = 2500, burnin = 500, typeNames = 'DA', 
              reductList = list(r = 8, N = 20) )
form <- as.formula( ~ temp*deficit + I(temp^2) + I(deficit^2) )
out  <- gjam(form, xdata = xdata, ydata = treeYdata, modelList = ml)
specNames <- colnames(treeYdata)
specColor <- rep('black',ncol(treeYdata))
specColor[ c(grep('quer',specNames),grep('cary',specNames)) ] <- '#d95f02'
specColor[ c(grep('acer',specNames),grep('frax',specNames)) ] <- '#1b9e77'
specColor[ c(grep('abie',specNames),grep('pice',specNames)) ] <- '#377eb8'

gjamPlot( output = out, plotPars = list(GRIDPLOTS=T, specColor = specColor) )

Additional information on variable types and their treatment in gjam is included later in this document and in the other gjam vignettes.

heterogeneous effort

Here is an example with discrete abundance data, now with heterogeneous sample effort. Heterogeneous effort applies where vegetation plots have different areas, animal survey data have variable search times, or catch returns from fishing vessels have different gear and/or trawling times. Here I simulate a list containing the columns and the effort that applies to those columns, shown for 50 observations:

S  <- 5                             
n  <- 50
ef <- list( columns = 1:S, values = round(runif(n,.5,5),1) )
f  <- gjamSimData(n, S, typeNames = 'DA', effort = ef)
ef

If ef$values consists of a length-n vector, then gjam assumes each value applies to all columns in the corresponding row specified in the vector ef$columns. This is the case shown above and would apply when effort is plot area, search time, sample volume, and so forth. Alternatively, values can be supplied as a matrix, having the same dimensions as ydata. As a matrix, ef$values can have elements that differ by observation and species. For example, camera trap data detect large animals at greater distances than small animals (column differences), and each camera can be deployed for different lengths of time (row differences). For simulation purposes gjamSimData only accepts a vector. However, for fitting with gjam effort$values can be supplied as a matrix with as many columns as are listed in effort$columns.

Because observations are discrete the continuous latent variables \(w_{is}\) are censored. Unlike the previous continuous example, observations \(y_{is}\) now assume only discrete values:

plot( f$w, f$y, cex = .2, xlab = 'Latent w scale', ylab = 'Observed y' )

The large scatter reflects the variable effort represented by each observation. Incorporating the effort scale gives this plot:

plot(f$w*ef$values, f$y, cex = .2, xlab = 'Latent w time effort', ylab = 'Observed y')

The heterogeneous effort affects the weight of each observation in model fitting. The effort is entered in modelList. Increase the number of iterations and look at plots:

S   <- 10                             
n   <- 1500
ef  <- list( columns = 1:S, values = round(runif(n,.5,5),1) )
f   <- gjamSimData(n, S, typeNames = 'DA', effort = ef)
ml  <- list(ng = 1000, burnin = 250, typeNames = f$typeNames, effort = ef)
out <- gjam(f$formula, f$xdata, f$ydata, modelList = ml)
pl  <- list( trueValues = f$trueValues )
gjamPlot(output = out, plotPars = pl)

Use summary(out) to see a summary of effort.

specifying censored intervals

To analyze data that are censored at specific intervals, I specify the censored values and the intervals to which those values apply. In the help page for gjamCensorY I discuss a vector of values and a corresponding matrix of upper and lower bounds in two rows and one column for each element in values. For example, if an observer stops counting wildebeest or sea lions at, say, 100 animals, I can set this as a censored interval:

# assumes there is a data matrix ydata
upper <- 100
intv  <- matrix(c(upper,Inf),2)
rownames(intv) <- c('lower','upper')
tmp <- gjamCensorY(values = upper, intervals = intv, y = f$ydata, type='DA')

There are two objects returned by gjamCensor. The list $censor holds two lists, the $columns, indicating which columns in $y are censored, and the $partition, giving the values and bounds for the partition. Because I did not specify whichcol in my call to gjamCensorY, censoring defaults to all columns in ydata. The object $y matches the mode and dimensions of the input y (a matrix or data.frame). This new version of the response matrix replaces any values in ydata that fall within a censored interval with the censored value. This feature is useful when observers are inconsistent on the assignment of intervals or where an analyst distrusts the precision reported in data. For example, if counts range to thousands, but it is known that counts beyond 100 are inaccurate, then all counts above 100 are censored and, thus, uncertain.

Censoring can be applied differently to each response variable. For example, chemical constitents reported as concentrations in a sample, sometimes reach non-zero minimum values, taken as detection limits for that instrument and that constituent (detection limits can differ for each constituent). In this case, there is a list for each column in ydata. I can generate this list with a loop, where the censored interval for each column j spans from -Inf to min(ydata[,j]),

miny   <- apply(f$ydata, 2, min, na.rm=T)     #minimum value by column
censor <-  numeric(0)     
p      <- matrix(0, 3, dimnames=list(c("values","lower","upper"), NULL))

for(j in 1:ncol(f$ydata)){
  p[1:3] <- c(miny[j], -Inf, miny[j])
  jlist  <- list("columns" = j, "partition" = p)
  censor <- c(censor, list(jlist))
  names(censor)[j] <- 'CA'
}

This list censor can be passed directly to gjam in the list modelList.

prior distribution on coefficients

Informative prior distributions in regression models are rare, perhaps partly because it is hard to assign both magnitude (e.g., large or small prior mean value) and weight (large or small prior variance) without obscuring the relative contributions of prior and data. Prior distributions for regression coefficients are typically Gaussian, having support on \((-\infty, \infty)\). In many cases, the sign of the effect is known, but the magnitude is not. Ad hoc experimentation with prior mean and variance can, at best, only insure that ‘most’ of the posterior distribution is positive or negative. Yet, it can be important to use prior information, especially in the multivariate setting, where covariances between responses can result in estimates where the sign of a coefficient effect makes no sense.

The knowledge of the ‘direction’ of the effect can be readily implmented with uniform priors truncated at zero, having the advantage that the posterior distribution that preserves the shape of the likelihood, but is restricted to positive or negative values (Clark et al. 2013).

The prior distribution for \(\mathbf{B}\) is either non-informative (if unspecified) or truncated by limits provided in the list betaPrior. The betaPrior list contains the two matrices lo and hi. The rows of these matrices have rownames that match explanatory variables in the formula and colnames in xdata. In the example that follows I fit a model for FIA data to winter temperature temp, climatic deficit, and their interaction. For this example, I use a prior distribution having positive effects of warm winters and negative effects of climate deficit. This prior is set up with the function gjamPriorTemplate.

The prior distribution is also the place to exclude specific predictor/response combinations, by setting them equal to NA in lo or hi. Here is an example of informative priors on some coefficients:

xdata   <- forestTraits$xdata
formula <- as.formula(~ temp*deficit)
snames  <- colnames(treeYdata)

# warm winter
hot <- c("liquStyr","liriTuli","pinuEchi","pinuElli","pinuPalu","pinuTaed",
         "querImbr","querLaur","querLyra","querMich","querMueh","querNigr",
         "querPhel","querVirg")  # arbitrary spp, positive winter temp
nh <- length(hot)
lo  <- vector('list', nh)
names(lo) <- paste('temp',hot,sep='_')
for(j in 1:nh)lo[[j]] <- 0

# humid climate (negative deficit)
humid <- c("abieBals", "betuAlle", "querNigr", "querPhel")  #again, arbitrary
nh <- length(humid)
hi  <- vector('list', nh)
names(hi) <- paste('deficit',humid,sep='_')
for(j in 1:nh)hi[[j]] <- 0
                      
bp <- gjamPriorTemplate(formula, xdata, ydata = treeYdata, lo = lo, hi=hi)
rl <- list(N = 10, r = 5) 
ml <- list(ng=1000, burnin=200, typeNames = 'CA', betaPrior = bp, reductList=rl)
out <- gjam(formula, xdata, ydata = treeYdata, modelList = ml)

sc  <- rep('grey',ncol(treeYdata))
sc[snames %in% hot] <- '#ff7f00'      # highlight informative priors
sc[snames %in% humid] <- '#1f78b4'
pl  <- list( GRIDPLOTS = T, specColor = sc)
gjamPlot(output = out, plotPars = pl)

The combination of lo and hi set the limits for posterior draws from the truncated multivariate normal distribution. The help page for gjamPriorTemplate provides an example with informative priors specified for individual predictor-response combinations.

species-factor levels that do not occur

There can be times when specific species-predictor combinations could be omitted. For example, some species may never occur on some soil types. If I want to estimate \(\beta\) conditioned on some elements being zero, I could do the following:

formula <- as.formula(~ temp + soil)

# find species-soil type combinations that occur less than 5 times
y0 <- treeYdata
y0[ y0 > 0 ] <- 1
soil <- rep( as.character(xdata$soil), ncol(y0) )
soil <- paste('soil', soil, sep='') # soil is a factor, so full name begins with variable name
spec <- rep( colnames(y0), each = nrow(y0) )
specBySoil <- tapply( as.vector(y0), list(spec, soil), sum )
wna  <- which( specBySoil < 5, arr.ind = T )                   # only if at least 5 observations

# assign NA to excluded coefficients using species-variable names
nh <- nrow(wna)
lo  <- vector('list', nh)
names(lo) <- paste( rownames(specBySoil)[wna[,1]], colnames(specBySoil)[wna[,2]], sep='_' )
hi <- lo
for(j in 1:nh)lo[[j]] <- NA
hi <- lo

# setup prior and fit model
bp  <- gjamPriorTemplate(formula, xdata, ydata = treeYdata, lo = lo, hi=hi)
rl  <- list(N = 10, r = 5) 
ml  <- list(ng=1000, burnin=200, typeNames = 'CA', betaPrior = bp, reductList=rl)
out <- gjam(formula, xdata, ydata = treeYdata, modelList = ml)

The values in out$parameters$betaMu will show these omitted values as NA.

sample effort in count data

Discrete count ('DA') data have effort associated with search time, plot area, and so forth. When the effort values are measured in units that result in large differences between data types sampled with different efforts, model fitting deteriorates. The modeling scale is \(W = Y/E\), where \(Y\) are counts, and \(E\) is effort. Units that make \(E\) large, can make \(W\) vanishingly small. When counts per unit effort span different orders of magnitude for different data types, consider changing the units on effort to bring them more in alignment with one another. Model predictions for ground beetles (pitfall traps) and small mammals (live traps) in our NEON analysis improved by shifting from trap-days to trap-months.

sample effort in composition data

Composition count ('CC') data have heterogenous effort due to different numbers of counts for each sample. For example, in microbiome data, the number of reads per sample can range from \(10^{2}\) to \(10^{6}\). Typically, the number of reads does not depend on total abundance. It is generally agreed that only relative differences are important. gjam knows that the effort in CC data is the total count for the sample, so effort does not need to be specified. Here is an example with simulated data:

f   <- gjamSimData( S = 8, typeNames = 'CC' )
ml  <- list( ng = 2000, burnin = 500, typeNames = f$typeNames )
out <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )
gjamPlot( output = out, plotPars = list( trueValues = f$trueValues) )

For comparison, here is an example with fractional composition, where there is no effort:

f   <- gjamSimData( S = 10, typeNames = 'FC' )
ml  <- list( ng = 2000, burnin = 500, typeNames = f$typeNames )
out <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )
gjamPlot( output = out, plotPars = list( trueValues = f$trueValues ) )

the partition in ordinal data

Ordinal count ('OC') data are collected where abundance must be evaluated rapidly or precise measurements are difficult. Because there is no absolute scale the partition must be inferred. Here is an example with 10 species:

f   <- gjamSimData( typeNames = 'OC' ) 
ml  <- list( ng = 2000, burnin = 500, typeNames = f$typeNames )
out <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )

A simple plot of the posterior mean values of cutMu shows the estimates with true values from simulation:

par( mfrow=c(1,1), bty = 'n' )
keep <- strsplit(colnames(out$parameters$cutMu),'C-') #only saved columns
keep <- matrix(as.numeric(unlist(keep)), ncol = 2, byrow = T)[,2]
plot( f$trueValues$cuts[,keep], out$parameters$cutMu, xlab = 'True partition',
      ylab = 'Estimated partition')

Here are plots:

pl  <- list(trueValues = f$trueValues)
gjamPlot(output = out, plotPars = pl)

categorical data

Categorical data have levels within groups. The levels are unordered. The columns in ydata that hold categorical responses are declared as typeNames = "CAT". In observation vector \(\mathbf{y}_{i}\) there are elements for one less than the number of factor levels. Suppose that observations are obtained on attributes of individual plants, each plant being an observation. The group leaf type might have four levels broadleaf decidious bd, needleleaf decidious nd, broadleaf evergreen be, and needleaf evergreen ne. A second group xylem anatomy might have three levels diffuse porous dp, ring porous rp, and tracheid tr. In both cases I assign the last class to be a reference class, other. Ten rows of the response matrix data might look like this:

##    leaf xylem
## 1    bd    dp
## 2    be    dp
## 3    bd    rp
## 4    nd other
## 5    be    rp
## 6    bd    rp
## 7 other    rp

This data.frame ydata becomes this response matrix y:

##      leaf_bd leaf_nd leaf_be leaf_other xylem_dp xylem_rp xylem_other
## [1,]       1       0       0          0        1        0           0
## [2,]       0       0       1          0        1        0           0
## [3,]       1       0       0          0        0        1           0
## [4,]       0       1       0          0        0        0           1
## [5,]       0       0       1          0        0        1           0
## [6,]       1       0       0          0        0        1           0
## [7,]       0       0       0          1        0        1           0

gjam expands the two groups into four and three columns in y, respectively. As for composition data there is one redundant column for each group. Here is an example with simulated data, having two categorical groups:

types <- c( 'CAT', 'CAT' )
f     <- gjamSimData( n = 3000, S = length(types), typeNames = types )
ml    <- list( ng = 1500, burnin = 500, typeNames = f$typeNames, PREDICTX = F )
out   <- gjam( f$formula, xdata = f$xdata, ydata = f$ydata, modelList = ml )
gjamPlot( out, plotPars = list(trueValues = f$trueValues, PLOTALLY = T) )

combinations of data types

One of the advantages of gjam is that it combines data of many types. Here is an example showing joint analysis of 12 species represented by five data types, specified by column:

types <- c('OC','OC','OC','OC','CC','CC','CC','CC','CC','CA','CA','PA','PA' ) 
f     <- gjamSimData( S = length(types), Q = 5, typeNames = types )
ml    <- list( ng = 2000, burnin = 500, typeNames = f$typeNames )
out   <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )
tmp   <- data.frame( f$typeNames, out$inputs$classBySpec[,1:10] )
print( tmp )

I have displayed the first 10 columns of classBySpec from the inputs of out, with their typeNames. The ordinal count ('OC') data occupy lower intervals. The width of each interval in OC data depends on the estimate of the partition in cutMu.

The composition count ('CC') data occupy a broader range of intervals. Because CC data are only relative, there is information on only \(S - 1\) species. One species is selected as other. The other class can be a collection of rare species (Clark et al. 2017).

Both continuous abundance ('CA') and presence-absence ('PA') data have two intervals. For CA data only the first interval is censored, the zeros (see above). For PA data both interval are censored; it is a multivariate probit.

Here are some plots for analysis of this model:

gjamPlot( output = out, plotPars = list(trueValues = f$trueValues, GRIDPLOTS = T) )

random effects

In addition to random effects on species used in dimension reduction, gjam allows for random groups. Examples could be observer effects for bird point counts or plot effects, where there are multiple observations per plot. Just like fixed effects, random effects should have replication. In other words, if I want to estimate an observer effect, I should have multiple observations for each observer. For plot effects, I want replication within plots.

To specify random groups, I provide the name of the column in xdata that holds the group indicator. This can be entered as an integer, a factor level, or a character variable.

modelList$random <- 'columnNameInXdata'

This column will be the basis for the random groups. Any rare groups (less than a few observations) will be aggregated into a single rareGroups category, again, to insure replication.

In the gjam output, here are objects related to random effects in output$parameters:

missing data, out-of-sample prediction

gjam identifies missing values in xdata and ydata and models them as part of the posterior distribution. These are located at $missing$xmiss and $missing$ymiss. The estimates for missing \(\mathbf{X}\) are $missing$xmissMu and $missing$xmissSe. The estimates for missing \(\mathbf{Y}\) are $missing$ymissMu and $missing$ymissSe.

To simulate missing data use nmiss to indicate the number of missing values. The actual value will be less than nmiss:

f <- gjamSimData(typeNames = 'OC', nmiss = 20)
which(is.na(f$xdata), arr.ind = T)

Note that missing values are assumed to occur in random rows and columns, but not in column one, which is the intercept and known. No further action is needed for model fitting, as gjam knows to treat these as missing data.

Out-of-sample prediction of \(\mathbf{Y}\) is not part of the posterior distribution. For model fitting, holdouts can be specified randomly in modelList with holdoutN (the number of plots to be held out at random) or with holdoutIndex (observation numbers, i.e., row numbers in x and y). The latter can be useful when a comparison of predictions is desired for different models using the same plots as holdouts.

When observations are held out, gjam provides out-of-sample prediction for both x[holdoutIndex,] and y[holdoutIndex,]. The holdouts are not included in the fitting of \(\boldsymbol{B}\),\(\boldsymbol{\Sigma}\), or \(\mathcal{P}\). For prediction of y[holdoutIndex,], the values of x[holdoutIndex,] are known, and sampling for w[holdoutIndex,] is done with multivariate normal distribution, without censoring. This is done because the censoring depends on y[holdoutIndex,], which taken to be unknown. This sample of w[holdoutIndex,] becomes a prediction for y[holdoutIndex,] using the partition (Figure 1a).

For inverse prediction of x[holdoutIndex,] the values of y[holdoutIndex,] are known. This represents the situation where a sample of the community is available, and the investigator would like to predict the environment of origin.

Here is an example with simulated data having both missing values and holdout observations:

f   <- gjamSimData( typeNames = 'CA', nmiss = 20 )
ml  <- list( ng = 2000, burnin = 500, typeNames = f$typeNames, holdoutN = 50 )
out <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )

par( mfrow=c(1,3), bty = 'n' )
xMu  <- out$prediction$xpredMu        # inverse prediction of x
xSd  <- out$prediction$xpredSd
yMu  <- out$prediction$ypredMu        # predicted y
hold <- out$modelList$holdoutIndex    # holdout observations (rows)

plot(out$inputs$xUnstand[hold,-1],xMu[hold,-1], cex=.2, xlab='True', ylab='Predictive mean')
title('holdouts in x')
abline( 0, 1, lty=2 )
plot(out$inputs$y[hold,], yMu[hold,], cex=.2, xlab='True', ylab='')
title('holdouts in y')
abline( 0, 1, lty=2 )

xmiss   <- out$missing$xmiss                              # locations of missing x
xmissMu <- out$missing$xmissMu
xmissSe <- out$missing$xmissSe
xmean   <- apply(f$xdata,2,mean,na.rm=T)[xmiss[,2]] # column means for missing values
plot(xmean, xmissMu, xlab='Variable mean', ylab='Missing estimate') #posterior estimates
segments(xmean, xmissMu - 1.96*xmissSe, xmean, xmissMu + 1.96*xmissSe)          #approx 95% CI
title('missing x')

Note that there are no ‘true’ values in the simulation for missing x–the last graph at right just shows estimates relative to mean values for respective variables.

prediction with heterogenous effort

Out-of-sample prediction can not only be done by holding out samples in gjam. It can also be done post-fitting, with the function gjamPredict. In this second case, a prediction grid is passed together with the fitted object generated by gjam. Here is an example with counts (DA) and continuous data (CA). I simulate, fit, and predict these data with heterogeneous sample effort:

sc   <- 3                                            # no. CA responses
sd   <- 10                                           # no. DA responses
tn   <- c( rep('CA',sc),rep('DA',sd) )               # combine CA and DA obs
S    <- length(tn)
n    <- 500
emat <- matrix( runif(n,.5,5), n, sd )                # simulated DA effort
eff  <- list( columns = c((sc+1):S), values = emat )
f    <- gjamSimData( n = n, typeNames = tn, effort = eff )
ml   <- list( ng = 2000, burnin = 500, typeNames = f$typeNames, effort = f$effort )
out  <- gjam( f$formula, f$xdata, f$ydata, modelList = ml )
gjamPredict( out, y2plot = colnames(f$ydata)[tn == 'DA'] ) # predict DA data

The prediction plot fits the data well, because it assumes the same effort. However, I might wish to predict data with a standard level of effort, say ‘1’. This new effort is taken in the same units as was used to fit the data, e.g., plot area, time observed, and so on. I use the same design matrix, but specify this new effort. Here I first predict the data with the actual effort, followed by the new effort of 1

cols <- c( '#1b9e77','#d95f02' )
new  <- list( xdata = f$xdata, effort=eff, nsim = 500 ) # effort unchanged 
p1   <- gjamPredict( output = out, newdata = new )

plot(f$y[,tn == 'DA'], p1$sdList$yMu[,tn == 'DA'], 
     xlab = 'Observed', ylab = 'Predicted', cex=.2, col = cols[1] )
abline( 0, 1, lty = 2 )

new$effort$values <- eff$values*0 + 1                 # predict for effort = 1
p2 <- gjamPredict( output = out, newdata = new )

points( f$y[,tn == 'DA'], p2$sdList$yMu[,tn == 'DA'], col= cols[2], cex=.2 )
legend( 'topleft', c('Actual effort', 'Effort = 1'), text.col = cols, bty='n' )

The orange dots show what the model would predict had effort on all observations been equal to 1.

conditional prediction

gjam can predict a subset of columns in y conditional on other columns using the function gjamPredict. Here is an example using the model fitted in the previous section. Consider model prediction in the case where the second plant species is absent, and the first species is at its mean value. In other words, for these plant species abundances, what is the effect of model predictors? I compare it to predictions where I first condition on the observed values for the first two plant species. I do not specify a new version of xdata, but rather include the columns of y to condition on:

cols <- c( '#1b9e77','#d95f02','#e7298a' )

condCols <- 1:3
pCols    <- c(1:S)[-condCols]

new <- list(ydataCond = f$y[,condCols], nsim=200)         # cond on observed CA data
p1  <- gjamPredict(output = out, newdata = new)$sdList$yMu[,pCols]

new    <- list(ydataCond = f$y[,condCols]*0, nsim=200)    # spp 1, 2 absent
p2     <- gjamPredict(output = out, newdata = new)$sdList$yMu[,pCols]

obs <- f$y[,pCols]
plot( jitter(obs), p2, xlab='Observed', ylab = 'Predicted', cex=.1, 
      ylim=range(c(p1,p2)), col = cols[1] )
points( jitter(obs), out$prediction$ypredMu[,pCols], col=cols[2], cex=.1)
points( jitter(obs), p1, col=cols[3], cex=.1)
abline( 0, 1, lty=2 )
legend( 'topleft', 
        c('Conditioned on absent spp 1:3', 'Unconditional', 'Conditioned on observed spp 1:3'), 
        text.col = cols, bty='n')

In the first case, I held the values for columns 1 and 2 at the values observed. This conditioning on observed values changes the predictions of other variables. In the second case, I conditioned on the specific values of y mentioned above. Both differ from the unconditional prediction.

When there are large covariances in the estimates of \(\Sigma\) the conditional predictions can differ dramatically from anything observed. In fact, if I condition on values of y that are well outside the data predictions will make no sense at all. Conversely, if covariances in \(\Sigma\) are near zero conditional distributions will not look much different from unconditional predictions. With dimension reduction we have only a crude estimate of \(\Sigma\), so conditional prediction was be judged accordingly.

Here is a second example, where I ask how knowledge of some species affects ability to predict others.

library(repmis)
d <- "https://github.com/jimclarkatduke/gjam/blob/master/forestTraits.RData?raw=True"
source_data(d)

xdata <- forestTraits$xdata                            # n X Q
ydata <- gjamReZero( forestTraits$treesDeZero )        # n X S
ydata <- ydata[,colnames(ydata) != 'other']
ydata <- gjamTrimY( ydata, 200, OTHER=F )$y

Here is a model fitted to the forest "DA" data:

ml      <- list( ng = 2500, burnin = 1000, typeNames = 'DA', 
                 PREDICTX = F )
formula <- as.formula(~ temp + deficit*moisture)
output  <- gjam( formula, xdata, ydata, modelList = ml )

# prediction
newdata <- list( nsim=100 )
tmp     <- gjamPredict( output, newdata=newdata )
full    <- tmp$sdList$yMu

Here I condition on half of the species and predict the other half:

cols    <- c( '#1b9e77','#d95f02' )
cnames  <- sample( colnames(ydata), S/2)
wc      <- match(cnames, colnames(ydata))
newdata <- list(ydataCond = ydata[,cnames], nsim=100)
tmp     <- gjamPredict(output, newdata = newdata)
condy   <- tmp$sdList$yMu

plot(ydata[,-wc], full[,-wc], ylim = c(range(ydata)), 
     xlab = 'Observed', ylab = 'Predicted', col = cols[1], cex = .3)
abline( 0, 1, lty=2 )
points( ydata[,-wc], condy[,-wc], col = cols[2], cex = .3)
legend( 'topleft', c('Unconditional', 'Conditional on half of species'), 
        text.col = cols[1:2], bty='n' )

Again, the conditional prediction can differ from the unconditional one, due to the covariances between species.

conditional parameters

The parameters used for conditional prediction can be obtained with the function gjamConditionalParameters. The coeffficients come from a partition of the response vector as \(\mathbf{w}' = [\mathbf{u}, \mathbf{v}]\) with coefficient and covariance matrices

\[ \mathbf{B} = \begin{bmatrix} \mathbf{B}_u \\ \mathbf{B}_v \\ \end{bmatrix}, \Sigma = \begin{bmatrix} \Sigma_{uu} & \Sigma_{uv} \\ \Sigma_{vu} & \Sigma_{vv} \\ \end{bmatrix} \]

Write the conditional distribution for responses in \(\mathbf{u}\) conditioned on values set for the others in \(\mathbf{v}\),

\[ \begin{aligned} \mathbf{u}|\mathbf{v} &\sim MVN \left( \mathbf{C}\mathbf{x} + \mathbf{A}\mathbf{v}, \mathbf{P} \right) \\ \mathbf{A} &= \Sigma_{uv} \Sigma^{-1}_{vv} \\ \mathbf{C} &= \mathbf{B}_u - \mathbf{A}\mathbf{B}_v \\ \mathbf{P} &= \Sigma_{uu} - \mathbf{A}\Sigma_{vu} \end{aligned} \] gjamConditionalParameters takes the fitted gjam object and returns these three matrices, each as a matrix of posterior means and standard errors and as a table with 95% credible intervals. The vector conditionOn holds the column names in ydata for the variables that I want to condition on, i.e., in the subvector \(\mathbf{v}\),

gjamConditionalParameters( output, conditionOn = c('acerSacc','betuAlle','faguGran','querRubr') )

Matrix \(\mathbf{A}\) is like a regression matrix on species. It is a \(S_u \times S_v\) matrix holding the effects of each response in conditionOn (\(\mathbf{v}\)) on the other columns (\(\mathbf{u}\)).

Matrix \(\mathbf{C}\) are regression coefficients on predictors in \(\mathbf{x}\), after extracting the contribution from other species. It is a \(S_u \times Q\) matrix.

Matrix \(\mathbf{P}\) is the \(S_u \times S_u\) conditional covariance.

presence-only data with effort

If I have a model for effort, then incidence data can have a likelihood, i.e., a probability assigned to observations that are aggregated into groups of known effort. I cannot model absence for the aggregate data unless I know how much effort was devoted to searching for it. Effort is widely known to have large spatio-temporal and taxonomic bias.

If I know the effort for a region, even in terms of a model (e.g., distance from universities and museums, from rivers or trails, numbers of a species already in collections), I can treat aggregate data as counts. If I do not know effort, out of desperation I might use the total numbers of all species counted in a region as a crude index of effort. The help page for function gjamPoints2Grid provides examples to aggregate incidence data with (x, y) locations to a lattice.

If effort is known I can supply a prediction grid predGrid for the known effort map and aggregate incidence to that map. I can then model the data are 'DA' data, specifying effort as in the example above.

If effort is unknown, I can model the data as composition data, 'CC'. Again, this is a desperate measure, because there are many reasons why even the total for all species at a lattice point might not represent relative effort.

big-S composition data: fungal endophytes

Microbiome data are often big-S, small-n; with thousands of response variables, columns in \(\mathbf{Y}\) (e.g., OTUs). They are also composition count ('CC') data, discrete counts, but not related to absolute abundance; they are meaningful in a relative sense. Because data only inform about relative abundance, there is information for only \(S - 1\) species. If there are thousands of OTUs, most of which are rare and thus not explained by the model, consider aggregating the many rare types into a single other class.

Fungal endophytes were sequenced on host tree seedlings (Hersh et al. 2016). In the data set fungEnd there is a compressed version of responses yDeZero containing OTU counts, a data.frame xdata containing predictors, and status, a vector of host responses, 0 for morbid, 1 for no signs of morbidity. Several histograms show the overwhelming numbers of zeros. Here we extract the data, stored in de-zeroed format, and generate some plots:

library(repmis)
d <- "https://github.com/jimclarkatduke/gjam/blob/master/fungEnd.RData?raw=True"
source_data(d)

xdata  <- fungEnd$xdata
otu    <- gjamReZero(fungEnd$yDeZero)
status <- fungEnd$status

par( mfrow=c(1,3), bty='n' )
hist( status, main='Host condition (morbid = 0)', ylab = 'Host obs' )
hist( otu, nclass=100, ylab = "Reads", main = "OTU's" )
nobs <- gjamTrimY( otu, minObs = 1, OTHER = F )$nobs
hist( nobs/ncol(otu), nclass=100, xlab = 'Fraction of observations', 
      ylab = 'Frequency', main='Incidence' )

The graph on the left shows that two thirds of host seedlings are in the morbid condition. The middle graph shows that an individual OTU may be counted > 40,000 times in an obseration, but the vast majority of observations (96%) are zero. The graph on the right shows the frequency for the fraction of observations in which each OTU occurs. The most common OTU occurs in 23% of observations, but most occur only once or twice.

The model will provide no information on the rarest taxa. Here we trim otu to include only OTUs that occur in > 100 observations. The rarest OTUs are aggregated into the last column of y with the column name other:

tmp <- gjamTrimY(otu, minObs = 100)
y   <- tmp$y
dim(fungEnd$y)               # all OTUs
dim(y)                       # trimmed data
tail(colnames(y))            # 'other' class added

The full response matrix includes the OTU composition counts and the host status in column 1:

ydata <- cbind(status, y)     # host status is also a response
S     <- ncol(ydata)
typeNames    <- rep('CC',S)   # composition count data
typeNames[1] <- 'PA'          # binary host status 

The interactions in the model involve two factors poly (two levels, polyculture vs monoculture) and host (eight factor levels, one for each host species). I assign acerRubr as the reference class for host,

xdata$host <- relevel(xdata$host,'acerRubr')

For this example we specify up to \(N = 20\) clusters with \(r = 3\) columns each. Here is an analysis of host seedling and polyculture effect on combined host morbidity status and the microbiome composition:

ml <- list( ng = 2000, burnin = 500, typeNames = typeNames, 
            reductList = list(r = 8, N = 15) )
output <- gjam(~ host*poly, xdata, ydata, modelList = ml)

Here is output:

S <- ncol(ydata)
specColor     <- rep('blue',S)
specColor[1]  <- '#b2182b'                 # highlight host status
fit <- gjamPlot( output, plotPars = list(specColor = specColor, GRIDPLOTS = T,
                                         SIGONLY = T) )
fit$eComs[1:5,]

Check the chains for convergence. Those coefficients that do not converge are poorly represented in factorial factor combinations.

Again, the low dimensional version of covariance \(\boldsymbol{\Sigma}\) is expected to perform best when there is structure in the data. The responses in matrix \(\mathbf{E}\), returned in fit$ematrix, classify OTUs in three main groups, contained in fit$eComs.

interactions and indirect effects

A plot of main effects, interactions, and indirect effects is used in this example to show contributions to host status, the response variable status. The effects on host status of on responses is available as a table with standard errors and credible intervals:

beta <- output$parameters$betaTable
ws   <- grep( 'status_',rownames(beta) )  # find coefficients for status
beta[ws,]

Following the intercept are rows showing main effects. These are followed by interaction terms. A quick visual of coefficients having credible intervals that exclude zero is here:

ws <- which(beta$sig95 == '*')
beta[ws,]

with - and + indicating negative and postive values.

Here I use the function gamIIE to create the object fit1, a list of these main, interaction, and indirect effects. I specify not to include the response variable other as an indirect effect on status, because we want to focus on the effects of microbes that have been assigned to known taxonomic groups.

I specify the values for main effects that are involved in the interactions between poly and host. Each factor has one less column in the design matrix x than factor levels. poly has two classes in xdata, one each for monoculture and polyculture, so there is one poly column in x. host has eight species in xdata, so there are seven columns in x. Recall that we assigned acerRubr to be the reference class, so there is no column for it in x.

The vector of predictor values xvector passed to gjamIIE has the same elements and names as columns in x. For this reason it is easiest to simply assign it a row in x, then change the values. The only values that influence interactions are those that are involved in interaction terms, as specified in formula.

For the following plots I copy the first row of x. In the first are main effects of all predictors on status. In the second plot is interactions with poly set to 1. In the third plot are indirect effects of the microbes:

x <- output$inputs$xUnstand
xvector <- x[1,]*0
names(xvector)  <- colnames(x)

xvector['hostfraxAmer'] <- 1
xvector['polypoly'] <- 1
fit1 <- gjamIIE(output, xvector, omitY = 'other')

par(mfrow=c(1,3), bty='n', oma = c(1,2,0,0), 
    mar = c(3,3,3,1), tcl = -0.5, mgp = c(3,1,0))
gjamIIEplot(fit1, response = 'status', effectMu = 'direct', 
            effectSd = 'direct', legLoc = 'bottomright', ylim=c(-10,10))
title('Direct effect by host')

gjamIIEplot(fit1, response = 'status', effectMu = 'int', effectSd = 'int',
            legLoc = 'topright', ylim=c(-10,10))
title('Interactions with polyculture')

gjamIIEplot(fit1, response = 'status', effectMu = 'ind', effectSd = 'ind',
            legLoc = 'topright', ylim=c(-5,5))
title('Indirect effect of microbiome')

The plot at left is the direct effect, which includes both the main effects plus interactions and plotted relative to the mean over all hosts. The interaction contribution at center is the effect of each host, when grown in polyculture (poly = ref1) and of polyculture when the host = 'fraxAmer.

The indirect effects bring with them the main effects and interaction effects on each microbial taxon. In this example the indirect effects are noisy, showing large 95% intervals.

trait response model (TRM)

There are \(n\) observations of \(M\) traits to be explained by \(Q - 1\) predictors in design matrix \(\mathbf{X}\). The Trait Response Model (TRM) in Clark 2016 is

\[\mathbf{w}_{i} \sim MVN(\mathbf{A}'\mathbf{x}_{i},\Omega)\]

where \(\mathbf{w}_{i}\) is a length-\(M\) vector of expected CWMM values, \(\mathbf{A}\) is the \(Q \times M\) matrix of coefficients, and \(\Omega\) is the \(M \times M\) residual covariance (diagram below). After describing the setup and model fitting I show how gjam summarizes the estimates and predictions.

Trait response model showing the sizes of matrices for a sample containing n observations, M traits, and Q predictors.

library(repmis)
d <- "https://github.com/jimclarkatduke/gjam/blob/master/forestTraits.RData?raw=True"
source_data(d)

xdata <- forestTraits$xdata                          # n X Q
types <- forestTraits$traitTypes                     # 12 trait types 
sbyt  <- forestTraits$specByTrait                    # S X 12
pbys  <- gjamReZero(forestTraits$treesDeZero)        # n X S
pbys  <- gjamTrimY(pbys,5)$y                         # at least 5 plots
sbyt  <- sbyt[match(colnames(pbys),rownames(sbyt)),] # trait matrix matches ydata
identical(rownames(sbyt),colnames(pbys))

table(sbyt$leaf)      # four levels

table(sbyt$xylem)     # diffuse/tracheid vs ring-porous

table(sbyt$repro)     # two levels

tmp         <- gjamSpec2Trait(pbys, sbyt, types)
tTypes      <- tmp$traitTypes                  # M = 15 values
u           <- tmp$plotByCWM                   # n X M
censor      <- tmp$censor                      # (0, 1) censoring, two-level CAT's
specByTrait <- tmp$specByTrait                 # S X M
M           <- ncol(u)
n           <- nrow(u)
cbind(colnames(u),tTypes)                      # M trait names and types

censorList    <- gjamCensorY(values = c(0,1), intervals = cbind( c(-Inf,0),c(1,Inf) ), 
                             y = u, whichcol = c(13:14))$censor

xdata$soil <- relevel(xdata$soil,'reference')

ml  <- list(ng = 3000, burnin = 500, typeNames = tTypes, holdoutN = 20,
            censor=censor, notStandard = c('u1','u2','u3'))
out <- gjam(~ temp + stdage + moisture*deficit + deficit*soil, 
                 xdata = xdata, ydata = u, modelList = ml)
tnames <- colnames(u)
sc <- rep('black', M)                                  # highlight types
wo <- which(tnames %in% c("leafN","leafP","SLA") )     # foliar traits
wf <- grep("leaf",tnames)                              # leaf habit
wc <- which(tnames %in% c("woodSG","diffuse","ring") ) # wood anatomy

sc[wc] <- 'brown'
sc[wf] <- 'darkblue'
sc[wo] <- 'darkgreen'

pl  <- list(GRIDPLOTS = TRUE, PLOTALLY = T, specColor = sc)
gjamPlot(output = out, plotPars = pl)
summary(out)

out$parameters$betaMu         # Q by M coefficient matrix alpha
out$parameters$betaStandXmu   # Q by M standardized for X
out$parameters$betaStandXWmu  # (Q-F) by M standardized for W/X, centered factors
out$parameters$betaTable      # QM by stats posterior summary
out$parameters$sigMu          # M by M covariance matrix omega
out$parameters$sigSe          # M by M covariance std errors

xdrydry <- xwetdry  <- out$inputs$xUnstand[1,]
xdrydry['moisture'] <- xdrydry['deficit'] <- -1
xwetdry['moisture'] <- 1
xwetdry['deficit']  <- -1

par(mfrow=c(2,2), bty='n', mar=c(1,3,1,1), oma = c(0,0,0,0), 
    mar = c(3,2,2,1), tcl = -0.5, mgp = c(3,1,0), family='')

fit1 <- gjamIIE(output = out, xvector = xdrydry)
fit2 <- gjamIIE(output = out, xvector = xwetdry)

gjamIIEplot(fit1, response = 'leafbroaddeciduous', 
            effectMu = c('main','int'), effectSd = c('main','int'), 
            legLoc = 'bottomleft', ylim=c(-.31,.3))
title('deciduous')
gjamIIEplot(fit1, response = 'leafneedleevergreen', 
            effectMu = c('main','int'), effectSd = c('main','int'), 
            legLoc = 'bottomleft', ylim=c(-.3,.3))
title('evergreen')

gjamIIEplot(fit2, response = 'leafbroaddeciduous', 
            effectMu = c('main','int'), effectSd = c('main','int'), 
            legLoc = 'bottomleft', ylim=c(-.3,.3))
gjamIIEplot(fit2, response = 'leafneedleevergreen', 
            effectMu = c('main','int'), effectSd = c('main','int'), 
            legLoc = 'bottomleft', ylim=c(-.3,.3))

xvector <- out$inputs$xUnstand[1,]
par(mfrow=c(2,1), bty='n', mar=c(1,1,1,1), oma = c(0,0,0,0), 
    mar = c(3,2,2,1), tcl = -0.5, mgp = c(3,1,0), family='')

omitY <- colnames(u)[colnames(u) != 'leafbroaddeciduous'] # omit all but deciduous

fit <- gjamIIE(out, xvector)
gjamIIEplot(fit, response = 'leafP', effectMu = c('main','ind'), 
            effectSd = c('main','ind'), legLoc = 'topright', ylim=c(-.6,.6))
title('foliar P')
gjamIIEplot(fit, response = 'leafN', effectMu = c('main','ind'), 
            effectSd = c('main','ind'), legLoc = 'bottomright', ylim=c(-.6,.6))
title('foliar N')

predictive trait model (PTM)

The PTM models species abundance data, then predicts traits. This approach has a number of advantages over TRM discussed in Clark 2016. The response is the \(n \times S\) matrix \(\mathbf{Y}\), which could be counts, biomass, and so forth.

tl  <- list(plotByTrait = u, traitTypes = tTypes, specByTrait = specByTrait)
rl  <- list(r = 8, N = 25)
ml  <- list(ng = 2000, burnin = 500, typeNames = 'CC', holdoutN = 20,
                  traitList = tl, reductList = rl)
out <- gjam(~ temp + stdage + deficit*soil, xdata = xdata, 
                     ydata = pbys, modelList = ml)

S  <- nrow(specByTrait)
sc <- rep('black', S)
wr <- which(specByTrait[,'ring'] == 1)                  # ring porous
wb <- which(specByTrait[,'leafneedleevergreen'] == 1)   # needle-leaf evergreen
wd <- which(specByTrait[,'leafneedledeciduous'] == 1)   # needle-leaf deciduous
ws <- which(specByTrait[,'shade'] >= 4)                 # shade tolerant
sc[wr] <- '#8c510a'
sc[ws] <- '#3288bd'
sc[wb] <- '#003c30'
sc[wd] <- '#80cdc1'

M  <- ncol(specByTrait)
tc <- rep('black', M)
names(tc) <- colnames(specByTrait)
tc[ 'ring' ] <- '#8c510a'
tc[ 'leafneedleevergreen' ] <- '#3288bd'
tc[ 'leafneedledeciduous' ] <- '#003c30'
tc[ 'shade' ] <- '#80cdc1'

par(family = '')
pl  <- list(GRIDPLOTS=T,
                  specColor = sc, traitColor = tc, ncluster = 5) 
fit <- gjamPlot(output = out, pl)

out$parameters$betaTraitXMu     # Q by M coefficient matrix alpha, standardized for X
out$parameters$betaTraitXWmu    # Q by M standardized for X/W
out$parameters$betaTraitXTable  # QM by stats posterior summary
out$parameters$betaTraitXWTable # QM by stats summary for X/W
out$parameters$varTraitMu       # M by M trait residual covariance, standardized for X
out$parameters$varTraitTable    # M^2 by stats summary, standardized for X/W

out$prediction$tMu[1:5,]     # n by M predictive means
out$prediction$tSe[1:5,]     # n by M predictive std errors

model summary

An observation consists of environmental variables and species attributes, \(\lbrace \mathbf{x}_{i}, \mathbf{y}_{i}\rbrace\), \(i = 1,..., n\). The vector \(\mathbf{x}_{i}\) contains predictors \(x_{iq}: q = 1,..., Q\). The vector \(\mathbf{y}_{i}\) contains attributes (responses), such as species abundance, presence-absence, and so forth, \(y_{is}: s = 1,..., S\). The effort \(E_{is}\) invested to obtain the observation of response \(s\) at location \(i\) can affect the observation. The combinations of continuous and discrete measurements in observed \(\mathbf{y}_{i}\) motivate the three elements of gjam.

A length-\(S\) vector \(\mathbf{w}_{i}\in{\Re}^S\) represents response \(\mathbf{y}_i\) in continuous space. This continuous space allows for the dependence structure with a covariance matrix. An element \(w_{is}\) can be known (e.g., continuous response \(y_{is}\)) or unknown (e.g., discrete responses).

A length-\(S\) vector of integers \(\mathbf{z}_{i}\) represents \(\mathbf{y}_i\) in discrete space. Each observed \(y_{is}\) is assigned to an interval \(z_{is} \in \{0,...,K_{is}\}\). The number of intervals \(K_{is}\) can differ between observations and between species, because each species can be observed in different ways.

The partition of continuous space at points \(p_{is,z} \in{\mathcal{P}}\) defines discrete intervals \(z_{is}\). Two values \((p_{is,k}, p_{is,k+1}]\) bound the \(k^{th}\) interval of \(s\) in observation \(i\). Intervals are contiguous and provide support over the real line \((-\infty, \infty)\). For discrete observations, \(k\) is a censored interval, and \(w_{is}\) is a latent variable. The set of censored intervals is \(\mathcal{C}\). The partition set \(\mathcal{P}\) can include both known (discrete counts, including composition data) and unknown (ordinal, categorical) points.

An observation \(y\) maps to \(w\),

\[y_{is} = \left \{ \begin{matrix} \ w_{is} & continuous\\ \ z_{is}, & w_{is} \in (p_{z_{is}}, p_{z_{is} + 1}] & discrete \end{matrix} \right.\]

Effort \(E_{is}\) affects the partition for discrete data. For the simple case where there is no error in the assignment of discrete intervals, \(\mathbf{z}_i\) is known, and the model for \(\mathbf{w}_i\) is

\[\mathbf{w}_i|\mathbf{x}_i, \mathbf{y}_i, \mathbf{E}_i \sim MVN(\boldsymbol{\mu}_i,\boldsymbol{\Sigma}) \times \prod_{s=1}^S\mathcal{I}_{is}\] \[\boldsymbol{\mu}_i = \mathbf{B}'\mathbf{x}_i\] \[\mathcal{I}_{is} = \prod_{k \in \mathcal{C}}I_{is,k}^{I(y_{is} = k)} (1 - I_{is,k})^{I(y_{is} \neq k)}\]

where \(I_{is} =I(p_{z_{is}} < w_{is} < p_{z_{is} + 1}]\), \(\mathcal{C}\) is the set of discrete intervals, \(\mathbf{B}\) is a \(Q \times S\) matrix of coefficients, and \(\boldsymbol{\Sigma}\) is a \(S \times S\) covariance matrix. There is a correlation matrix associated with \(\boldsymbol{\Sigma}\),

\[\mathbf{R}_{s,s'} = \frac{\boldsymbol{\Sigma}_{s,s'}}{\sqrt{\boldsymbol{\Sigma}_{s,s} \boldsymbol{\Sigma}_{s',s'}}}\]

For presence-absence (binary) data, \(\mathbf{p}_{is} = (-\infty, 0, \infty)\). This is equivalent to Chib and Greenberg’s (2008) model, which could be written \(\mathcal{I}_{is} = I(w_{is} > 0)^{y_{is}}I(w_{is} \leqslant 0)^{1 - y_{is}}\).

For a continous variable with point mass at zero, continuous abundance, this is a multivariate Tobit model, with \(\mathcal{I}_{is} = I(w_{is} = y_{is})^{I(y_{is} > 0)}I(w_{is} \leqslant 0)^{I(y_{is} = 0)}\). This is the same partition used for the probit model, the difference being that the positive values in the Tobit are uncensored.

Categorical responses fit within the same framework. Each categorical response occupies as many columns in \(\mathbf{Y}\) as there are independent levels in response \(s\), levels being \(k = 1,..., K_{s}-1\). For example, if randomly sampled plots are scored by one of five cover types, then there are four columns in \(\mathbf{Y}\) for the response \(s\). The four columns can have at most one \(1\). If all four columns are \(0\), then the reference level is observed. The observed level has the largest value of \(w_{is,k}\) (Table 1). This is similar to Zhang et al.’s (2008) model for categorical data.

For ordinal counts gjam is Lawrence et al.’s (2008) model having the partition \(\mathbf{p}_{is} = (-\infty, 0, p_{is,2}, p_{is,3},..., p_{is,K}, \infty)\), where all but the first two and the last elements must be inferred. The partition must be inferred, because the ordinal scale is only relative.

Like categorical data, composition data have one reference class. For this and other discrete count data, the partition for observation \(i\) can be defined to account for sample effort (next section).

more on parameter dimensions

Unlike a univariate model that has one \(Y\) per observation or multivariate models where all \(Y\)’s have the same dimension, gjam has \(Y\)’s on multiple dimensions. So there are two sets of dimensions to consider, the dimensions for the \(X\)’s and those for the \(Y\)’s. To avoid more notation I just refer to the dimension of a coefficient in the unstandardized coefficient matrix \(\mathbf{B}_u\) as \(W/X\) (Table 2). Except for responses that have no dimension (presence-absence, ordinal, categorical) \(W\) has the same dimension as \(Y\) (or \(Y/E\), if effort is specified). gjam returns unstandardized coefficients in tables parameters$betaMuUn, parameters$betaSeUn, and in MCMC chains$bgibbsUn, because they are not prone to be misinterpreted by a user who has supplied unstandardized predictors in xdata. I refer to the coefficient matrix as \(\mathbf{B}_u\) and corresponding unstandardized design matrix as \(\mathbf{X}_u\). Each row of $chains$bgibbsUn is one draw from the \(Q \times S\) matrix \(\mathbf{B}_u\).

Models are commonly fitted to covariates \(X\) that are ‘standardized’ for mean and variance. Standardization can stabilize posterior simulation. It is desirable when coefficients are needed in standard deviations. Inside gjam design matrix \(\mathbf{X}\), and thus \(\mathbf{B}\), are standardized, thus having dimension \(W\), not \(W/X\). Of course, for variables in xdata supplied in standarized form, \(\mathbf{B} = \mathbf{B}_u\). See Algorithm summary. Variables returned by gjam, including MCMC chains in $chains$bgibbs and posterior means and standard errors $parameters$betaMu and $parameters$betaSe are standardized.

The third set of scales comes from the correlation scale for the \(W\)’s. The correlation scale can be useful when considering responses that have different dimensions. In addition to \(\mathbf{B}_u\) I provide parameters on the correlation scale. This correlation scale is \(\mathbf{B}_r = \mathbf{B} \mathbf{D}^{-1/2}\), where \(\mathbf{D} = diag(\boldsymbol{\Sigma})\). If \(\mathbf{X}\) is standardized, \(\mathbf{B}_r\) is dimensionless. The MCMC chains in $chains$bFacGibbs and the estimates in $parameterTable$fBetaMu and $parameterTable$fBetaSd are standardized for \(X\) (standard deviation scale) and for \(W\) (correlation scale). They are dimensionless.

For sensitivity over all species and for comparisons between predictors I provide standardized sensitivity in a length-\(Q\) vector \(\mathbf{f}\). The sensitivity matrix to \(X\) across the full model is given by

\[\mathbf{F} = \mathbf{B}\boldsymbol{\Sigma}^{-1} \mathbf{B}'\]

Note that \(\mathbf{F}\) takes \(\mathbf{B}_r\), not \(\mathbf{B}\), and not \(\mathbf{B}_u\). The sensitivity matrix \(\mathbf{F}\) is $parameters$fmatrix. The sensitivity vector \(\mathbf{f} = diag( \mathbf{F} )\) is vector $parameters$fMu. Details are given in Clark et al. (2017).

more on factors in \(\mathbf{X}\)

The coefficient matrix \(\boldsymbol{\tilde{\mathbf{B}}}\) is useful when there are factors in the model. Factor treatment in gjam follows the convention where a reference level is taken as the overall intercept, and remaining coefficients are added to the intercept. The reference level can be controlled with the R function relevel. This approach makes sense in the ANOVA context, where an experiment has a control level to which other treatment levels are to be compared. A ‘significant’ level is different from the reference (e.g., control), but we are not told about its relationship to other levels. The coefficients in $parameters$betaMu and $parameters$betaSd are reported this way. Should it be needed, the contrasts matrix for this design is returned as a list for all factors in the model as $inputs$factorBeta$contrast and as a single contrasts matrix for the entire model as $inputs$factorBeta$eCont.

This standard structure is not the best way to compare factors in many ecological data sets, where a factor might represent soil type, cover type, fishing gear, and so on. In all of these cases there is no ‘control’ treatment. Here it is more useful to know how all levels relate to the mean taken across factor levels.

To provide more informative comparisons across factor levels and species I introduce a \(Q_1 \times S\) recontrast matrix that translates \(\mathbf{B}\), with intercept, to all factor levels, without intercept. Consider a three-level factor with levels a, b, c, the first being the reference class. There is a matrix \(\mathbf{G}\)

##   intercept  b  c
## a         1 -1 -1
## b         1  1  0
## c         1  0  1

and a matrix \(\mathbf{H}\),

##            a b c
## intercept  1 0 0
## b         -1 1 0
## c         -1 0 1

With \(\mathbf{L'} = \mathbf{G}^{-1}\), the recontrasted coefficients are

\[\mathbf{\tilde{B}} = \mathbf{L}\mathbf{B}\]

The rows of \(\mathbf{\tilde{B}}\) correspond to all factor levels. The intercept does not appear, because it has been distributed across factor levels. The corresponding design matrix is

\[\mathbf{\tilde{X}} = \mathbf{X}\mathbf{H}\]

If there are multiple factors then \(Q_1 > Q\), because the intercept expands to the reference classes for each factor. \(\mathbf{L}\) is provided as $inputs$factorBeta$lCont.

With factors, the sensitivity matrix reported in $parameters$fmatrixis

\[\mathbf{F} = \mathbf{\tilde{B}}\boldsymbol{\Sigma}^{-1} \mathbf{\tilde{B}'}\]

The response matrix in $parameters$ematrix is

\[\mathbf{E} = \mathbf{\tilde{B}} \mathbf{\tilde{V}} \mathbf{\tilde{B}'}\]

where \(\mathbf{\tilde{V}} = cov \left(\mathbf{X} \mathbf{H} \right)\).