# doc ----
#' Phylogenetic Generalized Linear Mixed Model for Community Data
#'
#' This function performs Generalized Linear Mixed Models for binary, count,
#' and continuous data, estimating regression coefficients with
#' approximate standard errors. It is specifically designed for community data
#' in which species occur within multiple sites (locations).
#' A Bayesian version of PGLMM uses the package \code{INLA},
#' which is not available on CRAN yet. If you wish to use this option,
#' you must first install \code{INLA} from \url{https://www.r-inla.org/} by running
#' \code{install.packages('INLA', repos='https://www.math.ntnu.no/inla/R/stable')} in R.
#'
#'
#' For Gaussian data, \code{pglmm} analyzes the phylogenetic linear mixed model
#'
#' \deqn{Y = \beta_0 + \beta_1x + b_0 + b_1x}{Y = beta_0 + beta_1x + b_0 + b_1x}
#' \deqn{b_0 ~ Gaussian(0, \sigma_0^2I_{sp})}{b_0 ~ Gaussian(0, sigma_0^2I_(sp))}
#' \deqn{b_1 ~ Gaussian(0, \sigma_0^2V_{sp})}{b_0 ~ Gaussian(0, sigma_0^2V_(sp))}
#' \deqn{\eta ~ Gaussian(0,\sigma^2)}{e ~ Gaussian(0,sigma^2)}
#'
#' where \eqn{\beta_0}{beta_0} and \eqn{\beta_1}{beta_1} are fixed
#' effects, and \eqn{V_{sp}}{V_(sp)} is a variance-covariance matrix
#' derived from a phylogeny (typically under the assumption of
#' Brownian motion evolution). Here, the variation in the mean
#' (intercept) for each species is given by the random effect
#' \eqn{b_0}{b_0} that is assumed to be independent among
#' species. Variation in species' responses to predictor variable
#' \eqn{x}{x} is given by a random effect \eqn{b_0}{b_0} that is
#' assumed to depend on the phylogenetic relatedness among species
#' given by \eqn{V_{sp}}{V_(sp)}; if species are closely related,
#' their specific responses to \eqn{x}{x} will be similar. This
#' particular model would be specified as
#'
#' \code{z <- pglmm(Y ~ X + (1|sp__), data = data, family = "gaussian", cov_ranef = list(sp = phy))}
#'
#' Or you can prepare the random terms manually (not recommended for simple models but may be necessary for complex models):
#'
#' \code{re.1 <- list(1, sp = dat$sp, covar = diag(nspp))}
#'
#' \code{re.2 <- list(dat$X, sp = dat$sp, covar = Vsp)}
#'
#' \code{z <- pglmm(Y ~ X, data = data, family = "gaussian", random.effects = list(re.1, re.2))}
#'
#' The covariance matrix covar is standardized to have its determinant
#' equal to 1. This in effect standardizes the interpretation of the
#' scalar \eqn{\sigma^2}{sigma^2}. Although mathematically this is
#' not required, it is a very good idea to standardize the predictor
#' (independent) variables to have mean 0 and variance 1. This will
#' make the function more robust and improve the interpretation of the
#' regression coefficients. For categorical (factor) predictor
#' variables, you will need to construct 0-1 dummy variables, and
#' these should not be standardized (for obvious reasons).
#'
#' For binary generalized linear mixed models (\code{family =
#' 'binomial'}), the function estimates parameters for the model of
#' the form, for example,
#'
#' \deqn{y = \beta_0 + \beta_1x + b_0 + b_1x}{y = beta_0 + beta_1x + b_0 + b_1x}
#' \deqn{Y = logit^{-1}(y)}{Y = logit^(-1)(y)}
#' \deqn{b_0 ~ Gaussian(0, \sigma_0^2I_{sp})}{b_0 ~ Gaussian(0, sigma_0^2I_(sp))}
#' \deqn{b_1 ~ Gaussian(0, \sigma_0^2V_{sp})}{b_0 ~ Gaussian(0, sigma_0^2V_(sp))}
#'
#' where \eqn{\beta_0}{beta_0} and \eqn{\beta_1}{beta_1} are fixed
#' effects, and \eqn{V_{sp}}{V_(sp)} is a variance-covariance matrix
#' derived from a phylogeny (typically under the assumption of
#' Brownian motion evolution).
#'
#' \code{z <- pglmm(Y ~ X + (1|sp__), data = data, family = "binomial", cov_ranef = list(sp = phy))}
#'
#' As with the linear mixed model, it is a very good idea to
#' standardize the predictor (independent) variables to have mean 0
#' and variance 1. This will make the function more robust and improve
#' the interpretation of the regression coefficients.
#'
#' @param formula A two-sided linear formula object describing the
#' mixed effects of the model.
#'
#' To specify that a random term should have phylogenetic covariance matrix along
#' with non-phylogenetic one, add \code{__} (two underscores) at the end of the group variable;
#' e.g., \code{+ (1 | sp__)} will construct two random terms,
#' one with phylogenetic covariance matrix and another with non-phylogenetic (identity) matrix.
#' In contrast, \code{__} in the nested terms (below) will only create a phylogenetic covariance matrix.
#' Nested random terms have the general form \code{(1|sp__@site__)} which represents
#' phylogenetically related species nested within correlated sites.
#' This form can be used for bipartite questions. For example, species could be
#' phylogenetically related pollinators and sites could be phylogenetically related plants, leading to
#' the random effect `(1|insects__@plants__)`. If more than one phylogeny is used, remember to add
#' all to the argument `cov_ranef = list(insects = insect_phylo, plants = plant_phylo)`. Phylogenetic correlations can
#' be dropped by removing the \code{__} underscores. Thus, the form \code{(1|sp@site__)} excludes the phylogenetic
#' correlations among species, while the form \code{(1|sp__@site)} excludes the correlations among sites.
#'
#' Note that correlated random terms are not allowed. For example,
#' \code{(x|g)} will be the same as \code{(0 + x|g)} in the \code{lme4::lmer} syntax. However,
#' \code{(x1 + x2|g)} won't work, so instead use \code{(x1|g) + (x2|g)}.
#' @param data A \code{\link{data.frame}} containing the variables named in formula.
#' @param family Either "gaussian" for a Linear Mixed Model, or
#' "binomial" or "poisson" for Generalized Linear Mixed Models.
#' "family" should be specified as a character string (i.e., quoted). For binomial and
#' Poisson data, we use the canonical logit and log link functions, respectively.
#' Binomial data can be either presence/absence, or a two-column array of 'successes' and 'failures'.
#' For both binomial and Poisson data, we add an observation-level
#' random term by default via \code{add.obs.re = TRUE}. If \code{bayes = TRUE} there are
#' two additional families available: "zeroinflated.binomial", and "zeroinflated.poisson",
#' which add a zero inflation parameter; this parameter gives the probability that the response is
#' a zero. The rest of the parameters of the model then reflect the "non-zero" part part
#' of the model. Note that "zeroinflated.binomial" only makes sense for success/failure
#' response data.
#' @param cov_ranef A named list of covariance matrices of random terms. The names should be the
#' group variables that are used as random terms with specified covariance matrices
#' (without the two underscores, e.g. \code{list(sp = tree1, site = tree2)}). The actual object
#' can be either a phylogeny with class "phylo" or a prepared covariance matrix. If it is a phylogeny,
#' `pglmm` will prune it and then convert it to a covariance matrix assuming Brownian motion evolution.
#' `pglmm` will also standardize all covariance matrices to have determinant of one. Group variables
#' will be converted to factors and all covariance matrices will be rearranged so that rows and
#' columns are in the same order as the levels of their corresponding group variables.
#' @param random.effects Optional pre-build list of random effects. If \code{NULL} (the default),
#' the function \code{\link{prep_dat_pglmm}} will prepare the random effects for you from the information
#' in \code{formula}, \code{data}, and \code{cov_ranef}. \code{random.effect} allows a list of
#' pre-generated random effects terms to increase flexibility; for example, this makes it
#' possible to construct models with both phylogenetic correlation and spatio-temporal autocorrelation.
#' In preparing \code{random.effect}, make sure that the orders of rows and columns of
#' covariance matrices in the list are the same as their corresponding group variables
#' in the data. Also, this should be _a list of lists_, e.g.
#' `random.effects = list(re1 = list(matrix_a), re2 = list(1, sp = sp, covar = Vsp))`.
#' @param REML Whether REML or ML is used for model fitting the random effects. Ignored if
#' \code{bayes = TRUE}.
#' @param optimizer nelder-mead-nlopt (default), bobyqa, Nelder-Mead, or subplex.
#' Nelder-Mead is from the stats package and the other optimizers are from the nloptr package.
#' Ignored if \code{bayes = TRUE}.
#' @param repulsion When there are nested random terms specified, \code{repulsion = FALSE} tests
#' for phylogenetic underdispersion while \code{repulsion = FALSE} tests for overdispersion.
#' This argument is a logical vector of length either 1 or >1.
#' If its length is 1, then all covariance matrices in nested terms will be either
#' inverted (overdispersion) or not. If its length is >1, then you can select
#' which covariance matrix in the nested terms to be inverted. Make sure to get
#' the length right: for all the terms with \code{@}, count the number of "__"
#' to determine the length of repulsion. For example, \code{sp__@site} and \code{sp@site__}
#' will each require one element of \code{repulsion}, while \code{sp__@site__} will take two
#' elements (repulsion for sp and repulsion for site). Therefore, if your nested terms are
#' \code{(1|sp__@site) + (1|sp@site__) + (1|sp__@site__)}, then you should set the
#' repulsion to be something like \code{c(TRUE, FALSE, TRUE, TRUE)} (length of 4).
#' @param add.obs.re Whether to add an observation-level random term for binomial or Poisson
#' distributions. Normally it would be a good idea to add this to account for overdispersion,
#' so \code{add.obs.re = TRUE} by default.
#' @param verbose If \code{TRUE}, the model deviance and running
#' estimates of \code{s2} and \code{B} are plotted each iteration
#' during optimization.
#' @param cpp Whether to use C++ function for optim. Default is TRUE. Ignored if \code{bayes = TRUE}.
#' @param bayes Whether to fit a Bayesian version of the PGLMM using \code{r-inla}.
#' @param s2.init An array of initial estimates of s2 for each random
#' effect that scales the variance. If s2.init is not provided for
#' \code{family="gaussian"}, these are estimated using \code{\link{lm}} assuming
#' no phylogenetic signal. A better approach might be to run \code{link[lme4:lmer]{lmer}}
#' and use the output random effects for \code{s2.init}. If \code{s2.init} is not
#' provided for \code{family = "binomial"}, these are set to 0.25.
#' @param B.init Initial estimates of \eqn{B}{B}, a matrix containing
#' regression coefficients in the model for the fixed effects. This
#' matrix must have \code{dim(B.init) = c(p + 1, 1)}, where \code{p} is the
#' number of predictor (independent) variables; the first element of
#' \code{B} corresponds to the intercept, and the remaining elements
#' correspond in order to the predictor (independent) variables in the
#' formula. If \code{B.init} is not provided, these are estimated
#' using \code{\link{lm}} or \code{\link{glm}} assuming no phylogenetic signal.
#' A better approach might be to run \code{\link[lme4:lmer]{lmer}} and use the
#' output fixed effects for \code{B.init}. When \code{bayes = TRUE}, initial values are estimated
#' using the maximum likelihood fit unless \code{ML.init = FALSE}, in
#' which case the default \code{INLA} initial values will be used.
#' @param reltol A control parameter dictating the relative tolerance
#' for convergence in the optimization; see \code{\link{optim}}.
#' @param maxit A control parameter dictating the maximum number of
#' iterations in the optimization; see \code{\link{optim}}.
#' @param tol.pql A control parameter dictating the tolerance for
#' convergence in the PQL estimates of the mean components of the
#' GLMM. Ignored if \code{family = "gaussian"} or \code{bayes = TRUE}.
#' @param maxit.pql A control parameter dictating the maximum number
#' of iterations in the PQL estimates of the mean components of the
#' GLMM. Ignored if \code{family = "gaussian"} or \code{bayes = TRUE}.
#' @param marginal.summ Summary statistic to use for the estimate of coefficients when
#' doing a Bayesian PGLMM (when \code{bayes = TRUE}). Options are: "mean",
#' "median", or "mode", referring to different characterizations of the central
#' tendency of the Bayesian posterior marginal distributions. Ignored if \code{bayes = FALSE}.
#' @param calc.DIC Should the Deviance Information Criterion be calculated and returned
#' when doing a Bayesian PGLMM? Ignored if \code{bayes = FALSE}.
#' @param calc.WAIC Should the WAIC be calculated and returned
#' when doing a Bayesian PGLMM? Ignored if \code{bayes = FALSE}.
#' @param prior Which type of default prior should be used by \code{pglmm}?
#' Only used if \code{bayes = TRUE}. There are currently four options:
#' "inla.default", which uses the default \code{INLA} priors; "pc.prior.auto", which uses a
#' complexity penalizing prior (as described in
#' \href{https://arxiv.org/abs/1403.4630v3}{Simpson et al. (2017)}) designed to automatically
#' choose good parameters (only available for gaussian and binomial responses); "pc.prior", which
#' allows the user to set custom parameters on the "pc.prior" prior, using the \code{prior_alpha}
#' and \code{prior_mu} parameters (Run \code{INLA::inla.doc("pc.prec")} for details on these
#' parameters); and "uninformative", which sets a very uninformative prior
#' (nearly uniform) by using a very flat exponential distribution. The last option is generally
#' not recommended but may in some cases give estimates closer to the maximum likelihood estimates.
#' "pc.prior.auto" is only implemented for \code{family = "gaussian"} and \code{family = "binomial"}
#' currently.
#' @param prior_alpha Only used if \code{bayes = TRUE} and \code{prior = "pc.prior"}, in
#' which case it sets the alpha parameter of \code{INLA}'s complexity penalizing prior for the
#' random effects. The prior is an exponential distribution where prob(sd > mu) = alpha,
#' where sd is the standard deviation of the random effect.
#' @param prior_mu Only used if \code{bayes = TRUE} and \code{prior = "pc.prior"}, in
#' which case it sets the mu parameter of \code{INLA}'s complexity penalizing prior for the
#' random effects. The prior is an exponential distribution where prob(sd > mu) = alpha,
#' where sd is the standard deviation of the random effect.
#' @param ML.init Only relevant if \code{bayes = TRUE}. Should maximum
#' likelihood estimates be calculated and used as initial values for
#' the Bayesian model fit? Sometimes this can be helpful, but it may not help; thus,
#' we set the default to \code{FALSE}. Also, it
#' does not work with the zero-inflated families.
#' @param tree A phylogeny for column sp, with "phylo" class, or a covariance matrix for sp.
#' Make sure to have all species in the matrix; if the matrix is not standardized,
#' (i.e., det(tree) != 1), `pglmm` will try to standardize it for you.
#' No longer used: keep here for compatibility.
#' @param tree_site A second phylogeny for "site". This is required only if the
#' site column contains species instead of sites. This can be used for bipartitie
#' questions; tree_site can also be a covariance matrix. Make sure to have all sites
#' in the matrix; if the matrix is not standardized (i.e., det(tree_site) != 1),
#' pglmm` will try to standardize it for you. No longer used: keep here for compatibility.
#' @param sp No longer used: keep here for compatibility.
#' @param site No longer used: keep here for compatibility.
#' @param bayes_options Additional options to pass to INLA for if \code{bayes = TRUE}. A named list where the names
#' correspond to parameters in the \code{inla} function. One special option is \code{diagonal}: if an element in
#' the options list is names \code{diagonal} this tells \code{INLA} to add its value to the diagonal of the random effects
#' precision matrices. This can help with numerical stability if the model is ill-conditioned (if you get a lot of warnings,
#' try setting this to \code{list(diagonal = 1e-4)}).
#' @param bayes_nested_matrix_as_list For `bayes = TRUE`, prepare the nested terms as a list of length of 4 as the old way?
#' @return An object (list) of class \code{communityPGLMM} with the following elements:
#' \item{formula}{the formula for fixed effects}
#' \item{formula_original}{the formula for both fixed effects and random effects}
#' \item{data}{the dataset}
#' \item{family}{\code{gaussian}, \code{binomial}, or \code{poisson} depending on the model fit}
#' \item{random.effects}{the list of random effects}
#' \item{B}{estimates of the regression coefficients}
#' \item{B.se}{approximate standard errors of the fixed effects regression coefficients.
#' This is set to NULL if \code{bayes = TRUE}.}
#' \item{B.ci}{approximate Bayesian credible interval of the fixed effects regression coefficients.
#' This is set to NULL if \code{bayes = FALSE}}
#' \item{B.cov}{approximate covariance matrix for the fixed effects regression coefficients}
#' \item{B.zscore}{approximate Z scores for the fixed effects regression coefficients. This is set to NULL if \code{bayes = TRUE}}
#' \item{B.pvalue}{approximate tests for the fixed effects regression coefficients being different from zero. This is set to NULL if \code{bayes = TRUE}}
#' \item{ss}{standard deviations of the random effects for the covariance matrix \eqn{\sigma^2V}{sigma^2 V} for each random effect in order. For the linear mixed model, the residual variance is listed last.}
#' \item{s2r}{random effects variances for non-nested random effects}
#' \item{s2n}{random effects variances for nested random effects}
#' \item{s2resid}{for linear mixed models, the residual variance}
#' \item{s2r.ci}{Bayesian credible interval for random effects variances for non-nested random effects.
#' This is set to NULL if \code{bayes = FALSE}}
#' \item{s2n.ci}{Bayesian credible interval for random effects variances for nested random effects.
#' This is set to NULL if \code{bayes = FALSE}}
#' \item{s2resid.ci}{Bayesian credible interval for linear mixed models, the residual variance.
#' This is set to NULL if \code{bayes = FALSE}}
#' \item{logLik}{for linear mixed models, the log-likelihood for either the restricted likelihood (\code{REML=TRUE}) or the overall likelihood (\code{REML=FALSE}). This is set to NULL for generalized linear mixed models. If \code{bayes = TRUE}, this is the marginal log-likelihood}
#' \item{AIC}{for linear mixed models, the AIC for either the restricted likelihood (\code{REML = TRUE}) or the overall likelihood (\code{REML = FALSE}). This is set to NULL for generalised linear mixed models}
#' \item{BIC}{for linear mixed models, the BIC for either the restricted likelihood (\code{REML = TRUE}) or the overall likelihood (\code{REML = FALSE}). This is set to NULL for generalised linear mixed models}
#' \item{DIC}{for Bayesian PGLMM, this is the Deviance Information Criterion metric of model fit. This is set to NULL if \code{bayes = FALSE}.}
#' \item{REML}{whether or not REML is used (\code{TRUE} or \code{FALSE}).}
#' \item{bayes}{whether or not a Bayesian model was fit.}
#' \item{marginal.summ}{The specified summary statistic used to summarize the Bayesian marginal distributions.
#' Only present if \code{bayes = TRUE}}
#' \item{s2.init}{the user-provided initial estimates of \code{s2}}
#' \item{B.init}{the user-provided initial estimates of \code{B}}
#' \item{Y}{the response (dependent) variable returned in matrix form}
#' \item{X}{the predictor (independent) variables returned in matrix form (including 1s in the first column)}
#' \item{H}{the residuals. For linear mixed models, this does not account for random terms,
# i.e. it is similar to \code{Y - predict(merMod, re.form = NA)} for models fitted with lme4.
#' To get residuals after accounting for both fixed and random terms, use \code{residuals()}.
#' For the generalized linear mixed model, these are the predicted residuals in the
#' logit -1 space.}
#' \item{iV}{the inverse of the covariance matrix for the entire system (of dimension (`nsp` * `nsite`)
#' by (`nsp` * `nsite`)). This is NULL if \code{bayes = TRUE}.}
#' \item{mu}{predicted mean values for the generalized linear mixed model (i.e., similar to \code{fitted(merMod)}).
#' Set to NULL for linear mixed models, for which we can use [fitted()].}
#' \item{nested}{matrices used to construct the nested design matrix. This is set to NULL if \code{bayes = TRUE}}
#' \item{Zt}{the design matrix for random effects. This is set to NULL if \code{bayes = TRUE}}
#' \item{St}{diagonal matrix that maps the random effects variances onto the design matrix}
#' \item{convcode}{the convergence code provided by \code{\link{optim}}. This is set to NULL if \code{bayes = TRUE}}
#' \item{niter}{number of iterations performed by \code{\link{optim}}. This is set to NULL if \code{bayes = TRUE}}
#' \item{inla.model}{Model object fit by underlying \code{inla} function. Only returned
#' if \code{bayes = TRUE}}
#' @author Anthony R. Ives, Daijiang Li, Russell Dinnage
#' @references Ives, A. R. and M. R. Helmus. 2011. Generalized linear
#' mixed models for phylogenetic analyses of community
#' structure. Ecological Monographs 81:511-525.
#' @references Ives A. R. 2018. Mixed and phylogenetic models: a conceptual introduction to correlated data.
#' https://leanpub.com/correlateddata.
#' @references Rafferty, N. E., and A. R. Ives. 2013. Phylogenetic
#' trait-based analyses of ecological networks. Ecology 94:2321-2333.
#' @references Simpson, Daniel, et al. 2017. Penalising model component complexity:
#' A principled, practical approach to constructing priors.
#' Statistical science 32(1): 1-28.
#' @references Li, D., Ives, A. R., & Waller, D. M. 2017.
#' Can functional traits account for phylogenetic signal in community composition?
#' New Phytologist, 214(2), 607-618.
#' @rdname pglmm
#' @export
#' @examples
#' ## Structure of examples:
#' # First, a (brief) description of model types, and how they are specified
#' # - these are *not* to be run 'as-is'; they show how models should be organised
#' # Second, a run-through of how to simulate, and then analyse, data
#' # - these *are* to be run 'as-is'; they show how to format and work with data
#'
#' \donttest{
#' #############################################
#' ### Brief summary of models and their use ###
#' #############################################
#' ## Model structures from Ives & Helmus (2011)
#' if(FALSE){
#' # dat = data set for regression (note: must have a column "sp" and a column "site")
#' # phy = phylogeny of class "phylo"
#' # repulsion = to test phylogenetic repulsion or not
#'
#' # Model 1 (Eq. 1)
#' z <- pglmm(freq ~ sp + (1|site) + (1|sp__@site), data = dat, family = "binomial",
#' cov_ranef = list(sp = phy), REML = TRUE, verbose = TRUE, s2.init = .1)
#'
#' # Model 2 (Eq. 2)
#' z <- pglmm(freq ~ sp + X + (1|site) + (X|sp__), data = dat, family = "binomial",
#' cov_ranef = list(sp = phy), REML = TRUE, verbose = TRUE, s2.init = .1)
#'
#' # Model 3 (Eq. 3)
#' z <- pglmm(freq ~ sp*X + (1|site) + (1|sp__@site), data = dat, family = "binomial",
#' cov_ranef = list(sp = phy), REML = TRUE, verbose = TRUE, s2.init = .1)
#'
#' ## Model structure from Rafferty & Ives (2013) (Eq. 3)
#' # dat = data set
#' # phyPol = phylogeny for pollinators (pol)
#' # phyPlt = phylogeny for plants (plt)
#'
#' z <- pglmm(freq ~ pol * X + (1|pol__) + (1|plt__) + (1|pol__@plt) +
#' (1|pol@plt__) + (1|pol__@plt__),
#' data = dat, family = "binomial",
#' cov_ranef = list(pol = phyPol, plt = phyPlt),
#' REML = TRUE, verbose = TRUE, s2.init = .1)
#' }
#'
#' #####################################################
#' ### Detailed analysis showing covariance matrices ###
#' #####################################################
#'
#' # This is the example from section 4.3 in Ives, A. R. (2018) Mixed
#' # and phylogenetic models: a conceptual introduction to correlated data.
#'
#' library(ape)
#' library(mvtnorm)
#'
#' # Investigating covariance matrices for different types of model structure
#' nspp <- 6
#' nsite <- 4
#'
#' # Simulate a phylogeny that has a lot of phylogenetic signal (power = 1.3)
#' phy <- compute.brlen(rtree(n = nspp), method = "Grafen", power = 1.3)
#'
#' # Simulate species means
#' sd.sp <- 1
#' mean.sp <- rTraitCont(phy, model = "BM", sigma=sd.sp^2)
#'
#' # Replicate values of mean.sp over sites
#' Y.sp <- rep(mean.sp, times=nsite)
#'
#' # Simulate site means
#' sd.site <- 1
#' mean.site <- rnorm(nsite, sd=sd.site)
#'
#' # Replicate values of mean.site over sp
#' Y.site <- rep(mean.site, each=nspp)
#'
#' # Compute a covariance matrix for phylogenetic attraction
#' sd.attract <- 1
#' Vphy <- vcv(phy)
#'
#' # Standardize the phylogenetic covariance matrix to have determinant = 1.
#' # (For an explanation of this standardization, see subsection 4.3.1 in Ives (2018))
#' Vphy <- Vphy/(det(Vphy)^(1/nspp))
#'
#' # Construct the overall covariance matrix for phylogenetic attraction.
#' # (For an explanation of Kronecker products, see subsection 4.3.1 in the book)
#' V <- kronecker(diag(nrow = nsite, ncol = nsite), Vphy)
#' Y.attract <- array(t(rmvnorm(n = 1, sigma = sd.attract^2*V)))
#'
#' # Simulate residual errors
#' sd.e <- 1
#' Y.e <- rnorm(nspp*nsite, sd = sd.e)
#'
#' # Construct the dataset
#' d <- data.frame(sp = rep(phy$tip.label, times = nsite),
#' site = rep(1:nsite, each = nspp))
#'
#' # Simulate abundance data
#' d$Y <- Y.sp + Y.site + Y.attract + Y.e
#'
#' # Analyze the model
#' pglmm(Y ~ 1 + (1|sp__) + (1|site) + (1|sp__@site), data = d, cov_ranef = list(sp = phy))
#'
#' # Display random effects: the function `pglmm_plot_ranef()` does what
#' # the name implies. You can set `show.image = TRUE` and `show.sim.image = TRUE`
#' # to see the matrices and simulations.
#' re <- pglmm_plot_ranef(Y ~ 1 + (1|sp__) + (1|site) + (1|sp__@site), data = d,
#' cov_ranef = list(sp = phy), show.image = FALSE,
#' show.sim.image = FALSE)
#'
#' #################################################
#' ### Example of a bipartite phylogenetic model ###
#' #################################################
#'
#' # Investigating covariance matrices for different types of model structure
#' nspp <- 20
#' nsite <- 15
#'
#' # Simulate a phylogeny that has a lot of phylogenetic signal (power = 1.3)
#' phy.sp <- compute.brlen(rtree(n = nspp), method = "Grafen", power = 1.3)
#' phy.site <- compute.brlen(rtree(n = nsite), method = "Grafen", power = 1.3)
#'
#' # Simulate species means
#' mean.sp <- rTraitCont(phy.sp, model = "BM", sigma = 1)
#'
#' # Replicate values of mean.sp over sites
#' Y.sp <- rep(mean.sp, times = nsite)
#'
#' # Simulate site means
#' mean.site <- rTraitCont(phy.site, model = "BM", sigma = 1)
#'
#' # Replicate values of mean.site over sp
#' Y.site <- rep(mean.site, each = nspp)
#'
#' # Generate covariance matrix for phylogenetic attraction among species
#' sd.sp.attract <- 1
#' Vphy.sp <- vcv(phy.sp)
#' Vphy.sp <- Vphy.sp/(det(Vphy.sp)^(1/nspp))
#' V.sp <- kronecker(diag(nrow = nsite, ncol = nsite), Vphy.sp)
#' Y.sp.attract <- array(t(rmvnorm(n = 1, sigma = sd.sp.attract^2*V.sp)))
#' # Generate covariance matrix for phylogenetic attraction among sites
#' sd.site.attract <- 1
#' Vphy.site <- vcv(phy.site)
#' Vphy.site <- Vphy.site/(det(Vphy.site)^(1/nsite))
#' V.site <- kronecker(Vphy.site, diag(nrow = nspp, ncol = nspp))
#' Y.site.attract <- array(t(rmvnorm(n = 1, sigma = sd.site.attract^2*V.site)))
#'
#' # Generate covariance matrix for phylogenetic attraction of species:site interaction
#' sd.sp.site.attract <- 1
#' V.sp.site <- kronecker(Vphy.site, Vphy.sp)
#' Y.sp.site.attract <- array(t(rmvnorm(n = 1, sigma = sd.sp.site.attract^2*V.sp.site)))
#'
#' # Simulate residual error
#' sd.e <- 0.5
#' Y.e <- rnorm(nspp*nsite, sd = sd.e)
#'
#' # Construct the dataset
#' d <- data.frame(sp = rep(phy.sp$tip.label, times = nsite),
#' site = rep(phy.site$tip.label, each = nspp))
#'
#' # Simulate abundance data
#' d$Y <- Y.sp + Y.site + Y.sp.attract + Y.site.attract + Y.sp.site.attract + Y.e
#'
#' # Plot random effects covariance matrices and then add phylogenies
#' # Note that, if show.image and show.sim are not specified, pglmm_plot_ranef() shows
#' # the covariance matrices if nspp * nsite < 200 and shows simulations
#' # if nspp * nsite > 100
#' re <- pglmm_plot_ranef(Y ~ 1 + (1|sp__) + (1|site__) + (1|sp__@site) +
#' (1|sp@site__) + (1|sp__@site__),
#' data=d, cov_ranef = list(sp = phy.sp, site = phy.site))
#'
#' # This flips the phylogeny to match to covariance matrices
#' rot.phy.site <- phy.site
#' for(i in (nsite+1):(nsite+Nnode(phy.site)))
#' rot.phy.site <- rotate(rot.phy.site, node = i)
#'
#' plot(phy.sp, main = "Species", direction = "upward")
#' plot(rot.phy.site, main = "Site")
#'
#' # Analyze the simulated data and compute a P-value for the (1|sp__@site__)
#' # random effect using a LRT. It is often better to fit the reduced model before
#' # the full model, because it s numerically easier to fit the reduced model,
#' # and then the parameter estimates from the reduced model can be given to the
#' # full model. In this case, I have used the estimates of the random effects
#' # from the reduce model, mod.r$ss, as the initial estimates for the same
#' # parameters in the full model in the statement s2.init=c(mod.r$ss, 0.01)^2.
#' # The final 0.01 is for the last random effect in the full model, (1|sp__@site__).
#' # Note also that the output of the random effects from communityPGLMM(), mod.r$ss,
#' # are the standard deviations, so they have to be squared for use as initial
#' # values of variances in mod.f.
#'
#' mod.r <- pglmm(Y ~ 1 + (1|sp__) + (1|site__) + (1|sp__@site) + (1|sp@site__),
#' data = d, cov_ranef = list(sp = phy.sp, site = phy.site))
#' mod.f <- pglmm(Y ~ 1 + (1|sp__) + (1|site__) + (1|sp__@site) + (1|sp@site__) +
#' (1|sp__@site__), data = d,
#' cov_ranef = list(sp = phy.sp, site = phy.site),
#' s2.init = c(mod.r$ss, 0.01)^2)
#' mod.f
#' pvalue <- pchisq(2*(mod.f$logLik - mod.r$logLik), df = 1, lower.tail = FALSE)
#' pvalue
#' }
pglmm <- function(formula, data = NULL, family = "gaussian", cov_ranef = NULL,
random.effects = NULL, REML = TRUE,
optimizer = c("nelder-mead-nlopt", "bobyqa", "Nelder-Mead", "subplex"),
repulsion = FALSE, add.obs.re = TRUE, verbose = FALSE,
cpp = TRUE, bayes = FALSE,
s2.init = NULL, B.init = NULL, reltol = 10^-6,
maxit = 500, tol.pql = 10^-6, maxit.pql = 200,
marginal.summ = "mean", calc.DIC = TRUE, calc.WAIC = TRUE, prior = "inla.default",
prior_alpha = 0.1, prior_mu = 1, ML.init = FALSE,
tree = NULL, tree_site = NULL, sp = NULL, site = NULL, bayes_options = NULL,
bayes_nested_matrix_as_list = FALSE
) {
optimizer = match.arg(optimizer)
if ((family %nin% c("gaussian", "binomial", "poisson")) & (bayes == FALSE)){
stop("\nSorry, but only binomial, poisson and gaussian options are available for
pglmm at this time")
}
if(bayes) {
if (!requireNamespace("INLA", quietly = TRUE)) {
stop("To run pglmm with bayes = TRUE, you need to install the packages 'INLA'. \
Please run in your R terminal:\
install.packages('INLA', repos='https://inla.r-inla-download.org/R/stable')")
}
if ((family %nin% c("gaussian", "binomial", "poisson", "zeroinflated.binomial", "zeroinflated.poisson"))){
stop("\nSorry, but only binomial (binary or success/failure), poisson (count), and gaussian options
are available for Bayesian pglmm at this time")
}
}
if(family %in% c("binomial", "poisson") & !is.null(tree)){
if(("phylo" %in% class(tree)) & !ape::is.ultrametric(tree)){
warning("The tree is not ultrametric, which will likely give misleading results for family=binomial and poisson models.")
}
}
data = as.data.frame(data) # in case of tibbles
fm_original = formula
prep_re = if(is.null(random.effects)) TRUE else FALSE
if(prep_re) {
# to make old code work ...
if(is.null(cov_ranef) & any(grepl("__", all.vars(formula)))){
if(!is.null(tree) | !is.null(tree_site))
warning("arguments `tree` and `tree_site` are deprecated; please use `cov_ranef` instead.",
call. = FALSE)
if(!is.null(tree) & is.null(tree_site)) cov_ranef = list(sp = tree) # column name must be sp
if(is.null(tree) & !is.null(tree_site)) cov_ranef = list(site = tree_site) # column name must be site
if(!is.null(tree) & !is.null(tree_site)) cov_ranef = list(sp = tree, site = tree_site)
}
dat_prepared = prep_dat_pglmm(formula, data, cov_ranef, repulsion, prep_re, family, add.obs.re, bayes, bayes_nested_matrix_as_list)
formula = dat_prepared$formula
random.effects = dat_prepared$random.effects
cov_ranef_updated = dat_prepared$cov_ranef_updated
} else {
formula = lme4::nobars(formula)
for(i in 1:length(random.effects)){
if(length(random.effects[[i]]) >= 3){
if(inherits(random.effects[[i]][[3]], c("matrix", "Matrix")) &
!is.null(rownames(random.effects[[i]][[3]]))){
if(!all(rownames(random.effects[[i]][[3]]) == colnames(random.effects[[i]][[3]])))
stop("the row and column names of cov matrix in random.effects[", i, "] not in the same order")
if(!all(rownames(random.effects[[i]][[3]]) == levels(random.effects[[i]][[2]]))){
warning("the row/column names of cov matrix in random.effects[", i,
"] not in the same order of its grouping variable, reordering now",
immediate. = TRUE, call. = FALSE)
if(length(setdiff(levels(random.effects[[i]][[2]]), rownames(random.effects[[i]][[3]]))))
stop("some levels in the grouping variable of random.effects[", i,
"] not in the cov matrix")
random.effects[[i]][[3]] =
random.effects[[i]][[3]][levels(random.effects[[i]][[2]]), levels(random.effects[[i]][[2]])]
}
}
}
}
}
# initial values for bayesian analysis: binomial and gaussian
if(bayes & ML.init & (family %in% c("binomial", "gaussian", "poisson"))) {
if (family == "gaussian") {
ML.init.z <- try(communityPGLMM.gaussian(formula = formula, data = data,
sp = sp, site = site,
random.effects = random.effects, REML = REML,
s2.init = s2.init, B.init = B.init,
reltol = reltol, maxit = maxit,
verbose = verbose, cpp = cpp, optimizer = optimizer))
if(!inherits(ML.init.z, "try-error")){
s2.init <- c(ML.init.z$s2r, ML.init.z$s2n, ML.init.z$s2resid)
B.init <- ML.init.z$B[ , 1, drop = TRUE]
} else {
warning("Initial model fitting with maximum likelihood approach failed.", immediate. = TRUE)
}
}
if (family %in% c("binomial", "poisson")) {# this may take too long if dataset is large...
if (is.null(s2.init)) s2.init <- 0.25
ML.init.z <- try(communityPGLMM.glmm(formula = formula, data = data,
sp = sp, site = site, family = family,
random.effects = random.effects, REML = REML,
s2.init = s2.init, B.init = B.init, reltol = reltol,
maxit = maxit, tol.pql = tol.pql, maxit.pql = maxit.pql,
verbose = verbose, cpp = cpp, optimizer = optimizer))
if(!inherits(ML.init.z, "try-error")){
s2.init <- c(ML.init.z$s2r, ML.init.z$s2n)
B.init <- ML.init.z$B[ , 1, drop = TRUE]
} else {
warning("Initial model fitting with maximum likelihood approach failed.", immediate. = TRUE)
}
}
}
if(bayes & ML.init & (family %nin% c("binomial", "gaussian", "poisson"))) {
warning('ML.init option is only available for binomial, poisson and gaussian families. You will have to
specify initial values manually if you think the default are problematic.')
}
if(bayes) {
z <- communityPGLMM.bayes(formula = formula, data = data, family = family,
sp = sp, site = site,
random.effects = random.effects,
s2.init = s2.init, B.init = B.init,
verbose = verbose,
marginal.summ = marginal.summ, calc.DIC = calc.DIC, calc.WAIC = calc.WAIC,
prior = prior,
prior_alpha = prior_alpha,
prior_mu = prior_mu,
bayes_options = bayes_options)
} else {# max likelihood
if (family == "gaussian") {
z <- communityPGLMM.gaussian(formula = formula, data = data,
sp = sp, site = site,
random.effects = random.effects, REML = REML,
s2.init = s2.init, B.init = B.init,
reltol = reltol, maxit = maxit,
verbose = verbose, cpp = cpp, optimizer = optimizer)
}
if (family %in% c("binomial", "poisson")) {
if (is.null(s2.init)) s2.init <- 0.25
z <- communityPGLMM.glmm(formula = formula, data = data, family = family,
sp = sp, site = site,
random.effects = random.effects, REML = REML,
s2.init = s2.init, B.init = B.init, reltol = reltol,
maxit = maxit, tol.pql = tol.pql, maxit.pql = maxit.pql,
verbose = verbose, cpp = cpp, optimizer = optimizer)
}
}
z$formula_original = fm_original
z$cov_ranef = if(is.null(cov_ranef)) NA else cov_ranef_updated
# # add names for ss
# if(!is.null(names(random.effects))){
# re.names = names(random.effects)[c(
# which(sapply(random.effects, length) %nin% c(1, 2, 4)), # non-nested terms
# which(sapply(random.effects, length) %in% c(1, 2, 4)) # nested terms
# )]
# if(family == "gaussian") re.names <- c(re.names, "residual")
# if((!bayes) & family != "gaussian") names(z$ss) = re.names
# }
return(z)
}
communityPGLMM.gaussian <- function(formula, data = list(), family = "gaussian",
sp = NULL, site = NULL, random.effects = list(),
REML = TRUE, s2.init = NULL, B.init = NULL,
reltol = 10^-8, maxit = 500, verbose = FALSE,
cpp = TRUE, optimizer = "bobyqa") {
dm = get_design_matrix(formula, data, random.effects, na.action = NULL)
X = dm$X; Y = dm$Y; St = dm$St; Zt = dm$Zt; nested = dm$nested
p <- ncol(X)
n <- nrow(X)
q <- length(random.effects)
# Compute initial estimates assuming no phylogeny if not provided
if (!is.null(B.init) & length(B.init) != p) {
warning("B.init not correct length, so computed B.init using glm()")
}
if ((is.null(B.init) | (!is.null(B.init) & length(B.init) != p)) & !is.null(s2.init)) {
B.init <- t(matrix(lm(formula = formula, data = data)$coefficients, ncol = p))
}
if (!is.null(B.init) & is.null(s2.init)) {
s2.init <- var(lm(formula = formula, data = data)$residuals)/q
}
if ((is.null(B.init) | (!is.null(B.init) & length(B.init) != p)) & is.null(s2.init)) {
B.init <- t(matrix(lm(formula = formula, data = data)$coefficients, ncol = p))
s2.init <- var(lm(formula = formula, data = data)$residuals)/q
}
B <- B.init
s <- as.vector(array(s2.init^0.5, dim = c(1, q)))
if(cpp){
if(is.null(St)) St = as(matrix(0, 0, 0), "dgTMatrix")
if(is.null(Zt)) Zt = as(matrix(0, 0, 0), "dgTMatrix")
out_res = pglmm_gaussian_internal_cpp(par = s, X, Y, Zt, St, nested, REML,
verbose, optimizer, maxit,
reltol, q, n, p, pi)
logLik = out_res$logLik
out = out_res$out
row.names(out$B) = colnames(X)
out$s2r = as.vector(out$s2r)
convcode = out_res$convcode
niter = out_res$niter[,1]
} else {# R version
if(optimizer == "Nelder-Mead"){
if (q > 1) {
opt <- optim(fn = pglmm_gaussian_LL_calc, par = s, X = X, Y = Y, Zt = Zt, St = St,
nested = nested, REML = REML, verbose = verbose,
method = "Nelder-Mead", control = list(maxit = maxit, reltol = reltol))
} else {
opt <- optim(fn = pglmm_gaussian_LL_calc, par = s, X = X, Y = Y, Zt = Zt, St = St,
nested = nested, REML = REML, verbose = verbose,
method = "L-BFGS-B", control = list(maxit = maxit))
}
} else {
# opts for nloptr
if (optimizer == "bobyqa") nlopt_algor = "NLOPT_LN_BOBYQA"
if (optimizer == "nelder-mead-nlopt") nlopt_algor = "NLOPT_LN_NELDERMEAD"
if (optimizer == "subplex") nlopt_algor = "NLOPT_LN_SBPLX"
opts <- list("algorithm" = nlopt_algor, "ftol_rel" = reltol, "ftol_abs" = reltol,
"xtol_rel" = 0.0001, "maxeval" = maxit)
S0 <- nloptr::nloptr(x0 = s, eval_f = pglmm_gaussian_LL_calc, opts = opts,
X = X, Y = Y, Zt = Zt, St = St, nested = nested,
REML = REML, verbose = verbose, optim_ll = TRUE)
opt = list(par = S0$solution, value = S0$objective, counts = S0$iterations,
convergence = S0$status, message = S0$message)
}
convcode = opt$convergence
niter = opt$counts
par_opt <- abs(Re(opt$par))
LL <- opt$value
out = pglmm_gaussian_LL_calc(par_opt, X, Y, Zt, St, nested, REML, verbose, optim_ll = FALSE)
out$B.cov = as.matrix(out$B.cov)
if (REML == TRUE) {
logLik <- as.numeric(-0.5 * (n - p) * log(2 * pi) + 0.5 * determinant(t(X) %*% X)$modulus[1] - LL)
} else {
logLik <- as.numeric(-0.5 * n * log(2 * pi) - LL)
}
}
# add names to variance estimates to make sure they are in order
if(!is.null(names(random.effects))){
re.len <- sapply(random.effects, length)
if(length(out$sr)) names(out$sr) <- names(out$s2r) <- names(random.effects)[re.len == 3]
if(length(out$sn)) names(out$sn) <- names(out$s2n) <- names(random.effects)[re.len %in% c(1, 4)]
names(out$s2resid) <- "residual"
}
ss <- c(out$sr, out$sn, out$s2resid^0.5)
B.zscore <- out$B/out$B.se
B.pvalue <- 2 * pnorm(abs(B.zscore), lower.tail = FALSE)
k <- p + q + 1
AIC <- -2 * logLik + 2 * k
BIC <- -2 * logLik + k * (log(n) - log(pi))
results <- list(formula = formula, data = data, family = family, random.effects = random.effects,
B = out$B, B.se = out$B.se, B.cov = out$B.cov, B.zscore = B.zscore,
B.pvalue = B.pvalue, ss = ss, s2n = out$s2n, s2r = out$s2r,
s2resid = out$s2resid, logLik = logLik, AIC = AIC, BIC = BIC,
REML = REML, bayes = FALSE, s2.init = s2.init, B.init = B.init, Y = Y, X = X, H = out$H,
iV = as.matrix(out$iV), mu = NULL, nested = nested, Zt = Zt, St = St,
convcode = convcode, niter = niter)
class(results) <- c("communityPGLMM", "pglmm")
results
}
communityPGLMM.glmm <- function(formula, data = list(), family = "binomial",
sp = NULL, site = NULL, random.effects = list(),
REML = TRUE, s2.init = 0.05, B.init = NULL,
reltol = 10^-5, maxit = 40, tol.pql = 10^-6,
maxit.pql = 200, verbose = FALSE, cpp = TRUE,
optimizer = "bobyqa") {
dm = get_design_matrix(formula, data, random.effects, na.action = NULL)
X = dm$X; Y = dm$Y; size = dm$size; St = dm$St; Zt = dm$Zt; nested = dm$nested
p <- ncol(X)
n <- nrow(X)
q <- length(random.effects)
# Compute initial estimates assuming no phylogeny if not provided
if (!is.null(B.init) & length(B.init) != p) {
warning("B.init not correct length, so computed B.init using glm()")
}
if ((is.null(B.init) | (!is.null(B.init) & length(B.init) != p))) {
B.init <- t(matrix(glm(formula = formula, data = data, family = family, na.action = na.omit)$coefficients, ncol = p))
} else {
B.init <- matrix(B.init, ncol = 1)
}
ss <- as.vector(array(s2.init^0.5, dim = c(1, q)))
if(cpp){
if(is.null(St)) St = as(matrix(0, 0, 0), "dgTMatrix")
if(is.null(Zt)) Zt = as(matrix(0, 0, 0), "dgTMatrix")
internal_res = pglmm_internal_cpp(X = X, Y = Y, Zt = Zt, St = St,
nested = nested, REML = REML, verbose = verbose,
n = n, p = p, q = q, maxit = maxit,
reltol = reltol, tol_pql = tol.pql,
maxit_pql = maxit.pql, optimizer = optimizer,
B_init = B.init, ss = ss,
family = family, totalSize = size)
B = internal_res$B
row.names(B) = colnames(X)
ss = internal_res$ss[,1]
iV = as(internal_res$iV, "dgCMatrix")
mu = internal_res$mu
row.names(mu) = 1:nrow(mu)
H = internal_res$H
convcode = internal_res$convcode
niter = internal_res$niter[, 1]
LL = internal_res$LL
} else {
B <- B.init
b <- matrix(0, nrow = n)
beta <- rbind(B, b)
if(family == "binomial") mu <- exp(X %*% B)/(1 + exp(X %*% B))
if(family == "poisson") mu <- exp(X %*% B)
XX <- cbind(X, diag(1, nrow = n, ncol = n))
est.ss <- ss
est.B <- B
oldest.ss <- 10^6
oldest.B <- matrix(10^6, nrow = length(est.B))
iteration <- 0
# exitflag <- 0; rcondflag <- 0
while (((t(est.ss - oldest.ss) %*% (est.ss - oldest.ss) > tol.pql^2) |
(t(est.B - oldest.B) %*% (est.B - oldest.B) > tol.pql^2)) &
(iteration <= maxit.pql)) {
iteration <- iteration + 1
oldest.ss <- est.ss
oldest.B <- est.B
est.B.m <- B
oldest.B.m <- matrix(10^6, nrow = length(est.B))
iteration.m <- 0
# mean component
while ((t(est.B.m - oldest.B.m) %*% (est.B.m - oldest.B.m) > tol.pql^2) &
(iteration.m <= maxit.pql)) {
iteration.m <- iteration.m + 1
oldest.B.m <- est.B.m
iV <- pglmm.iV.logdetV(par = ss, Zt = Zt, St = St, mu = mu, nested = nested, logdet = FALSE, family = family, size = size)$iV
if(family == "binomial") Z <- X %*% B + b + (Y/size - mu)/(mu * (1 - mu))
if(family == "poisson") Z <- X %*% B + b + (Y - mu)/mu
denom <- t(X) %*% iV %*% X
num <- t(X) %*% iV %*% Z
B <- solve(denom, matrix(num))
B <- as.matrix(B)
V = pglmm.V(par = ss, Zt = Zt, St = St, mu = mu, nested = nested, family = family, size = size)
if(family == "binomial") iW <- diag(as.vector(1/(size * mu * (1 - mu))))
if(family == "poisson") iW <- diag(as.vector(1/mu))
C <- V - iW
b <- C %*% iV %*% (Z - X %*% B)
beta <- rbind(B, matrix(b))
if(family == "binomial") mu <- exp(XX %*% beta)/(1 + exp(XX %*% beta))
if(family == "poisson") mu <- exp(XX %*% beta)
est.B.m <- B
if (verbose == TRUE) show(c(iteration, B))
# cat("mean part:", iteration.m, t(B), "\n")
# cat(" denom", as.matrix(denom), "\n")
# cat(" num", as.matrix(num), "\n")
if (any(is.nan(B))) {
stop("Estimation of B failed. Check for lack of variation in Y. You could try with a smaller s2.init, but this might not help.")
}
}
# variance component
if(family == "binomial") Z <- X %*% B + b + (Y/size - mu)/(mu * (1 - mu))
if(family == "poisson") Z <- X %*% B + b + (Y - mu)/mu
H <- matrix(Z - X %*% B)
if(optimizer == "Nelder-Mead"){
if (q > 1) {
opt <- optim(fn = pglmm.LL, par = ss, H = H, X = X, Zt = Zt, St = St,
mu = mu, nested = nested, family = family, size = size, REML = REML, verbose = verbose,
method = "Nelder-Mead", control = list(maxit = maxit, reltol = reltol))
} else {
opt <- optim(fn = pglmm.LL, par = ss, H = H, X = X, Zt = Zt, St = St,
mu = mu, nested = nested, family = family, size = size, REML = REML, verbose = verbose,
method = "L-BFGS-B", control = list(maxit = maxit))
}
} else {
if (optimizer == "bobyqa") nlopt_algor = "NLOPT_LN_BOBYQA"
if (optimizer == "nelder-mead-nlopt") nlopt_algor = "NLOPT_LN_NELDERMEAD"
if (optimizer == "subplex") nlopt_algor = "NLOPT_LN_SBPLX"
opts <- list("algorithm" = nlopt_algor, "ftol_rel" = reltol, "ftol_abs" = reltol,
"xtol_rel" = 0.0001, "maxeval" = maxit)
S0 <- nloptr::nloptr(x0 = ss, eval_f = pglmm.LL, opts = opts,
H = H, X = X, Zt = Zt, St = St, mu = mu,
nested = nested, family = family, size = size, REML = REML, verbose = verbose)
opt = list(par = S0$solution, value = S0$objective, counts = S0$iterations,
convergence = S0$status, message = S0$message)
}
ss <- abs(opt$par)
LL <- opt$value
convcode = opt$convergence
niter = opt$counts
if(verbose) cat("var part:", iteration, LL, ss, "\n")
est.ss <- ss
est.B <- B
}
row.names(B) = colnames(X)
}
# Extract parameters
q.nonNested = dm$q.nonNested
if (q.nonNested > 0) {
sr <- ss[1:q.nonNested]
} else {
sr <- NULL
}
q.Nested = dm$q.Nested
if (q.Nested > 0) {
sn <- ss[(q.nonNested + 1):(q.nonNested + q.Nested)]
} else {
sn <- NULL
}
# add names to variance estimates to make sure they are in order
if(!is.null(names(random.effects))){
re.len <- sapply(random.effects, length)
if(length(sr)) names(sr) <- names(random.effects)[re.len == 3]
if(length(sn)) names(sn) <- names(random.effects)[re.len %in% c(1, 4)]
names(ss) <- c(names(sr), names(sn))
}
s2r <- sr^2
s2n <- sn^2
if (family == 'binomial') {
mu_hat <- exp(X %*% B) / (1 + exp(X %*% B))
if (any(size > 1)) {
logLik.glm <-
sum(Y * log(mu_hat) + (size - Y) * log(1 - mu_hat) +
log(factorial(size) / (factorial(Y) * factorial(size - Y))))
} else{
logLik.glm <- sum(Y * log(mu_hat) + (1 - Y) * log(1 - mu_hat))
}
} else{
mu_hat <- exp(X %*% B)
logLik.glm <- sum(-mu_hat + Y * log(mu_hat) - log(factorial(Y)))
}
if(!is.null(St) && all(dim(St) == 0)) St <- NULL
logLik <- logLik.glm +
as.numeric(-LL + pglmm.LL(0 * ss, H = H, X = X, Zt = Zt, St = St, mu = mu,
nested = nested, REML = REML, family = family,
size = size, verbose = verbose))
k <- p + q + 1
AIC <- -2 * logLik + 2 * k
BIC <- -2 * logLik + k * (log(n) - log(pi))
B.cov <- solve(t(X) %*% iV %*% X)
B.se <- as.matrix(diag(B.cov))^0.5
B.zscore <- B/B.se
B.pvalue <- 2 * pnorm(abs(B/B.se), lower.tail = FALSE)
results <- list(formula = formula, data = data, family = family, random.effects = random.effects,
B = B, B.se = B.se, B.cov = B.cov, B.zscore = B.zscore, B.pvalue = B.pvalue,
ss = ss, s2n = s2n, s2r = s2r, s2resid = NULL, logLik = logLik, AIC = AIC,
BIC = BIC, REML = REML, bayes = FALSE, s2.init = s2.init, B.init = B.init, Y = Y, size = size, X = X,
H = as.matrix(H), iV = iV, mu = mu, nested = nested, Zt = Zt, St = St,
convcode = convcode, niter = niter)
class(results) <- c("communityPGLMM", "pglmm")
return(results)
}
communityPGLMM.bayes <- function(formula, data = list(), family = "gaussian",
sp = NULL, site = NULL, random.effects = list(),
s2.init = NULL, B.init = NULL,
verbose = FALSE,
marginal.summ = "mean", calc.DIC = FALSE, calc.WAIC = FALSE,
prior = "inla.default",
prior_alpha = 0.1, prior_mu = 1, bayes_options = NULL) {
mf <- model.frame(formula = formula, data = data, na.action = NULL)
X <- model.matrix(attr(mf, "terms"), data = mf)
Y <- model.response(mf)
p <- ncol(X)
n <- nrow(X)
q <- length(random.effects)
if(is.matrix(Y) && ncol(Y) == 2){ # success, fails for binomial data
Ntrials <- rowSums(Y) # total trials
Y <- Y[, 1] # success
# update formula
left_side = all.vars(update(formula, .~0))[1]
formula_bayes <- formula
formula_bayes[[2]] <- as.name(left_side)
# formula_bayes = as.formula(gsub(pattern = "^(cbind[(].*[)])",
# replacement = left_side, x = deparse(formula)))
} else {
formula_bayes = formula
Ntrials = NULL
}
base_family <- gsub("zeroinflated.", "", family, fixed = TRUE)
if(family == "zeroinflated.binomial") {
family <- "zeroinflatedbinomial1"
}
if(family == "zeroinflated.poisson") {
family <- "zeroinflatedpoisson1"
}
if(family == "gaussian") q <- q + 1
# Compute initial estimates assuming no phylogeny if not provided
if(family == "gaussian") {
if (!is.null(B.init) & length(B.init) != p) {
warning("B.init not correct length, so computed B.init using glm()")
}
if ((is.null(B.init) | (!is.null(B.init) & length(B.init) != p)) & !is.null(s2.init)) {
B.init <- lm(formula = formula, data = data)$coefficients
}
if (!is.null(B.init) & is.null(s2.init)) {
s2.init <- rep(var(lm(formula = formula, data = data)$residuals)/q, q)
}
if ((is.null(B.init) | (!is.null(B.init) & length(B.init) != p)) & is.null(s2.init)) {
B.init <- lm(formula = formula, data = data)$coefficients
s2.init <- rep(var(lm(formula = formula, data = data)$residuals)/q, q)
}
} else {
if (!is.null(B.init) & length(B.init) != p) {
warning("B.init not correct length, so computed B.init using glm()")
}
glm_bayes = glm(formula = formula, data = data, family = base_family, na.action = na.omit)
if ((is.null(B.init) | (!is.null(B.init) & length(B.init) != p))) {
B.init <- t(matrix(glm_bayes$coefficients, ncol = p))
} else {
B.init <- matrix(B.init, ncol = 1)
}
if (is.null(s2.init)) {
s2.init <- rep(var(glm_bayes$residuals)/q, q)
}
}
#B <- B.init
#s <- as.vector(array(s2.init^0.5, dim = c(1, q)))
s2.init <- log(1/s2.init)
if(family == "gaussian") {
resid.init <- s2.init[q]
s2.init <- s2.init[-q]
}
if(prior == "pc.prior.auto") {
if(family == "gaussian") {
lmod <- lm(formula, data)
sdres <- sd(residuals(lmod))
pcprior <- list(prec = list(prior = "pc.prec", param = c(3 * sdres, 0.01)))
} else {
if(family == "binomial") {
# lmod <- glm(formula, data = data, family = family)
# sdres <- sd(lmod$y - lmod$fitted.values)
pcprior <- list(prec = list(prior = "pc.prec", param = c(1, 0.1)))
} else {
warning("pc.prior.auto not yet implemented for this family. switching to default INLA prior...")
prior <- "inla.default"
}
}
}
if(prior == "pc.prior") {
pcprior <- list(prec = list(prior = "pc.prec", param = c(prior_mu, prior_alpha)))
}
if(prior == "uninformative") {
pcprior <- list(prec = list(prior = "pc.prec", param = c(100, 0.99)))
## very flat prior, generally not recommended!
}
# contruct INLA formula
inla_formula <- Reduce(paste, deparse(formula_bayes))
inla_effects <- vector("list", length = length(random.effects))
if(is.null(names(random.effects))){
names(inla_effects) <- as.character(1:length(inla_effects))
} else { # assign names so that the inla.model has names in random terms
names(inla_effects) <- names(random.effects)
}
inla_Cmat <- inla_weights <- inla_reps <- inla_effects
for(i in seq_along(random.effects)) {
if(length(random.effects[[i]]) == 1) {
# nested term: 1|sp@site, 1|sp__@site, 1|sp@site__, 1|sp__@site__
inla_effects[[i]] <- 1:nrow(data)
inla_Cmat[[i]] <- solve(random.effects[[i]][[1]])
} else if(length(random.effects[[i]]) == 2) {
# nested term: x|sp@site, x|sp__@site, x|sp@site__, x|sp__@site__
inla_effects[[i]] <- 1:nrow(data)
inla_weights[[i]] <- random.effects[[i]][[1]]
inla_Cmat[[i]] <- solve(random.effects[[i]][[2]])
} else if(length(random.effects[[i]]) == 3) {
# non-nested term: e.g. 1|sp__, x|sp__
inla_effects[[i]] <- as.numeric(random.effects[[i]][[2]])
inla_Cmat[[i]] <- solve(random.effects[[i]][[3]])
inla_weights[[i]] <- random.effects[[i]][[1]]
} else { # nested term: 1|sp@site, 1|sp__@site, 1|sp@site__, 1|sp__@site__
inla_effects[[i]] <- as.numeric(random.effects[[i]][[2]])
inla_Cmat[[i]] <- solve(random.effects[[i]][[3]])
inla_weights[[i]] <- random.effects[[i]][[1]]
inla_reps[[i]] <- as.numeric(random.effects[[i]][[4]])
}
}
if(!is.null(bayes_options$diagonal)) {
diagonal <- bayes_options$diagonal
bayes_options <- bayes_options[-which(names(bayes_options) == "diagonal")]
if(length(bayes_options) == 0) {
bayes_options <- NULL
}
} else {
diagonal <- NULL
}
f_form = vector(mode = "character", length = length(random.effects))
for(i in seq_along(random.effects)) {
if(length(random.effects[[i]]) == 3) { # non-nested term
if(length(random.effects[[i]][[1]]) == 1) { # 1 | sp, 1 | sp__, 1 | site, 1 | site__
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
} else { # x | sp, x | sp__, x | site, x | site__
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], inla_weights[[", i, "]], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
}
} else { # nested term 1 | sp__@site, etc.
if(length(random.effects[[i]]) == 4) {
if(length(random.effects[[i]][[1]]) == 1) {
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], replicate = inla_reps[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
} else {
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], replicate = inla_reps[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
}
} else { # length of 1 or 2: specified as a matrix (1|sp__@site) or list of 2 (x|sp__@site)
if(length(random.effects[[i]]) == 1) { # (1|sp__@site) etc.
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
} else {
if(length(random.effects[[i]]) == 2) { # (x|sp__@site) etc.
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], inla_weights[[", i, "]], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
} else { # other lengths? just in case ...
f_form[i] <- paste0("f(inla_effects[['", names(inla_effects)[i], "']], model = 'generic0', constr = FALSE, Cmatrix = inla_Cmat[[", i, "]], initial = s2.init[", i, "], diagonal = diagonal)")
}
}
}
}
}
if(prior != "inla.default"){
f_form = unname(sapply(f_form, function(x){
gsub(pattern = "[)]$", replacement = ", hyper = pcprior)", x)
}))
}
f_form = paste(f_form, collapse = " + ")
inla_formula <- paste(inla_formula, f_form, sep = " + ")
if(calc.DIC) {
if(calc.WAIC) {
control.compute <- list(dic = TRUE, waic = TRUE)
} else {
control.compute <- list(dic = TRUE)
}
} else {
if(calc.WAIC) {
control.compute <- list(waic = TRUE)
} else {
control.compute <- list()
}
}
argus <- c(list(formula = as.formula(inla_formula),
data = data,
verbose = verbose,
family = family),
bayes_options)
if(is.null(argus$control.fixed)) {
argus$control.fixed = list(prec.intercept = 0.0001, correlation.matrix = TRUE)
} else {
if(is.null(argus$control.fixed$prec.intercept)) {
argus$control.fixed$prec.intercept <- 0.0001
}
if(is.null(argus$control.fixed$correlation.matrix)) {
argus$control.fixed$correlation.matrix <- TRUE
}
}
if(is.null(argus$control.compute$dic)) {
argus$control.compute$dic <- calc.DIC
}
if(is.null(argus$control.compute$waic)) {
argus$control.compute$waic <- calc.WAIC
}
if(is.null(argus$control.compute$config)) {
argus$control.compute$config <- TRUE
}
if(compareVersion(as.character(packageVersion("INLA")), "21.07.10-1") == 1) {
if(is.null(argus$control.compute$return.marginals.predictor)) {
argus$control.compute$return.marginals.predictor <- TRUE
}
}
if(is.null(argus$control.predictor)) {
argus$control.predictor = list(compute = TRUE, link = 1)
} else {
if(is.null(argus$control.predictor$compute)) {
argus$control.predictor$compute <- TRUE
}
if(is.null(argus$control.predictor$link)) {
argus$control.predictor$link <- 1
}
}
if(family == "gaussian") {
if(is.null(argus$control.family)) {
argus$control.family = list(hyper = list(prec = list(initial = resid.init)))
} else {
if(is.null(argus$control.family$hyper)) {
argus$control.family$hyper <- list(prec = list(initial = resid.init))
}
}
# out <- INLA::inla(as.formula(inla_formula), data = data,
# verbose = verbose,
# control.family = list(hyper = list(prec = list(initial = resid.init))),
# control.fixed = list(prec.intercept = 0.0001, correlation.matrix = TRUE),
# control.compute = control.compute,
# control.predictor = list(compute = TRUE))
} else {
argus$Ntrials <- Ntrials
}
out <- do.call(INLA::inla, argus)
# out <- INLA::inla(as.formula(inla_formula), data = data,
# verbose = verbose,
# family = family,
# control.fixed = list(prec.intercept = 0.0001, correlation.matrix = TRUE),
# control.compute = control.compute,
# control.predictor=list(compute = TRUE),
# Ntrials = Ntrials)
#summary(out)
#print(out$summary.fitted.values)
if(calc.DIC) {
DIC <- out$dic$dic
} else {
DIC <- NULL
}
if(calc.WAIC){
WAIC <- out$waic$waic
} else {
WAIC <- NULL
}
if(marginal.summ == "median") marginal.summ <- "0.5quant"
nested <- sapply(random.effects, length) %in% c(1, 2, 4)
variances <- 1/out$summary.hyperpar[ , marginal.summ]
# names(variances) <- gsub("Precision for ", "", rownames(out$summary.hyperpar))
variances.ci <- 1/out$summary.hyperpar[ , c("0.975quant", "0.025quant")]
if(family == "gaussian") {
resid_var <- variances[1]
names(resid_var) <- "residual"
variances <- variances[-1]
resid_var.ci <- variances.ci[1, ]
variances.ci <- variances.ci[-1, ]
} else {
resid_var <- NULL
resid_var.ci <- NULL
}
if(grepl("zeroinflated", family)) {
zeroinlated_param <- 1/variances[1]
variances <- variances[-1]
zeroinflated_param.ci <- rev(1/variances.ci[1, ])
variances.ci <- variances.ci[-1, ]
} else {
zeroinlated_param <- NULL
zeroinflated_param.ci <- NULL
}
if(!is.null(names(random.effects))){
names(variances) <- names(random.effects)
row.names(variances.ci) <- paste("Precision for", names(random.effects))
}
std.vars <- variances^0.5
ss <- c(std.vars[!nested], std.vars[nested], resid_var)
if(marginal.summ == "median") marginal.summ <- "0.5quant"
B <- out$summary.fixed[ , marginal.summ]
H <- Y - out$summary.fitted.values[ , marginal.summ, drop = TRUE]
mu <- out$summary.fitted.values[ , marginal.summ, drop = FALSE]
#H <- NULL
results <- list(formula = formula, data = data, family = family, random.effects = random.effects,
B = B, B.se = NULL,
B.ci = out$summary.fixed[ , c("0.025quant", "0.975quant")],
B.cov = out$misc$lincomb.derived.correlation.matrix, B.zscore = NULL,
B.pvalue = NULL, ss = ss, s2n = variances[nested], s2r = variances[!nested],
s2resid = resid_var, zi = zeroinlated_param, s2n.ci = variances.ci[nested, ],
s2r.ci = variances.ci[!nested, ], s2resid.ci = resid_var.ci,
zi.ci = zeroinflated_param.ci,
logLik = out$mlik[1, 1], AIC = NULL, BIC = NULL, DIC = DIC, WAIC = WAIC,
REML = NULL, bayes = TRUE, marginal.summ = marginal.summ,
s2.init = s2.init, B.init = B.init, Y = Y, X = X, H = H,
iV = NULL, mu = mu, nested = nested, Zt = NULL, St = NULL,
convcode = NULL, niter = NULL, inla.model = out)
class(results) <- c("communityPGLMM", "pglmm")
results
}
#' @export
#' @rdname pglmm
communityPGLMM <- pglmm # to be compatible with old code
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.