Version:	1.1-13
Date:	2025-02-11
Title:	Vector Generalized Linear and Additive Models
Author:	Thomas Yee [aut, cre], Cleve Moler [ctb] (LINPACK routines in src)
Maintainer:	Thomas Yee <t.yee@auckland.ac.nz>
Depends:	R (≥ 4.1.0), methods, stats, stats4, splines
Suggests:	VGAMextra, MASS, mgcv
Enhances:	VGAMdata
Description:	An implementation of about 6 major classes of statistical regression models. The central algorithm is Fisher scoring and iterative reweighted least squares. At the heart of this package are the vector generalized linear and additive model (VGLM/VGAM) classes. VGLMs can be loosely thought of as multivariate GLMs. VGAMs are data-driven VGLMs that use smoothing. The book "Vector Generalized Linear and Additive Models: With an Implementation in R" (Yee, 2015) <doi:10.1007/978-1-4939-2818-7> gives details of the statistical framework and the package. Currently only fixed-effects models are implemented. Many (100+) models and distributions are estimated by maximum likelihood estimation (MLE) or penalized MLE. The other classes are RR-VGLMs (reduced-rank VGLMs), quadratic RR-VGLMs, doubly constrained RR-VGLMs, quadratic RR-VGLMs, reduced-rank VGAMs, RCIMs (row-column interaction models)—these classes perform constrained and unconstrained quadratic ordination (CQO/UQO) models in ecology, as well as constrained additive ordination (CAO). Hauck-Donner effect detection is implemented. Note that these functions are subject to change; see the NEWS and ChangeLog files for latest changes.
License:	GPL-3
URL:	https://www.stat.auckland.ac.nz/~yee/VGAM/
NeedsCompilation:	yes
BuildVignettes:	yes
LazyLoad:	yes
LazyData:	yes
Packaged:	2025-02-11 08:14:16 UTC; tyee001
Repository:	CRAN
Date/Publication:	2025-02-12 16:30:02 UTC

Vector Generalized Linear and Additive Models and Other Associated Models

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! The WSDM is an extension of the HDE (wsdm).

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, t.yee@auckland.ac.nz, with contributions from Victor Miranda and several graduate students over the years, especially Xiangjie (Albert) Xue and Chanatda Somchit.

Maintainer: Thomas Yee t.yee@auckland.ac.nz.

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, 568–588.

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, wsdm, 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) 

The A1A2A3 Blood Group System

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

## Not run: 
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)))

## End(Not run)

The AA-Aa-aa Blood Group System

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)

The AB-Ab-aB-ab Blood Group System

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)

The ABO Blood Group System

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)))

Ordinal Regression with Adjacent Categories Probabilities

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)

Add or Drop All Possible Single Terms to/from a Model

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)

Akaike's Information Criterion

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)

Arcsine–Logit Link Mixtures

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)

Altered, Inflated, Truncated and Deflated Values in GAITD Regression

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, 568–588.

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)

Binomial Logistic Regression by Asymmetric Maximum Likelihood Estimation

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)

Exponential Regression by Asymmetric Maximum Likelihood Estimation

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)

Asymmetric Least Squares Quantile Regression

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)

Poisson Regression by Asymmetric Maximum Likelihood Estimation

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)

Analysis of Deviance for Vector Generalized Linear Model Fits

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")

Autoregressive Process with Order-1 Family Function

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)

Computation of the Exact EIM of an Order-1 Autoregressive Process

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!

Arcsine Link Function

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)

Auckland University Undergraduate Counts Data

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.

Details

Each student is cross-classified by their colleges (Science and Engineering have been combined) and the socio-economic status (SES) of their fathers (1 = highest, down to 4 = lowest).

Source

Dr Tony Morrison.

References

Wild, C. J. and Seber, G. A. F. (2000). Chance Encounters: A First Course in Data Analysis and Inference, New York: Wiley.

Examples

auuc
## Not run: 
round(fitted(grc(auuc)))
round(fitted(grc(auuc, Rank = 2)))

## End(Not run)

Auxiliary Function for the Positive Bernoulli Family Function with Time Effects

Description

Returns behavioural effects indicator variables from a capture history matrix.

Usage

aux.posbernoulli.t(y, check.y = FALSE, rename = TRUE, name = "bei")

Arguments

Details

This function can help fit certain capture–recapture models (commonly known as M_{tb} or M_{tbh} (no prefix h means it is an intercept-only model) in the literature). See posbernoulli.t for details.

Value

A list with the following components.

cap.hist1

A matrix the same dimension as y. In any particular row there are 0s up to the first capture. Then there are 1s thereafter.

cap1

A vector specifying which time occasion the animal was first captured.

y0i

Number of noncaptures before the first capture.

yr0i

Number of noncaptures after the first capture.

yr1i

Number of recaptures after the first capture.

See Also

posbernoulli.t, deermice.

Examples

# Fit a M_tbh model to the deermice data:
(pdata <- aux.posbernoulli.t(with(deermice,
                                  cbind(y1, y2, y3, y4, y5, y6))))

deermice <- data.frame(deermice,
                    bei = 0,  # Add this
                    pdata$cap.hist1)  # Incorporate these
head(deermice)  # Augmented with behavioural effect indicator variables
tail(deermice)

Data on Back Pain Prognosis, from Anderson (1984)

Description

Data from a study of patients suffering from back pain. Prognostic variables were recorded at presentation and progress was categorised three weeks after treatment.

Usage

data(backPain)

Format

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

x2

length of previous attack.

x3

pain change.

x4

lordosis.

pain

an ordered factor describing the progress of each patient with levels worse < same < slight.improvement < moderate.improvement < marked.improvement < complete.relief.

Source

http://ideas.repec.org/c/boc/bocode/s419001.html

The data set and this help file was copied from gnm so that a vignette in VGAM could be run; the analysis is described in Yee (2010).

The data frame backPain2 is a modification of backPain where the variables have been renamed (x1 becomes x2, x2 becomes x3, x3 becomes x4) and converted into factors.

References

Anderson, J. A. (1984). Regression and Ordered Categorical Variables. J. R. Statist. Soc. B, 46(1), 1-30.

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

Examples

summary(backPain)
summary(backPain2)

Bacon and Eggs Data

Description

Purchasing of bacon and eggs.

Usage

data(beggs)

Format

Data frame of a two way table.

b0, b1, b2, b3, b4

The b refers to bacon. The number of times bacon was purchased was 0, 1, 2, 3, or 4.

e0, e1, e2, e3, e4

The e refers to eggs. The number of times eggs was purchased was 0, 1, 2, 3, or 4.

Details

The data is from Information Resources, Inc., a consumer panel based in a large US city [see Bell and Lattin (1998) for further details]. Starting in June 1991, the purchases in the bacon and fresh eggs product categories for a sample of 548 households over four consecutive store trips was tracked. Only those grocery shopping trips with a total basket value of at least five dollars was considered. For each household, the total number of bacon purchases in their four eligible shopping trips and the total number of egg purchases (usually a package of eggs) for the same trips, were counted.

Source

Bell, D. R. and Lattin, J. M. (1998) Shopping Behavior and Consumer Preference for Store Price Format: Why ‘Large Basket’ Shoppers Prefer EDLP. Marketing Science, 17, 66–88.

References

Danaher, P. J. and Hardie, B. G. S. (2005). Bacon with Your Eggs? Applications of a New Bivariate Beta-Binomial Distribution. American Statistician, 59(4), 282–286.

See Also

rrvglm, rcim, grc.

Examples

beggs
colSums(beggs)
rowSums(beggs)

The Bell Series of Integers

Description

Returns the values of the Bell series.

Usage

bell(n)

Arguments

Details

The Bell numbers emerge from a series expansion of \exp(e^x - 1) for real x. The first few values are B_{0}=1, B_{1}=1, B_{2}=2, B_{3}=5, B_{4}=15. The series increases quickly so that overflow occurs when its argument is more than 218.

Value

This function returns B_{n}.

Author(s)

T. W. Yee

References

Bell, E. T. (1934). Exponential polynomials. Ann. Math., 35, 258–277.

Bell, E. T. (1934). Exponential numbers. Amer. Math. Monthly, 41, 411–419.

See Also

bellff, rbell.

Examples

 ## Not run: 
plot(0:10, bell(0:10), log = "y", type = "h", col = "blue")

## End(Not run)

Benford's Distribution

Description

Density, distribution function, quantile function, and random generation for Benford's distribution.

Usage

dbenf(x, ndigits = 1, log = FALSE)
pbenf(q, ndigits = 1, lower.tail = TRUE, log.p = FALSE)
qbenf(p, ndigits = 1, lower.tail = TRUE, log.p = FALSE)
rbenf(n, ndigits = 1)

Arguments

Details

Benford's Law (aka the significant-digit law) is the empirical observation that in many naturally occuring tables of numerical data, the leading significant (nonzero) digit is not uniformly distributed in \{1,2,\ldots,9\}. Instead, the leading significant digit (=D, say) obeys the law

P(D=d) = \log_{10} \left( 1 + \frac1d \right)

for d=1,\ldots,9. This means the probability the first significant digit is 1 is approximately 0.301, etc.

Benford's Law was apparently first discovered in 1881 by astronomer/mathematician S. Newcombe. It started by the observation that the pages of a book of logarithms were dirtiest at the beginning and progressively cleaner throughout. In 1938, a General Electric physicist called F. Benford rediscovered the law on this same observation. Over several years he collected data from different sources as different as atomic weights, baseball statistics, numerical data from Reader's Digest, and drainage areas of rivers.

Applications of Benford's Law has been as diverse as to the area of fraud detection in accounting and the design computers.

Benford's distribution has been called “a” logarithmic distribution; see logff.

Value

dbenf gives the density, pbenf gives the distribution function, and qbenf gives the quantile function, and rbenf generates random deviates.

Author(s)

T. W. Yee and Kai Huang

References

Benford, F. (1938). The Law of Anomalous Numbers. Proceedings of the American Philosophical Society, 78, 551–572.

Newcomb, S. (1881). Note on the Frequency of Use of the Different Digits in Natural Numbers. American Journal of Mathematics, 4, 39–40.

Examples

dbenf(x <- c(0:10, NA, NaN, -Inf, Inf))
pbenf(x)

## Not run: 
xx <- 1:9
barplot(dbenf(xx), col = "lightblue", xlab = "Leading digit",
        ylab = "Probability", names.arg = as.character(xx),
        main = "Benford's distribution", las = 1)

hist(rbenf(1000), border = "blue", prob = TRUE,
     main = "1000 random variates from Benford's distribution",
     xlab = "Leading digit", sub="Red is the true probability",
     breaks = 0:9 + 0.5, ylim = c(0, 0.35), xlim = c(0, 10.0))
lines(xx, dbenf(xx), col = "red", type = "h")
points(xx, dbenf(xx), col = "red")

## End(Not run)

The Benini Distribution

Description

Density, distribution function, quantile function and random generation for the Benini distribution with parameter shape.

Usage

dbenini(x, y0, shape, log = FALSE)
pbenini(q, y0, shape, lower.tail = TRUE, log.p = FALSE)
qbenini(p, y0, shape, lower.tail = TRUE, log.p = FALSE)
rbenini(n, y0, shape)

Arguments

Details

See benini1, the VGAM family function for estimating the parameter s by maximum likelihood estimation, for the formula of the probability density function and other details.

Value

dbenini gives the density, pbenini gives the distribution function, qbenini gives the quantile function, and rbenini generates random deviates.

Author(s)

T. W. Yee and Kai Huang

References

Kleiber, C. and Kotz, S. (2003). Statistical Size Distributions in Economics and Actuarial Sciences, Hoboken, NJ, USA: Wiley-Interscience.

See Also

benini1.

Examples

## Not run: 
y0 <- 1; shape <- exp(1)
xx <- seq(0.0, 4, len = 101)
plot(xx, dbenini(xx, y0 = y0, shape = shape), col = "blue",
     main = "Blue is density, orange is the CDF", type = "l",
     sub = "Purple lines are the 10,20,...,90 percentiles",
     ylim = 0:1, las = 1, ylab = "", xlab = "x")
abline(h = 0, col = "blue", lty = 2)
lines(xx, pbenini(xx, y0 = y0, shape = shape), col = "orange")
probs <- seq(0.1, 0.9, by = 0.1)
Q <- qbenini(probs, y0 = y0, shape = shape)
lines(Q, dbenini(Q, y0 = y0, shape = shape),
      col = "purple", lty = 3, type = "h")
pbenini(Q, y0 = y0, shape = shape) - probs  # Should be all zero

## End(Not run)

Benini Distribution Family Function

Description

Estimating the 1-parameter Benini distribution by maximum likelihood estimation.

Usage

benini1(y0 = stop("argument 'y0' must be specified"),
        lshape = "loglink", ishape = NULL, imethod = 1,
        zero = NULL, parallel = FALSE,
        type.fitted = c("percentiles", "Qlink"),
        percentiles = 50)

Arguments

Details

The Benini distribution has a probability density function that can be written

f(y) = 2 s \exp(-s[(\log(y/y_0))^2]) \log(y/y_0) / y

for 0 < y_0 < y, and shape s > 0. The cumulative distribution function for Y is

F(y) = 1 - \exp(-s[(\log(y/y_0))^2]).

Here, Newton-Raphson and Fisher scoring coincide. The median of Y is now returned as the fitted values, by default. This VGAM family function can handle a multiple responses, which is inputted as a matrix.

On fitting, the extra slot has a component called y0 which contains the value of the y0 argument.

Value

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

Note

Yet to do: the 2-parameter Benini distribution estimates another shape parameter a too. Hence, the code may change in the future.

Author(s)

T. W. Yee

References

Kleiber, C. and Kotz, S. (2003). Statistical Size Distributions in Economics and Actuarial Sciences, Hoboken, NJ, USA: Wiley-Interscience.

See Also

Benini.

Examples

