Package 'VGAM'

Description

VGAM provides functions for fitting vector generalized linear and additive models (VGLMs and VGAMs), and associated models (Reduced-rank VGLMs or RR-VGLMs, Doubly constrained RR-VGLMs (DRR-VGLMs), Quadratic RR-VGLMs, Reduced-rank VGAMs). This package fits many models and distributions by maximum likelihood estimation (MLE) or penalized MLE, under this statistical framework. Also fits constrained ordination models in ecology such as constrained quadratic ordination (CQO).

Details

This package centers on the iteratively reweighted least squares (IRLS) algorithm. Other key words include Fisher scoring, additive models, reduced-rank regression, penalized likelihood, and constrained ordination. The central modelling functions are vglm, vgam, rrvglm, rcim, cqo, cao. Function vglm operates very similarly to glm but is much more general, and many methods functions such as coef and predict are available. The package uses S4 (see methods-package).

Some notable companion packages: (1) VGAMdata mainly contains data sets useful for illustrating VGAM. Some of the big ones were initially from VGAM. Recently, some older VGAM family functions have been shifted into this package. (2) VGAMextra written by Victor Miranda has some additional VGAM family and link functions, with a bent towards time series models. (3) svyVGAM provides design-based inference, e.g., to survey sampling settings. This is because the weights argument of vglm can be assigned any positive values including survey weights.

Compared to other similar packages, such as gamlss and mgcv, VGAM has more models implemented (150+ of them) and they are not restricted to a location-scale-shape framework or (largely) the 1-parameter exponential family. The general statistical framework behind it all, once grasped, makes regression modelling unified. Some features of the package are: (i) many family functions handle multiple responses; (ii) reduced-rank regression is available by operating on latent variables (optimal linear combinations of the explanatory variables); (iii) basic automatic smoothing parameter selection is implemented for VGAMs (sm.os and sm.ps with a call to magic), although it has to be refined; (iv) smart prediction allows correct prediction of nested terms in the formula provided smart functions are used.

The GLM and GAM classes are special cases of VGLMs and VGAMs. The VGLM/VGAM framework is intended to be very general so that it encompasses as many distributions and models as possible. VGLMs are limited only by the assumption that the regression coefficients enter through a set of linear predictors. The VGLM class is very large and encompasses a wide range of multivariate response types and models, e.g., it includes univariate and multivariate distributions, categorical data analysis, extreme values, correlated binary data, quantile and expectile regression, time series problems. Potentially, it can handle generalized estimating equations, survival analysis, bioassay data and nonlinear least-squares problems.

Crudely, VGAMs are to VGLMs what GAMs are to GLMs. Two types of VGAMs are implemented: 1st-generation VGAMs with s use vector backfitting, while 2nd-generation VGAMs with sm.os and sm.ps use O-splines and P-splines so have a direct solution (hence avoids backfitting) and have automatic smoothing parameter selection. The former is older and is based on Yee and Wild (1996). The latter is more modern (Yee, Somchit and Wild, 2024) but it requires a reasonably large number of observations to work well because it is based on optimizing over a predictive criterion rather than using a Bayesian approach.

An important feature of the framework is that of constraint matrices. They apportion the regression coefficients according to each explanatory variable. For example, since each parameter has a link function applied to it to turn it into a linear or additive predictor, does a covariate have an equal effect on each parameter? Or no effect? Arguments such as zero, parallel and exchangeable, are merely easy ways to have them constructed internally. Users may input them explicitly using the constraint argument, and CM.symm0 etc. can make this easier.

Another important feature is implemented by xij. It allows different linear/additive predictors to have a different values of the same explanatory variable, e.g., multinomial for the conditional logit model and the like.

VGLMs with dimension reduction form the class of RR-VGLMs. This is achieved by reduced rank regression. Here, a subset of the constraint matrices are estimated rather than being known and prespecified. Optimal linear combinations of the explanatory variables are taken (creating latent variables) which are used for fitting a VGLM. Thus the regression can be thought of as being in two stages. The class of DRR-VGLMs provides further structure to RR-VGLMs by allowing constraint matrices to be specified for each column of A and row of C. Thus the reduced rank regression can be fitted with greater control.

This package is the first to check for the Hauck-Donner effect (HDE) in regression models; see hdeff. This is an aberration of the Wald statistics when the parameter estimates are too close to the boundary of the parameter space. When present the p-value of a regression coefficient is biased upwards so that a highly significant variable might be deemed nonsignificant. Thus the HDE can create havoc for variable selection!

Somewhat related to the previous paragraph, hypothesis testing using the likelihood ratio test, Rao's score test (Lagrange multiplier test) and (modified) Wald's test are all available; see summaryvglm. For all regression coefficients of a model, taken one at a time, all three methods require further IRLS iterations to obtain new values of the other regression coefficients after one of the coefficients has had its value set (usually to 0). Hence the computation load is overall significant.

For a complete list of this package, use library(help = "VGAM"). New VGAM family functions are continually being written and added to the package.

Warning

This package is undergoing continual development and improvement, therefore users should treat many things as subject to change. This includes the family function names, argument names, many of the internals, moving some functions to VGAMdata, the use of link functions, and slot names. For example, many link functions were renamed in 2019 so that they all end in "link", e.g., loglink() instead of loge(). Some future pain can be avoided by using good programming techniques, e.g., using extractor functions such as coef(), weights(), vcov(), predict(). Although changes are now less frequent, please expect changes in all aspects of the package. See the NEWS file for a list of changes from version to version.

Author(s)

Thomas W. Yee, [email protected], with contributions from Victor Miranda and several graduate students over the years, especially Xiangjie (Albert) Xue and Chanatda Somchit.

Maintainer: Thomas Yee [email protected].

References

Yee, T. W. (2015). Vector Generalized Linear and Additive Models: With an Implementation in R. New York, USA: Springer.

Yee, T. W. and Hastie, T. J. (2003). Reduced-rank vector generalized linear models. Statistical Modelling, 3, 15–41.

Yee, T. W. and Stephenson, A. G. (2007). Vector generalized linear and additive extreme value models. Extremes, 10, 1–19.

Yee, T. W. and Wild, C. J. (1996). Vector generalized additive models. Journal of the Royal Statistical Society, Series B, Methodological, 58, 481–493.

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.

Yee, T. W. (2008). The VGAM Package. R News, 8, 28–39.

Yee, T. W. (2010). The VGAM package for categorical data analysis. Journal of Statistical Software, 32, 1–34. doi:10.18637/jss.v032.i10.

Yee, T. W. (2014). Reduced-rank vector generalized linear models with two linear predictors. Computational Statistics and Data Analysis, 71, 889–902.

Yee, T. W. and Ma, C. (2024). Generally altered, inflated, truncated and deflated regression. Statistical Science, 39 (in press).

Yee, T. W. (2022). On the Hauck-Donner effect in Wald tests: Detection, tipping points and parameter space characterization, Journal of the American Statistical Association, 117, 1763–1774. doi:10.1080/01621459.2021.1886936.

Yee, T. W. and Somchit, C. and Wild, C. J. (2024). Penalized vector generalized additive models. Manuscript in preparation.

The website for the VGAM package and book is https://www.stat.auckland.ac.nz/~yee/. There are some resources there, especially as relating to my book and new features added to VGAM.

Some useful background reference for the package include:

Chambers, J. and Hastie, T. (1991). Statistical Models in S. Wadsworth & Brooks/Cole.

Green, P. J. and Silverman, B. W. (1994). Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach. Chapman and Hall.

Hastie, T. J. and Tibshirani, R. J. (1990). Generalized Additive Models. Chapman and Hall.

See Also

vglm, vgam, rrvglm, rcim, cqo, TypicalVGAMfamilyFunction, CommonVGAMffArguments, Links, hdeff, glm, lm, https://CRAN.R-project.org/package=VGAM.

Examples

