VGAMextra (version 0.0-6)

weibullRff: Distribution--specified quantile regression: 2--parameter Weibull Distribution

Description

Estimates the 2--parameter Weibull distribution by maximum likelihood. An extension of weibullR from VGAM. Weibull quantile regression and Weibull--mean modelling are also handled via the first linear predictor.

Usage

     weibullRff(link1 = c("loglink", "weibullMlink", "weibullQlink")[1],
                lshape = "loglink", percentile = 50,
                imu = NULL, iscale = NULL, ishape = NULL,
                lss = TRUE, nrfs = 1, probs.y = c(0.2, 0.5, 0.8),
                imethod = 1, zero = "shape")

Value

An object of class "vglm". See vglm-class for full details.

Arguments

link1

Link function for the first linear predictor. Default is link1 = "loglink", mimicking weibullR. The other options are the 2--parameter weibullQlink, applied to the Weibull quantile function, and the 2--parameter weibullMlink, applied to the Weibull mean function. See below for more details.

percentile

Numeric. A vector with the percentiles of interest, between 0 and 100. Used only in Weibull quantile regression, that is, when link1 = "weibullQlink".

lshape,imu, iscale, ishape, lss, nrfs, probs.y, imethod

Same as weibullR.

zero

Specifies the parameters to be modelled as intercept--only. Further details below.

See CommonVGAMffArguments.

Author

V. Miranda and Thomas W. Yee.

Details

weibullRff is a modified version of weibullR adapted to handle weibullQlink and weibullMlink, two 2-parameter linear predictors that model the Weibull mean and quantiles respectively. The underlying density is the ordinary scale(\(\beta\)) & shape(\(\alpha\)) Weibull density (see weibullR).

The second linear predictor is always \(\eta_2 = \log~\alpha\). The argument link1 handles the first linear predictor.

** Mimicking weibullR **

The default is link1 = "loglink", i.e., \(\eta_1 = \log~\beta = \log~scale\), and \(\eta_2 = \log~\alpha = \log~shape\), as with weibullR. The mean (\(\mu\)) is returned as the fitted value.

** Weibull quantile regression **

For Weibull quantile regression set link1 = "weibullQlink" and enter a numeric vector of percentiles of interest via percentile. See examples.

NOTE: Enter the response using Q.reg. See example below. The Weibull quantiles are returned as the fitted values.

** Weibull-mean modelling **

For Weibull-mean modelling (viz. mean time to failure) set link1 = "weibullMlink". The mean (\(\mu\)) is returned as the fitted value.

References

Miranda & Yee (2021) Two--Parameter Link Functions, With Application to Negative Binomial, Weibull and Quantile Regression. In preparation.

See Also

Q.reg, weibullQlink, weibullMlink, weibullR, gamma, CommonVGAMffArguments.

Examples

if (FALSE) {
set.seed(18121)
nn <- 300
x2 <- sort(runif(nn, 0, 3))  # Predictor/covariate.
bb <- exp(1.1 + 0.2 * x2)    # Scale parameter as function of x2.
aa <- exp(1.0 - 0.35 * x2)     # Shape parameter as function of x2.
mymu <- bb * gamma(1 + 1/aa)  # The Weibull mean.

## Use weibullMlink to generate appropriate scale parameter.
newbb <- weibullMlink(theta = log(mymu), shape = aa, inverse = TRUE, deriv = 0)

## A single random response
wdata <- data.frame(y1 = rweibull(nn, shape = aa, scale = newbb), x2 = x2)

# Plotting the data / Histogram
plot(y1  ~ x2, xlim = c(0, 3.1), ylim = c(-1, 35),
     pch = 20, data = wdata, col = "black", 
     main = "Weibull Quantile regression~ x2")
abline(h = 0, v = 0, col = "grey", lty = "dashed")
with(wdata, hist(y1, col = "red", breaks = 15))

## Weibull regression - percentile = c(25, 50, 75)
## Note the use of Q.reg.
fit1 <- vglm(Q.reg(y1, length.arg = 3) ~ x2, 
             weibullRff(link1 = "weibullQlink", zero = NULL,
                                 percentile = c(25, 50, 75)), 
             trace = TRUE, data = wdata)
head(fitted(fit1))
summary(fit1)
my.coef3Q <- coef(fit1, mat = TRUE)

### Proportion of data below the estimated 25% Quantile line.
100 * (1 - (sum(wdat$y1 >= fitted(fit2)[, 1]) / nn))  # Around 25%
### Proportion of data below the estimated 50% Quantile line.
100 * (1 - (sum(wdat$y1 >= fitted(fit2)[, 2]) / nn))   # Around 50%
### Proportion of data below the estimated 75% Quantile line.
100 * (1 - ( sum(wdat$y1 >= fitted(fit2)[, 3]) / nn ))   # Around 75%

## The quantile plots ##
my.coef3Q <- coef(fit2, matrix = TRUE)
with(wdat, lines(x2, exp(my.coef3Q[1, 1] + my.coef3Q[2, 1] * x2), 
                    col = "red", lty = "dotted", lwd = 4))
with(wdat, lines(x2, exp(my.coef3Q[1, 3] + my.coef3Q[2, 3] * x2), 
                 col = "orange", lty = "dotted", lwd = 4))
with(wdat, lines(x2, exp(my.coef3Q[1, 5] + my.coef3Q[2, 5] * x2), 
                 col = "blue", lty = "dotted", lwd = 4))

## Adding the 'mean' or expected Weibull regression line.
fit2 <- vglm(y1 ~ x2, 
             weibullRff(link1 = "weibullMlink", zero = NULL), 
             trace = TRUE, data= wdat)
my.coef3Q <- coef(fit2, mat = TRUE)
with(wdat, lines(x2, exp(my.coef3Q[1, 1] + my.coef3Q[2, 1] * x2), 
                 col = "yellow", lty = "dashed", lwd = 3))


legend("topleft", c("25h Perc", "50th Perc", "Mean", "75th Perc"),
       col = c("red", "orange", "cyan", "blue"),
       lty = c("dashed", "dashed", "dashed", "dashed"), lwd = rep(4, 4))
     }