y0 <- 1; nn <- 3000
bdata <- data.frame(y  = rbenini(nn, y0 = y0, shape = exp(2)))
fit <- vglm(y ~ 1, benini1(y0 = y0), data = bdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
fit@extra$y0
c(head(fitted(fit), 1), with(bdata, median(y)))  # Should be equal

The Beta-Binomial Distribution

Description

Density, distribution function, and random generation for the beta-binomial distribution and the inflated beta-binomial distribution.

Usage

dbetabinom(x, size, prob, rho = 0, log = FALSE)
pbetabinom(q, size, prob, rho = 0, log.p = FALSE)
rbetabinom(n, size, prob, rho = 0)
dbetabinom.ab(x, size, shape1, shape2, log = FALSE,
              Inf.shape = exp(20), limit.prob = 0.5)
pbetabinom.ab(q, size, shape1, shape2, limit.prob = 0.5,
              log.p = FALSE)
rbetabinom.ab(n, size, shape1, shape2, limit.prob = 0.5,
              .dontuse.prob = NULL)
dzoibetabinom(x, size, prob, rho = 0, pstr0 = 0, pstrsize = 0,
              log = FALSE)
pzoibetabinom(q, size, prob, rho, pstr0 = 0, pstrsize = 0,
              lower.tail = TRUE, log.p = FALSE)
rzoibetabinom(n, size, prob, rho = 0, pstr0 = 0, pstrsize = 0)
dzoibetabinom.ab(x, size, shape1, shape2, pstr0 = 0, pstrsize = 0,
                 log = FALSE)
pzoibetabinom.ab(q, size, shape1, shape2, pstr0 = 0, pstrsize = 0,
              lower.tail = TRUE, log.p = FALSE)
rzoibetabinom.ab(n, size, shape1, shape2, pstr0 = 0, pstrsize = 0)

Arguments

Details

The beta-binomial distribution is a binomial distribution whose probability of success is not a constant but it is generated from a beta distribution with parameters shape1 and shape2. Note that the mean of this beta distribution is mu = shape1/(shape1+shape2), which therefore is the mean or the probability of success.

See betabinomial and betabinomialff, the VGAM family functions for estimating the parameters, for the formula of the probability density function and other details.

For the inflated beta-binomial distribution, the probability mass function is

P(Y = y) = (1 - pstr0 - pstrsize) \times BB(y) + pstr0 \times I[y = 0] + pstrsize \times I[y = size]

where BB(y) is the probability mass function of the beta-binomial distribution with the same shape parameters (pbetabinom.ab), pstr0 is the inflated probability at 0 and pstrsize is the inflated probability at 1. The default values of pstr0 and pstrsize mean that these functions behave like the ordinary Betabinom when only the essential arguments are inputted.

Value

dbetabinom and dbetabinom.ab give the density, pbetabinom and pbetabinom.ab give the distribution function, and rbetabinom and rbetabinom.ab generate random deviates.

dzoibetabinom and dzoibetabinom.ab give the inflated density, pzoibetabinom and pzoibetabinom.ab give the inflated distribution function, and rzoibetabinom and rzoibetabinom.ab generate random inflated deviates.

Warning

Setting rho = 1 is not recommended, however the code may be modified in the future to handle this special case.

Note

pzoibetabinom, pzoibetabinom.ab, pbetabinom and pbetabinom.ab can be particularly slow. The functions here ending in .ab are called from those functions which don't. The simple transformations \mu=\alpha / (\alpha + \beta) and \rho=1/(1 + \alpha + \beta) are used, where \alpha and \beta are the two shape parameters.

Author(s)

T. W. Yee and Xiangjie Xue

See Also

Extbetabinom, betabinomial, betabinomialff, Zoabeta, Beta.

Examples

set.seed(1); rbetabinom(10, 100, prob = 0.5)
set.seed(1);     rbinom(10, 100, prob = 0.5)  # The same as rho = 0

## Not run:  N <- 9; xx <- 0:N; s1 <- 2; s2 <- 3
dy <- dbetabinom.ab(xx, size = N, shape1 = s1, shape2 = s2)
barplot(rbind(dy, dbinom(xx, size = N, prob = s1 / (s1+s2))),
        beside = TRUE, col = c("blue","green"), las = 1,
        main = paste("Beta-binomial (size=",N,", shape1=", s1,
                   ", shape2=", s2, ") (blue) vs\n",
        " Binomial(size=", N, ", prob=", s1/(s1+s2), ") (green)",
                     sep = ""),
        names.arg = as.character(xx), cex.main = 0.8)
sum(dy * xx)  # Check expected values are equal
sum(dbinom(xx, size = N, prob = s1 / (s1+s2)) * xx)
# Should be all 0:
cumsum(dy) - pbetabinom.ab(xx, N, shape1 = s1, shape2 = s2)

y <- rbetabinom.ab(n = 1e4, size = N, shape1 = s1, shape2 = s2)
ty <- table(y)
barplot(rbind(dy, ty / sum(ty)),
        beside = TRUE, col = c("blue", "orange"), las = 1,
        main = paste("Beta-binomial (size=", N, ", shape1=", s1,
                     ", shape2=", s2, ") (blue) vs\n",
        " Random generated beta-binomial(size=", N, ", prob=",
        s1/(s1+s2), ") (orange)", sep = ""), cex.main = 0.8,
        names.arg = as.character(xx))

N <- 1e5; size <- 20; pstr0 <- 0.2; pstrsize <- 0.2
kk <- rzoibetabinom.ab(N, size, s1, s2, pstr0, pstrsize)
hist(kk, probability = TRUE, border = "blue", ylim = c(0, 0.25),
     main = "Blue/green = inflated; orange = ordinary beta-binomial",
     breaks = -0.5 : (size + 0.5))
sum(kk == 0) / N  # Proportion of 0
sum(kk == size) / N  # Proportion of size
lines(0 : size,
      dbetabinom.ab(0 : size, size, s1, s2), col = "orange")
lines(0 : size, col = "green", type = "b",
      dzoibetabinom.ab(0 : size, size, s1, s2, pstr0, pstrsize))

## End(Not run)

Beta-binomial Distribution Family Function

Description

Fits a beta-binomial distribution by maximum likelihood estimation. The two parameters here are the mean and correlation coefficient.

Usage

betabinomial(lmu = "logitlink", lrho = "logitlink",
   irho = NULL, imethod = 1,
   ishrinkage = 0.95, nsimEIM = NULL, zero = "rho")

Arguments

Details

There are several parameterizations of the beta-binomial distribution. This family function directly models the mean and correlation parameter, i.e., the probability of success. The model can be written T|P=p \sim Binomial(N,p) where P has a beta distribution with shape parameters \alpha and \beta. Here, N is the number of trials (e.g., litter size), T=NY is the number of successes, and p is the probability of a success (e.g., a malformation). That is, Y is the proportion of successes. Like binomialff, the fitted values are the estimated probability of success (i.e., E[Y] and not E[T]) and the prior weights N are attached separately on the object in a slot.

The probability function is

P(T=t) = {N \choose t} \frac{Be(\alpha+t, \beta+N-t)} {Be(\alpha, \beta)}

where t=0,1,\ldots,N, and Be is the beta function with shape parameters \alpha and \beta. Recall Y = T/N is the real response being modelled.

The default model is \eta_1 = logit(\mu) and \eta_2 = logit(\rho) because both parameters lie between 0 and 1. The mean (of Y) is p=\mu=\alpha/(\alpha+\beta) and the variance (of Y) is \mu(1-\mu)(1+(N-1)\rho)/N. Here, the correlation \rho is given by 1/(1 + \alpha + \beta) and is the correlation between the N individuals within a litter. A litter effect is typically reflected by a positive value of \rho. It is known as the over-dispersion parameter.

This family function uses Fisher scoring. Elements of the second-order expected derivatives with respect to \alpha and \beta are computed numerically, which may fail for large \alpha, \beta, N or else take a long time.

Value

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

Suppose fit is a fitted beta-binomial model. Then depvar(fit) are the sample proportions y, fitted(fit) returns estimates of E(Y), and weights(fit, type = "prior") returns the number of trials N.

Warning

If the estimated rho parameter is close to 0 then a good solution is to use extbetabinomial. Or you could try lrho = "rhobitlink".

This family function is prone to numerical difficulties due to the expected information matrices not being positive-definite or ill-conditioned over some regions of the parameter space. If problems occur try setting irho to some numerical value, nsimEIM = 100, say, or else use etastart argument of vglm, etc.

Note

This function processes the input in the same way as binomialff. But it does not handle the case N=1 very well because there are two parameters to estimate, not one, for each row of the input. Cases where N=1 can be omitted via the subset argument of vglm.

The extended beta-binomial distribution of Prentice (1986) implemented by extbetabinomial is the preferred VGAM family function for BBD regression.

Author(s)

T. W. Yee

References

Moore, D. F. and Tsiatis, A. (1991). Robust estimation of the variance in moment methods for extra-binomial and extra-Poisson variation. Biometrics, 47, 383–401.

See Also

extbetabinomial, betabinomialff, betabinomial.rho, Betabinom, binomialff, betaff, dirmultinomial, log1plink, cloglink, lirat, simulate.vlm.

Examples

# Example 1
bdata <- data.frame(N = 10, mu = 0.5, rho = 0.8)
bdata <- transform(bdata,
            y = rbetabinom(100, size = N, prob = mu, rho = rho))
fit <- vglm(cbind(y, N-y) ~ 1, betabinomial, bdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
head(cbind(depvar(fit), weights(fit, type = "prior")))


# Example 2
fit <- vglm(cbind(R, N-R) ~ 1, betabinomial, lirat,
            trace = TRUE, subset = N > 1)
coef(fit, matrix = TRUE)
Coef(fit)
t(fitted(fit))
t(depvar(fit))
t(weights(fit, type = "prior"))


# Example 3, which is more complicated
lirat <- transform(lirat, fgrp = factor(grp))
summary(lirat)  # Only 5 litters in group 3
fit2 <- vglm(cbind(R, N-R) ~ fgrp + hb, betabinomial(zero = 2),
             data = lirat, trace = TRUE, subset = N > 1)
coef(fit2, matrix = TRUE)
## Not run:  with(lirat, plot(hb[N > 1], fit2@misc$rho,
         xlab = "Hemoglobin", ylab = "Estimated rho",
         pch = as.character(grp[N > 1]), col = grp[N > 1])) 
## End(Not run)
## Not run:   # cf. Figure 3 of Moore and Tsiatis (1991)
with(lirat, plot(hb, R / N, pch = as.character(grp), col = grp,
         xlab = "Hemoglobin level", ylab = "Proportion Dead",
         main = "Fitted values (lines)", las = 1))
smalldf <- with(lirat, lirat[N > 1, ])
for (gp in 1:4) {
  xx <- with(smalldf, hb[grp == gp])
  yy <- with(smalldf, fitted(fit2)[grp == gp])
  ooo <- order(xx)
  lines(xx[ooo], yy[ooo], col = gp, lwd = 2)
} 
## End(Not run)

Beta-binomial Distribution Family Function (with known rho)

Description

Fits a beta-binomial distribution by maximum likelihood estimation, where the correlation coefficient rho is inputted. The parameter estimated is the mean.

Usage

betabinomial.rho(lmu = "logitlink", imethod = 1, ishrinkage = 0.95)

Arguments

Details

This family function conducts a logistic-like regression where the correlation parameter \rho of a betabinomial distribution is inputted by the user. The family function is somewhat like a simplified betabinomial. The argument form2 (see vglm) is used to input the \rho values, which must lie in [0, 1].

The default model has \eta_1 = logit(\mu).

Value

Same as betabinomial.

Author(s)

T. W. Yee

See Also

betabinomial, extbetabinomial, betabinomialff, Betabinom, vglm, binomialff, betaff.

Examples

## Not run: 
# Example 1
nn <- 10000; NN <- 100
bdata <- data.frame(N = NN, x2 = rnorm(nn),
                    x3 = rnorm(nn))
bdata <-
    transform(bdata,
              mu1  = logitlink(-0.5, inverse = TRUE),
              rho1 = logitlink( 0.5, inverse = TRUE),
              mu2  = logitlink(-0.5 + x2, inverse = TRUE),
              rho2 = logitlink(-0.5 + x3, inverse = TRUE))
bdata <- transform(bdata,
         y1 = rbetabinom(nn, size = N, prob = mu1, rho = rho1),
         y2 = rbetabinom(nn, size = N, prob = mu2, rho = rho2))
fit1 <- vglm(cbind(y1, N - y1) ~ 1, betabinomial.rho,
         form2 = ~ rho1, crit = "c", bdata, trace = TRUE)
coef(fit1, matrix = TRUE)
head(fit1@extra$rho)
max(abs(fitted(fit1) - with(bdata, mu1)))  # Should be 0

# Example 2
fit2 <- vglm(cbind(y2, N - y2) ~ x2, form2 = ~ rho2,
         betabinomial.rho, crit = "c",
         bdata, trace = TRUE)
coef(fit2, matrix = TRUE)
max(abs(fit2@extra$rho - with(bdata, rho2)))  # Should be 0
max(abs(fitted(fit2) - with(bdata, mu2)))  # Should be 0

## End(Not run)

Beta-binomial Distribution Family Function

Description

Fits a beta-binomial distribution by maximum likelihood estimation. The two parameters here are the shape parameters of the underlying beta distribution.

Usage

betabinomialff(lshape1 = "loglink", lshape2 = "loglink",
   ishape1 = 1, ishape2 = NULL, imethod = 1, ishrinkage = 0.95,
   nsimEIM = NULL, zero = NULL)

Arguments

Details

There are several parameterizations of the beta-binomial distribution. This family function directly models the two shape parameters of the associated beta distribution rather than the probability of success (however, see Note below). The model can be written T|P=p \sim Binomial(N,p) where P has a beta distribution with shape parameters \alpha and \beta. Here, N is the number of trials (e.g., litter size), T=NY is the number of successes, and p is the probability of a success (e.g., a malformation). That is, Y is the proportion of successes. Like binomialff, the fitted values are the estimated probability of success (i.e., E[Y] and not E[T]) and the prior weights N are attached separately on the object in a slot.

The probability function is

P(T=t) = {N \choose t} \frac{B(\alpha+t, \beta+N-t)} {B(\alpha, \beta)}

where t=0,1,\ldots,N, and B is the beta function with shape parameters \alpha and \beta. Recall Y = T/N is the real response being modelled.

The default model is \eta_1 = \log(\alpha) and \eta_2 = \log(\beta) because both parameters are positive. The mean (of Y) is p=\mu=\alpha/(\alpha+\beta) and the variance (of Y) is \mu(1-\mu)(1+(N-1)\rho)/N. Here, the correlation \rho is given by 1/(1 + \alpha + \beta) and is the correlation between the N individuals within a litter. A litter effect is typically reflected by a positive value of \rho. It is known as the over-dispersion parameter.

This family function uses Fisher scoring. The two diagonal elements of the second-order expected derivatives with respect to \alpha and \beta are computed numerically, which may fail for large \alpha, \beta, N or else take a long time.

Value

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

Suppose fit is a fitted beta-binomial model. Then fit@y (better: depvar(fit)) contains the sample proportions y, fitted(fit) returns estimates of E(Y), and weights(fit, type = "prior") returns the number of trials N.

Warning

This family function is prone to numerical difficulties due to the expected information matrices not being positive-definite or ill-conditioned over some regions of the parameter space. If problems occur try setting ishape1 to be some other positive value, using ishape2 and/or setting zero = 2.

This family function may be renamed in the future. See the warnings in betabinomial.

Note

This function processes the input in the same way as binomialff. But it does not handle the case N=1 very well because there are two parameters to estimate, not one, for each row of the input. Cases where N=1 can be omitted via the subset argument of vglm.

Although the two linear/additive predictors given above are in terms of \alpha and \beta, basic algebra shows that the default amounts to fitting a logit link to the probability of success; subtracting the second linear/additive predictor from the first gives that logistic regression linear/additive predictor. That is, logit(p) = \eta_1 - \eta_2. This is illustated in one of the examples below.

The extended beta-binomial distribution of Prentice (1986) implemented by extbetabinomial is the preferred VGAM family function for BBD regression.

Author(s)

T. W. Yee

References

Moore, D. F. and Tsiatis, A. (1991). Robust estimation of the variance in moment methods for extra-binomial and extra-Poisson variation. Biometrics, 47, 383–401.

Prentice, R. L. (1986). Binary regression using an extended beta-binomial distribution, with discussion of correlation induced by covariate measurement errors. Journal of the American Statistical Association, 81, 321–327.

See Also

extbetabinomial, betabinomial, Betabinom, binomialff, betaff, dirmultinomial, lirat, simulate.vlm.

Examples

# Example 1
N <- 10; s1 <- exp(1); s2 <- exp(2)
y <- rbetabinom.ab(n = 100, size = N, shape1 = s1, shape2 = s2)
fit <- vglm(cbind(y, N-y) ~ 1, betabinomialff, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
head(fit@misc$rho)  # The correlation parameter
head(cbind(depvar(fit), weights(fit, type = "prior")))


# Example 2
fit <- vglm(cbind(R, N-R) ~ 1, betabinomialff, data = lirat,
            trace = TRUE, subset = N > 1)
coef(fit, matrix = TRUE)
Coef(fit)
fit@misc$rho  # The correlation parameter
t(fitted(fit))
t(depvar(fit))
t(weights(fit, type = "prior"))
# A "loglink" link for the 2 shape params is a logistic regression:
all.equal(c(fitted(fit)),
          as.vector(logitlink(predict(fit)[, 1] -
                          predict(fit)[, 2], inverse = TRUE)))


# Example 3, which is more complicated
lirat <- transform(lirat, fgrp = factor(grp))
summary(lirat)  # Only 5 litters in group 3
fit2 <- vglm(cbind(R, N-R) ~ fgrp + hb, betabinomialff(zero = 2),
           data = lirat, trace = TRUE, subset = N > 1)
coef(fit2, matrix = TRUE)
coef(fit2, matrix = TRUE)[, 1] -
coef(fit2, matrix = TRUE)[, 2]  # logitlink(p)
## Not run:  with(lirat, plot(hb[N > 1], fit2@misc$rho,
   xlab = "Hemoglobin", ylab = "Estimated rho",
   pch = as.character(grp[N > 1]), col = grp[N > 1])) 
## End(Not run)
## Not run:   # cf. Figure 3 of Moore and Tsiatis (1991)
with(lirat, plot(hb, R / N, pch = as.character(grp), col = grp,
   xlab = "Hemoglobin level", ylab = "Proportion Dead", las = 1,
   main = "Fitted values (lines)"))

smalldf <- with(lirat, lirat[N > 1, ])
for (gp in 1:4) {
  xx <- with(smalldf, hb[grp == gp])
  yy <- with(smalldf, fitted(fit2)[grp == gp])
  ooo <- order(xx)
  lines(xx[ooo], yy[ooo], col = gp, lwd = 2)
} 
## End(Not run)

The Two-parameter Beta Distribution Family Function

Description

Estimation of the mean and precision parameters of the beta distribution.

Usage

betaff(A = 0, B = 1, lmu = "logitlink", lphi = "loglink",
       imu = NULL, iphi = NULL,
       gprobs.y = ppoints(8), gphi  = exp(-3:5)/4, zero = NULL)

Arguments

Details

The two-parameter beta distribution can be written f(y) =

(y-A)^{\mu_1 \phi-1} \times (B-y)^{(1-\mu_1) \phi-1} / [beta(\mu_1 \phi,(1-\mu_1) \phi) \times (B-A)^{\phi-1}]

for A < y < B, and beta(.,.) is the beta function (see beta). The parameter \mu_1 satisfies \mu_1 = (\mu - A) / (B-A) where \mu is the mean of Y. That is, \mu_1 is the mean of of a standard beta distribution: E(Y) = A + (B-A) \times \mu_1, and these are the fitted values of the object. Also, \phi is positive and A < \mu < B. Here, the limits A and B are known.

Another parameterization of the beta distribution involving the raw shape parameters is implemented in betaR.

For general A and B, the variance of Y is (B-A)^2 \times \mu_1 \times (1-\mu_1) / (1+\phi). Then \phi can be interpreted as a precision parameter in the sense that, for fixed \mu, the larger the value of \phi, the smaller the variance of Y. Also, \mu_1 = shape1/(shape1+shape2) and \phi = shape1+shape2. Fisher scoring is implemented.

Value

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

Note

The response must have values in the interval (A, B). The user currently needs to manually choose lmu to match the input of arguments A and B, e.g., with extlogitlink; see the example below.

Author(s)

Thomas W. Yee

References

Ferrari, S. L. P. and Francisco C.-N. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31, 799–815.

See Also

betaR, Beta, dzoabeta, genbetaII, betaII, betabinomialff, betageometric, betaprime, rbetageom, rbetanorm, kumar, extlogitlink, simulate.vlm.

Examples

bdata <- data.frame(y = rbeta(nn <- 1000, shape1 = exp(0),
                              shape2 = exp(1)))
fit1 <- vglm(y ~ 1, betaff, data = bdata, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1)  # Useful for intercept-only models

# General A and B, and with a covariate
bdata <- transform(bdata, x2 = runif(nn))
bdata <- transform(bdata, mu = logitlink(0.5 - x2, inverse = TRUE),
                          prec = exp(3.0 + x2))  # prec == phi
bdata <- transform(bdata, shape2 = prec * (1 - mu),
                          shape1 = mu * prec)
bdata <- transform(bdata,
                   y = rbeta(nn, shape1 = shape1, shape2 = shape2))
bdata <- transform(bdata, Y = 5 + 8 * y)  # From 5--13, not 0--1
fit <- vglm(Y ~ x2, data = bdata, trace = TRUE,
   betaff(A = 5, B = 13, lmu = "extlogitlink(min = 5, max = 13)"))
coef(fit, matrix = TRUE)

The Beta-Geometric Distribution

Description

Density, distribution function, and random generation for the beta-geometric distribution.

Usage

dbetageom(x, shape1, shape2, log = FALSE)
pbetageom(q, shape1, shape2, log.p = FALSE)
rbetageom(n, shape1, shape2)

Arguments

Details

The beta-geometric distribution is a geometric distribution whose probability of success is not a constant but it is generated from a beta distribution with parameters shape1 and shape2. Note that the mean of this beta distribution is shape1/(shape1+shape2), which therefore is the mean of the probability of success.

Value

dbetageom gives the density, pbetageom gives the distribution function, and rbetageom generates random deviates.

Note

pbetageom can be particularly slow.

Author(s)

T. W. Yee

See Also

geometric, betaff, Beta.

Examples

## Not run: 
shape1 <- 1; shape2 <- 2; y <- 0:30
proby <- dbetageom(y, shape1, shape2, log = FALSE)
plot(y, proby, type = "h", col = "blue", ylab = "P[Y=y]", main = paste0(
     "Y ~ Beta-geometric(shape1=", shape1,", shape2=", shape2, ")"))
sum(proby)

## End(Not run)

Beta-geometric Distribution Family Function

Description

Maximum likelihood estimation for the beta-geometric distribution.

Usage

betageometric(lprob = "logitlink", lshape = "loglink",
    iprob = NULL,    ishape = 0.1,
    moreSummation = c(2, 100), tolerance = 1.0e-10, zero = NULL)

Arguments

Details

A random variable Y has a 2-parameter beta-geometric distribution if P(Y=y) = p (1-p)^y for y=0,1,2,\ldots where p are generated from a standard beta distribution with shape parameters shape1 and shape2. The parameterization here is to focus on the parameters p and \phi = 1/(shape1+shape2), where \phi is shape. The default link functions for these ensure that the appropriate range of the parameters is maintained. The mean of Y is E(Y) = shape2 / (shape1-1) = (1-p) / (p-\phi) if shape1 > 1, and if so, then this is returned as the fitted values.