# Example 1; proportional odds model
pneumo <- transform(pneumo, let = log(exposure.time))
(fit1 <- vglm(cbind(normal, mild, severe) ~ let, propodds, data = pneumo))
depvar(fit1)  # Better than using fit1@y; dependent variable (response)
weights(fit1, type = "prior")  # Number of observations
coef(fit1, matrix = TRUE)      # p.179, in McCullagh and Nelder (1989)
constraints(fit1)              # Constraint matrices
summary(fit1)  # HDE could affect these results
summary(fit1, lrt0 = TRUE, score0 = TRUE, wald0 = TRUE)  # No HDE
hdeff(fit1)  # Check for any Hauck-Donner effect

# Example 2; zero-inflated Poisson model
zdata <- data.frame(x2 = runif(nn <- 2000))
zdata <- transform(zdata, pstr0  = logitlink(-0.5 + 1*x2, inverse = TRUE),
                          lambda = loglink(  0.5 + 2*x2, inverse = TRUE))
zdata <- transform(zdata, y = rzipois(nn, lambda, pstr0 = pstr0))
with(zdata, table(y))
fit2 <- vglm(y ~ x2, zipoisson, data = zdata, trace = TRUE)
coef(fit2, matrix = TRUE)  # These should agree with the above values


# Example 3; fit a two species GAM simultaneously
fit3 <- vgam(cbind(agaaus, kniexc) ~ s(altitude, df = c(2, 3)),
             binomialff(multiple.responses = TRUE), data = hunua)
coef(fit3, matrix = TRUE)   # Not really interpretable
## Not run:  plot(fit3, se = TRUE, overlay = TRUE, lcol = 3:4, scol = 3:4)

ooo <- with(hunua, order(altitude))
with(hunua,  matplot(altitude[ooo], fitted(fit3)[ooo, ], type = "l",
     lwd = 2, col = 3:4,
     xlab = "Altitude (m)", ylab = "Probability of presence", las = 1,
     main = "Two plant species' response curves", ylim = c(0, 0.8)))
with(hunua, rug(altitude)) 
## End(Not run)


# Example 4; LMS quantile regression
fit4 <- vgam(BMI ~ s(age, df = c(4, 2)), lms.bcn(zero = 1),
             data = bmi.nz, trace = TRUE)
head(predict(fit4))
head(fitted(fit4))
head(bmi.nz)  # Person 1 is near the lower quartile among people his age
head(cdf(fit4))

## Not run:  par(mfrow = c(1,1), bty = "l", mar = c(5,4,4,3)+0.1, xpd=TRUE)
qtplot(fit4, percentiles = c(5,50,90,99), main = "Quantiles", las = 1,
       xlim = c(15, 90), ylab = "BMI", lwd=2, lcol=4)  # Quantile plot

ygrid <- seq(15, 43, len = 100)  # BMI ranges
par(mfrow = c(1, 1), lwd = 2)  # Density plot
aa <- deplot(fit4, x0 = 20, y = ygrid, xlab = "BMI", col = "black",
    main = "Density functions at Age=20 (black), 42 (red) and 55 (blue)")
aa
aa <- deplot(fit4, x0 = 42, y = ygrid, add = TRUE, llty = 2, col = "red")
aa <- deplot(fit4, x0 = 55, y = ygrid, add = TRUE, llty = 4, col = "blue",
            Attach = TRUE)
aa@post$deplot  # Contains density function values

## End(Not run)


# Example 5; GEV distribution for extremes
(fit5 <- vglm(maxtemp ~ 1, gevff, data = oxtemp, trace = TRUE))
head(fitted(fit5))
coef(fit5, matrix = TRUE)
Coef(fit5)
vcov(fit5)
vcov(fit5, untransform = TRUE)
sqrt(diag(vcov(fit5)))  # Approximate standard errors
## Not run:  rlplot(fit5)

Description

Estimates the three independent parameters of the the A1A2A3 blood group system.

Usage

A1A2A3(link = "logitlink", inbreeding = FALSE, ip1 = NULL, ip2 = NULL, iF = NULL)

Arguments

Details

The parameters p1 and p2 are probabilities, so that p3=1-p1-p2 is the third probability. The parameter f is the third independent parameter if inbreeding = TRUE. If inbreeding = FALSE then $f = 0$ and Hardy-Weinberg Equilibrium (HWE) is assumed.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Note

The input can be a 6-column matrix of counts, with columns corresponding to A1A1, A1A2, A1A3, A2A2, A2A3, A3A3 (in order). Alternatively, the input can be a 6-column matrix of proportions (so each row adds to 1) and the weights argument is used to specify the total number of counts for each row.

Author(s)

T. W. Yee

References

Lange, K. (2002). Mathematical and Statistical Methods for Genetic Analysis, 2nd ed. New York: Springer-Verlag.

See Also

AA.Aa.aa, AB.Ab.aB.ab, ABO, MNSs.

Examples

ymat <- cbind(108, 196, 429, 143, 513, 559)
fit <- vglm(ymat ~ 1, A1A2A3(link = probitlink), trace = TRUE, crit = "coef")
fit <- vglm(ymat ~ 1, A1A2A3(link = logitlink, ip1 = 0.3, ip2 = 0.3, iF = 0.02),
            trace = TRUE, crit = "coef")
Coef(fit)  # Estimated p1 and p2
rbind(ymat, sum(ymat) * fitted(fit))
sqrt(diag(vcov(fit)))

Description

Estimates the parameter of the AA-Aa-aa blood group system, with or without Hardy Weinberg equilibrium.

Usage

AA.Aa.aa(linkp = "logitlink", linkf = "logitlink", inbreeding = FALSE,
         ipA = NULL, ifp = NULL, zero = NULL)

Arguments

Details

This one or two parameter model involves a probability called pA. The probability of getting a count in the first column of the input (an AA) is pA*pA. When inbreeding = TRUE, an additional parameter f is used. If inbreeding = FALSE then $f = 0$ and Hardy-Weinberg Equilibrium (HWE) is assumed. The EIM is used if inbreeding = FALSE.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Warning

Setting inbreeding = FALSE makes estimation difficult with non-intercept-only models. Currently, this code seems to work with intercept-only models.

Note

The input can be a 3-column matrix of counts, where the columns are AA, Ab and aa (in order). Alternatively, the input can be a 3-column matrix of proportions (so each row adds to 1) and the weights argument is used to specify the total number of counts for each row.

Author(s)

T. W. Yee

References

Weir, B. S. (1996). Genetic Data Analysis II: Methods for Discrete Population Genetic Data, Sunderland, MA: Sinauer Associates, Inc.

See Also

AB.Ab.aB.ab, ABO, A1A2A3, MNSs.

Examples

y <- cbind(53, 95, 38)
fit1 <- vglm(y ~ 1, AA.Aa.aa, trace = TRUE)
fit2 <- vglm(y ~ 1, AA.Aa.aa(inbreeding = TRUE), trace = TRUE)
rbind(y, sum(y) * fitted(fit1))
Coef(fit1)  # Estimated pA
Coef(fit2)  # Estimated pA and f
summary(fit1)

Description

Estimates the parameter of the AB-Ab-aB-ab blood group system.

Usage

AB.Ab.aB.ab(link = "logitlink", init.p = NULL)

Arguments

Details

This one parameter model involves a probability called p.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Note

The input can be a 4-column matrix of counts, where the columns are AB, Ab, aB and ab (in order). Alternatively, the input can be a 4-column matrix of proportions (so each row adds to 1) and the weights argument is used to specify the total number of counts for each row.

Author(s)

T. W. Yee

References

Lange, K. (2002). Mathematical and Statistical Methods for Genetic Analysis, 2nd ed. New York: Springer-Verlag.

See Also

AA.Aa.aa, ABO, A1A2A3, MNSs.

Examples

ymat <- cbind(AB=1997, Ab=906, aB=904, ab=32)  # Data from Fisher (1925)
fit <- vglm(ymat ~ 1, AB.Ab.aB.ab(link = "identitylink"), trace = TRUE)
fit <- vglm(ymat ~ 1, AB.Ab.aB.ab, trace = TRUE)
rbind(ymat, sum(ymat)*fitted(fit))
Coef(fit)  # Estimated p
p <- sqrt(4*(fitted(fit)[, 4]))
p*p
summary(fit)

Description

