invgamma2mr: 2 - parameter Inverse Gamma Distribution

Description

Estimates the 2-parameter Inverse Gamma distribution by maximum likelihood estimation.

Usage

invgamma2mr(lmu      = "loglink", 
              lshape   = logofflink(offset = -2), 
              parallel = FALSE, 
              ishape   = NULL, 
              imethod  = 1, 
              zero     = "shape")

Value

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

and vgam.

Arguments

lmu, lshape: Link functions applied to the (positives) mu and shape parameters (called $\mu$ and $a$ respectively), according to gamma2. See CommonVGAMffArguments for further information.
parallel: Same as gamma2. Details at CommonVGAMffArguments.
ishape: Optional initial value for shape, same as gamma2
imethod: Same as gamma2.
zero: Numeric or character vector. Position or name(s) of the parameters/linear predictors to be modeled as intercept--only. Default is "shape". Details at CommonVGAMffArguments.

Warning

Note that zero can be a numeric or a character vector specifying the position of the names (partially or not) of the linear predictor modeled as intercept only. In this family function such names are

c("mu", "shape").

Numeric values can be entered as usual. See CommonVGAMffArguments for further details.

Author

Victor Miranda and T. W. Yee

Details

The Gamma distribution and the Inverse Gamma distribution are related as follows:Let X be a random variable distributed as $Gamma (a, \beta)$, where $a > 0$ denotes the shape parameter and $\beta > 0$ is the scale paramater. Then $Y = 1/X$ is an Inverse Gamma random variable with parameters scale = $a$ and shape = $1/\beta$.

The Inverse Gamma density function is given by

$$f(y;\mu, a) = \frac{(a - 1)^{a} \mu^{a}}{\Gamma(a)}y^{-a- 1} \ e^{-\mu(a - 1)/y},$$

for $\mu > 0$, $a > 0$ and $y > 0$. Here, $\Gamma(\cdot)$ is the gamma function, as in gamma. The mean of Y is $\mu=\mu$ (returned as the fitted values) with variance $\sigma^2 = \mu^2 / (a - 2)$ if $a > 2$, else is infinite. Thus, the link function for the shape parameter is logloglink. Then, by default, the two linear/additive predictors are $\eta_1=\log(\mu)$, and $\eta_2=\log(a)$, i.e in the VGLM context, $\eta = (log(\mu), loglog(a)$

This VGAM family function handles multiple reponses by implementing Fisher scoring and unlike gamma2, the working-weight matrices are not diagonal. The Inverse Gamma distribution is right-skewed and either for small values of $a$ (plus modest $\mu$) or very large values of $\mu$ (plus moderate $a > 2$), the density has values too close to zero.

References

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

Examples

Run this code

#------------------------------------------------------------------------#
# Essentially fitting a 2-parameter inverse gamma distribution
# with 2 responses.

set.seed(101)
y1 = rinvgamma(n = 500, scale = exp(2.0), shape = exp(2.0))
y2 = rinvgamma(n = 500, scale = exp(2.5), shape = exp(2.5))
gdata <- data.frame(y1, y2)

fit1 <- vglm(cbind(y1, y2) ~ 1, 
            family = invgamma2mr(zero = NULL, 
            
                                 # OPTIONAL INITIAL VALUE
                                 # ishape = exp(2),
                                 
                                 imethod = 1),
            data = gdata, trace = TRUE)

Coef(fit1)
c(Coef(fit1), log(mean(gdata$y1)), log(mean(gdata$y2)))
summary(fit1)
vcov(fit1, untransform = TRUE)

#------------------------------------------------------------------------#
# An example including one covariate.
# Note that the x2 affects the shape parameter, which implies that both,
# 'mu' and 'shape' are affected.
# Consequently, zero must be set as NULL !

x2    <- runif(1000)
gdata <- data.frame(y3 = rinvgamma(n = 1000, 
                                   scale = exp(2.0), 
                                   shape = exp(2.0 + x2)))

fit2 <- vglm(y3 ~ x2, 
            family = invgamma2mr(lshape = "loglink", zero = NULL), 
            data = gdata, trace = TRUE)

coef(fit2, matrix = TRUE)
summary(fit2)
vcov(fit2)

Run the code above in your browser using DataLab