qrrvglm.control: Control function for QRR-VGLMs (CQO)

Description

Algorithmic constants and parameters for a constrained quadratic ordination (CQO), by fitting a quadratic reduced-rank vector generalized linear model (QRR-VGLM), are set using this function. It is the control function for cqo.

Usage

qrrvglm.control(Rank = 1,
                Bestof = if(length(Cinit)) 1 else 10,
                checkwz=TRUE,
                Cinit = NULL,
                Crow1positive = TRUE,
                epsilon = 1.0e-06,
                EqualTolerances = ITolerances,
                Etamat.colmax = 10,
                FastAlgorithm = TRUE,
                GradientFunction=TRUE,
                Hstep = 0.001,
                isdlv = rep(c(2, 1, rep(0.5, len=Rank)), len=Rank),
                iKvector = 0.1,
                iShape = 0.1,
                ITolerances = TRUE,
                maxitl = 40,
                method.init = 1,
                Maxit.optim = 250,
                MUXfactor = rep(7, length=Rank),
                Norrr = ~ 1,
                optim.maxit = 20,
                Parscale = if(ITolerances) 0.001 else 1.0,
                SD.Cinit = 0.02,
                SmallNo = 5.0e-13, 
                trace = TRUE,
                Use.Init.Poisson.QO=TRUE, 
                wzepsilon = .Machine$double.eps^0.75, ...)

Arguments

Rank

The numerical rank $R$ of the model, i.e., the number of ordination axes. Must be an element from the set {1,2,...,min($M$,$p_2$)} where the vector of explanatory variables $x$ is partitioned into ($x_1$,$x_2$), which is of dimension $p_1

Bestof

Integer. The best of Bestof models fitted is returned. This argument helps guard against local solutions by (hopefully) finding the global solution from many fits. The argument has value 1 if an initial value for $C$ is inputted

checkwz

logical indicating whether the diagonal elements of the working weight matrices should be checked whether they are sufficiently positive, i.e., greater than wzepsilon. If not, any values less than wzepsilon are replac

Cinit

Optional initial $C$ matrix, which must be a $p_2$ by $R$ matrix. The default is to apply .Init.Poisson.QO() to obtain initial values.

Crow1positive

Logical vector of length Rank (recycled if necessary): are the elements of the first row of $C$ positive? For example, if Rank is 4, then specifying

Crow1positive=c(FALSE,
      TRUE)

will force $C[1,1]$

epsilon

Positive numeric. Used to test for convergence for GLMs fitted in FORTRAN. Larger values mean a loosening of the convergence criterion. If an error code of 3 is reported, try increasing this value.

EqualTolerances

Logical indicating whether each (quadratic) predictor will have equal tolerances. Setting EqualTolerances=TRUE can help avoid numerical problems, especially with binary data. Note that the estimated (common) tolerance matri

Etamat.colmax

Positive integer, no smaller than Rank. Controls the amount of memory used by .Init.Poisson.QO(). It is the maximum number of columns allowed for the pseudo-response and its weights. In general, the larger the value

FastAlgorithm

Logical. Whether a new fast algorithm is to be used. The fast algorithm results in a large speed increases compared to Yee (2004). Some details of the fast algorithm are found in Appendix A of Yee (2006). Setting FastAlgorithm=FALSE

GradientFunction

Logical. Whether optim's argument gr
   is used or not, i.e., to compute gradient values.  Used only if
   FastAlgorithm is TRUE.  The default value is usually
   f

Hstep

Positive value. Used as the step size in the finite difference
   approximation to the derivatives by optim.

isdlv

Initial standard deviations for the latent variables (site scores).
   Numeric, positive and of length $R$ (recycled if necessary).
   This argument is used only if ITolerances=TRUE.  Used by
   .Init.Poisson.QO() to obtain initi

iKvector, iShape

Numeric, recycled to length $S$ if necessary.
   Initial values used for estimating the positive $k$ and
   $\lambda$ parameters of the negative binomial and
   2-parameter gamma distributions respectively. For further information
   see

ITolerances

Logical. If TRUE then the (common) tolerance matrix is the
   $R$ by $R$ identity matrix by definition.  Note that having
   ITolerances=TRUE implies EqualTolerances=TRUE, but
   not vice versa.  Internally, the quad

maxitl

Maximum number of times the optimizer is called or restarted.
    Most users should ignore this argument.

method.init

Method of initialization. A positive integer 1 or 2 or 3 etc.
    depending on the VGAM family function.
    Currently it is used for negbinomial and

Maxit.optim

Positive integer.  Number of iterations given to the function
    optim at each of the optim.maxit
    iterations.

MUXfactor

Multiplication factor for detecting large offset values.  Numeric,
   positive and of length $R$ (recycled if necessary).  This argument
   is used only if ITolerances=TRUE. Offsets are $-0.5$
   multiplied by the sum of the squares of all $R

optim.maxit

Positive integer.  Number of times optim
    is invoked.  At iteration i, the ith value of
    Maxit.optim is fed into opti

Norrr