Estimates the two independent parameters of the the ABO blood group system.

Usage

ABO(link.pA = "logitlink", link.pB = "logitlink", ipA = NULL, ipB = NULL,
    ipO = NULL, zero = NULL)

Arguments

Details

The parameters pA and pB are probabilities, so that pO=1-pA-pB is the third probability. The probabilities pA and pB correspond to A and B respectively, so that pO is the probability for O. It is easier to make use of initial values for pO than for pB. In documentation elsewhere I sometimes use pA=p, pB=q, pO=r.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Note

The input can be a 4-column matrix of counts, where the columns are A, B, AB, O (in order). Alternatively, the input can be a 4-column matrix of proportions (so each row adds to 1) and the weights argument is used to specify the total number of counts for each row.

Author(s)

T. W. Yee

References

Lange, K. (2002). Mathematical and Statistical Methods for Genetic Analysis, 2nd ed. New York: Springer-Verlag.

See Also

AA.Aa.aa, AB.Ab.aB.ab, A1A2A3, MNSs.

Examples

ymat <- cbind(A = 725, B = 258, AB = 72, O = 1073)  # Order matters, not the name
fit <- vglm(ymat ~ 1, ABO(link.pA = "identitylink",
                          link.pB = "identitylink"), trace = TRUE,
            crit = "coef")
coef(fit, matrix = TRUE)
Coef(fit)  # Estimated pA and pB
rbind(ymat, sum(ymat) * fitted(fit))
sqrt(diag(vcov(fit)))

Description

Fits an adjacent categories regression model to an ordered (preferably) factor response.

Usage

acat(link = "loglink", parallel = FALSE, reverse = FALSE,
     zero = NULL, ynames = FALSE, Thresh = NULL, Trev = reverse,
     Tref = if (Trev) "M" else 1, whitespace = FALSE)

Arguments

Details

In this help file the response $Y$ is assumed to be a factor with ordered values $1,2,\ldots,M+1$, so that $M$ is the number of linear/additive predictors $\eta_j$. By default, the log link is used because the ratio of two probabilities is positive.

Internally, deriv3 is called to perform symbolic differentiation and consequently this family function will struggle if $M$ becomes too large. If this occurs, try combining levels so that $M$ is effectively reduced. One idea is to aggregate levels with the fewest observations in them first.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, rrvglm and vgam.

Warning

No check is made to verify that the response is ordinal if the response is a matrix; see ordered.

Note

The response should be either a matrix of counts (with row sums that are all positive), or an ordered factor. In both cases, the y slot returned by vglm/vgam/rrvglm is the matrix of counts.

For a nominal (unordered) factor response, the multinomial logit model (multinomial) is more appropriate.

Here is an example of the usage of the parallel argument. If there are covariates x1, x2 and x3, then parallel = TRUE ~ x1 + x2 -1 and parallel = FALSE ~ x3 are equivalent. This would constrain the regression coefficients for x1 and x2 to be equal; those of the intercepts and x3 would be different.

Author(s)

Thomas W. Yee

References

Agresti, A. (2013). Categorical Data Analysis, 3rd ed. Hoboken, NJ, USA: Wiley.
Tutz, G. (2012). Regression for Categorical Data, Cambridge: Cambridge University Press.
Yee, T. W. (2010). The VGAM package for categorical data analysis. Journal of Statistical Software, 32, 1–34. doi:10.18637/jss.v032.i10.

See Also

cumulative, cratio, sratio, multinomial, CM.equid, CommonVGAMffArguments, margeff, pneumo, budworm, deriv3.

Examples

pneumo <- transform(pneumo, let = log(exposure.time))
(fit <- vglm(cbind(normal, mild, severe) ~ let, acat, pneumo))
coef(fit, matrix = TRUE)
constraints(fit)
model.matrix(fit)

Description

Compute all the single terms in the scope argument that can be added to or dropped from the model, fit those models and compute a table of the changes in fit.

Usage

## S3 method for class 'vglm'
add1(object, scope, test = c("none", "LRT"), k = 2, ...)
## S3 method for class 'vglm'
drop1(object, scope, test = c("none", "LRT"), k = 2, ...)

Arguments

Details

These functions are a direct adaptation of add1.glm and drop1.glm for vglm-class objects. For drop1 methods, a missing scope is taken to be all terms in the model. The hierarchy is respected when considering terms to be added or dropped: all main effects contained in a second-order interaction must remain, and so on. In a scope formula . means ‘what is already there’.

Compared to add1.glm and drop1.glm these functions are simpler, e.g., there is no Cp, F and Rao (score) tests, x and scale arguments. Most models do not have a deviance, however twice the log-likelihood differences are used to test the significance of terms.

The default output table gives AIC, defined as minus twice log likelihood plus $2p$ where $p$ is the rank of the model (the number of effective parameters). This is only defined up to an additive constant (like log-likelihoods).

Value

An object of class "anova" summarizing the differences in fit between the models.

Warning

In general, the same warnings in add1.glm and drop1.glm apply here. Furthermore, these functions have not been rigorously tested for all models, so treat the results cautiously and please report any bugs.

Care is needed to check that the constraint matrices of added terms are correct. Also, if object is of the form vglm(..., constraints = list(x1 = cm1, x2 = cm2)) then add1.vglm may fail because the constraints argument needs to have the constaint matrices for all terms.

Note

Most VGAM family functions do not compute a deviance, but instead the likelihood function is evaluated at the MLE. Hence a column name "Deviance" only appears for a few models; and almost always there is a column labelled "logLik".

See Also

step4vglm, vglm, extractAIC.vglm, trim.constraints, anova.vglm, backPain2, update.

Examples

data("backPain2", package = "VGAM")
summary(backPain2)
fit1 <- vglm(pain ~ x2 + x3 + x4, propodds, data = backPain2)
coef(fit1)
add1(fit1, scope = ~ x2 * x3 * x4, test = "LRT")
drop1(fit1, test = "LRT")
fit2 <- vglm(pain ~ x2 * x3 * x4, propodds, data = backPain2)
drop1(fit2)

Description

Calculates the Akaike information criterion for a fitted model object for which a log-likelihood value has been obtained.

Usage

AICvlm(object, ..., corrected = FALSE, k = 2)
   AICvgam(object, ..., k = 2)
 AICrrvglm(object, ..., k = 2)
AICdrrvglm(object, ..., k = 2)
AICqrrvglm(object, ..., k = 2)
 AICrrvgam(object, ..., k = 2)

Arguments

Details

