Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/family.bivariate.R

Estimate the four parameters of the Freund (1961) bivariate extension of the exponential distribution by maximum likelihood estimation.

1 2 3 |

`la, lap, lb, lbp` |
Link functions applied to the (positive)
parameters |

`ia, iap, ib, ibp` |
Initial value for the four parameters respectively. The default is to estimate them all internally. |

`independent` |
Logical. If |

`zero` |
A vector specifying which
linear/additive predictors are modelled as intercepts only.
The values can be from the set {1,2,3,4}.
The default is none of them.
See |

This model represents one type of bivariate extension of the exponential
distribution that is applicable to certain problems, in particular,
to two-component systems which can function if one of the components
has failed. For example, engine failures in two-engine planes, paired
organs such as peoples' eyes, ears and kidneys.
Suppose *y1* and *y2* are random variables
representing the lifetimes of two components *A* and *B*
in a two component system.
The dependence between *y1* and *y2*
is essentially such that the failure of the *B* component
changes the parameter of the exponential life distribution
of the *A* component from *alpha* to
*alpha'*, while the failure of the *A* component
changes the parameter of the exponential life distribution
of the *B* component from *beta* to
*beta'*.

The joint probability density function is given by

*
f(y1,y2) = alpha * beta' * exp(-beta' * y2 -
(alpha+beta-beta') * y1) *

for *0 < y1 < y2*, and

*
f(y1,y2) = beta * alpha' * exp(-alpha' * y1 -
(alpha+beta-alpha') * y2) *

for *0 < y2 < y1*.
Here, all four parameters are positive, as well as the responses
*y1* and *y2*.
Under this model, the probability that component *A*
is the first to fail is
*alpha/(alpha+beta)*.
The time to the first failure is distributed as an
exponential distribution with rate
*alpha+beta*. Furthermore, the
distribution of the time from first failure to failure
of the other component is a mixture of
Exponential(*alpha'*) and
Exponential(*beta'*) with proportions
*beta/(alpha+beta)*
and *alpha/(alpha+beta)*
respectively.

The marginal distributions are, in general, not exponential.
By default, the linear/additive predictors are
*eta1=log(alpha)*,
*eta2=log(alpha')*,
*eta3=log(beta)*,
*eta4=log(beta')*.

A special case is when *alpha=alpha'*
and *beta'=beta'*, which means that
*y1* and *y2* are independent, and
both have an ordinary exponential distribution with means
*1/alpha* and *1/beta*
respectively.

Fisher scoring is used, and the initial values correspond to the MLEs of an intercept model. Consequently, convergence may take only one iteration.

An object of class `"vglmff"`

(see `vglmff-class`

).
The object is used by modelling functions such as `vglm`

and `vgam`

.

To estimate all four parameters, it is necessary to have some
data where *y1<y2* and *y2<y1*.

The response must be a two-column matrix, with columns
*y1* and *y2*.
Currently, the fitted value is a matrix with two columns; the
first column has values
*(alpha'+beta)/(alpha' * (alpha+beta))*
for the mean of *y1*,
while the second column has values
*(beta'+alpha)/(beta' * (alpha+beta))*
for the mean of *y2*.
The variance of *y1* is

*
[(alpha')^2 + 2 * alpha * beta + beta^2]/
[(alpha')^2 * (alpha + beta)^2], *

the variance of *y2* is

*
[(beta')^2 + 2 * alpha * beta + alpha^2]/
[(beta')^2 * (alpha + beta)^2], *

the covariance of *y1* and *y2* is

*
[alpha' * beta' - alpha * beta]/
[alpha' * beta' * (alpha + beta)^2]. *

T. W. Yee

Freund, J. E. (1961).
A bivariate extension of the exponential distribution.
*Journal of the American Statistical Association*,
**56**, 971–977.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | ```
fdata <- data.frame(y1 = rexp(nn <- 1000, rate = exp(1)))
fdata <- transform(fdata, y2 = rexp(nn, rate = exp(2)))
fit1 <- vglm(cbind(y1, y2) ~ 1, fam = freund61, data = fdata, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1)
vcov(fit1)
head(fitted(fit1))
summary(fit1)
# y1 and y2 are independent, so fit an independence model
fit2 <- vglm(cbind(y1, y2) ~ 1, freund61(indep = TRUE),
data = fdata, trace = TRUE)
coef(fit2, matrix = TRUE)
constraints(fit2)
pchisq(2 * (logLik(fit1) - logLik(fit2)), # p-value
df = df.residual(fit2) - df.residual(fit1), lower.tail = FALSE)
lrtest(fit1, fit2) # Better alternative
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.