twophase {msm}R Documentation

Coxian phase-type distribution with two phases

Description

Density, distribution, quantile functions and other utilities for the Coxian phase-type distribution with two phases.

Usage

d2phase(x, l1, mu1, mu2, log = FALSE)

p2phase(q, l1, mu1, mu2, lower.tail = TRUE, log.p = FALSE)

q2phase(p, l1, mu1, mu2, lower.tail = TRUE, log.p = FALSE)

r2phase(n, l1, mu1, mu2)

h2phase(x, l1, mu1, mu2, log = FALSE)

Arguments

x, q

vector of quantiles.

l1

Intensity for transition between phase 1 and phase 2.

mu1

Intensity for transition from phase 1 to exit.

mu2

Intensity for transition from phase 2 to exit.

log

logical; if TRUE, return log density or log hazard.

lower.tail

logical; if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x].

log.p

logical; if TRUE, probabilities p are given as log(p).

p

vector of probabilities.

n

number of observations. If length(n) > 1, the length is taken to be the number required.

Details

This is the distribution of the time to reach state 3 in a continuous-time Markov model with three states and transitions permitted from state 1 to state 2 (with intensity λ1\lambda_1) state 1 to state 3 (intensity μ1\mu_1) and state 2 to state 3 (intensity μ2\mu_2). States 1 and 2 are the two "phases" and state 3 is the "exit" state.

The density is

f(tλ1,μ1)=e(λ1+μ1)t(μ1+(λ1+μ1)λ1t)f(t | \lambda_1, \mu_1) = e^{-(\lambda_1+\mu_1)t}(\mu_1 + (\lambda_1+\mu_1)\lambda_1 t)

if λ1+μ1=μ2\lambda_1 + \mu_1 = \mu_2, and

f(tλ1,μ1,μ2)=(λ1+μ1)e(λ1+μ1)t(μ2μ1)+μ2λ1eμ2tλ1+μ1μ2f(t | \lambda_1, \mu_1, \mu_2) = \frac{(\lambda_1+\mu_1)e^{-(\lambda_1+\mu_1)t}(\mu_2-\mu_1) + \mu_2\lambda_1e^{-\mu_2t}}{\lambda_1+\mu_1-\mu_2}

otherwise. The distribution function is

F(tλ1,μ1)=1e(λ1+μ1)t(1+λ1t)F(t | \lambda_1, \mu_1) = 1 - e^{-(\lambda_1+\mu_1) t} (1 + \lambda_1 t)

if λ1+μ1=μ2\lambda_1 + \mu_1 = \mu_2, and

F(tλ1,μ1,μ2)=1e(λ1+μ1)t(μ2μ1)+λ1eμ2tλ1+μ1μ2F(t | \lambda_1, \mu_1, \mu_2) = 1 - \frac{e^{-(\lambda_1 + \mu_1)t} (\mu_2 - \mu_1) + \lambda_1 e^{-\mu_2 t}}{ \lambda_1 + \mu_1 - \mu_2}

otherwise. Quantiles are calculated by numerically inverting the distribution function.

The mean is (1+λ1/μ2)/(λ1+μ1)(1 + \lambda_1/\mu_2) / (\lambda_1 + \mu_1).

The variance is (2+2λ1(λ1+μ1+μ2)/μ22(1+λ1/μ2)2)/(λ1+μ1)2(2 + 2\lambda_1(\lambda_1+\mu_1+ \mu_2)/\mu_2^2 - (1 + \lambda_1/\mu_2)^2)/(\lambda_1+\mu_1)^2.

If μ1=μ2\mu_1=\mu_2 it reduces to an exponential distribution with rate μ1\mu_1, and the parameter λ1\lambda_1 is redundant. Or also if λ1=0\lambda_1=0.

The hazard at x=0x=0 is μ1\mu_1, and smoothly increasing if μ1<μ2\mu_1<\mu_2. If λ1+μ1μ2\lambda_1 + \mu_1 \geq \mu_2 it increases to an asymptote of μ2\mu_2, and if λ1+μ1μ2\lambda_1 + \mu_1 \leq \mu_2 it increases to an asymptote of λ1+μ1\lambda_1 + \mu_1. The hazard is decreasing if μ1>μ2\mu_1>\mu_2, to an asymptote of μ2\mu_2.

Value

d2phase gives the density, p2phase gives the distribution function, q2phase gives the quantile function, r2phase generates random deviates, and h2phase gives the hazard.

Alternative parameterisation

An individual following this distribution can be seen as coming from a mixture of two populations:

1) "short stayers" whose mean sojourn time is M1=M_1 = 1/(λ1+μ1) 1/(\lambda_1+\mu_1) and sojourn distribution is exponential with rate λ1+μ1\lambda_1 + \mu_1.

2) "long stayers" whose mean sojourn time M2=M_2 = 1/(λ1+μ1)+1/μ2 1/(\lambda_1+\mu_1) + 1/\mu_2 and sojourn distribution is the sum of two exponentials with rate λ1+\lambda_1 + μ1 \mu_1 and μ2\mu_2 respectively. The individual is a "long stayer" with probability p=λ1/(λ1+μ1)p=\lambda_1/(\lambda_1 + \mu_1).

Thus a two-phase distribution can be more intuitively parameterised by the short and long stay means M1<M2M_1 < M_2 and the long stay probability pp. Given these parameters, the transition intensities are λ1=p/M1\lambda_1=p/M_1, μ1=(1p)/M1\mu_1=(1-p)/M_1, and μ2=1/(M2M1)\mu_2=1/(M_2-M_1). This can be useful for choosing intuitively reasonable initial values for procedures to fit these models to data.

The hazard is increasing at least if M2<2M1M_2 < 2M_1, and also only if (M22M1)/(M2M1)<p(M_2 - 2M_1)/(M_2 - M_1) < p.

For increasing hazards with λ1+μ1μ2\lambda_1 + \mu_1 \leq \mu_2, the maximum hazard ratio between any time tt and time 0 is 1/(1p)1/(1-p).

For increasing hazards with λ1+μ1μ2\lambda_1 + \mu_1 \geq \mu_2, the maximum hazard ratio is M1/((1p)(M2M_1/((1-p)(M_2 - M1)) M_1)). This is the minimum hazard ratio for decreasing hazards.

Author(s)

C. H. Jackson chris.jackson@mrc-bsu.cam.ac.uk

References

C. Dutang, V. Goulet and M. Pigeon (2008). actuar: An R Package for Actuarial Science. Journal of Statistical Software, vol. 25, no. 7, 1-37. URL http://www.jstatsoft.org/v25/i07


[Package msm version 1.7.1 Index]