The following formula is used for VGLMs: $-2 \mbox{log-likelihood} + k n_{par}$, where $n_{par}$ represents the number of parameters in the fitted model, and $k = 2$ for the usual AIC. One could assign $k = \log(n)$ ($n$ the number of observations) for the so-called BIC or SBC (Schwarz's Bayesian criterion). This is the function AICvlm().

This code relies on the log-likelihood being defined, and computed, for the object. When comparing fitted objects, the smaller the AIC, the better the fit. The log-likelihood and hence the AIC is only defined up to an additive constant.

Any estimated scale parameter (in GLM parlance) is used as one parameter.

For VGAMs and CAO the nonlinear effective degrees of freedom for each smoothed component is used. This formula is heuristic. These are the functions AICvgam() and AICcao().

The finite sample correction is usually recommended when the sample size is small or when the number of parameters is large. When the sample size is large their difference tends to be negligible. The correction is described in Hurvich and Tsai (1989), and is based on a (univariate) linear model with normally distributed errors.

Value

Returns a numeric value with the corresponding AIC (or BIC, or ..., depending on k).

Warning

This code has not been double-checked. The general applicability of AIC for the VGLM/VGAM classes has not been developed fully. In particular, AIC should not be run on some VGAM family functions because of violation of certain regularity conditions, etc.

Note

AIC has not been defined for QRR-VGLMs, yet.

Using AIC to compare posbinomial models with, e.g., posbernoulli.tb models, requires posbinomial(omit.constant = TRUE). See posbinomial for an example. A warning is given if it suspects a wrong omit.constant value was used.

Where defined, AICc(...) is the same as AIC(..., corrected = TRUE).

Author(s)

T. W. Yee.

References

Hurvich, C. M. and Tsai, C.-L. (1989). Regression and time series model selection in small samples, Biometrika, 76, 297–307.

See Also

VGLMs are described in vglm-class; VGAMs are described in vgam-class; RR-VGLMs are described in rrvglm-class; AIC, BICvlm, TICvlm, drop1.vglm, extractAIC.vglm.

Examples

pneumo <- transform(pneumo, let = log(exposure.time))
(fit1 <- vglm(cbind(normal, mild, severe) ~ let,
              cumulative(parallel = TRUE, reverse = TRUE), data = pneumo))
coef(fit1, matrix = TRUE)
AIC(fit1)
AICc(fit1)  # Quick way
AIC(fit1, corrected = TRUE)  # Slow way
(fit2 <- vglm(cbind(normal, mild, severe) ~ let,
              cumulative(parallel = FALSE, reverse = TRUE), data = pneumo))
coef(fit2, matrix = TRUE)
AIC(fit2)
AICc(fit2)
AIC(fit2, corrected = TRUE)

Description

Computes some arcsine–logit mixture link transformations, including their inverse and the first few derivatives.

Usage

alogitlink(theta, bvalue = NULL, taumix.logit = 1,
     tol = 1e-13, nmax = 99, inverse = FALSE, deriv = 0,
     short = TRUE, tag = FALSE, c10 = c(4, -pi))
lcalogitlink(theta, bvalue = NULL, pmix.logit = 0.01,
     tol = 1e-13, nmax = 99, inverse = FALSE, deriv = 0,
     short = TRUE, tag = FALSE, c10 = c(4, -pi))

Arguments

Details

lcalogitlink is a linear combination (LC) of asinlink and logitlink.

Value

The following holds for the LC variant. For deriv >= 0, (1 - pmix.logit) * asinlink(p, deriv = deriv) + pmix.logit * logitlink(p, deriv = deriv) when inverse = FALSE, and if inverse = TRUE then a nonlinear equation is solved for the probability, given eta. For deriv = 1, then the function returns d eta / d theta as a function of theta if inverse = FALSE, else if inverse = TRUE then it returns the reciprocal.

Warning

The default values for taumix.logit and pmix.logit may change in the future. The name and order of the arguments may change too.

Author(s)

Thomas W. Yee

References

Hauck, J. W. W. and A. Donner (1977). Wald's test as applied to hypotheses in logit analysis. Journal of the American Statistical Association, 72, 851–853.

See Also

asinlink, logitlink, Links, probitlink, clogloglink, cauchitlink, binomialff, sloglink, hdeff, https://www.cia.gov/index.html.

Examples

p <- seq(0.01, 0.99, length= 10)
alogitlink(p)
max(abs(alogitlink(alogitlink(p), inv = TRUE) - p))  # 0?

## Not run: 
par(mfrow = c(2, 2), lwd = (mylwd <- 2))
y <- seq(-4, 4, length = 100)
p <- seq(0.01, 0.99, by = 0.01)

for (d in 0:1) {
  matplot(p, cbind(logitlink(p, deriv = d), probitlink(p, deriv = d)),
          type = "n", col = "blue", ylab = "transformation",
          las = 1, main = if (d == 0) "Some probability link functions"
          else "First derivative")
  lines(p,   logitlink(p, deriv = d), col = "green")
  lines(p,  probitlink(p, deriv = d), col = "blue")
  lines(p, clogloglink(p, deriv = d), col = "tan")
  lines(p,  alogitlink(p, deriv = d), col = "red3")
  if (d == 0) {
    abline(v = 0.5, h = 0, lty = "dashed")
    legend(0, 4.5, c("logitlink", "probitlink", "clogloglink",
           "alogitlink"), lwd = mylwd,
           col = c("green", "blue", "tan", "red3"))
  } else
    abline(v = 0.5, lwd = 0.5, col = "gray")
}

for (d in 0) {
  matplot(y, cbind( logitlink(y, deriv = d, inverse = TRUE),
                   probitlink(y, deriv = d, inverse = TRUE)),
          type  = "n", col = "blue", xlab = "transformation", ylab = "p",
          main = if (d == 0) "Some inverse probability link functions"
          else "First derivative", las=1)
  lines(y,   logitlink(y, deriv = d, inverse = TRUE), col = "green")
  lines(y,  probitlink(y, deriv = d, inverse = TRUE), col = "blue")
  lines(y, clogloglink(y, deriv = d, inverse = TRUE), col = "tan")
  lines(y,  alogitlink(y, deriv = d, inverse = TRUE), col = "red3")
  if (d == 0) {
      abline(h = 0.5, v = 0, lwd = 0.5, col = "gray")
      legend(-4, 1, c("logitlink", "probitlink", "clogloglink",
             "alogitlink"), lwd = mylwd,
             col = c("green", "blue", "tan", "red3"))
  }
}
par(lwd = 1)

## End(Not run)

Description

Return the altered, inflated, truncated and deflated values in a GAITD regression object, else test whether the model is altered, inflated, truncated or deflated.

Usage

altered(object, ...)
inflated(object, ...)
truncated(object, ...)
is.altered(object, ...)
is.deflated(object, ...)
is.inflated(object, ...)
is.truncated(object, ...)

Arguments

Details

Yee and Ma (2023) propose GAITD regression where values from four (or seven since there are parametric and nonparametric forms) disjoint sets are referred to as special. These extractor functions return one set each; they are the alter, inflate, truncate, deflate (and sometimes max.support) arguments from the family function.

Value

Returns one type of ‘special’ sets associated with GAITD regression. This is a vector, else a list for truncation. All three sets are returned by specialsvglm.

Warning

Some of these functions are subject to change. Only family functions beginning with "gaitd" will work with these functions, hence zipoisson fits will return FALSE or empty values.

References

Yee, T. W. and Ma, C. (2024). Generally altered, inflated, truncated and deflated regression. Statistical Science, 39 (in press).

See Also

vglm, vglm-class, specialsvglm, gaitdpoisson, gaitdlog, gaitdzeta, Gaitdpois.

Examples

## Not run: 
abdata <- data.frame(y = 0:7, w = c(182, 41, 12, 2, 2, 0, 0, 1))
fit1 <- vglm(y ~ 1, gaitdpoisson(a.mix = 0),
             data = abdata, weight = w, subset = w > 0)
specials(fit1)  # All three sets
altered(fit1)  # Subject to change
inflated(fit1)  # Subject to change
truncated(fit1)  # Subject to change
is.altered(fit1)
is.inflated(fit1)
is.truncated(fit1)  
## End(Not run)

Description

Binomial quantile regression estimated by maximizing an asymmetric likelihood function.

Usage

amlbinomial(w.aml = 1, parallel = FALSE, digw = 4, link = "logitlink")

Arguments

Details

The general methodology behind this VGAM family function is given in Efron (1992) and full details can be obtained there. This model is essentially a logistic regression model (see binomialff) but the usual deviance is replaced by an asymmetric squared error loss function; it is multiplied by $w.aml$ for positive residuals. The solution is the set of regression coefficients that minimize the sum of these deviance-type values over the data set, weighted by the weights argument (so that it can contain frequencies). Newton-Raphson estimation is used here.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Warning

If w.aml has more than one value then the value returned by deviance is the sum of all the (weighted) deviances taken over all the w.aml values. See Equation (1.6) of Efron (1992).

Note

On fitting, the extra slot has list components "w.aml" and "percentile". The latter is the percent of observations below the “w-regression plane”, which is the fitted values. Also, the individual deviance values corresponding to each element of the argument w.aml is stored in the extra slot.

For amlbinomial objects, methods functions for the generic functions qtplot and cdf have not been written yet.

See amlpoisson about comments on the jargon, e.g., expectiles etc.

In this documentation the word quantile can often be interchangeably replaced by expectile (things are informal here).

Author(s)

Thomas W. Yee

References

Efron, B. (1992). Poisson overdispersion estimates based on the method of asymmetric maximum likelihood. Journal of the American Statistical Association, 87, 98–107.

See Also

amlpoisson, amlexponential, amlnormal, extlogF1, alaplace1, denorm.

Examples

# Example: binomial data with lots of trials per observation
set.seed(1234)
sizevec <- rep(100, length = (nn <- 200))
mydat <- data.frame(x = sort(runif(nn)))
mydat <- transform(mydat,
                   prob = logitlink(-0 + 2.5*x + x^2, inverse = TRUE))
mydat <- transform(mydat, y = rbinom(nn, size = sizevec, prob = prob))
(fit <- vgam(cbind(y, sizevec - y) ~ s(x, df = 3),
             amlbinomial(w = c(0.01, 0.2, 1, 5, 60)),
             mydat, trace = TRUE))
fit@extra

## Not run: 
par(mfrow = c(1,2))
# Quantile plot
with(mydat, plot(x, jitter(y), col = "blue", las = 1, main =
     paste(paste(round(fit@extra$percentile, digits = 1), collapse = ", "),
           "percentile-expectile curves")))
with(mydat, matlines(x, 100 * fitted(fit), lwd = 2, col = "blue", lty=1))

# Compare the fitted expectiles with the quantiles
with(mydat, plot(x, jitter(y), col = "blue", las = 1, main =
     paste(paste(round(fit@extra$percentile, digits = 1), collapse = ", "),
           "percentile curves are red")))
with(mydat, matlines(x, 100 * fitted(fit), lwd = 2, col = "blue", lty = 1))

for (ii in fit@extra$percentile)
    with(mydat, matlines(x, 100 *
         qbinom(p = ii/100, size = sizevec, prob = prob) / sizevec,
                  col = "red", lwd = 2, lty = 1))

## End(Not run)

Description

Exponential expectile regression estimated by maximizing an asymmetric likelihood function.

Usage

amlexponential(w.aml = 1, parallel = FALSE, imethod = 1, digw = 4,
               link = "loglink")

Arguments

Details

The general methodology behind this VGAM family function is given in Efron (1992) and full details can be obtained there. This model is essentially an exponential regression model (see exponential) but the usual deviance is replaced by an asymmetric squared error loss function; it is multiplied by $w.aml$ for positive residuals. The solution is the set of regression coefficients that minimize the sum of these deviance-type values over the data set, weighted by the weights argument (so that it can contain frequencies). Newton-Raphson estimation is used here.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Warning

Note that the link argument of exponential and amlexponential are currently different: one is the rate parameter and the other is the mean (expectile) parameter.

If w.aml has more than one value then the value returned by deviance is the sum of all the (weighted) deviances taken over all the w.aml values. See Equation (1.6) of Efron (1992).

Note

On fitting, the extra slot has list components "w.aml" and "percentile". The latter is the percent of observations below the “w-regression plane”, which is the fitted values. Also, the individual deviance values corresponding to each element of the argument w.aml is stored in the extra slot.

For amlexponential objects, methods functions for the generic functions qtplot and cdf have not been written yet.

See amlpoisson about comments on the jargon, e.g., expectiles etc.

In this documentation the word quantile can often be interchangeably replaced by expectile (things are informal here).

Author(s)

Thomas W. Yee

References

Efron, B. (1992). Poisson overdispersion estimates based on the method of asymmetric maximum likelihood. Journal of the American Statistical Association, 87, 98–107.

See Also

exponential, amlbinomial, amlpoisson, amlnormal, extlogF1, alaplace1, lms.bcg, deexp.

Examples

nn <- 2000
mydat <- data.frame(x = seq(0, 1, length = nn))
mydat <- transform(mydat,
                   mu = loglink(-0 + 1.5*x + 0.2*x^2, inverse = TRUE))
mydat <- transform(mydat, mu = loglink(0 - sin(8*x), inverse = TRUE))
mydat <- transform(mydat,  y = rexp(nn, rate = 1/mu))
(fit <- vgam(y ~ s(x, df=5), amlexponential(w=c(0.001, 0.1, 0.5, 5, 60)),
             mydat, trace = TRUE))
fit@extra

## Not run:  # These plots are against the sqrt scale (to increase clarity)
par(mfrow = c(1,2))
# Quantile plot
with(mydat, plot(x, sqrt(y), col = "blue", las = 1, main =
     paste(paste(round(fit@extra$percentile, digits = 1), collapse=", "),
           "percentile-expectile curves")))
with(mydat, matlines(x, sqrt(fitted(fit)), lwd = 2, col = "blue", lty=1))

# Compare the fitted expectiles with the quantiles
with(mydat, plot(x, sqrt(y), col = "blue", las = 1, main =
     paste(paste(round(fit@extra$percentile, digits = 1), collapse=", "),
           "percentile curves are orange")))
with(mydat, matlines(x, sqrt(fitted(fit)), lwd = 2, col = "blue", lty=1))

for (ii in fit@extra$percentile)
  with(mydat, matlines(x, sqrt(qexp(p = ii/100, rate = 1/mu)),
                       col = "orange")) 
## End(Not run)

Description

Asymmetric least squares, a special case of maximizing an asymmetric likelihood function of a normal distribution. This allows for expectile/quantile regression using asymmetric least squares error loss.

Usage

amlnormal(w.aml = 1, parallel = FALSE, lexpectile = "identitylink",
          iexpectile = NULL, imethod = 1, digw = 4)

Arguments

Details

This is an implementation of Efron (1991) and full details can be obtained there. Equation numbers below refer to that article. The model is essentially a linear model (see lm), however, the asymmetric squared error loss function for a residual $r$ is $r^2$ if $r \leq 0$ and $w r^2$ if $r > 0$. The solution is the set of regression coefficients that minimize the sum of these over the data set, weighted by the weights argument (so that it can contain frequencies). Newton-Raphson estimation is used here.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Note

On fitting, the extra slot has list components "w.aml" and "percentile". The latter is the percent of observations below the “w-regression plane”, which is the fitted values.

One difficulty is finding the w.aml value giving a specified percentile. One solution is to fit the model within a root finding function such as uniroot; see the example below.

For amlnormal objects, methods functions for the generic functions qtplot and cdf have not been written yet.

See the note in amlpoisson on the jargon, including expectiles and regression quantiles.

The deviance slot computes the total asymmetric squared error loss (2.5). If w.aml has more than one value then the value returned by the slot is the sum taken over all the w.aml values.

This VGAM family function could well be renamed amlnormal() instead, given the other function names amlpoisson, amlbinomial, etc.

In this documentation the word quantile can often be interchangeably replaced by expectile (things are informal here).

Author(s)

Thomas W. Yee

References

Efron, B. (1991). Regression percentiles using asymmetric squared error loss. Statistica Sinica, 1, 93–125.

See Also

amlpoisson, amlbinomial, amlexponential, bmi.nz, extlogF1, alaplace1, denorm, lms.bcn and similar variants are alternative methods for quantile regression.

Examples

## Not run: 
# Example 1
ooo <- with(bmi.nz, order(age))
bmi.nz <- bmi.nz[ooo, ]  # Sort by age
(fit <- vglm(BMI ~ sm.bs(age), amlnormal(w.aml = 0.1), bmi.nz))
fit@extra  # Gives the w value and the percentile
coef(fit, matrix = TRUE)

# Quantile plot
with(bmi.nz, plot(age, BMI, col = "blue", main =
     paste(round(fit@extra$percentile, digits = 1),
           "expectile-percentile curve")))
with(bmi.nz, lines(age, c(fitted(fit)), col = "black"))

# Example 2
# Find the w values that give the 25, 50 and 75 percentiles
find.w <- function(w, percentile = 50) {
  fit2 <- vglm(BMI ~ sm.bs(age), amlnormal(w = w), data = bmi.nz)
  fit2@extra$percentile - percentile
}
# Quantile plot
with(bmi.nz, plot(age, BMI, col = "blue", las = 1, main =
     "25, 50 and 75 expectile-percentile curves"))
for (myp in c(25, 50, 75)) {
# Note: uniroot() can only find one root at a time
  bestw <- uniroot(f = find.w, interval = c(1/10^4, 10^4),
                   percentile = myp)
  fit2 <- vglm(BMI ~ sm.bs(age), amlnormal(w = bestw$root), bmi.nz)
  with(bmi.nz, lines(age, c(fitted(fit2)), col = "orange"))
}

# Example 3; this is Example 1 but with smoothing splines and
# a vector w and a parallelism assumption.
ooo <- with(bmi.nz, order(age))
bmi.nz <- bmi.nz[ooo, ]  # Sort by age
fit3 <- vgam(BMI ~ s(age, df = 4), data = bmi.nz, trace = TRUE,
             amlnormal(w = c(0.1, 1, 10), parallel = TRUE))
fit3@extra  # The w values, percentiles and weighted deviances

# The linear components of the fit; not for human consumption:
coef(fit3, matrix = TRUE)

# Quantile plot
with(bmi.nz, plot(age, BMI, col="blue", main =
  paste(paste(round(fit3@extra$percentile, digits = 1), collapse = ", "),
        "expectile-percentile curves")))
with(bmi.nz, matlines(age, fitted(fit3), col = 1:fit3@extra$M, lwd = 2))
with(bmi.nz, lines(age, c(fitted(fit )), col = "black"))  # For comparison

## End(Not run)

Description

Poisson quantile regression estimated by maximizing an asymmetric likelihood function.

Usage

amlpoisson(w.aml = 1, parallel = FALSE, imethod = 1, digw = 4,
           link = "loglink")

Arguments

Details

This method was proposed by Efron (1992) and full details can be obtained there. The model is essentially a Poisson regression model (see poissonff) but the usual deviance is replaced by an asymmetric squared error loss function; it is multiplied by $w.aml$ for positive residuals. The solution is the set of regression coefficients that minimize the sum of these deviance-type values over the data set, weighted by the weights argument (so that it can contain frequencies). Newton-Raphson estimation is used here.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Warning

If w.aml has more than one value then the value returned by deviance is the sum of all the (weighted) deviances taken over all the w.aml values. See Equation (1.6) of Efron (1992).

Note

On fitting, the extra slot has list components "w.aml" and "percentile". The latter is the percent of observations below the “w-regression plane”, which is the fitted values. Also, the individual deviance values corresponding to each element of the argument w.aml is stored in the extra slot.

For amlpoisson objects, methods functions for the generic functions qtplot and cdf have not been written yet.

About the jargon, Newey and Powell (1987) used the name expectiles for regression surfaces obtained by asymmetric least squares. This was deliberate so as to distinguish them from the original regression quantiles of Koenker and Bassett (1978). Efron (1991) and Efron (1992) use the general name regression percentile to apply to all forms of asymmetric fitting. Although the asymmetric maximum likelihood method very nearly gives regression percentiles in the strictest sense for the normal and Poisson cases, the phrase quantile regression is used loosely in this VGAM documentation.

In this documentation the word quantile can often be interchangeably replaced by expectile (things are informal here).

Author(s)

Thomas W. Yee

References

Efron, B. (1991). Regression percentiles using asymmetric squared error loss. Statistica Sinica, 1, 93–125.

Efron, B. (1992). Poisson overdispersion estimates based on the method of asymmetric maximum likelihood. Journal of the American Statistical Association, 87, 98–107.

Koenker, R. and Bassett, G. (1978). Regression quantiles. Econometrica, 46, 33–50.

Newey, W. K. and Powell, J. L. (1987). Asymmetric least squares estimation and testing. Econometrica, 55, 819–847.

See Also

amlnormal, amlbinomial, extlogF1, alaplace1.

Examples

set.seed(1234)
mydat <- data.frame(x = sort(runif(nn <- 200)))
mydat <- transform(mydat, y = rpois(nn, exp(0 - sin(8*x))))
(fit <- vgam(y ~ s(x), fam = amlpoisson(w.aml = c(0.02, 0.2, 1, 5, 50)),
             mydat, trace = TRUE))
fit@extra

## Not run: 
# Quantile plot
with(mydat, plot(x, jitter(y), col = "blue", las = 1, main =
     paste(paste(round(fit@extra$percentile, digits = 1), collapse = ", "),
           "percentile-expectile curves")))
with(mydat, matlines(x, fitted(fit), lwd = 2)) 
## End(Not run)

Description

Compute an analysis of deviance table for one or more vector generalized linear model fits.

Usage

## S3 method for class 'vglm'
anova(object, ..., type = c("II", "I", "III", 2, 1, 3),
     test = c("LRT", "none"), trydev = TRUE, silent = TRUE)

Arguments

Details

anova.vglm is intended to be similar to anova.glm so specifying a single object and type = 1 gives a sequential analysis of deviance table for that fit. By analysis of deviance, it is meant loosely that if the deviance of the model is not defined or implemented, then twice the difference between the log-likelihoods of two nested models remains asymptotically chi-squared distributed with degrees of freedom equal to the difference in the number of parameters of the two models. Of course, the usual regularity conditions are assumed to hold. For Type I, the analysis of deviance table has the reductions in the residual deviance as each term of the formula is added in turn are given in as the rows of a table, plus the residual deviances themselves. Type I or sequential tests (as in anova.glm). are computationally the easiest of the three methods. For this, the order of the terms is important, and the each term is added sequentially from first to last.

The Anova() function in car allows for testing Type II and Type III (SAS jargon) hypothesis tests, although the definitions used are not precisely that of SAS. As car notes, Type I rarely test interesting hypotheses in unbalanced designs. Type III enter each term last, keeping all the other terms in the model.

Type II tests, according to SAS, add the term after all other terms have been added to the model except terms that contain the effect being tested; an effect is contained in another effect if it can be derived by deleting variables from the latter effect. Type II tests are currently the default.

As in anova.glm, but not as Anova.glm() in car, if more than one object is specified, then the table has a row for the residual degrees of freedom and deviance for each model. For all but the first model, the change in degrees of freedom and deviance is also given. (This only makes statistical sense if the models are nested.) It is conventional to list the models from smallest to largest, but this is up to the user. It is necessary to have type = 1 with more than one objects are specified.

See anova.glm for more details and warnings. The VGAM package now implements full likelihood models only, therefore no dispersion parameters are estimated.

Value

An object of class "anova" inheriting from class "data.frame".

Warning

See anova.glm. Several VGAM family functions implement distributions which do not satisfying the usual regularity conditions needed for the LRT to work. No checking or warning is given for these.

As car says, be careful of Type III tests because they violate marginality. Type II tests (the default) do not have this problem.

Note

It is possible for this function to stop when type = 2 or 3, e.g., anova(vglm(cans ~ myfactor, poissonff, data = boxcar)) where myfactor is a factor.

The code was adapted directly from anova.glm and Anova.glm() in car by T. W. Yee. Hence the Type II and Type III tests do not correspond precisely with the SAS definition.

See Also

anova.glm, stat.anova, stats:::print.anova, Anova.glm() in car if car is installed, vglm, lrtest, add1.vglm, drop1.vglm, lrt.stat.vlm, score.stat.vlm, wald.stat.vlm, backPain2, update.

Examples

# Example 1: a proportional odds model fitted to pneumo.
set.seed(1)
pneumo <- transform(pneumo, let = log(exposure.time), x3 = runif(8))
fit1 <- vglm(cbind(normal, mild, severe) ~ let     , propodds, pneumo)
fit2 <- vglm(cbind(normal, mild, severe) ~ let + x3, propodds, pneumo)
fit3 <- vglm(cbind(normal, mild, severe) ~ let + x3, cumulative, pneumo)
anova(fit1, fit2, fit3, type = 1)  # Remember to specify 'type'!!
anova(fit2)
anova(fit2, type = "I")
anova(fit2, type = "III")

# Example 2: a proportional odds model fitted to backPain2.
data("backPain2", package = "VGAM")
summary(backPain2)
fitlogit <- vglm(pain ~ x2 * x3 * x4, propodds, data = backPain2)
coef(fitlogit)
anova(fitlogit)
anova(fitlogit, type = "I")
anova(fitlogit, type = "III")

Description

Maximum likelihood estimation of the three-parameter AR-1 model

Usage

AR1(ldrift = "identitylink", lsd  = "loglink", lvar = "loglink", lrho = "rhobitlink",
    idrift  = NULL, isd  = NULL, ivar = NULL, irho = NULL, imethod = 1,
    ishrinkage = 0.95, type.likelihood = c("exact", "conditional"),
    type.EIM  = c("exact", "approximate"), var.arg = FALSE, nodrift = FALSE,
    print.EIM = FALSE, zero = c(if (var.arg) "var" else "sd", "rho"))

Arguments

Details

The AR-1 model implemented here has

$Y_1 \sim N(\mu, \sigma^2 / (1-\rho^2)),$

and

$Y_i = \mu^* + \rho Y_{i-1} + e_i,$

where the $e_i$ are i.i.d. Normal(0, sd = $\sigma$) random variates.

Here are a few notes: (1). A test for weak stationarity might be to verify whether $1/\rho$ lies outside the unit circle. (2). The mean of all the $Y_i$ is $\mu^* /(1-\rho)$ and these are returned as the fitted values. (3). The correlation of all the $Y_i$ with $Y_{i-1}$ is $\rho$. (4). The default link function ensures that $-1 < \rho < 1$.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Warning

Monitoring convergence is urged, i.e., set trace = TRUE.

Moreover, if the exact EIMs are used, set print.EIM = TRUE to compare the computed exact to the approximate EIM.

Under the VGLM/VGAM approach, parameters can be modelled in terms of covariates. Particularly, if the standard deviation of the white noise is modelled in this way, then type.EIM = "exact" may certainly lead to unstable results. The reason is that white noise is a stationary process, and consequently, its variance must remain as a constant. Consequently, the use of variates to model this parameter contradicts the assumption of stationary random components to compute the exact EIMs proposed by Porat and Friedlander (1987).

To prevent convergence issues in such cases, this family function internally verifies whether the variance of the white noise remains as a constant at each Fisher scoring iteration. If this assumption is violated and type.EIM = "exact" is set, then AR1 automatically shifts to type.EIM = "approximate". Also, a warning is accordingly displayed.

Note

Multiple responses are handled. The mean is returned as the fitted values.

Author(s)

Victor Miranda (exact method) and Thomas W. Yee (approximate method).

References

Porat, B. and Friedlander, B. (1987). The Exact Cramer-Rao Bond for Gaussian Autoregressive Processes. IEEE Transactions on Aerospace and Electronic Systems, AES-23(4), 537–542.

See Also

AR1EIM, vglm.control, dAR1, arima.sim.

Examples

## Not run: 
### Example 1: using  arima.sim() to generate a 0-mean stationary time series.
nn <- 500
tsdata <- data.frame(x2 =  runif(nn))
ar.coef.1 <- rhobitlink(-1.55, inverse = TRUE)  # Approx -0.65
ar.coef.2 <- rhobitlink( 1.0, inverse = TRUE)   # Approx  0.50
set.seed(1)
tsdata  <- transform(tsdata,
              index = 1:nn,
              TS1 = arima.sim(nn, model = list(ar = ar.coef.1),
                              sd = exp(1.5)),
              TS2 = arima.sim(nn, model = list(ar = ar.coef.2),
                              sd = exp(1.0 + 1.5 * x2)))

### An autoregressive intercept--only model.   ###
### Using the exact EIM, and "nodrift = TRUE"  ###
fit1a <- vglm(TS1 ~ 1, data = tsdata, trace = TRUE,
              AR1(var.arg = FALSE, nodrift = TRUE,
                  type.EIM = "exact",
                  print.EIM = FALSE),
              crit = "coefficients")
Coef(fit1a)
summary(fit1a)

### Two responses. Here, the white noise standard deviation of TS2   ###
### is modelled in terms of 'x2'. Also, 'type.EIM = exact'.  ###
fit1b <- vglm(cbind(TS1, TS2) ~ x2,
              AR1(zero = NULL, nodrift = TRUE,
                  var.arg = FALSE,
                  type.EIM = "exact"),
              constraints = list("(Intercept)" = diag(4),
                                 "x2" = rbind(0, 0, 1, 0)),
              data = tsdata, trace = TRUE, crit = "coefficients")
coef(fit1b, matrix = TRUE)
summary(fit1b)

### Example 2: another stationary time series
nn     <- 500
my.rho <- rhobitlink(1.0, inverse = TRUE)
my.mu  <- 1.0
my.sd  <- exp(1)
tsdata  <- data.frame(index = 1:nn, TS3 = runif(nn))

set.seed(2)
for (ii in 2:nn)
  tsdata$TS3[ii] <- my.mu/(1 - my.rho) +
                    my.rho * tsdata$TS3[ii-1] + rnorm(1, sd = my.sd)
tsdata <- tsdata[-(1:ceiling(nn/5)), ]  # Remove the burn-in data:

### Fitting an AR(1). The exact EIMs are used.
fit2a <- vglm(TS3 ~ 1, AR1(type.likelihood = "exact",  # "conditional",
                                type.EIM = "exact"),
              data = tsdata, trace = TRUE, crit = "coefficients")

Coef(fit2a)
summary(fit2a)      # SEs are useful to know

Coef(fit2a)["rho"]    # Estimate of rho, for intercept-only models
my.rho                # The 'truth' (rho)
Coef(fit2a)["drift"]  # Estimate of drift, for intercept-only models
my.mu /(1 - my.rho)   # The 'truth' (drift)

## End(Not run)

Description

Computation of the exact Expected Information Matrix of the Autoregressive process of order-$1$ (AR($1$)) with Gaussian white noise and stationary random components.

Usage

AR1EIM(x = NULL, var.arg = NULL, p.drift = NULL,
       WNsd = NULL, ARcoeff1 = NULL, eps.porat = 1e-2)

Arguments

Details

This function implements the algorithm of Porat and Friedlander (1986) to recursively compute the exact expected information matrix (EIM) of Gaussian time series with stationary random components.

By default, when the VGLM/VGAM family function AR1 is used to fit an AR($1$) model via vglm, Fisher scoring is executed using the approximate EIM for the AR process. However, this model can also be fitted using the exact EIMs computed by AR1EIM.

Given $N$ consecutive data points, ${y_{0}, y_{1}, \ldots, y_{N - 1} }$ with probability density $f(\boldsymbol{y})$, the Porat and Friedlander algorithm calculates the EIMs $[J_{n-1}(\boldsymbol{\theta})]$, for all $1 \leq n \leq N$. This is done based on the Levinson-Durbin algorithm for computing the orthogonal polynomials of a Toeplitz matrix. In particular, for the AR($1$) model, the vector of parameters to be estimated under the VGAM/VGLM approach is

$\boldsymbol{\eta} = (\mu^{*}, \log(\sigma^2), rhobit(\rho)),$

where $\sigma^2$ is the variance of the white noise and $mu^{*}$ is the drift parameter (See AR1 for further details on this).

Consequently, for each observation $n = 1, \ldots, N$, the EIM, $J_{n}(\boldsymbol{\theta})$, has dimension $3 \times 3$, where the diagonal elements are:

$J_{[n, 1, 1]} = E[ -\partial^2 \log f(\boldsymbol{y}) / \partial ( \mu^{*} )^2 ],$

$J_{[n, 2, 2]} = E[ -\partial^2 \log f(\boldsymbol{y}) / \partial (\sigma^2)^2 ],$

and

$J_{[n, 3, 3]} = E[ -\partial^2 \log f(\boldsymbol{y}) / \partial ( \rho )^2 ].$

As for the off-diagonal elements, one has the usual entries, i.e.,

$J_{[n, 1, 2]} = J_{[n, 2, 1]} = E[ -\partial^2 \log f(\boldsymbol{y}) / \partial \sigma^2 \partial \rho],$

etc.

If var.arg = FALSE, then $\sigma$ instead of $\sigma^2$ is estimated. Therefore, $J_{[n, 2, 2]}$, $J_{[n, 1, 2]}$, etc., are correspondingly replaced.

Once these expected values are internally computed, they are returned in an array of dimension $N \times 1 \times 6$, of the form

$J[, 1, ] = [ J_{[ , 1, 1]}, J_{[ , 2, 2]}, J_{[ , 3, 3]}, J_{[ , 1, 2]}, J_{[, 2, 3]}, J_{[ , 1, 3]} ].$

AR1EIM handles multiple time series, say $NOS$. If this happens, then it accordingly returns an array of dimension $N \times NOS \times 6$. Here, $J[, k, ]$, for $k = 1, \ldots, NOS$, is a matrix of dimension $N \times 6$, which stores the EIMs for the $k^{th}$th response, as above, i.e.,

$J[, k, ] = [ J_{[ , 1, 1]}, J_{[ , 2, 2]}, J_{[ , 3, 3]}, \ldots ],$

the bandwith form, as per required by AR1.

Value

An array of dimension $N \times NOS \times 6$, as above.

This array stores the EIMs calculated from the joint density as a function of

$\boldsymbol{\theta} = (\mu^*, \sigma^2, \rho).$

Nevertheless, note that, under the VGAM/VGLM approach, the EIMs must be correspondingly calculated in terms of the linear predictors, $\boldsymbol{\eta}$.

Asymptotic behaviour of the algorithm

For large enough $n$, the EIMs, $J_n(\boldsymbol{\theta})$, become approximately linear in $n$. That is, for some $n_0$,

$J_n(\boldsymbol{\theta}) \equiv J_{n_0}(\boldsymbol{\theta}) + (n - n_0) \bar{J}(\boldsymbol{\theta}),~~~~~~(**)$

where $\bar{J}(\boldsymbol{\theta})$ is a constant matrix.

This relationsihip is internally considered if a proper value of $n_0$ is determined. Different ways can be adopted to find $n_0$. In AR1EIM, this is done by checking the difference between the internally estimated variances and the entered ones at WNsd. If this difference is less than eps.porat at some iteration, say at iteration $n_0$, then AR1EIM takes $\bar{J}(\boldsymbol{\theta})$ as the last computed increment of $J_n(\boldsymbol{\theta})$, and extraplotates $J_k(\boldsymbol{\theta})$, for all $k \geq n_0$ using $(*)$. Else, the algorithm will complete the iterations for $1 \leq n \leq N$.

Finally, note that the rate of convergence reasonably decreases if the asymptotic relationship $(*)$ is used to compute $J_k(\boldsymbol{\theta})$, $k \geq n_0$. Normally, the number of operations involved on this algorithm is proportional to $N^2$.

See Porat and Friedlander (1986) for full details on the asymptotic behaviour of the algorithm.

Warning

Arguments WNsd, and ARcoeff1 are matrices of dimension $N \times NOS$. Else, these arguments are accordingly recycled.

Note

For simplicity, one can assume that the time series analyzed has a 0-mean. Consequently, where the family function AR1 calls AR1EIM to compute the EIMs, the argument p.drift is internally set to zero-vector, whereas x is centered by subtracting its mean value.

Author(s)

V. Miranda and T. W. Yee.

References

Porat, B. and Friedlander, B. (1986). Computation of the Exact Information Matrix of Gaussian Time Series with Stationary Random Components. IEEE Transactions on Acoustics, Speech, and Signal Processing, 54(1), 118–130.

See Also

AR1.

Examples

set.seed(1)
  nn <- 500
  ARcoeff1 <- c(0.3, 0.25)        # Will be recycled.
  WNsd     <- c(exp(1), exp(1.5)) # Will be recycled.
  p.drift  <- c(0, 0)             # Zero-mean gaussian time series.

  ### Generate two (zero-mean) AR(1) processes ###
  ts1 <- p.drift[1]/(1 - ARcoeff1[1]) +
                   arima.sim(model = list(ar = ARcoeff1[1]), n = nn,
                   sd = WNsd[1])
  ts2 <- p.drift[2]/(1 - ARcoeff1[2]) +
                   arima.sim(model = list(ar = ARcoeff1[2]), n = nn,
                   sd = WNsd[2])

  ARdata <- matrix(cbind(ts1, ts2), ncol = 2)


  ### Compute the exact EIMs: TWO responses. ###
  ExactEIM <- AR1EIM(x = ARdata, var.arg = FALSE, p.drift = p.drift,
                           WNsd = WNsd, ARcoeff1 = ARcoeff1)

  ### For response 1:
  head(ExactEIM[, 1 ,])      # NOTICE THAT THIS IS A (nn x 6) MATRIX!

  ### For response 2:
  head(ExactEIM[, 2 ,])      # NOTICE THAT THIS IS A (nn x 6) MATRIX!

Description

Computes the arcsine link, including its inverse and the first few derivatives.

Usage

asinlink(theta, bvalue = NULL, inverse = FALSE,
   deriv = 0, short = TRUE, tag = FALSE, c10 = c(4, -pi))

Arguments

Details

Function alogitlink gives some motivation for this link. However, the problem with this link is that it is bounded by default between (-pi, pi) so that it can be unsuitable for regression. This link is a scaled and centred CDF of the arcsine distribution. The centring is chosen so that asinlink(0.5) is 0, and the scaling is chosen so that asinlink(0.5, deriv = 1) and logitlink(0.5, deriv = 1) are equal (the value 4 actually), hence this link will operate similar to the logitlink when close to 0.5.

Value

Similar to logitlink but using different formulas.

Warning

It is possible that the scaling might change in the future.

Author(s)

Thomas W. Yee

See Also

logitlink, alogitlink, Links, probitlink, clogloglink, cauchitlink, binomialff, sloglink, hdeff.

Examples

p <- seq(0.01, 0.99, length= 10)
asinlink(p)
max(abs(asinlink(asinlink(p), inv = TRUE) - p))  # 0?

## Not run: 
par(mfrow = c(2, 2), lwd = (mylwd <- 2))
y <- seq(-4, 4, length = 100)
p <- seq(0.01, 0.99, by = 0.01)

for (d in 0:1) {
  matplot(p, cbind(logitlink(p, deriv = d), probitlink(p, deriv = d)),
          type = "n", col = "blue", ylab = "transformation",
          log = ifelse(d == 1, "y", ""),
          las = 1, main = if (d == 0) "Some probability link functions"
          else "First derivative")
  lines(p,   logitlink(p, deriv = d), col = "green")
  lines(p,  probitlink(p, deriv = d), col = "blue")
  lines(p, clogloglink(p, deriv = d), col = "tan")
  lines(p,    asinlink(p, deriv = d), col = "red3")
  if (d == 0) {
    abline(v = 0.5, h = 0, lty = "dashed")
    legend(0, 4.5, c("logitlink", "probitlink", "clogloglink",
           "asinlink"), lwd = mylwd,
           col = c("green", "blue", "tan", "red3"))
  } else
    abline(v = 0.5, lwd = 0.5, col = "gray")
}

for (d in 0) {
  matplot(y, cbind( logitlink(y, deriv = d, inverse = TRUE),
                   probitlink(y, deriv = d, inverse = TRUE)),
          type  = "n", col = "blue", xlab = "transformation", ylab = "p",
          main = if (d == 0) "Some inverse probability link functions"
          else "First derivative", las=1)
  lines(y,   logitlink(y, deriv = d, inverse = TRUE), col = "green")
  lines(y,  probitlink(y, deriv = d, inverse = TRUE), col = "blue")
  lines(y, clogloglink(y, deriv = d, inverse = TRUE), col = "tan")
  lines(y,    asinlink(y, deriv = d, inverse = TRUE), col = "red3")
  if (d == 0) {
      abline(h = 0.5, v = 0, lwd = 0.5, col = "gray")
      legend(-4, 1, c("logitlink", "probitlink", "clogloglink",
             "asinlink"), lwd = mylwd,
             col = c("green", "blue", "tan", "red3"))
  }
}
par(lwd = 1)

## End(Not run)

Description

Undergraduate student enrolments at the University of Auckland in 1990.

Usage

data(auuc)

Format

A data frame with 4 observations on the following 5 variables.

Commerce

a numeric vector of counts.

Arts

a numeric vector of counts.

SciEng

a numeric vector of counts.

Law

a numeric vector of counts.

Medicine

a numeric vector of counts.