Formula giving terms that are not to be included in the
    reduced-rank regression (or formation of the latent variables),
    i.e., those belong to $x_1$.
    Those variables which do not make up the latent variable (reduced-rank
    regression

Parscale

Numerical and positive-valued vector of length $C$
   (recycled if necessary).
   Passed into optim(..., control=list(parscale=Parscale));
   the elements of $C$ become $C$ / Parscale.
   Setting ITolerances=TRUE res

SD.Cinit

Standard deviation of the initial values for the elements
      of $C$.
      These are normally distributed with mean zero.
      This argument is used only if Use.Init.Poisson.QO = FALSE
      and $C$ is not inputted using Cinit

trace

Logical indicating if output should be produced for
    each iteration. The default is TRUE because the
    calculations are numerically intensive, meaning it may take
    a long time, so that the user might think the computer has
    locked

SmallNo

Positive numeric between .Machine$double.eps and 0.0001. 
  Used to avoid under- or over-flow in the IRLS algorithm. 
  Used only if FastAlgorithm is TRUE.

Use.Init.Poisson.QO

Logical. If TRUE then the function .Init.Poisson.QO() is
    used to obtain initial values for the canonical coefficients $C$.
    If FALSE then random numbers are used instead.

wzepsilon

Small positive number used to test whether the diagonals of the working
  weight matrices are sufficiently positive.

...

Ignored at present.

`Value`

A list with components matching the input names.

`Warning`

The default value of Bestof is a bare minimum for many datasets,
therefore it will be necessary to increase its value to increase the
chances of obtaining the global solution.

`Details`

Recall that the central formula for CQO is
  $$\eta = B_1^T x_1 + A \nu +
               \sum_{m=1}^M (\nu^T D_m \nu) e_m$$
  where $x_1$ is a vector (usually just a 1 for an intercept),
  $x_2$ is a vector of environmental variables, $\nu=C^T
  x_2$ is a $R$-vector of latent variables, $e_m$ is
  a vector of 0s but with a 1 in the $m$th position.
    QRR-VGLMs are an extension of RR-VGLMs and allow for maximum
    likelihood solutions to constrained quadratic ordination (CQO) models.
   Having ITolerances=TRUE means all the tolerance matrices
   are the order-$R$ identity matrix, i.e., it forces
   bell-shaped curves/surfaces on all species. This results in a
   more difficult optimization problem (especially for 2-parameter
   models such as the negative binomial and gamma) because of overflow
   errors and it appears there are more local solutions.  To help avoid
   the overflow errors, scaling $C$ by the factor Parscale
   can help enormously.  Even better, scaling $C$ by specifying
   isdlv is more understandable to humans.  If failure to
   converge occurs, try adjusting Parscale, or better, setting
   EqualTolerances=TRUE (and hope that the estimated tolerance
   matrix is positive-definite).  To fit an equal-tolerances model, it
   is firstly best to try setting ITolerances=TRUE and varying
   isdlv and/or MUXfactor if it fails to converge.
   If it still fails to converge after many attempts, try setting
   EqualTolerances=TRUE, however this will usually be a lot slower
   because it requires a lot more memory.
   With a $R>1$ model, the latent variables are always uncorrelated,
   i.e., the variance-covariance matrix of the site scores is a diagonal
   matrix.
   If setting EqualTolerances=TRUE is used and the common
   estimated tolerance matrix is positive-definite then that model is
   effectively the same as the ITolerances=TRUE model (the two are
   transformations of each other). In general, ITolerances=TRUE
   is numerically more unstable and presents a more difficult problem
   to optimize; the arguments isdlv and/or MUXfactor often
   must be assigned some good value(s) (possibly found by trial and error)
   in order for convergence to occur.  Setting ITolerances=TRUE
   forces a bell-shaped curve or surface onto all the species data,
   therefore this option should be used with deliberation.  If unsuitable,
   the resulting fit may be very misleading.  Usually it is a good idea
   for the user to set EqualTolerances=FALSE to see which species
   appear to have a bell-shaped curve or surface. Improvements to the
   fit can often be achieved using transformations, e.g., nitrogen
   concentration to log nitrogen concentration.
   Fitting a CAO model (see cao) first is a good idea for
   pre-examining the data and checking whether it is appropriate to fit
   a CQO model.

`References`

Yee, T. W. (2004)
A new technique for maximum-likelihood
canonical Gaussian ordination.
Ecological Monographs,
74, 685--701.
Yee, T. W. (2006)
Constrained additive ordination.
Ecology,
87, 203--213.

`See Also`

cqo,
  rcqo,
  Coef.qrrvglm,
  Coef.qrrvglm-class,  optim,
  binomialff,
  poissonff,
  negbinomial,
  gamma2,
  gaussianff.

`Examples`

Run this code# Poisson CQO with equal tolerances
data(hspider)
set.seed(111)  # This leads to the global solution
hspider[,1:6]=scale(hspider[,1:6]) # Good idea when ITolerances = TRUE
p1 = cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi,
               Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~
         WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
         ITolerances = TRUE, 
         fam = quasipoissonff, data = hspider)
sort(p1@misc$deviance.Bestof) # A history of all the iterations

(isdlv = sd(lv(p1))) # should be approx isdlv
 
# Refit the model with better initial values
set.seed(111)  # This leads to the global solution
p1 = cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, 
               Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~
         WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
         ITolerances = TRUE, isdlv = isdlv,   # Note the use of isdlv here
         fam = quasipoissonff, data = hspider)
sort(p1@misc$deviance.Bestof) # A history of all the iterations

# Negative binomial CQO; smallest deviance is about 275.389
set.seed(111)  # This leads to the global solution
nb1 = cqo(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi, 
                Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~
          WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
          ITol = FALSE, EqualTol = TRUE, # A good idea for negbinomial
          fam = negbinomial, data = hspider)
sort(nb1@misc$deviance.Bestof) # A history of all the iterations
summary(nb1)
lvplot(nb1, lcol=1:12, y=TRUE, pcol=1:12)
Run the code above in your browser using DataLab