The geometric distribution is a special case of the beta-geometric distribution with \phi=0 (see geometric). However, fitting data from a geometric distribution may result in numerical problems because the estimate of \log(\phi) will 'converge' to -Inf.

Value

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

Note

The first iteration may be very slow; if practical, it is best for the weights argument of vglm etc. to be used rather than inputting a very long vector as the response, i.e., vglm(y ~ 1, ..., weights = wts) is to be preferred over vglm(rep(y, wts) ~ 1, ...). If convergence problems occur try inputting some values of argument ishape.

If an intercept-only model is fitted then the misc slot of the fitted object has list components shape1 and shape2.

Author(s)

T. W. Yee

References

Paul, S. R. (2005). Testing goodness of fit of the geometric distribution: an application to human fecundability data. Journal of Modern Applied Statistical Methods, 4, 425–433.

See Also

geometric, betaff, rbetageom.

Examples

## Not run: 
bdata <- data.frame(y = 0:11,
                    wts = c(227,123,72,42,21,31,11,14,6,4,7,28))
fitb <- vglm(y ~ 1, betageometric, bdata, weight = wts, trace = TRUE)
fitg <- vglm(y ~ 1,     geometric, bdata, weight = wts, trace = TRUE)
coef(fitb, matrix = TRUE)
Coef(fitb)
sqrt(diag(vcov(fitb, untransform = TRUE)))
fitb@misc$shape1
fitb@misc$shape2
# Very strong evidence of a beta-geometric:
pchisq(2 * (logLik(fitb) - logLik(fitg)), df = 1, lower.tail = FALSE)

## End(Not run)

Beta Distribution of the Second Kind

Description

Maximum likelihood estimation of the 3-parameter beta II distribution.

Usage

betaII(lscale = "loglink", lshape2.p = "loglink",
       lshape3.q = "loglink", iscale = NULL, ishape2.p = NULL,
       ishape3.q = NULL, imethod = 1,
       gscale = exp(-5:5), gshape2.p = exp(-5:5),
       gshape3.q = seq(0.75, 4, by = 0.25),
       probs.y = c(0.25, 0.5, 0.75), zero = "shape")

Arguments

Details

The 3-parameter beta II is the 4-parameter generalized beta II distribution with shape parameter a=1. It is also known as the Pearson VI distribution. Other distributions which are special cases of the 3-parameter beta II include the Lomax (p=1) and inverse Lomax (q=1). More details can be found in Kleiber and Kotz (2003).

The beta II distribution has density

f(y) = y^{p-1} / [b^p B(p,q) \{1 + y/b\}^{p+q}]

for b > 0, p > 0, q > 0, y \geq 0. Here, b is the scale parameter scale, and the others are shape parameters. The mean is

E(Y) = b \, \Gamma(p + 1) \, \Gamma(q - 1) / (\Gamma(p) \, \Gamma(q))

provided q > 1; these are returned as the fitted values. This family function handles multiple responses.

Value

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

Note

See the notes in genbetaII.

Author(s)

T. W. Yee

References

Kleiber, C. and Kotz, S. (2003). Statistical Size Distributions in Economics and Actuarial Sciences, Hoboken, NJ, USA: Wiley-Interscience.

See Also

betaff, genbetaII, dagum, sinmad, fisk, inv.lomax, lomax, paralogistic, inv.paralogistic.

Examples

bdata <- data.frame(y = rsinmad(2000, shape1.a = 1,
         shape3.q = exp(2), scale = exp(1)))  # Not genuine data!
# fit <- vglm(y ~ 1, betaII, data = bdata, trace = TRUE)
fit <- vglm(y ~ 1, betaII(ishape2.p = 0.7, ishape3.q = 0.7),
            data = bdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
summary(fit)

The Beta-Normal Distribution

Description

Density, distribution function, quantile function and random generation for the univariate beta-normal distribution.

Usage

dbetanorm(x, shape1, shape2, mean = 0, sd = 1, log = FALSE)
pbetanorm(q, shape1, shape2, mean = 0, sd = 1,
          lower.tail = TRUE, log.p = FALSE)
qbetanorm(p, shape1, shape2, mean = 0, sd = 1,
          lower.tail = TRUE, log.p = FALSE)
rbetanorm(n, shape1, shape2, mean = 0, sd = 1)

Arguments

Details

The function betauninormal, the VGAM family function for estimating the parameters, has not yet been written.

Value

dbetanorm gives the density, pbetanorm gives the distribution function, qbetanorm gives the quantile function, and rbetanorm generates random deviates.

Author(s)

T. W. Yee

References

Gupta, A. K. and Nadarajah, S. (2004). Handbook of Beta Distribution and Its Applications, pp.146–152. New York: Marcel Dekker.

Examples

## Not run: 
shape1 <- 0.1; shape2 <- 4; m <- 1
x <- seq(-10, 2, len = 501)
plot(x, dbetanorm(x, shape1, shape2, m = m), type = "l",
     ylim = 0:1, las = 1,
     ylab = paste0("betanorm(",shape1,", ",shape2,", m=",m, ", sd=1)"),
     main = "Blue is density, orange is the CDF",
     sub = "Gray lines are the 10,20,...,90 percentiles", col = "blue")
lines(x, pbetanorm(x, shape1, shape2, m = m), col = "orange")
abline(h = 0, col = "black")
probs <- seq(0.1, 0.9, by = 0.1)
Q <- qbetanorm(probs, shape1, shape2, m = m)
lines(Q, dbetanorm(Q, shape1, shape2, m = m),
      col = "gray50", lty = 2, type = "h")
lines(Q, pbetanorm(Q, shape1, shape2, m = m),
      col = "gray50", lty = 2, type = "h")
abline(h = probs, col = "gray50", lty = 2)
pbetanorm(Q, shape1, shape2, m = m) - probs  # Should be all 0

## End(Not run)

The Beta-Prime Distribution

Description

Estimation of the two shape parameters of the beta-prime distribution by maximum likelihood estimation.

Usage

betaprime(lshape = "loglink", ishape1 = 2, ishape2 = NULL,
          zero = NULL)

Arguments

Details

The beta-prime distribution is given by

f(y) = y^{shape1-1} (1+y)^{-shape1-shape2} / B(shape1,shape2)

for y > 0. The shape parameters are positive, and here, B is the beta function. The mean of Y is shape1 / (shape2-1) provided shape2>1; these are returned as the fitted values.

If Y has a Beta(shape1,shape2) distribution then Y/(1-Y) and (1-Y)/Y have a Betaprime(shape1,shape2) and Betaprime(shape2,shape1) distribution respectively. Also, if Y_1 has a gamma(shape1) distribution and Y_2 has a gamma(shape2) distribution then Y_1/Y_2 has a Betaprime(shape1,shape2) distribution.

Value

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

Note

The response must have positive values only.

The beta-prime distribution is also known as the beta distribution of the second kind or the inverted beta distribution.

Author(s)

Thomas W. Yee

References

Johnson, N. L. and Kotz, S. and Balakrishnan, N. (1995). Chapter 25 of: Continuous Univariate Distributions, 2nd edition, Volume 2, New York: Wiley.

See Also

betaff, Beta.

Examples

nn <- 1000
bdata <- data.frame(shape1 = exp(1), shape2 = exp(3))
bdata <- transform(bdata, yb = rbeta(nn, shape1, shape2))
bdata <- transform(bdata, y1 = (1-yb) /    yb,
                          y2 =    yb  / (1-yb),
                          y3 = rgamma(nn, exp(3)) / rgamma(nn, exp(2)))

fit1 <- vglm(y1 ~ 1, betaprime, data = bdata, trace = TRUE)
coef(fit1, matrix = TRUE)

fit2 <- vglm(y2 ~ 1, betaprime, data = bdata, trace = TRUE)
coef(fit2, matrix = TRUE)

fit3 <- vglm(y3 ~ 1, betaprime, data = bdata, trace = TRUE)
coef(fit3, matrix = TRUE)

# Compare the fitted values
with(bdata, mean(y3))
head(fitted(fit3))
Coef(fit3)  # Useful for intercept-only models

The Two-parameter Beta Distribution Family Function

Description

Estimation of the shape parameters of the two-parameter beta distribution.

Usage

betaR(lshape1 = "loglink", lshape2 = "loglink",
      i1 = NULL, i2 = NULL, trim = 0.05,
      A = 0, B = 1, parallel = FALSE, zero = NULL)

Arguments

Details

The two-parameter beta distribution is given by f(y) =

(y-A)^{shape1-1} \times (B-y)^{shape2-1} / [Beta(shape1,shape2) \times (B-A)^{shape1+shape2-1}]

for A < y < B, and Beta(.,.) is the beta function (see beta). The shape parameters are positive, and here, the limits A and B are known. The mean of Y is E(Y) = A + (B-A) \times shape1 / (shape1 + shape2), and these are the fitted values of the object.

For the standard beta distribution the variance of Y is shape1 \times shape2 / [(1+shape1+shape2) \times (shape1+shape2)^2]. If \sigma^2= 1 / (1+shape1+shape2) then the variance of Y can be written \sigma^2 \mu (1-\mu) where \mu=shape1 / (shape1 + shape2) is the mean of Y.

Another parameterization of the beta distribution involving the mean and a precision parameter is implemented in betaff.

Value

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

Note

The response must have values in the interval (A, B). VGAM 0.7-4 and prior called this function betaff.

Author(s)

Thomas W. Yee

References

Johnson, N. L. and Kotz, S. and Balakrishnan, N. (1995). Chapter 25 of: Continuous Univariate Distributions, 2nd edition, Volume 2, New York: Wiley.

Gupta, A. K. and Nadarajah, S. (2004). Handbook of Beta Distribution and Its Applications, New York: Marcel Dekker.

See Also

betaff, Beta, genbetaII, betaII, betabinomialff, betageometric, betaprime, rbetageom, rbetanorm, kumar, simulate.vlm.

Examples

bdata <- data.frame(y = rbeta(1000, shape1 = exp(0), shape2 = exp(1)))
fit <- vglm(y ~ 1, betaR(lshape1 = "identitylink",
            lshape2 = "identitylink"), bdata, trace = TRUE, crit = "coef")
fit <- vglm(y ~ 1, betaR, data = bdata, trace = TRUE, crit = "coef")
coef(fit, matrix = TRUE)
Coef(fit)  # Useful for intercept-only models

bdata <- transform(bdata, Y = 5 + 8 * y)  # From 5 to 13, not 0 to 1
fit <- vglm(Y ~ 1, betaR(A = 5, B = 13), data = bdata, trace = TRUE)
Coef(fit)
c(meanY = with(bdata, mean(Y)), head(fitted(fit),2))

Ali-Mikhail-Haq Distribution Family Function

Description

Estimate the association parameter of Ali-Mikhail-Haq's bivariate distribution by maximum likelihood estimation.

Usage

biamhcop(lapar = "rhobitlink", iapar = NULL, imethod = 1,
         nsimEIM = 250)

Arguments

Details

The cumulative distribution function is

P(Y_1 \leq y_1, Y_2 \leq y_2) = y_1 y_2 / ( 1 - \alpha (1 - y_1) (1 - y_2) )

for -1 < \alpha < 1. The support of the function is the unit square. The marginal distributions are the standard uniform distributions. When \alpha = 0 the random variables are independent. This is an Archimedean copula.

Value

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

Note

The response must be a two-column matrix. Currently, the fitted value is a matrix with two columns and values equal to 0.5. This is because each marginal distribution corresponds to a standard uniform distribution.

Author(s)

T. W. Yee and C. S. Chee

References

Balakrishnan, N. and Lai, C.-D. (2009). Continuous Bivariate Distributions, 2nd ed. New York: Springer.

See Also

rbiamhcop, bifgmcop, bigumbelIexp, rbilogis, simulate.vlm.

Examples

ymat <- rbiamhcop(1000, apar = rhobitlink(2, inverse = TRUE))
fit <- vglm(ymat ~ 1, biamhcop, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)

Ali-Mikhail-Haq Bivariate Distribution

Description

Density, distribution function, and random generation for the (one parameter) bivariate Ali-Mikhail-Haq distribution.

Usage

dbiamhcop(x1, x2, apar, log = FALSE)
pbiamhcop(q1, q2, apar)
rbiamhcop(n, apar)

Arguments

Details

See biamhcop, the VGAM family functions for estimating the parameter by maximum likelihood estimation, for the formula of the cumulative distribution function and other details.

Value

dbiamhcop gives the density, pbiamhcop gives the distribution function, and rbiamhcop generates random deviates (a two-column matrix).

Author(s)

T. W. Yee and C. S. Chee

See Also

biamhcop.

Examples

 x <- seq(0, 1, len = (N <- 101)); apar <- 0.7
ox <- expand.grid(x, x)
zedd <- dbiamhcop(ox[, 1], ox[, 2], apar = apar)
## Not run: 
contour(x, x, matrix(zedd, N, N), col = "blue")
zedd <- pbiamhcop(ox[, 1], ox[, 2], apar = apar)
contour(x, x, matrix(zedd, N, N), col = "blue")

plot(r <- rbiamhcop(n = 1000, apar = apar), col = "blue")
par(mfrow = c(1, 2))
hist(r[, 1])  # Should be uniform
hist(r[, 2])  # Should be uniform

## End(Not run)

Clayton Copula (Bivariate) Family Function

Description

Estimate the correlation parameter of the (bivariate) Clayton copula distribution by maximum likelihood estimation.

Usage

biclaytoncop(lapar = "loglink", iapar = NULL, imethod = 1,
             parallel = FALSE, zero = NULL)

Arguments

Details

The cumulative distribution function is

P(u_1, u_2;\alpha) = (u_1^{-\alpha} + u_2^{-\alpha}-1)^{-1/\alpha}

for 0 \leq \alpha . Here, \alpha is the association parameter. The support of the function is the interior of the unit square; however, values of 0 and/or 1 are not allowed (currently). The marginal distributions are the standard uniform distributions. When \alpha = 0 the random variables are independent.

This VGAM family function can handle multiple responses, for example, a six-column matrix where the first 2 columns is the first out of three responses, the next 2 columns being the next response, etc.

Value

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

Note

The response matrix must have a multiple of two-columns. Currently, the fitted value is a matrix with the same number of columns and values equal to 0.5. This is because each marginal distribution corresponds to a standard uniform distribution.

This VGAM family function is fragile; each response must be in the interior of the unit square.

Author(s)

R. Feyter and T. W. Yee

References

Clayton, D. (1982). A model for association in bivariate survival data. Journal of the Royal Statistical Society, Series B, Methodological, 44, 414–422.

Schepsmeier, U. and Stober, J. (2014). Derivatives and Fisher information of bivariate copulas. Statistical Papers 55, 525–542.

See Also

rbiclaytoncop, dbiclaytoncop, kendall.tau.

Examples

ymat <- rbiclaytoncop(n = (nn <- 1000), apar = exp(2))
bdata <- data.frame(y1 = ymat[, 1], y2 = ymat[, 2],
                    y3 = ymat[, 1], y4 = ymat[, 2], x2 = runif(nn))
summary(bdata)
## Not run:  plot(ymat, col = "blue") 
fit1 <-
  vglm(cbind(y1, y2, y3, y4) ~ 1,  # 2 responses, e.g., (y1,y2) is the 1st
       biclaytoncop, data = bdata,
       trace = TRUE, crit = "coef")  # Sometimes a good idea
coef(fit1, matrix = TRUE)
Coef(fit1)
head(fitted(fit1))
summary(fit1)

# Another example; apar is a function of x2
bdata <- transform(bdata, apar = exp(-0.5 + x2))
ymat <- rbiclaytoncop(n = nn, apar = with(bdata, apar))
bdata <- transform(bdata, y5 = ymat[, 1], y6 = ymat[, 2])
fit2 <- vgam(cbind(y5, y6) ~ s(x2), data = bdata,
             biclaytoncop(lapar = "loglink"), trace = TRUE)
## Not run: plot(fit2, lcol = "blue", scol = "orange", se = TRUE) 

Clayton Copula (Bivariate) Distribution

Description

Density and random generation for the (one parameter) bivariate Clayton copula distribution.

Usage

dbiclaytoncop(x1, x2, apar = 0, log = FALSE)
rbiclaytoncop(n, apar = 0)

Arguments

Details

See biclaytoncop, the VGAM family functions for estimating the parameter by maximum likelihood estimation, for the formula of the cumulative distribution function and other details.

Value

dbiclaytoncop gives the density at point (x1,x2), rbiclaytoncop generates random deviates (a two-column matrix).

Note

dbiclaytoncop() does not yet handle x1 = 0 and/or x2 = 0.

Author(s)

R. Feyter and T. W. Yee

References

Clayton, D. (1982). A model for association in bivariate survival data. Journal of the Royal Statistical Society, Series B, Methodological, 44, 414–422.

See Also

biclaytoncop, binormalcop, binormal.

Examples

## Not run:  edge <- 0.01  # A small positive value
N <- 101; x <- seq(edge, 1.0 - edge, len = N); Rho <- 0.7
ox <- expand.grid(x, x)
zedd <- dbiclaytoncop(ox[, 1], ox[, 2], apar = Rho, log = TRUE)
par(mfrow = c(1, 2))
contour(x, x, matrix(zedd, N, N), col = 4, labcex = 1.5, las = 1)
plot(rbiclaytoncop(1000, 2), col = 4, las = 1)  
## End(Not run)

Bayesian Information Criterion

Description

Calculates the Bayesian information criterion (BIC) for a fitted model object for which a log-likelihood value has been obtained.

Usage

BICvlm(object, ..., k = log(nobs(object)))

Arguments

Details

The so-called BIC or SBC (Schwarz's Bayesian criterion) can be computed by calling AICvlm with a different k argument. See AICvlm for information and caveats.

Value

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

Warning

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

Many VGAM family functions such as cumulative can have the number of observations absorbed into the prior weights argument (e.g., weights in vglm), either before or after fitting. Almost all VGAM family functions can have the number of observations defined by the weights argument, e.g., as an observed frequency. BIC simply uses the number of rows of the model matrix, say, as defining n, hence the user must be very careful of this possible error. Use at your own risk!!

Note

BIC, AIC and other ICs can have have many additive constants added to them. The important thing are the differences since the minimum value corresponds to the best model.

BIC has not been defined for QRR-VGLMs yet.

Author(s)

T. W. Yee.

See Also

AICvlm, VGLMs are described in vglm-class; VGAMs are described in vgam-class; RR-VGLMs are described in rrvglm-class; BIC, AIC.

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)
BIC(fit1)
(fit2 <- vglm(cbind(normal, mild, severe) ~ let,
              cumulative(parallel = FALSE, reverse = TRUE), data = pneumo))
coef(fit2, matrix = TRUE)
BIC(fit2)

Farlie-Gumbel-Morgenstern's Bivariate Distribution Family Function

Description

Estimate the association parameter of Farlie-Gumbel-Morgenstern's bivariate distribution by maximum likelihood estimation.

Usage

bifgmcop(lapar = "rhobitlink", iapar = NULL, imethod = 1)

Arguments

Details

The cumulative distribution function is

P(Y_1 \leq y_1, Y_2 \leq y_2) = y_1 y_2 ( 1 + \alpha (1 - y_1) (1 - y_2) )

for -1 < \alpha < 1. The support of the function is the unit square. The marginal distributions are the standard uniform distributions. When \alpha = 0 the random variables are independent.

Value

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

Note

The response must be a two-column matrix. Currently, the fitted value is a matrix with two columns and values equal to 0.5. This is because each marginal distribution corresponds to a standard uniform distribution.

Author(s)

T. W. Yee

References

Castillo, E., Hadi, A. S., Balakrishnan, N. and Sarabia, J. S. (2005). Extreme Value and Related Models with Applications in Engineering and Science, Hoboken, NJ, USA: Wiley-Interscience.

Smith, M. D. (2007). Invariance theorems for Fisher information. Communications in Statistics—Theory and Methods, 36(12), 2213–2222.

See Also

rbifgmcop, bifrankcop, bifgmexp, simulate.vlm.

Examples

ymat <- rbifgmcop(1000, apar = rhobitlink(3, inverse = TRUE))
## Not run: plot(ymat, col = "blue")
fit <- vglm(ymat ~ 1, fam = bifgmcop, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
head(fitted(fit))

Farlie-Gumbel-Morgenstern's Bivariate Distribution

Description

Density, distribution function, and random generation for the (one parameter) bivariate Farlie-Gumbel-Morgenstern's distribution.

Usage

dbifgmcop(x1, x2, apar, log = FALSE)
pbifgmcop(q1, q2, apar)
rbifgmcop(n, apar)

Arguments

Details

See bifgmcop, the VGAM family functions for estimating the parameter by maximum likelihood estimation, for the formula of the cumulative distribution function and other details.

Value

dbifgmcop gives the density, pbifgmcop gives the distribution function, and rbifgmcop generates random deviates (a two-column matrix).

Author(s)

T. W. Yee

See Also

bifgmcop.

Examples

## Not run:  N <- 101; x <- seq(0.0, 1.0, len = N); apar <- 0.7
ox <- expand.grid(x, x)
zedd <- dbifgmcop(ox[, 1], ox[, 2], apar = apar)
contour(x, x, matrix(zedd, N, N), col = "blue")
zedd <- pbifgmcop(ox[, 1], ox[, 2], apar = apar)
contour(x, x, matrix(zedd, N, N), col = "blue")

plot(r <- rbifgmcop(n = 3000, apar = apar), col = "blue")
par(mfrow = c(1, 2))
hist(r[, 1])  # Should be uniform
hist(r[, 2])  # Should be uniform

## End(Not run)

Bivariate Farlie-Gumbel-Morgenstern Exponential Distribution Family Function

Description

Estimate the association parameter of FGM bivariate exponential distribution by maximum likelihood estimation.

Usage

bifgmexp(lapar = "rhobitlink", iapar = NULL, tola0 = 0.01,
         imethod = 1)

Arguments

Details

The cumulative distribution function is

P(Y_1 \leq y_1, Y_2 \leq y_2) = e^{-y_1-y_2} ( 1 + \alpha [1 - e^{-y_1}] [1 - e^{-y_2}] ) + 1 - e^{-y_1} - e^{-y_2}

for \alpha between -1 and 1. The support of the function is for y_1>0 and y_2>0. The marginal distributions are an exponential distribution with unit mean. When \alpha = 0 then the random variables are independent, and this causes some problems in the estimation process since the distribution no longer depends on the parameter.

A variant of Newton-Raphson is used, which only seems to work for an intercept model. It is a very good idea to set trace = TRUE.

Value

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

Note

The response must be a two-column matrix. Currently, the fitted value is a matrix with two columns and values equal to 1. This is because each marginal distribution corresponds to a exponential distribution with unit mean.

This VGAM family function should be used with caution.

Author(s)

T. W. Yee

References

Castillo, E., Hadi, A. S., Balakrishnan, N. and Sarabia, J. S. (2005). Extreme Value and Related Models with Applications in Engineering and Science, Hoboken, NJ, USA: Wiley-Interscience.

See Also

bifgmcop, bigumbelIexp.

Examples

N <- 1000; mdata <- data.frame(y1 = rexp(N), y2 = rexp(N))
## Not run: plot(ymat)
fit <- vglm(cbind(y1, y2) ~ 1, bifgmexp, data = mdata, trace = TRUE)
fit <- vglm(cbind(y1, y2) ~ 1, bifgmexp, data = mdata, # May fail
            trace = TRUE, crit = "coef")
coef(fit, matrix = TRUE)
Coef(fit)
head(fitted(fit))

Frank's Bivariate Distribution Family Function

Description

Estimate the association parameter of Frank's bivariate distribution by maximum likelihood estimation.

Usage

bifrankcop(lapar = "loglink", iapar = 2, nsimEIM = 250)

Arguments

Details

The cumulative distribution function is

P(Y_1 \leq y_1, Y_2 \leq y_2) = H_{\alpha}(y_1,y_2) = \log_{\alpha} [1 + (\alpha^{y_1}-1)(\alpha^{y_2}-1)/ (\alpha-1)]

for \alpha \ne 1. Note the logarithm here is to base \alpha. The support of the function is the unit square.

When 0 < \alpha < 1 the probability density function h_{\alpha}(y_1,y_2) is symmetric with respect to the lines y_2=y_1 and y_2=1-y_1. When \alpha > 1 then h_{\alpha}(y_1,y_2) = h_{1/\alpha}(1-y_1,y_2).

\alpha=1 then H(y_1,y_2) = y_1 y_2, i.e., uniform on the unit square. As \alpha approaches 0 then H(y_1,y_2) = \min(y_1,y_2). As \alpha approaches infinity then H(y_1,y_2) = \max(0, y_1+y_2-1).

The default is to use Fisher scoring implemented using rbifrankcop. For intercept-only models an alternative is to set nsimEIM=NULL so that a variant of Newton-Raphson is used.

Value

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

Note

The response must be a two-column matrix. Currently, the fitted value is a matrix with two columns and values equal to a half. This is because the marginal distributions correspond to a standard uniform distribution.

Author(s)

T. W. Yee

References

Genest, C. (1987). Frank's family of bivariate distributions. Biometrika, 74, 549–555.

See Also

rbifrankcop, bifgmcop, simulate.vlm.

Examples

## Not run: 
ymat <- rbifrankcop(n = 2000, apar = exp(4))
plot(ymat, col = "blue")
fit <- vglm(ymat ~ 1, fam = bifrankcop, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
vcov(fit)
head(fitted(fit))
summary(fit)

## End(Not run)

Gumbel's Type I Bivariate Distribution Family Function

Description

Estimate the association parameter of Gumbel's Type I bivariate distribution by maximum likelihood estimation.

Usage

bigumbelIexp(lapar = "identitylink", iapar = NULL, imethod = 1)

Arguments

Details

The cumulative distribution function is

P(Y_1 \leq y_1, Y_2 \leq y_2) = e^{-y_1-y_2+\alpha y_1 y_2} + 1 - e^{-y_1} - e^{-y_2}

for real \alpha. The support of the function is for y_1>0 and y_2>0. The marginal distributions are an exponential distribution with unit mean.

A variant of Newton-Raphson is used, which only seems to work for an intercept model. It is a very good idea to set trace=TRUE.

Value

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

Note

The response must be a two-column matrix. Currently, the fitted value is a matrix with two columns and values equal to 1. This is because each marginal distribution corresponds to a exponential distribution with unit mean.

This VGAM family function should be used with caution.

Author(s)

T. W. Yee

References

Gumbel, E. J. (1960). Bivariate Exponential Distributions. Journal of the American Statistical Association, 55, 698–707.

See Also

bifgmexp.

Examples

nn <- 1000
gdata <- data.frame(y1 = rexp(nn), y2 = rexp(nn))
## Not run:  with(gdata, plot(cbind(y1, y2))) 
fit <- vglm(cbind(y1, y2) ~ 1, bigumbelIexp, gdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
head(fitted(fit))

Bivariate Logistic Distribution

Description

Density, distribution function, quantile function and random generation for the 4-parameter bivariate logistic distribution.

Usage

dbilogis(x1, x2, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1,
         log = FALSE)
pbilogis(q1, q2, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1)
rbilogis(n, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1)

Arguments

Details

See bilogis, the VGAM family function for estimating the four parameters by maximum likelihood estimation, for the formula of the cumulative distribution function and other details.

Value

dbilogis gives the density, pbilogis gives the distribution function, and rbilogis generates random deviates (a two-column matrix).

Note

Gumbel (1961) proposed two bivariate logistic distributions with logistic distribution marginals, which he called Type I and Type II. The Type I is this one. The Type II belongs to the Morgenstern type. The biamhcop distribution has, as a special case, this distribution, which is when the random variables are independent.

Author(s)

T. W. Yee

References

Gumbel, E. J. (1961). Bivariate logistic distributions. Journal of the American Statistical Association, 56, 335–349.

See Also

bilogistic, biamhcop.

Examples

## Not run:  par(mfrow = c(1, 3))
ymat <- rbilogis(n = 2000, loc1 = 5, loc2 = 7, scale2 = exp(1))
myxlim <- c(-2, 15); myylim <- c(-10, 30)
plot(ymat, xlim = myxlim, ylim = myylim)

N <- 100
x1 <- seq(myxlim[1], myxlim[2], len = N)
x2 <- seq(myylim[1], myylim[2], len = N)
ox <- expand.grid(x1, x2)
z <- dbilogis(ox[,1], ox[,2], loc1 = 5, loc2 = 7, scale2 = exp(1))
contour(x1, x2, matrix(z, N, N), main = "density")
z <- pbilogis(ox[,1], ox[,2], loc1 = 5, loc2 = 7, scale2 = exp(1))
contour(x1, x2, matrix(z, N, N), main = "cdf") 
## End(Not run)

Bivariate Logistic Distribution Family Function

Description

Estimates the four parameters of the bivariate logistic distribution by maximum likelihood estimation.

Usage

bilogistic(llocation = "identitylink", lscale = "loglink",
           iloc1 = NULL, iscale1 = NULL, iloc2 = NULL, iscale2 =
           NULL, imethod = 1, nsimEIM = 250, zero = NULL)

Arguments

Details

The four-parameter bivariate logistic distribution has a density that can be written as

f(y_1,y_2;l_1,s_1,l_2,s_2) = 2 \frac{\exp[-(y_1-l_1)/s_1 - (y_2-l_2)/s_2]}{ s_1 s_2 \left( 1 + \exp[-(y_1-l_1)/s_1] + \exp[-(y_2-l_2)/s_2] \right)^3}

where s_1>0 and s_2>0 are the scale parameters, and l_1 and l_2 are the location parameters. Each of the two responses are unbounded, i.e., -\infty<y_j<\infty. The mean of Y_1 is l_1 etc. The fitted values are returned in a 2-column matrix. The cumulative distribution function is

F(y_1,y_2;l_1,s_1,l_2,s_2) = \left( 1 + \exp[-(y_1-l_1)/s_1] + \exp[-(y_2-l_2)/s_2] \right)^{-1}

The marginal distribution of Y_1 is

P(Y_1 \leq y_1) = F(y_1;l_1,s_1) = \left( 1 + \exp[-(y_1-l_1)/s_1] \right)^{-1} .

By default, \eta_1=l_1, \eta_2=\log(s_1), \eta_3=l_2, \eta_4=\log(s_2) are the linear/additive predictors.

Value

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

Author(s)

T. W. Yee

References

Gumbel, E. J. (1961). Bivariate logistic distributions. Journal of the American Statistical Association, 56, 335–349.

Castillo, E., Hadi, A. S., Balakrishnan, N. and Sarabia, J. S. (2005). Extreme Value and Related Models with Applications in Engineering and Science, Hoboken, NJ, USA: Wiley-Interscience.

See Also

logistic, rbilogis.

Examples

## Not run: 
ymat <- rbilogis(n <- 50, loc1 = 5, loc2 = 7, scale2 = exp(1))
plot(ymat)
bfit <- vglm(ymat ~ 1, family = bilogistic, trace = TRUE)
coef(bfit, matrix = TRUE)
Coef(bfit)
head(fitted(bfit))
vcov(bfit)
head(weights(bfit, type = "work"))
summary(bfit)

## End(Not run)

Bivariate Binary Regression with an Odds Ratio (Family Function)

Description

Fits a Palmgren (bivariate odds-ratio model, or bivariate logistic regression) model to two binary responses. Actually, a bivariate logistic/probit/cloglog/cauchit model can be fitted. The odds ratio is used as a measure of dependency.

Usage

binom2.or(lmu = "logitlink", lmu1 = lmu, lmu2 = lmu, loratio = "loglink",
          imu1 = NULL, imu2 = NULL, ioratio = NULL, zero = "oratio",
          exchangeable = FALSE, tol = 0.001, more.robust = FALSE)

Arguments

Details

Also known informally as the Palmgren model, the bivariate logistic model is a full-likelihood based model defined as two logistic regressions plus log(oratio) = eta3 where eta3 is the third linear/additive predictor relating the odds ratio to explanatory variables. Explicitly, the default model is

logit[P(Y_j=1)] = \eta_j,\ \ \ j=1,2

for the marginals, and

\log[P(Y_{00}=1) P(Y_{11}=1) / (P(Y_{01}=1) P(Y_{10}=1))] = \eta_3,

specifies the dependency between the two responses. Here, the responses equal 1 for a success and a 0 for a failure, and the odds ratio is often written \psi=p_{00}p_{11}/(p_{10}p_{01}). The model is fitted by maximum likelihood estimation since the full likelihood is specified. The two binary responses are independent if and only if the odds ratio is unity, or equivalently, the log odds ratio is 0. Fisher scoring is implemented.

The default models \eta_3 as a single parameter only, i.e., an intercept-only model, but this can be circumvented by setting zero = NULL in order to model the odds ratio as a function of all the explanatory variables. The function binom2.or() can handle other probability link functions such as probitlink, clogloglink and cauchitlink links as well, so is quite general. In fact, the two marginal probabilities can each have a different link function. A similar model is the bivariate probit model (binom2.rho), which is based on a standard bivariate normal distribution, but the bivariate probit model is less interpretable and flexible.

The exchangeable argument should be used when the error structure is exchangeable, e.g., with eyes or ears data.

Value

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

When fitted, the fitted.values slot of the object contains the four joint probabilities, labelled as (Y_1,Y_2) = (0,0), (0,1), (1,0), (1,1), respectively. These estimated probabilities should be extracted with the fitted generic function.

Note

At present we call binom2.or families a bivariate odds-ratio model. The response should be either a 4-column matrix of counts (whose columns correspond to (Y_1,Y_2) = (0,0), (0,1), (1,0), (1,1) respectively), or a two-column matrix where each column has two distinct values, or a factor with four levels. The function rbinom2.or may be used to generate such data. Successful convergence requires at least one case of each of the four possible outcomes.

By default, a constant odds ratio is fitted because zero = 3. Set zero = NULL if you want the odds ratio to be modelled as a function of the explanatory variables; however, numerical problems are more likely to occur.

The argument lmu, which is actually redundant, is used for convenience and for upward compatibility: specifying lmu only means the link function will be applied to lmu1 and lmu2. Users who want a different link function for each of the two marginal probabilities should use the lmu1 and lmu2 arguments, and the argument lmu is then ignored. It doesn't make sense to specify exchangeable = TRUE and have different link functions for the two marginal probabilities.

Regarding Yee and Dirnbock (2009), the xij (see vglm.control) argument enables environmental variables with different values at the two time points to be entered into an exchangeable binom2.or model. See the author's webpage for sample code.

Author(s)

Thomas W. Yee

References

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models, 2nd ed. London: Chapman & Hall.

le Cessie, S. and van Houwelingen, J. C. (1994). Logistic regression for correlated binary data. Applied Statistics, 43, 95–108.

Palmgren, J. (1989). Regression Models for Bivariate Binary Responses. Technical Report no. 101, Department of Biostatistics, University of Washington, Seattle.

Yee, T. W. and Dirnbock, T. (2009). Models for analysing species' presence/absence data at two time points. Journal of Theoretical Biology, 259(4), 684–694.

See Also

rbinom2.or, binom2.rho, loglinb2, binom3.or, loglinb3, coalminers, binomialff, logitlink, probitlink, clogloglink, cauchitlink.

Examples

# Fit the model in Table 6.7 in McCullagh and Nelder (1989)
coalminers <- transform(coalminers, Age = (age - 42) / 5)
fit <- vglm(cbind(nBnW, nBW, BnW, BW) ~ Age,
            binom2.or(zero = NULL), data = coalminers)
fitted(fit)
summary(fit)
coef(fit, matrix = TRUE)
c(weights(fit, type = "prior")) * fitted(fit)  # Table 6.8

## Not run:  with(coalminers, matplot(Age, fitted(fit), type = "l", las = 1,
                         xlab = "(age - 42) / 5", lwd = 2))
with(coalminers, matpoints(Age, depvar(fit), col=1:4))
legend(x = -4, y = 0.5, lty = 1:4, col = 1:4, lwd = 2,
       legend = c("1 = (Breathlessness=0, Wheeze=0)",
                  "2 = (Breathlessness=0, Wheeze=1)",
                  "3 = (Breathlessness=1, Wheeze=0)",
                  "4 = (Breathlessness=1, Wheeze=1)")) 
## End(Not run)


# Another model: pet ownership
## Not run:  data(xs.nz, package = "VGAMdata")
# More homogeneous:
petdata <- subset(xs.nz, ethnicity == "European" & age < 70 &
                         sex == "M")
petdata <- na.omit(petdata[, c("cat", "dog", "age")])
summary(petdata)
with(petdata, table(cat, dog))  # Can compute the odds ratio

fit <- vgam(cbind((1-cat) * (1-dog), (1-cat) * dog,
                     cat  * (1-dog),    cat  * dog) ~ s(age, df = 5),
            binom2.or(zero =    3), data = petdata, trace = TRUE)
colSums(depvar(fit))
coef(fit, matrix = TRUE)

## End(Not run)

## Not run:  # Plot the estimated probabilities
ooo <- order(with(petdata, age))
matplot(with(petdata, age)[ooo], fitted(fit)[ooo, ], type = "l",
        xlab = "Age", ylab = "Probability", main = "Pet ownership",
        ylim = c(0, max(fitted(fit))), las = 1, lwd = 1.5)
legend("topleft", col=1:4, lty = 1:4, leg = c("no cat or dog ",
       "dog only", "cat only", "cat and dog"), lwd = 1.5) 
## End(Not run)

Bivariate Odds Ratio Model

Description

Density and random generation for a bivariate binary regression model using an odds ratio as the measure of dependency.

Usage

rbinom2.or(n, mu1, mu2 = if (exchangeable) mu1 else
    stop("'mu2' not specified"),
    oratio = 1, exchangeable = FALSE, tol = 0.001,
    twoCols = TRUE, colnames = if (twoCols) c("y1","y2")
    else c("00", "01", "10", "11"), ErrorCheck = TRUE)
dbinom2.or(mu1, mu2 = if (exchangeable) mu1 else
    stop("'mu2' not specified"), oratio = 1,
    exchangeable = FALSE, tol = 0.001,
    colnames = c("00", "01", "10", "11"), ErrorCheck = TRUE)

Arguments

Details

The function rbinom2.or generates data coming from a bivariate binary response model. The data might be fitted with the VGAM family function binom2.or.

The function dbinom2.or does not really compute the density (because that does not make sense here) but rather returns the four joint probabilities.

Value

The function rbinom2.or returns either a 2 or 4 column matrix of 1s and 0s, depending on the argument twoCols.

The function dbinom2.or returns a 4 column matrix of joint probabilities; each row adds up to unity.

Author(s)

T. W. Yee

See Also

binom2.or.

Examples

nn <- 1000  # Example 1
ymat <- rbinom2.or(nn, mu1 = logitlink(1, inv = TRUE),
                   oratio = exp(2), exch = TRUE)
(mytab <- table(ymat[, 1], ymat[, 2], dnn = c("Y1", "Y2")))
(myor <- mytab["0","0"] * mytab["1","1"] / (mytab["1","0"] *
         mytab["0","1"]))
fit <- vglm(ymat ~ 1, binom2.or(exch = TRUE))
coef(fit, matrix = TRUE)

bdata <- data.frame(x2 = sort(runif(nn)))  # Example 2
bdata <- transform(bdata,
           mu1 = logitlink(-2 + 4 * x2, inverse = TRUE),
           mu2 = logitlink(-1 + 3 * x2, inverse = TRUE))
dmat <- with(bdata, dbinom2.or(mu1 = mu1, mu2 = mu2,
                               oratio = exp(2)))
ymat <- with(bdata, rbinom2.or(n = nn, mu1 = mu1, mu2 = mu2,
                               oratio = exp(2)))
fit2 <- vglm(ymat ~ x2, binom2.or, data = bdata)
coef(fit2, matrix = TRUE)
## Not run: 
matplot(with(bdata, x2), dmat, lty = 1:4, col = 1:4,
        main = "Joint probabilities", ylim = 0:1, type = "l",
        ylab = "Probabilities", xlab = "x2", las = 1)
legend("top", lty = 1:4, col = 1:4,
       legend = c("1 = (y1=0, y2=0)", "2 = (y1=0, y2=1)",
                  "3 = (y1=1, y2=0)", "4 = (y1=1, y2=1)"))

## End(Not run)

Bivariate Probit Regression

Description

Fits a bivariate probit model to two binary responses. The correlation parameter rho is the measure of dependency.

Usage

binom2.rho(lmu = "probitlink", lrho = "rhobitlink",
           imu1 = NULL, imu2 = NULL,
           irho = NULL, imethod = 1, zero = "rho",
           exchangeable = FALSE, grho = seq(-0.95, 0.95, by = 0.05),
           nsimEIM = NULL)
binom2.Rho(rho = 0, imu1 = NULL, imu2 = NULL,
           exchangeable = FALSE, nsimEIM = NULL)

Arguments

Details

The bivariate probit model was one of the earliest regression models to handle two binary responses jointly. It has a probit link for each of the two marginal probabilities, and models the association between the responses by the \rho parameter of a standard bivariate normal distribution (with zero means and unit variances). One can think of the joint probabilities being \Phi(\eta_1,\eta_2;\rho) where \Phi is the cumulative distribution function of a standard bivariate normal distribution.

Explicitly, the default model is

probit[P(Y_j=1)] = \eta_j,\ \ \ j=1,2

for the marginals, and

rhobit[rho] = \eta_3.

The joint probability P(Y_1=1,Y_2=1)= \Phi(\eta_1,\eta_2;\rho), and from these the other three joint probabilities are easily computed. The model is fitted by maximum likelihood estimation since the full likelihood is specified. Fisher scoring is implemented.

The default models \eta_3 as a single parameter only, i.e., an intercept-only model for rho, but this can be circumvented by setting zero = NULL in order to model rho as a function of all the explanatory variables.

The bivariate probit model should not be confused with a bivariate logit model with a probit link (see binom2.or). The latter uses the odds ratio to quantify the association. Actually, the bivariate logit model is recommended over the bivariate probit model because the odds ratio is a more natural way of measuring the association between two binary responses.

Value

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

When fitted, the fitted.values slot of the object contains the four joint probabilities, labelled as (Y_1,Y_2) = (0,0), (0,1), (1,0), (1,1), respectively.

Note

See binom2.or about the form of input the response should have.

By default, a constant \rho is fitted because zero = "rho". Set zero = NULL if you want the \rho parameter to be modelled as a function of the explanatory variables. The value \rho lies in the interval (-1,1), therefore a rhobitlink link is default.

Converge problems can occur. If so, assign irho a range of values and monitor convergence (e.g., set trace = TRUE). Else try imethod. Practical experience shows that local solutions can occur, and that irho needs to be quite close to the (global) solution. Also, imu1 and imu2 may be used.

This help file is mainly about binom2.rho(). binom2.Rho() fits a bivariate probit model with known \rho. The inputted rho is saved in the misc slot of the fitted object, with rho as the component name.

In some econometrics applications (e.g., Freedman 2010, Freedman and Sekhon 2010) one response is used as an explanatory variable, e.g., a recursive binomial probit model. Such will not work here. Historically, the bivariate probit model was the first VGAM I ever wrote, based on Ashford and Sowden (1970). I don't think they ever thought of it either! Hence the criticisms raised go beyond the use of what was originally intended.

Author(s)

Thomas W. Yee

References

Ashford, J. R. and Sowden, R. R. (1970). Multi-variate probit analysis. Biometrics, 26, 535–546.

Freedman, D. A. (2010). Statistical Models and Causal Inference: a Dialogue with the Social Sciences, Cambridge: Cambridge University Press.

Freedman, D. A. and Sekhon, J. S. (2010). Endogeneity in probit response models. Political Analysis, 18, 138–150.

See Also

rbinom2.rho, rhobitlink, pbinorm, binom2.or, loglinb2, coalminers, binomialff, rhobitlink, fisherzlink.

Examples

coalminers <- transform(coalminers, Age = (age - 42) / 5)
fit <- vglm(cbind(nBnW, nBW, BnW, BW) ~ Age,
            binom2.rho, data = coalminers, trace = TRUE)
summary(fit)
coef(fit, matrix = TRUE)

Bivariate Probit Model

Description

Density and random generation for a bivariate probit model. The correlation parameter rho is the measure of dependency.

Usage

rbinom2.rho(n, mu1, mu2 = if (exchangeable) mu1 else
  stop("'mu2' not specified"), rho = 0,
  exchangeable = FALSE, twoCols = TRUE,
  colnames = if (twoCols) c("y1","y2") else
  c("00", "01", "10", "11"), ErrorCheck = TRUE)
dbinom2.rho(mu1, mu2 = if (exchangeable) mu1 else
  stop("'mu2' not specified"), rho = 0, exchangeable = FALSE,
  colnames = c("00", "01", "10", "11"), ErrorCheck = TRUE)

Arguments

Details

The function rbinom2.rho generates data coming from a bivariate probit model. The data might be fitted with the VGAM family function binom2.rho.

The function dbinom2.rho does not really compute the density (because that does not make sense here) but rather returns the four joint probabilities.

Value

The function rbinom2.rho returns either a 2 or 4 column matrix of 1s and 0s, depending on the argument twoCols.

The function dbinom2.rho returns a 4 column matrix of joint probabilities; each row adds up to unity.

Author(s)

T. W. Yee

See Also

binom2.rho.

Examples

(myrho <- rhobitlink(2, inverse = TRUE))  # Example 1
nn <- 2000
ymat <- rbinom2.rho(nn, mu1 = 0.8, rho = myrho, exch = TRUE)
(mytab <- table(ymat[, 1], ymat[, 2], dnn = c("Y1", "Y2")))
fit <- vglm(ymat ~ 1, binom2.rho(exch = TRUE))
coef(fit, matrix = TRUE)

bdata <- data.frame(x2 = sort(runif(nn)))  # Example 2
bdata <- transform(bdata, mu1 = probitlink(-2+4*x2, inv = TRUE),
                          mu2 = probitlink(-1+3*x2, inv = TRUE))
dmat <- with(bdata, dbinom2.rho(mu1, mu2, myrho))
ymat <- with(bdata, rbinom2.rho(nn, mu1, mu2, myrho))
fit2 <- vglm(ymat ~ x2, binom2.rho, data = bdata)
coef(fit2, matrix = TRUE)
## Not run:  matplot(with(bdata, x2), dmat, lty = 1:4, col = 1:4,
        type = "l", main = "Joint probabilities",
        ylim = 0:1, lwd = 2, ylab = "Probability")
legend(x = 0.25, y = 0.9, lty = 1:4, col = 1:4, lwd = 2,
       legend = c("1 = (y1=0, y2=0)", "2 = (y1=0, y2=1)",
                  "3 = (y1=1, y2=0)", "4 = (y1=1, y2=1)")) 
## End(Not run)

Trivariate Binary Regression with Three Odds Ratios (Family Function)

Description

Fits three Palmgren (bivariate odds-ratio model, or bivariate logistic regression) models simultaneously to three binary responses. The odds ratios are used to measure dependencies between the responses. Several options for the joint probabilities are available.

Usage

binom3.or(lmu = "logitlink", lmu1 = lmu, lmu2 = lmu,
  lmu3 = lmu, loratio = "loglink", zero = "oratio",
  exchangeable = FALSE, eq.or = FALSE, jpmethod =
  c("min", "mean", "median", "max", "1", "2", "3"),
  imu1 = NULL, imu2 = NULL, imu3 = NULL,
  ioratio12 = NULL, ioratio13 = NULL, ioratio23 = NULL,
  tol = 0.001, more.robust = FALSE)

Arguments

Details

This heuristic model is an extension of binom2.or for handling three binary responses. Rather than allowing something like vglm(cbind(y1,y2, y1,y3, y2,y3) ~ x2, binom2.or), which has three pairs of bivariate responses in the usual form of multiple responses allowed in VGAM, I have decided to write binom3.or which operates on cbind(y1, y2, y3) instead. This model thus uses three odds ratios to allow for dependencies between the pairs of responses. It is heuristic because the joint probability P(y_1=1,y_2=1,y_3=1)=p_{123} is computed by a number of conditional independence assumptions.

This trivariate logistic model has a fully specified likelihood. Explicitly, the default model is

logit\;P(Y_j=1)] = \eta_j,\ \ \ j=1,2

for the first two marginals,

logit\; P(Y_1=1) = \eta_4,

logit\; P(Y_3=1) = \eta_5,

logit\; P(Y_2=1) = \eta_7,

logit\; P(Y_3=1) = \eta_8,

and

\log \psi_{12} = \eta_3,

\log \psi_{13} = \eta_6,

\log \psi_{23} = \eta_9,

specifies the dependency between each possible pair of responses. Many details on such quantities are similar to binom2.or.

By default, all odds ratios are intercept-only. The exchangeable argument should be used when the error structure is exchangeable. However, there is a difference between full and partial exchangeability and setting exchangeable = TRUE results in the full version. The partial version would require the manual input of certain constraint matrices via constraints.

Value

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

When fitted, the fitted.values slot of the object is a matrix with successive columns equalling the eight joint probabilities, labelled as (Y_1,Y_2,Y_3) = (0,0,0), (0,0,1), (0,1,0), (0,1,1), (1,0,0), (1,0,1), (1,1,0), (1,1,1), respectively. These estimated probabilities should be extracted with the fitted generic function.

Warning

Because the parameter space of this model is restricted (see rbinom3.or), this family function is more limited than loglinb3. However, this model is probably more interpretable since the marginal probabilities and odds ratios are modelled by conventional link functions directly.

If the data is very sparse then convergence problems will occur. It is recommended that the sample size is several hundred at least. Opinion: anything less than n=100 is liable for failure. Setting trace = TRUE is urged.

Note

At present we call binom3.or a trivariate odds-ratio model (TOM). The response should be either a 8-column matrix of counts (whose columns correspond to (Y_1,Y_2,Y_3) ordered as above), or a three-column matrix where each column has two distinct values, or a factor with 8 levels. The function rbinom3.or may be used to generate such data.

Because some of the \eta_j are repeated, the constraint matrices have a special form in order to provide consistency.

By default, intercept-only odds ratios are fitted because zero = "oratio". Set zero = NULL for the odds ratios to be modelled as a function of the explanatory variables; however, numerical problems are more likely to occur.

The argument lmu, which is actually redundant, is used for convenience and for upward compatibility: specifying lmu only means the link function will be applied to lmu1, lmu2 and lmu3. Users who want a different link function for each of the marginal probabilities should use lmu1, lmu2 and lmu3, and the argument lmu is then ignored. It doesn't make sense to specify exchangeable = TRUE and have different link functions for the marginal probabilities.

See Also

rbinom3.or, binom2.or, loglinb3.

Examples

  set.seed(1)
## Not run: 
nn <- 1000  # Example 1
ymat <- rbinom3.or(nn, mu1 = logitlink(0.5, inv = TRUE),
                   oratio12 = exp(1), exch = TRUE)
fit1 <- vglm(ymat ~ 1, binom3.or(exc = TRUE), tra = TRUE)
coef(fit1, matrix = TRUE)
constraints(fit1)

bdata <- data.frame(x2 = sort(runif(nn)))  # Example 2
bdata <- transform(bdata,
       mu1 = logitlink(-1 + 1 * x2, inv = TRUE),
       mu2 = logitlink(-1 + 2 * x2, inv = TRUE),
       mu3 = logitlink( 2 - 1 * x2, inv = TRUE))
ymat2 <- with(bdata,
     rbinom3.or(nn, mu1, mu2, mu3, exp(0.25),
                oratio13 = exp(0.25), exp(0.25)))
fit2 <- vglm(ymat2 ~ x2, binom3.or(eq.or = TRUE),
             bdata, trace = TRUE)
coef(fit2, matrix = TRUE)

## End(Not run)

Trivariate Odds Ratio Model

Description

Density and random generation for a trivariate binary regression model using three odds ratios to measure dependencies.

Usage

dbinom3.or(mu1, mu2 = if (exchangeable) mu1 else
   stop("'mu2' not specified"), mu3 = if (exchangeable)
   mu3 else stop("'mu3' not specified"), oratio12 = 1, 
   oratio13 = 1, oratio23 = 1, exchangeable = FALSE,
   jpmethod = c("min", "mean", "median", "max", "1", "2", "3"),
   tol = 0.001, ErrorCheck = TRUE)
rbinom3.or(n, mu1, mu2 = if (exchangeable) mu1 else
   stop("'mu2' not specified"), mu3 = if (exchangeable) mu1
   else stop("'mu3' not specified"), oratio12 = 1,
   oratio13 = 1, oratio23 = 1, exchangeable = FALSE,
   jpmethod = c("min", "mean", "median", "max", "1", "2", "3"),
   threeCols = TRUE, tol = 0.001, ErrorCheck = TRUE)

Arguments

Details

The function dbinom3.or does not really compute the density (because that does not make sense here) but rather returns the eight joint probabilities if the parameters are in the parameter space. Simulations have shown that if all the marginal probabilities are uniformly distributed and all the odds ratios have a standard lognormal distribution (with joint independence) then about 31 percent of the parameter space is valid. With exchangeability, it is about 33 percent. This means that binom3.or has quite some severe limitations for general use.

The function rbinom3.or generates data coming from a trivariate binary response model. Valid data from this might be fitted with the VGAM family function binom3.or. Any invalid data (because the parameters are outside the parameter space) are NaNs.

Value

The function dbinom3.or returns a 8 column matrix of joint probabilities; each row adds up to unity if the parameters are in the parameter space. If not, then NaNs are returned.

The function rbinom3.or returns either a 3 or 8 column matrix of 1s and 0s, depending on the argument threeCols.

References

Yee, T. W. (2024). New regression methods for three or four binary responses. In preparation.

See Also

binom3.or, is.nan.

Examples

dbinom3.or(0.5, 0.5, 0.5, 1, 2, 2)
# Outside the parameter space:
dbinom3.or(0.9, 0.9, 0.9, 1, 2, 2)

## Not run:  nn <- 100000
for (Exch in c(TRUE, FALSE)) {
  zdata <- data.frame(orat12 = rlnorm(nn), p1 = runif(nn))
  zdata <- transform(zdata,
             orat13 = if (Exch) orat12 else rlnorm(nn),
             orat23 = if (Exch) orat12 else rlnorm(nn),
             p2 = if (Exch) p1 else runif(nn),
             p3 = if (Exch) p1 else runif(nn))

  mat1 <- with(zdata, dbinom3.or(p1, p2, p3,
               orat12, orat13, orat23, exch = Exch))
# Important statistic: Pr(in the parameter space) =
print(1 - nrow(na.omit(mat1)) / nrow(mat1))
}
round(head(mat1), 4)

## End(Not run)

Binomial Family Function

Description

Family function for fitting generalized linear models to binomial responses

Usage

binomialff(link = "logitlink", multiple.responses = FALSE,
       parallel = FALSE, zero = NULL, bred = FALSE, earg.link = FALSE)

Arguments

Details

This function is largely to mimic binomial, however there are some differences.

When used with cqo and cao, it may be preferable to use the clogloglink link.

Value

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

Warning

See the above note regarding bred.

The maximum likelihood estimate will not exist if the data is completely separable or quasi-completely separable. See Chapter 10 of Altman et al. (2004) for more details, and safeBinaryRegression and hdeff.vglm. Yet to do: add a sepcheck = TRUE, say, argument to further detect this problem and give an appropriate warning.

Note

If multiple.responses is FALSE (default) then the response can be of one of two formats: a factor (first level taken as failure), or a 2-column matrix (first column = successes) of counts. The argument weights in the modelling function can also be specified as any vector of positive values. In general, 1 means success and 0 means failure (to check, see the y slot of the fitted object). Note that a general vector of proportions of success is no longer accepted.

The notation M is used to denote the number of linear/additive predictors.

If multiple.responses is TRUE, then the matrix response can only be of one format: a matrix of 1's and 0's (1 = success).

Fisher scoring is used. This can sometimes fail to converge by oscillating between successive iterations (Ridout, 1990). See the example below.

Author(s)

Thomas W. Yee

References

McCullagh, P. and Nelder, J. A. (1989). Generalized Linear Models, 2nd ed. London: Chapman & Hall.

Altman, M. and Gill, J. and McDonald, M. P. (2004). Numerical Issues in Statistical Computing for the Social Scientist, Hoboken, NJ, USA: Wiley-Interscience.

Ridout, M. S. (1990). Non-convergence of Fisher's method of scoring—a simple example. GLIM Newsletter, 20(6).

See Also

hdeff.vglm, Links, wsdm, alogitlink, asinlink, N1binomial, rrvglm, cqo, cao, betabinomial, posbinomial, zibinomial, binom2.or, binom3.or, double.expbinomial, seq2binomial, amlbinomial, simplex, binomial, simulate.vlm, safeBinaryRegression, residualsvglm.

Examples

shunua <- hunua[sort.list(with(hunua, altitude)), ]  # Sort by altitude
fit <- vglm(agaaus ~ poly(altitude, 2), binomialff(link = clogloglink),
            data = shunua)
## Not run: 
plot(agaaus ~ jitter(altitude), shunua, ylab = "Pr(Agaaus = 1)",
     main = "Presence/absence of Agathis australis", col = 4, las = 1)
with(shunua, lines(altitude, fitted(fit), col = "orange", lwd = 2)) 
## End(Not run)

# Fit two species simultaneously
fit2 <- vgam(cbind(agaaus, kniexc) ~ s(altitude),
             binomialff(multiple.responses = TRUE), data = shunua)
## Not run: 
with(shunua, matplot(altitude, fitted(fit2), type = "l",
     main = "Two species response curves", las = 1)) 
## End(Not run)

# Shows that Fisher scoring can sometime fail. See Ridout (1990).
ridout <- data.frame(v = c(1000, 100, 10), r = c(4, 3, 3), n = rep(5, 3))
(ridout <- transform(ridout, logv = log(v)))
# The iterations oscillates between two local solutions:
glm.fail <- glm(r / n ~ offset(logv) + 1, weight = n,
               binomial(link = 'cloglog'), ridout, trace = TRUE)
coef(glm.fail)
# vglm()'s half-stepping ensures the MLE of -5.4007 is obtained:
vglm.ok <- vglm(cbind(r, n-r) ~ offset(logv) + 1,
               binomialff(link = clogloglink), ridout, trace = TRUE)
coef(vglm.ok)
wsdm(vglm.ok)

# Separable data
set.seed(123)
threshold <- 0
bdata <- data.frame(x2 = sort(rnorm(nn <- 100)))
bdata <- transform(bdata, y1 = ifelse(x2 < threshold, 0, 1))
fit <- vglm(y1 ~ x2, binomialff(bred = TRUE),
            data = bdata, criter = "coef", trace = TRUE)
coef(fit, matrix = TRUE)  # Finite!!
summary(fit, wsdm = TRUE, hdiff = 0.0001)

## Not run:  plot(depvar(fit) ~ x2, data = bdata, col = "blue", las = 1)
lines(fitted(fit) ~ x2, data = bdata, col = "orange")
abline(v = threshold, col = "gray", lty = "dashed") 
## End(Not run)

Bivariate Normal Distribution Cumulative Distribution Function

Description

Density, cumulative distribution function and random generation for the bivariate normal distribution distribution.

Usage

dbinorm(x1, x2, mean1 = 0, mean2 = 0, var1 = 1, var2 = 1, cov12 = 0,
        log = FALSE)
pbinorm(q1, q2, mean1 = 0, mean2 = 0, var1 = 1, var2 = 1, cov12 = 0)
rbinorm(n,      mean1 = 0, mean2 = 0, var1 = 1, var2 = 1, cov12 = 0)
 pnorm2(x1, x2, mean1 = 0, mean2 = 0, var1 = 1, var2 = 1, cov12 = 0)

Arguments

Details

The default arguments correspond to the standard bivariate normal distribution with correlation parameter \rho = 0. That is, two independent standard normal distributions. Let sd1 (say) be sqrt(var1) and written \sigma_1, etc. Then the general formula for the correlation coefficient is \rho = cov / (\sigma_1 \sigma_2) where cov is argument cov12. Thus if arguments var1 and var2 are left alone then cov12 can be inputted with \rho.

One can think of this function as an extension of pnorm to two dimensions, however note that the argument names have been changed for VGAM 0.9-1 onwards.

Value

dbinorm gives the density, pbinorm gives the cumulative distribution function, rbinorm generates random deviates (n by 2 matrix).

Warning

Being based on an approximation, the results of pbinorm() may be negative! Also, pnorm2() should be withdrawn soon; use pbinorm() instead because it is identical.

Note

For rbinorm(), if the ith variance-covariance matrix is not positive-definite then the ith row is all NAs.

References

pbinorm() is based on Donnelly (1973), the code was translated from FORTRAN to ratfor using struct, and then from ratfor to C manually. The function was originally called bivnor, and TWY only wrote a wrapper function.

Donnelly, T. G. (1973). Algorithm 462: Bivariate Normal Distribution. Communications of the ACM, 16, 638.

See Also

pnorm, binormal, uninormal.

Examples

yvec <- c(-5, -1.96, 0, 1.96, 5)
ymat <- expand.grid(yvec, yvec)
cbind(ymat, pbinorm(ymat[, 1], ymat[, 2]))

## Not run:  rhovec <- seq(-0.95, 0.95, by = 0.01)
plot(rhovec, pbinorm(0, 0, cov12 = rhovec),
     xlab = expression(rho), lwd = 2,
     type = "l", col = "blue", las = 1)
abline(v = 0, h = 0.25, col = "gray", lty = "dashed") 
## End(Not run)

Bivariate Normal Distribution Family Function

Description

Maximum likelihood estimation of the five parameters of a bivariate normal distribution.

Usage

binormal(lmean1 = "identitylink", lmean2 = "identitylink",
         lsd1   = "loglink",     lsd2   = "loglink",
         lrho   = "rhobitlink",
         imean1 = NULL,       imean2 = NULL,
         isd1   = NULL,       isd2   = NULL,
         irho   = NULL,       imethod = 1,
         eq.mean = FALSE,     eq.sd = FALSE,
         zero = c("sd", "rho"), rho.arg = NA)

Arguments

Details

For the bivariate normal distribution, this fits a linear model (LM) to the means, and by default, the other parameters are intercept-only. The response should be a two-column matrix. The correlation parameter is rho, which lies between -1 and 1 (thus the rhobitlink link is a reasonable choice). The fitted means are returned as the fitted values, which is in the form of a two-column matrix. Fisher scoring is implemented.

Value

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

Warning

This function may be renamed to normal2() or something like that at a later date.

Note

If both equal means and equal standard deviations are desired then use something like constraints = list("(Intercept)" = matrix(c(1,1,0,0,0, 0,0,1,1,0 ,0,0,0,0,1), 5, 3)) and maybe zero = NULL etc.

Author(s)

T. W. Yee

See Also

uninormal, trinormal, pbinorm, bistudentt, rhobitlink.

Examples

set.seed(123); nn <- 1000
bdata <- data.frame(x2 = runif(nn), x3 = runif(nn))
bdata <- transform(bdata, y1 = rnorm(nn, 1 + 2 * x2),
                          y2 = rnorm(nn, 3 + 4 * x2))
fit1 <- vglm(cbind(y1, y2) ~ x2,
             binormal(eq.sd = TRUE), bdata, trace = TRUE)
coef(fit1, matrix = TRUE)
constraints(fit1)
summary(fit1)

# Estimated P(Y1 <= y1, Y2 <= y2) under the fitted model
var1  <- loglink(2 * predict(fit1)[, "loglink(sd1)"], inv = TRUE)
var2  <- loglink(2 * predict(fit1)[, "loglink(sd2)"], inv = TRUE)
cov12 <- rhobitlink(predict(fit1)[, "rhobitlink(rho)"], inv = TRUE)
head(with(bdata, pbinorm(y1, y2,
                         mean1 = predict(fit1)[, "mean1"],
                         mean2 = predict(fit1)[, "mean2"],
                         var1 = var1, var2 = var2, cov12 = cov12)))

Gaussian Copula (Bivariate) Family Function

Description

Estimate the correlation parameter of the (bivariate) Gaussian copula distribution by maximum likelihood estimation.

Usage

binormalcop(lrho = "rhobitlink", irho = NULL,
            imethod = 1, parallel = FALSE, zero = NULL)

Arguments

Details

The cumulative distribution function is

P(Y_1 \leq y_1, Y_2 \leq y_2) = \Phi_2 ( \Phi^{-1}(y_1), \Phi^{-1}(y_2); \rho )

for -1 < \rho < 1, \Phi_2 is the cumulative distribution function of a standard bivariate normal (see pbinorm), and \Phi is the cumulative distribution function of a standard univariate normal (see pnorm).

The support of the function is the interior of the unit square; however, values of 0 and/or 1 are not allowed. The marginal distributions are the standard uniform distributions. When \rho = 0 the random variables are independent.

This VGAM family function can handle multiple responses, for example, a six-column matrix where the first 2 columns is the first out of three responses, the next 2 columns being the next response, etc.

Value

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

Note

The response matrix must have a multiple of two-columns. Currently, the fitted value is a matrix with the same number of columns and values equal to 0.5. This is because each marginal distribution corresponds to a standard uniform distribution.

This VGAM family function is fragile; each response must be in the interior of the unit square. Setting crit = "coef" is sometimes a good idea because inaccuracies in pbinorm might mean unnecessary half-stepping will occur near the solution.

Author(s)

T. W. Yee

References

Schepsmeier, U. and Stober, J. (2014). Derivatives and Fisher information of bivariate copulas. Statistical Papers 55, 525–542.

See Also

rbinormcop, rhobitlink, pnorm, kendall.tau.

Examples

nn <- 1000
ymat <- rbinormcop(nn, rho = rhobitlink(-0.9, inverse = TRUE))
bdata <- data.frame(y1 = ymat[, 1], y2 = ymat[, 2],
                    y3 = ymat[, 1], y4 = ymat[, 2],
                    x2 = runif(nn))
summary(bdata)
## Not run:  plot(ymat, col = "blue") 
fit1 <-  # 2 responses, e.g., (y1,y2) is the 1st
  vglm(cbind(y1, y2, y3, y4) ~ 1, fam = binormalcop,
       crit = "coef",  # Sometimes a good idea
       data = bdata, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1)
head(fitted(fit1))
summary(fit1)

# Another example; rho is a linear function of x2
bdata <- transform(bdata, rho = -0.5 + x2)
ymat <- rbinormcop(n = nn, rho = with(bdata, rho))
bdata <- transform(bdata, y5 = ymat[, 1], y6 = ymat[, 2])
fit2 <- vgam(cbind(y5, y6) ~ s(x2), data = bdata,
             binormalcop(lrho = "identitylink"), trace = TRUE)
## Not run: plot(fit2, lcol = "blue", scol = "orange", se = TRUE)

Gaussian Copula (Bivariate) Distribution

Description

Density, distribution function, and random generation for the (one parameter) bivariate Gaussian copula distribution.

Usage

dbinormcop(x1, x2, rho = 0, log = FALSE)
pbinormcop(q1, q2, rho = 0)
rbinormcop(n, rho = 0)

Arguments

Details

See binormalcop, the VGAM family functions for estimating the parameter by maximum likelihood estimation, for the formula of the cumulative distribution function and other details.

Value

dbinormcop gives the density, pbinormcop gives the distribution function, and rbinormcop generates random deviates (a two-column matrix).

Note

Yettodo: allow x1 and/or x2 to have values 1, and to allow any values for x1 and/or x2 to be outside the unit square.

Author(s)

T. W. Yee

See Also

binormalcop, binormal.

Examples

## Not run:  edge <- 0.01  # A small positive value
N <- 101; x <- seq(edge, 1.0 - edge, len = N); Rho <- 0.7
ox <- expand.grid(x, x)
zedd <- dbinormcop(ox[, 1], ox[, 2], rho = Rho, log = TRUE)
contour(x, x, matrix(zedd, N, N), col = "blue", labcex = 1.5)
zedd <- pbinormcop(ox[, 1], ox[, 2], rho = Rho)
contour(x, x, matrix(zedd, N, N), col = "blue", labcex = 1.5)

## End(Not run)

Plackett's Bivariate Copula

Description

Density, distribution function, and random generation for the (one parameter) bivariate Plackett copula.

Usage

dbiplackcop(x1, x2, oratio, log = FALSE)
pbiplackcop(q1, q2, oratio)
rbiplackcop(n, oratio)

Arguments

Details

See biplackettcop, the VGAM family functions for estimating the parameter by maximum likelihood estimation, for the formula of the cumulative distribution function and other details.

Value

dbiplackcop gives the density, pbiplackcop gives the distribution function, and rbiplackcop generates random deviates (a two-column matrix).

Author(s)

T. W. Yee

References

Mardia, K. V. (1967). Some contributions to contingency-type distributions. Biometrika, 54, 235–249.

See Also

biplackettcop, bifrankcop.

Examples

## Not run:  N <- 101; oratio <- exp(1)
x <- seq(0.0, 1.0, len = N)
ox <- expand.grid(x, x)
zedd <- dbiplackcop(ox[, 1], ox[, 2], oratio = oratio)
contour(x, x, matrix(zedd, N, N), col = "blue")
zedd <- pbiplackcop(ox[, 1], ox[, 2], oratio = oratio)
contour(x, x, matrix(zedd, N, N), col = "blue")

plot(rr <- rbiplackcop(n = 3000, oratio = oratio))
par(mfrow = c(1, 2))
hist(rr[, 1])  # Should be uniform
hist(rr[, 2])  # Should be uniform

## End(Not run)

Plackett's Bivariate Copula Family Function

Description

Estimate the association parameter of Plackett's bivariate distribution (copula) by maximum likelihood estimation.

Usage

biplackettcop(link = "loglink", ioratio = NULL, imethod = 1,
              nsimEIM = 200)

Arguments

Details

The defining equation is

\psi = H \times (1-y_1-y_2+H) / ((y_1-H) \times (y_2-H))

where P(Y_1 \leq y_1, Y_2 \leq y_2) = H_{\psi}(y_1,y_2) is the cumulative distribution function. The density function is h_{\psi}(y_1,y_2) =

\psi [1 + (\psi-1)(y_1 + y_2 - 2 y_1 y_2) ] / \left( [1 + (\psi-1)(y_1 + y_2) ]^2 - 4 \psi (\psi-1) y_1 y_2 \right)^{3/2}

for \psi > 0. Some writers call \psi the cross product ratio but it is called the odds ratio here. The support of the function is the unit square. The marginal distributions here are the standard uniform although it is commonly generalized to other distributions.

If \psi = 1 then h_{\psi}(y_1,y_2) = y_1 y_2, i.e., independence. As the odds ratio tends to infinity one has y_1=y_2. As the odds ratio tends to 0 one has y_2=1-y_1.

Fisher scoring is implemented using rbiplackcop. Convergence is often quite slow.

Value

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

Note

The response must be a two-column matrix. Currently, the fitted value is a 2-column matrix with 0.5 values because the marginal distributions correspond to a standard uniform distribution.

Author(s)

T. W. Yee

References

Plackett, R. L. (1965). A class of bivariate distributions. Journal of the American Statistical Association, 60, 516–522.

See Also

rbiplackcop, bifrankcop.

Examples

## Not run: 
ymat <- rbiplackcop(n = 2000, oratio = exp(2))
plot(ymat, col = "blue")
fit <- vglm(ymat ~ 1, fam = biplackettcop, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
vcov(fit)
head(fitted(fit))
summary(fit)

## End(Not run)

Biplot of Constrained Regression Models

Description

biplot is a generic function applied to RR-VGLMs and QRR-VGLMs etc. These apply to rank-1 and rank-2 models of these only. For RR-VGLMs these plot the second latent variable scores against the first latent variable scores.

Methods

x

The object from which the latent variables are extracted and/or plotted.

Note

See lvplot which is very much related to biplots.

Birnbaum-Saunders Regression Family Function

Description

Estimates the shape and scale parameters of the Birnbaum-Saunders distribution by maximum likelihood estimation.

Usage

bisa(lscale = "loglink", lshape = "loglink", iscale = 1,
     ishape = NULL, imethod = 1, zero = "shape", nowarning = FALSE)

Arguments

Details

The (two-parameter) Birnbaum-Saunders distribution has a cumulative distribution function that can be written as

F(y;a,b) = \Phi[ \xi(y/b)/a]

where \Phi(\cdot) is the cumulative distribution function of a standard normal (see pnorm), \xi(t) = \sqrt{t} - 1 / \sqrt{t}, y > 0, a>0 is the shape parameter, b>0 is the scale parameter. The mean of Y (which is the fitted value) is b(1 + a^2/2). and the variance is a^2 b^2 (1 + \frac{5}{4}a^2). By default, \eta_1 = \log(a) and \eta_2 = \log(b) for this family function.

Note that a and b are orthogonal, i.e., the Fisher information matrix is diagonal. This family function implements Fisher scoring, and it is unnecessary to compute any integrals numerically.

Value

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

Author(s)

T. W. Yee

References

Lemonte, A. J. and Cribari-Neto, F. and Vasconcellos, K. L. P. (2007). Improved statistical inference for the two-parameter Birnbaum-Saunders distribution. Computational Statistics & Data Analysis, 51, 4656–4681.

Birnbaum, Z. W. and Saunders, S. C. (1969). A new family of life distributions. Journal of Applied Probability, 6, 319–327.

Birnbaum, Z. W. and Saunders, S. C. (1969). Estimation for a family of life distributions with applications to fatigue. Journal of Applied Probability, 6, 328–347.

Engelhardt, M. and Bain, L. J. and Wright, F. T. (1981). Inferences on the parameters of the Birnbaum-Saunders fatigue life distribution based on maximum likelihood estimation. Technometrics, 23, 251–256.

Johnson, N. L. and Kotz, S. and Balakrishnan, N. (1995). Continuous Univariate Distributions, 2nd edition, Volume 2, New York: Wiley.

See Also

pbisa, inv.gaussianff, CommonVGAMffArguments.

Examples

bdata1 <- data.frame(x2 = runif(nn <- 1000))
bdata1 <- transform(bdata1, shape = exp(-0.5 + x2),
                            scale = exp(1.5))
bdata1 <- transform(bdata1, y = rbisa(nn, scale, shape))
fit1 <- vglm(y ~ x2, bisa(zero = 1), data = bdata1, trace = TRUE)
coef(fit1, matrix = TRUE)

## Not run: 
bdata2 <- data.frame(shape = exp(-0.5), scale = exp(0.5))
bdata2 <- transform(bdata2, y = rbisa(nn, scale, shape))
fit <- vglm(y ~ 1, bisa, data = bdata2, trace = TRUE)
with(bdata2, hist(y, prob = TRUE, ylim = c(0, 0.5),
                  col = "lightblue"))
coef(fit, matrix = TRUE)
with(bdata2, mean(y))
head(fitted(fit))
x <- with(bdata2, seq(0, max(y), len = 200))
lines(dbisa(x, Coef(fit)[1], Coef(fit)[2]) ~ x, data = bdata2,
      col = "orange", lwd = 2) 
## End(Not run)

The Birnbaum-Saunders Distribution

Description

Density, distribution function, and random generation for the Birnbaum-Saunders distribution.

Usage

dbisa(x, scale = 1, shape, log = FALSE)
pbisa(q, scale = 1, shape, lower.tail = TRUE, log.p = FALSE)
qbisa(p, scale = 1, shape, lower.tail = TRUE, log.p = FALSE)
rbisa(n, scale = 1, shape)

Arguments

Details

The Birnbaum-Saunders distribution is a distribution which is used in survival analysis. See bisa, the VGAM family function for estimating the parameters, for more details.

Value

dbisa gives the density, pbisa gives the distribution function, and qbisa gives the quantile function, and rbisa generates random deviates.

Author(s)

T. W. Yee and Kai Huang

See Also

bisa.

Examples

## Not run: 
x <- seq(0, 6, len = 400)
plot(x, dbisa(x, shape = 1), type = "l", col = "blue",
     ylab = "Density", lwd = 2, ylim = c(0,1.3), lty = 3,
     main = "X ~ Birnbaum-Saunders(shape, scale = 1)")
lines(x, dbisa(x, shape = 2), col = "orange", lty = 2, lwd = 2)
lines(x, dbisa(x, shape = 0.5), col = "green", lty = 1, lwd = 2)
legend(x = 3, y = 0.9, legend = paste("shape  = ",c(0.5, 1,2)),
       col = c("green","blue","orange"), lty = 1:3, lwd = 2)

shape <- 1; x <- seq(0.0, 4, len = 401)
plot(x, dbisa(x, shape = shape), type = "l", col = "blue",
     main = "Blue is density, orange is the CDF", las = 1,
     sub = "Red lines are the 10,20,...,90 percentiles",
     ylab = "", ylim = 0:1)
abline(h = 0, col = "blue", lty = 2)
lines(x, pbisa(x, shape = shape), col = "orange")
probs <- seq(0.1, 0.9, by = 0.1)
Q <- qbisa(probs, shape = shape)
lines(Q, dbisa(Q, shape = shape), col = "red", lty = 3, type = "h")
pbisa(Q, shape = shape) - probs  # Should be all zero
abline(h = probs, col = "red", lty = 3)
lines(Q, pbisa(Q, shape = shape), col = "red", lty = 3, type = "h")

## End(Not run)

Bivariate Student-t Family Function

Description

Estimate the degrees of freedom and correlation parameters of the (bivariate) Student-t distribution by maximum likelihood estimation.

Usage

bistudentt(ldf = "logloglink", lrho = "rhobitlink",
           idf = NULL, irho = NULL, imethod = 1,
           parallel = FALSE, zero = "rho")

Arguments

Details

The density function is

f(y_1, y_2; \nu, \rho) = \frac{1}{2\pi\sqrt{1-\rho^2}} (1 + (y_1^2 + y_2^2 - 2\rho y_1 y_2) / (\nu (1-\rho^2)))^{-(\nu+2)/2}

for -1 < \rho < 1, and real y_1 and y_2.

This VGAM family function can handle multiple responses, for example, a six-column matrix where the first 2 columns is the first out of three responses, the next 2 columns being the next response, etc.

Value

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

Warning

The working weight matrices have not been fully checked.

Note

The response matrix must have a multiple of two-columns. Currently, the fitted value is a matrix with the same number of columns and values equal to 0.0.

Author(s)

T. W. Yee, with help from Thibault Vatter.

References

Schepsmeier, U. and Stober, J. (2014). Derivatives and Fisher information of bivariate copulas. Statistical Papers 55, 525–542.

See Also

dbistudentt, binormal, pt.

Examples

nn <- 1000
mydof <- logloglink(1, inverse = TRUE)
ymat <- cbind(rt(nn, df = mydof), rt(nn, df = mydof))
bdata <- data.frame(y1 = ymat[, 1], y2 = ymat[, 2],
                    y3 = ymat[, 1], y4 = ymat[, 2],
                    x2 = runif(nn))
summary(bdata)
## Not run:  plot(ymat, col = "blue") 
fit1 <-    # 2 responses, e.g., (y1,y2) is the 1st
  vglm(cbind(y1, y2, y3, y4) ~ 1,
       bistudentt,  # crit = "coef",  # Sometimes a good idea
       data = bdata, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1)
head(fitted(fit1))
summary(fit1)

Bivariate Student-t Distribution Density Function

Description

Density for the bivariate Student-t distribution.

Usage

dbistudentt(x1, x2, df, rho = 0, log = FALSE)

Arguments

Details

One can think of this function as an extension of dt to two dimensions. See bistudentt for more information.

Value

dbistudentt gives the density.

See Also

bistudentt, dt.

Examples

## Not run:  N <- 101; x <- seq(-4, 4, len = N); Rho <- 0.7
mydf <- 10; ox <- expand.grid(x, x)
zedd <- dbistudentt(ox[, 1], ox[, 2], df = mydf,
                    rho = Rho, log = TRUE)
contour(x, x, matrix(zedd, N, N), col = "blue", labcex = 1.5)

## End(Not run)

Body Mass Index of New Zealand Adults Data

Description

The body mass indexes and ages from an approximate random sample of 700 New Zealand adults.

Usage

data(bmi.nz)

Format

A data frame with 700 observations on the following 2 variables.

age

a numeric vector; their age (years).

BMI

a numeric vector; their body mass indexes, which is their weight divided by the square of their height (kg / m^2).

Details

They are a random sample from the Fletcher Challenge/Auckland Heart and Health survey conducted in the early 1990s.

There are some outliers in the data set.

A variable gender would be useful, and may be added later.

Source

Formerly the Clinical Trials Research Unit, University of Auckland, New Zealand.

References

MacMahon, S., Norton, R., Jackson, R., Mackie, M. J., Cheng, A., Vander Hoorn, S., Milne, A., McCulloch, A. (1995) Fletcher Challenge-University of Auckland Heart & Health Study: design and baseline findings. New Zealand Medical Journal, 108, 499–502.

Examples

## Not run:  with(bmi.nz, plot(age, BMI, col = "blue"))
fit <- vgam(BMI ~ s(age, df = c(2, 4, 2)), lms.yjn,
            data = bmi.nz, trace = TRUE)
qtplot(fit, pcol = "blue", tcol = "brown", lcol = "brown") 
## End(Not run)

Borel-Tanner Distribution Family Function

Description

Estimates the parameter of a Borel-Tanner distribution by maximum likelihood estimation.

Usage

borel.tanner(Qsize = 1, link = "logitlink", imethod = 1)

Arguments

Details

The Borel-Tanner distribution (Tanner, 1953) describes the distribution of the total number of customers served before a queue vanishes given a single queue with random arrival times of customers (at a constant rate r per unit time, and each customer taking a constant time b to be served). Initially the queue has Q people and the first one starts to be served. The two parameters appear in the density only in the form of the product rb, therefore we use a=rb, say, to denote the single parameter to be estimated. The density function is

f(y;a) = \frac{ Q }{(y-Q)!} y^{y-Q-1} a^{y-Q} \exp(-ay)

where y=Q,Q+1,Q+2,\ldots. The case Q=1 corresponds to the Borel distribution (Borel, 1942). For the Q=1 case it is necessary for 0 < a < 1 for the distribution to be proper. The Borel distribution is a basic Lagrangian distribution of the first kind. The Borel-Tanner distribution is an Q-fold convolution of the Borel distribution.

The mean is Q/(1-a) (returned as the fitted values) and the variance is Q a / (1-a)^3. The distribution has a very long tail unless a is small. Fisher scoring is implemented.

Value

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

Author(s)

T. W. Yee

References

Tanner, J. C. (1953). A problem of interference between two queues. Biometrika, 40, 58–69.

Borel, E. (1942). Sur l'emploi du theoreme de Bernoulli pour faciliter le calcul d'une infinite de coefficients. Application au probleme de l'attente a un guichet. Comptes Rendus, Academie des Sciences, Paris, Series A, 214, 452–456.

Johnson N. L., Kemp, A. W. and Kotz S. (2005). Univariate Discrete Distributions, 3rd edition, p.328. Hoboken, New Jersey: Wiley.

Consul, P. C. and Famoye, F. (2006). Lagrangian Probability Distributions, Boston, MA, USA: Birkhauser.

See Also

rbort, poissonff, felix.

Examples

bdata <- data.frame(y = rbort(n <- 200))
fit <- vglm(y ~ 1, borel.tanner, bdata, trace = TRUE, crit = "c")
coef(fit, matrix = TRUE)
Coef(fit)
summary(fit)

The Borel-Tanner Distribution

Description

Density and random generation for the Borel-Tanner distribution.

Usage

dbort(x, Qsize = 1, a = 0.5, log = FALSE)
rbort(n, Qsize = 1, a = 0.5)

Arguments

Details

See borel.tanner, the VGAM family function for estimating the parameter, for the formula of the probability density function and other details.

Value

dbort gives the density, rbort generates random deviates.

Warning

Looping is used for rbort, therefore values of a close to 1 will result in long (or infinite!) computational times. The default value of a is subjective.

Author(s)

T. W. Yee

See Also

borel.tanner.

Examples

## Not run:  qsize <- 1; a <- 0.5; x <- qsize:(qsize+10)
plot(x, dbort(x, qsize, a), type = "h", las = 1, col = "blue",
     ylab = paste("fbort(qsize=", qsize, ", a=", a, ")"),
     log = "y", main = "Borel-Tanner density function") 
## End(Not run)

Bradley Terry Model

Description

Fits a Bradley Terry model (intercept-only model) by maximum likelihood estimation.

Usage

brat(refgp = "last", refvalue = 1, ialpha = 1)

Arguments

Details

The Bradley Terry model involves M+1 competitors who either win or lose against each other (no draws/ties allowed in this implementation–see bratt if there are ties). The probability that Competitor i beats Competitor j is \alpha_i / (\alpha_i+\alpha_j), where all the \alphas are positive. Loosely, the \alphas can be thought of as the competitors' ‘abilities’. For identifiability, one of the \alpha_i is set to a known value refvalue, e.g., 1. By default, this function chooses the last competitor to have this reference value. The data can be represented in the form of a M+1 by M+1 matrix of counts, where winners are the rows and losers are the columns. However, this is not the way the data should be inputted (see below).

Excluding the reference value/group, this function chooses \log(\alpha_j) as the M linear predictors. The log link ensures that the \alphas are positive.

The Bradley Terry model can be fitted by logistic regression, but this approach is not taken here. The Bradley Terry model can be fitted with covariates, e.g., a home advantage variable, but unfortunately, this lies outside the VGLM theoretical framework and therefore cannot be handled with this code.

Value

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

Warning

Presently, the residuals are wrong, and the prior weights are not handled correctly. Ideally, the total number of counts should be the prior weights, after the response has been converted to proportions. This would make it similar to family functions such as multinomial and binomialff.

Note

The function Brat is useful for coercing a M+1 by M+1 matrix of counts into a one-row matrix suitable for brat. Diagonal elements are skipped, and the usual S order of c(a.matrix) of elements is used. There should be no missing values apart from the diagonal elements of the square matrix. The matrix should have winners as the rows, and losers as the columns. In general, the response should be a 1-row matrix with M(M+1) columns.

Only an intercept model is recommended with brat. It doesn't make sense really to include covariates because of the limited VGLM framework.

Notationally, note that the VGAM family function brat has M+1 contestants, while bratt has M contestants.

Author(s)

T. W. Yee

References

Agresti, A. (2013). Categorical Data Analysis, 3rd ed. Hoboken, NJ, USA: Wiley.

Stigler, S. (1994). Citation patterns in the journals of statistics and probability. Statistical Science, 9, 94–108.

The BradleyTerry2 package has more comprehensive capabilities than this function.

See Also

bratt, Brat, multinomial, binomialff.

Examples

# Citation statistics: being cited is a 'win'; citing is a 'loss'
journal <- c("Biometrika", "Comm.Statist", "JASA", "JRSS-B")
mat <- matrix(c( NA, 33, 320, 284,
                730, NA, 813, 276,
                498, 68,  NA, 325,
                221, 17, 142,  NA), 4, 4)
dimnames(mat) <- list(winner = journal, loser = journal)
fit <- vglm(Brat(mat) ~ 1, brat(refgp = 1), trace = TRUE)
fit <- vglm(Brat(mat) ~ 1, brat(refgp = 1), trace = TRUE, crit = "coef")
summary(fit)
c(0, coef(fit))  # Log-abilities (in order of "journal")
c(1, Coef(fit))  # Abilities (in order of "journal")
fitted(fit)     # Probabilities of winning in awkward form
(check <- InverseBrat(fitted(fit)))  # Probabilities of winning
check + t(check)  # Should be 1's in the off-diagonals

Inputting Data to fit a Bradley Terry Model

Description

Takes in a square matrix of counts and outputs them in a form that is accessible to the brat and bratt family functions.

Usage

Brat(mat, ties = 0 * mat, string = c(">", "=="), whitespace = FALSE)

Arguments

Details

In the VGAM package it is necessary for each matrix to be represented as a single row of data by brat and bratt. Hence the non-diagonal elements of the M+1 by M+1 matrix are concatenated into M(M+1) values (no ties), while if there are ties, the non-diagonal elements of the M by M matrix are concatenated into M(M-1) values.

Value

A matrix with 1 row and either M(M+1) or M(M-1) columns.

Note

This is a data preprocessing function for brat and bratt.

Yet to do: merge InverseBrat into brat.

Author(s)

T. W. Yee

References

Agresti, A. (2013). Categorical Data Analysis, 3rd ed. Hoboken, NJ, USA: Wiley.

See Also

brat, bratt, InverseBrat.

Examples

journal <- c("Biometrika", "Comm Statist", "JASA", "JRSS-B")
mat <- matrix(c( NA, 33, 320, 284,   730, NA, 813, 276,
                498, 68,  NA, 325,   221, 17, 142, NA), 4, 4)
dimnames(mat) <- list(winner = journal, loser = journal)
Brat(mat)  # Less readable
Brat(mat, whitespace = TRUE)  # More readable
vglm(Brat(mat, whitespace = TRUE) ~ 1, brat, trace = TRUE)

Bradley Terry Model With Ties

Description

Fits a Bradley Terry model with ties (intercept-only model) by maximum likelihood estimation.

Usage

bratt(refgp = "last", refvalue = 1, ialpha = 1, i0 = 0.01)

Arguments

Details

There are several models that extend the ordinary Bradley Terry model to handle ties. This family function implements one of these models. It involves M competitors who either win or lose or tie against each other. (If there are no draws/ties then use brat). The probability that Competitor i beats Competitor j is \alpha_i / (\alpha_i+\alpha_j+\alpha_0), where all the \alphas are positive. The probability that Competitor i ties with Competitor j is \alpha_0 / (\alpha_i+\alpha_j+\alpha_0). Loosely, the \alphas can be thought of as the competitors' ‘abilities’, and \alpha_0 is an added parameter to model ties. For identifiability, one of the \alpha_i is set to a known value refvalue, e.g., 1. By default, this function chooses the last competitor to have this reference value. The data can be represented in the form of a M by M matrix of counts, where winners are the rows and losers are the columns. However, this is not the way the data should be inputted (see below).

Excluding the reference value/group, this function chooses \log(\alpha_j) as the first M-1 linear predictors. The log link ensures that the \alphas are positive. The last linear predictor is \log(\alpha_0).

The Bradley Terry model can be fitted with covariates, e.g., a home advantage variable, but unfortunately, this lies outside the VGLM theoretical framework and therefore cannot be handled with this code.

Value

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

Note

The function Brat is useful for coercing a M by M matrix of counts into a one-row matrix suitable for bratt. Diagonal elements are skipped, and the usual S order of c(a.matrix) of elements is used. There should be no missing values apart from the diagonal elements of the square matrix. The matrix should have winners as the rows, and losers as the columns. In general, the response should be a matrix with M(M-1) columns.

Also, a symmetric matrix of ties should be passed into Brat. The diagonal of this matrix should be all NAs.

Only an intercept model is recommended with bratt. It doesn't make sense really to include covariates because of the limited VGLM framework.

Notationally, note that the VGAM family function brat has M+1 contestants, while bratt has M contestants.

Author(s)

T. W. Yee

References

Torsney, B. (2004). Fitting Bradley Terry models using a multiplicative algorithm. In: Antoch, J. (ed.) Proceedings in Computational Statistics COMPSTAT 2004, Physica-Verlag: Heidelberg. Pages 513–526.

See Also

brat, Brat, binomialff.

Examples

# citation statistics: being cited is a 'win'; citing is a 'loss'
journal <- c("Biometrika", "Comm.Statist", "JASA", "JRSS-B")
mat <- matrix(c( NA, 33, 320, 284,
                730, NA, 813, 276,
                498, 68,  NA, 325,
                221, 17, 142,  NA), 4, 4)
dimnames(mat) <- list(winner = journal, loser = journal)

# Add some ties. This is fictitional data.
ties <- 5 + 0 * mat
ties[2, 1] <- ties[1,2] <- 9

# Now fit the model
fit <- vglm(Brat(mat, ties) ~ 1, bratt(refgp = 1), trace = TRUE,
            crit = "coef")

summary(fit)
c(0, coef(fit))  # Log-abilities (last is log(alpha0))
c(1, Coef(fit))  #     Abilities (last is alpha0)

fit@misc$alpha   # alpha_1,...,alpha_M
fit@misc$alpha0  # alpha_0

fitted(fit)  # Probabilities of winning and tying, in awkward form
predict(fit)
(check <- InverseBrat(fitted(fit)))    # Probabilities of winning
qprob <- attr(fitted(fit), "probtie")  # Probabilities of a tie
qprobmat <- InverseBrat(c(qprob), NCo = nrow(ties))  # Pr(tie)
check + t(check) + qprobmat  # Should be 1s in the off-diagonals

Western Spuce Budworm

Description

Counts of western spuce budworm (Choristoneura freemani) across seven developmental stages (five larval instars, pupae, and adults) on 12 sampling occasions.

Usage

data(budworm)

Format

A data frame with the following variables.

ddeg

Degree days.

total

Sum of stages 1–7.

stage1, stage2, stage3, stage4

Successive stages.

stage5, stage6, stage7

Successive stages.

Details

This data concerns the development of a defoliating moth widespread in western North America (i.e., north of Mexico). According to Boersch-Supan (2021), the insect passes hrough successive stages j=1,\ldots,r, delimited by r-1 moults. The data was originally used in a 1986 publication but has been corrected for two sampling occasions; the data appears in Candy (1990) and was analyzed in Boersch-Supan (2021). See the latter for more references.

Source

Candy, S. G. (1990). Biology of the mountain pinhole borer, Platypus subgranosus Scheld, in Tasmania. MA thesis, University of Tasmania, Australia. https://eprints.utas.edu.au/18864/.

References

Boersch-Supan, P. H. (2021). Modeling insect phenology using ordinal regression and continuation ratio models. ReScience C, 7.1, 1–14.

Examples

budworm
summary(budworm)

Model Calibrations

Description

calibrate is a generic function used to produce calibrations from various model fitting functions. The function invokes particular ‘methods’ which depend on the ‘class’ of the first argument.

Usage

calibrate(object, ...)

Arguments

Details

Given a regression model with explanatory variables X and response Y, calibration involves estimating X from Y using the regression model. It can be loosely thought of as the opposite of predict (which takes an X and returns a Y of some sort.) In general, the central algorithm is maximum likelihood calibration.

Value

In general, given a new response Y, some function of the explanatory variables X are returned. For example, for constrained ordination models such as CQO and CAO models, it is usually not possible to return X, so the latent variables are returned instead (they are linear combinations of the X). See the specific calibrate methods functions to see what they return.

Note

This function was not called predictx because of the inability of constrained ordination models to return X; they can only return the latent variable values (also known as site scores) instead.

Author(s)

T. W. Yee

References

ter Braak, C. J. F. and van Dam, H. (1989). Inferring pH from diatoms: a comparison of old and new calibration methods. Hydrobiologia, 178, 209–223.

See Also

predict, calibrate.rrvglm, calibrate.qrrvglm.

Examples

## Not run: 
hspider[, 1:6] <- scale(hspider[, 1:6])  # Stdzed environmental vars
set.seed(123)
pcao1 <- cao(cbind(Pardlugu, Pardmont, Pardnigr, Pardpull, Zoraspin) ~
         WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
         family = poissonff, data = hspider, Rank = 1, Bestof = 3,
         df1.nl = c(Zoraspin = 2, 1.9), Crow1positive = TRUE)

siteNos <- 1:2  # Calibrate these sites
cpcao1 <- calibrate(pcao1, trace = TRUE,
                    newdata = data.frame(depvar(pcao1)[siteNos, ],
                                         model.matrix(pcao1)[siteNos, ]))

# Graphically compare the actual site scores with their calibrated values
persp(pcao1, main = "Site scores: solid=actual, dashed=calibrated",
      label = TRUE, col = "blue", las = 1)
abline(v = latvar(pcao1)[siteNos], col = seq(siteNos))  # Actual scores
abline(v = cpcao1, lty = 2, col = seq(siteNos))  # Calibrated values

## End(Not run)

Calibration for Constrained Regression Models

Description

calibrate is a generic function applied to RR-VGLMs, QRR-VGLMs and RR-VGAMs, etc.

Methods

object

The object from which the calibration is performed.

Calibration for CQO and CAO models

Description

Performs maximum likelihood calibration for constrained quadratic and additive ordination models (CQO and CAO models are better known as QRR-VGLMs and RR-VGAMs respectively).

Usage

calibrate.qrrvglm(object, newdata = NULL,
    type = c("latvar", "predictors", "response", "vcov", "everything"),
    lr.confint = FALSE, cf.confint = FALSE,
    level = 0.95, initial.vals = NULL, ...)

Arguments

Details

Given a fitted regression CQO/CAO model, maximum likelihood calibration is theoretically easy and elegant. However, the method assumes that all the responses are independent, which is often not true in practice. More details and references are given in Yee (2018) and ch.6 of Yee (2015).

The function optim is used to search for the maximum likelihood solution. Good initial values are needed, and arguments in calibrate.qrrvglm.control allows the user some control over the choice of these.

Value

Several methods are implemented to obtain confidence intervals/regions for the calibration estimates. One method is when lr.confint = TRUE, then a 4-column matrix is returned with the confidence limits being the final 2 columns (if type = "everything" then the matrix is returned in the lr.confint list component). Another similar method is when cf.confint = TRUE. There may be some redundancy in whatever is returned. Other methods are returned by using type and they are described as follows.

The argument type determines what is returned. If type = "everything" then all the type values are returned in a list, with the following components. Each component has length nrow(newdata).

Warning

This function is computationally expensive. Setting trace = TRUE to get a running log can be a good idea. This function has been tested but not extensively.

Note

Despite the name of this function, CAO models are handled as well to a certain extent. Some combinations of parameters are not handled, e.g., lr.confint = TRUE only works for rank-1, type = "vcov" only works for binomialff and poissonff models with canonical links and noRRR = ~ 1, and higher-order rank models need eq.tolerances = TRUE or I.tolerances = TRUE as well. For rank-1 objects, lr.confint = TRUE is recommended above type = "vcov" in terms of accuracy and overall generality. For class "qrrvglm" objects it is necessary that all response' tolerance matrices are positive-definite which correspond to bell-shaped response curves/surfaces.

For binomialff and poissonff models the deviance slot is used for the optimization rather than the loglikelihood slot, therefore one can calibrate using real-valued responses. (If the loglikelihood slot were used then functions such as dpois would be used with log = TRUE and then would be restricted to feed in integer-valued response values.)

Maximum likelihood calibration for Gaussian logit regression models may be performed by rioja but this applies to a single environmental variable such as pH in data("SWAP", package = "rioja"). In VGAM calibrate() estimates values of the latent variable rather than individual explanatory variables, hence the setting is more on ordination.

Author(s)

T. W. Yee. Recent work on the standard errors by David Zucker and Sam Oman at HUJI is gratefully acknowledged—these are returned in the vcov component and provided inspiration for lr.confint and cf.confint. A joint publication is being prepared on this subject.

References

Abate, J. and Whitt, W. (1992). The Fourier-series method for inverting transforms of probability distributions. Queueing Systems, 10, 5–88.

ter Braak, C. J. F. (1995). Calibration. In: Data Analysis in Community and Landscape Ecology by Jongman, R. H. G., ter Braak, C. J. F. and van Tongeren, O. F. R. (Eds.) Cambridge University Press, Cambridge.

See Also

calibrate.qrrvglm.control, calibrate.rrvglm, calibrate, cqo, cao, optim, uniroot.

Examples

## Not run: 
hspider[, 1:6] <- scale(hspider[, 1:6])  # Stdze environmental variables
set.seed(123)
siteNos <- c(1, 5)  # Calibrate these sites
pet1 <- cqo(cbind(Pardlugu, Pardmont, Pardnigr, Pardpull, Zoraspin) ~
        WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
        trace = FALSE,
        data = hspider[-siteNos, ],  # Sites not in fitted model
        family = poissonff, I.toler = TRUE, Crow1positive = TRUE)
y0 <- hspider[siteNos, colnames(depvar(pet1))]  # Species counts
(cpet1 <- calibrate(pet1, trace = TRUE, newdata = data.frame(y0)))
(clrpet1 <- calibrate(pet1, lr.confint = TRUE, newdata = data.frame(y0)))
(ccfpet1 <- calibrate(pet1, cf.confint = TRUE, newdata = data.frame(y0)))
(cp1wald <- calibrate(pet1, newdata = y0, type = "everything"))

## End(Not run)

## Not run: 
# Graphically compare the actual site scores with their calibrated
# values. 95 percent likelihood-based confidence intervals in green.
persp(pet1, main = "Site scores: solid=actual, dashed=calibrated",
      label = TRUE, col = "gray50", las = 1)
# Actual site scores:
xvars <- rownames(concoef(pet1))  # Variables comprising the latvar
est.latvar <- as.matrix(hspider[siteNos, xvars]) %*% concoef(pet1)
abline(v = est.latvar, col = seq(siteNos))
abline(v = cpet1, lty = 2, col = seq(siteNos))  # Calibrated values
arrows(clrpet1[,  3], c(60, 60), clrpet1[,  4], c(60, 60),  # Add CIs
       length = 0.08, col = "orange", angle = 90, code = 3, lwd = 2)
arrows(ccfpet1[,  3], c(70, 70), ccfpet1[,  4], c(70, 70),  # Add CIs
       length = 0.08, col = "limegreen", angle = 90, code = 3, lwd = 2)
arrows(cp1wald$latvar - 1.96 * sqrt(cp1wald$vcov), c(65, 65),
       cp1wald$latvar + 1.96 * sqrt(cp1wald$vcov), c(65, 65),  # Wald CIs
       length = 0.08, col = "blue", angle = 90, code = 3, lwd = 2)
legend("topright", lwd = 2,
       leg = c("CF interval", "Wald  interval", "LR interval"),
       col = c("limegreen", "blue", "orange"), lty = 1)

## End(Not run)

Control Function for CQO/CAO Calibration

Description

Algorithmic constants and parameters for running calibrate.qrrvglm are set using this function.

Usage

calibrate.qrrvglm.control(object, trace = FALSE, method.optim = "BFGS",
    gridSize = ifelse(Rank == 1, 21, 9), varI.latvar = FALSE, ...)

Arguments

Details

Most CQO/CAO users will only need to make use of trace and gridSize. These arguments should be used inside their call to calibrate.qrrvglm, not this function directly.

Value

A list which with the following components.

Note

Despite the name of this function, CAO models are handled as well.

References

Yee, T. W. (2020). On constrained and unconstrained quadratic ordination. Manuscript in preparation.

See Also

calibrate.qrrvglm, Coef.qrrvglm.

Examples

## Not run:  hspider[, 1:6] <- scale(hspider[, 1:6])  # Needed for I.tol=TRUE
set.seed(123)
p1 <- cqo(cbind(Alopacce, Alopcune, Pardlugu, Pardnigr,
                Pardpull, Trocterr, Zoraspin) ~
          WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
          family = poissonff, data = hspider, I.tol = TRUE)
sort(deviance(p1, history = TRUE))  # A history of all the iterations
siteNos <- 3:4  # Calibrate these sites
cp1 <- calibrate(p1, trace = TRUE,
                 new = data.frame(depvar(p1)[siteNos, ]))

## End(Not run)
## Not run: 
# Graphically compare the actual site scores with their calibrated values
persp(p1, main = "Site scores: solid=actual, dashed=calibrated",
      label = TRUE, col = "blue", las = 1)
abline(v = latvar(p1)[siteNos], col = seq(siteNos))  # Actual site scores
abline(v = cp1, lty = 2, col = seq(siteNos))  # Calibrated values

## End(Not run)

Calibration for CLO models (RR-VGLMs)

Description

Performs maximum likelihood calibration for constrained linear ordination models (CLO models are better known as RR-VGLMs).

Usage

calibrate.rrvglm(object, newdata = NULL,
    type = c("latvar", "predictors", "response", "vcov", "everything"),
    lr.confint = FALSE, cf.confint = FALSE,
    level = 0.95, initial.vals = NULL, ...)

Arguments

Details

Given a fitted regression CLO model, maximum likelihood calibration is theoretically easy and elegant. However, the method assumes that all responses are independent. More details and references are given in Yee (2015).

Calibration requires grouped or non-sparse data as the response. For example, if the family function is multinomial then one cannot usually calibrate y0 if it is a vector of 0s except for one 1. Instead, the response vector should be from grouped data so that there are few 0s. Indeed, it is found empirically that the stereotype model (also known as a reduced-rank multinomial logit model) calibrates well only with grouped data, and if the response vector is all 0s except for one 1 then the MLE will probably be at -Inf or +Inf. As another example, if the family function is poissonff then y0 must not be a vector of all 0s; instead, the response vector should have few 0s ideally. In general, you can use simulation to see what type of data calibrates acceptably.

Internally, this function is a simplification of calibrate.qrrvglm and users should look at that function for details. Good initial values are needed, and a grid is constructed to obtain these. The function calibrate.rrvglm.control allows the user some control over the choice of these.

Value

See calibrate.qrrvglm. Of course, the quadratic term in the latent variables vanishes for RR-VGLMs, so the model is simpler.

Warning

See calibrate.qrrvglm.

Note

See calibrate.qrrvglm about, e.g., calibration using real-valued responses.

Author(s)

T. W. Yee

See Also

calibrate.qrrvglm, calibrate, rrvglm, weightsvglm, optim, uniroot.

Examples

## Not run:   # Example 1
nona.xs.nz <- na.omit(xs.nz)  # Overkill!! (Data in VGAMdata package)
nona.xs.nz$dmd     <- with(nona.xs.nz, round(drinkmaxday))
nona.xs.nz$feethr  <- with(nona.xs.nz, round(feethour))
nona.xs.nz$sleephr <- with(nona.xs.nz, round(sleep))
nona.xs.nz$beats   <- with(nona.xs.nz, round(pulse))

p2 <- rrvglm(cbind(dmd, feethr, sleephr, beats) ~ age + smokenow +
  depressed + embarrassed + fedup + hurt + miserable +  # 11 psychological
  nofriend + moody + nervous + tense + worry + worrier, # variables
  noRRR = ~ age + smokenow, trace = FALSE, poissonff, data = nona.xs.nz,
  Rank = 2)
cp2 <- calibrate(p2, newdata = head(nona.xs.nz, 9), trace = TRUE)
cp2

two.cases <- nona.xs.nz[1:2, ]  # Another calibration example
two.cases$dmd       <- c(4, 10)
two.cases$feethr    <- c(4, 7)
two.cases$sleephr   <- c(7, 8)
two.cases$beats     <- c(62, 71)
(cp2b <- calibrate(p2, newdata = two.cases))

# Example 2
p1 <- rrvglm(cbind(dmd, feethr, sleephr, beats) ~ age + smokenow +
  depressed + embarrassed + fedup + hurt + miserable +  # 11 psychological
  nofriend + moody + nervous + tense + worry + worrier, # variables
  noRRR = ~ age + smokenow, trace = FALSE, poissonff, data = nona.xs.nz,
  Rank = 1)
(cp1c <- calibrate(p1, newdata = two.cases, lr.confint = TRUE))

## End(Not run)

Control Function for CLO (RR-VGLM) Calibration

Description

Algorithmic constants and parameters for running calibrate.rrvglm are set using this function.

Usage

calibrate.rrvglm.control(object, trace = FALSE, method.optim = "BFGS",
    gridSize = ifelse(Rank == 1, 17, 9), ...)

Arguments

Details

Most CLO users will only need to make use of trace and gridSize. These arguments should be used inside their call to calibrate.rrvglm, not this function directly.

Value

Similar to calibrate.qrrvglm.control.

See Also

calibrate.rrvglm, Coef.rrvglm.

Fitting Constrained Additive Ordination (CAO)

Description

A constrained additive ordination (CAO) model is fitted using the reduced-rank vector generalized additive model (RR-VGAM) framework.

Usage

cao(formula, family = stop("argument 'family' needs to be assigned"),
    data = list(),
    weights = NULL, subset = NULL, na.action = na.fail,
    etastart = NULL, mustart = NULL, coefstart = NULL,
    control = cao.control(...), offset = NULL,
    method = "cao.fit", model = FALSE, x.arg = TRUE, y.arg = TRUE,
    contrasts = NULL, constraints = NULL,
    extra = NULL, qr.arg = FALSE, smart = TRUE, ...)

Arguments

Details

The arguments of cao are a mixture of those from vgam and cqo, but with some extras in cao.control. Currently, not all of the arguments work properly.

CAO can be loosely be thought of as the result of fitting generalized additive models (GAMs) to several responses (e.g., species) against a very small number of latent variables. Each latent variable is a linear combination of the explanatory variables; the coefficients C (called C below) are called constrained coefficients or canonical coefficients, and are interpreted as weights or loadings. The C are estimated by maximum likelihood estimation. It is often a good idea to apply scale to each explanatory variable first.

For each response (e.g., species), each latent variable is smoothed by a cubic smoothing spline, thus CAO is data-driven. If each smooth were a quadratic then CAO would simplify to constrained quadratic ordination (CQO; formerly called canonical Gaussian ordination or CGO). If each smooth were linear then CAO would simplify to constrained linear ordination (CLO). CLO can theoretically be fitted with cao by specifying df1.nl=0, however it is more efficient to use rrvglm.

Currently, only Rank=1 is implemented, and only noRRR = ~1 models are handled.

With binomial data, the default formula is

logit(P[Y_s=1]) = \eta_s = f_s(\nu), \ \ \ s=1,2,\ldots,S

where x_2 is a vector of environmental variables, and \nu=C^T x_2 is a R-vector of latent variables. The \eta_s is an additive predictor for species s, and it models the probabilities of presence as an additive model on the logit scale. The matrix C is estimated from the data, as well as the smooth functions f_s. The argument noRRR = ~ 1 specifies that the vector x_1, defined for RR-VGLMs and QRR-VGLMs, is simply a 1 for an intercept. Here, the intercept in the model is absorbed into the functions. A clogloglink link may be preferable over a logitlink link.

With Poisson count data, the formula is

\log(E[Y_s]) = \eta_s = f_s(\nu)

which models the mean response as an additive models on the log scale.

The fitted latent variables (site scores) are scaled to have unit variance. The concept of a tolerance is undefined for CAO models, but the optimums and maximums are defined. The generic functions Max and Opt should work for CAO objects, but note that if the maximum occurs at the boundary then Max will return a NA. Inference for CAO models is currently undeveloped.

Value

An object of class "cao" (this may change to "rrvgam" in the future). Several generic functions can be applied to the object, e.g., Coef, concoef, lvplot, summary.

Warning

CAO is very costly to compute. With version 0.7-8 it took 28 minutes on a fast machine. I hope to look at ways of speeding things up in the future.

Use set.seed just prior to calling cao() to make your results reproducible. The reason for this is finding the optimal CAO model presents a difficult optimization problem, partly because the log-likelihood function contains many local solutions. To obtain the (global) solution the user is advised to try many initial values. This can be done by setting Bestof some appropriate value (see cao.control). Trying many initial values becomes progressively more important as the nonlinear degrees of freedom of the smooths increase.

Note

CAO models are computationally expensive, therefore setting trace = TRUE is a good idea, as well as running it on a simple random sample of the data set instead.

Sometimes the IRLS algorithm does not converge within the FORTRAN code. This results in warnings being issued. In particular, if an error code of 3 is issued, then this indicates the IRLS algorithm has not converged. One possible remedy is to increase or decrease the nonlinear degrees of freedom so that the curves become more or less flexible, respectively.