A single cluster-level frailty w drives both processes:
lambda_R(t|w) = w lambda_0R(t) and
lambda_T(t|w) = w^alpha lambda_0T(t), with
w ~ Gamma(1/theta, 1/theta) so that E[w] = 1 and Var[w] =
theta. alpha is the association parameter: 0 makes the terminal
event independent of the recurrent process, 1 is the ordinary shared
frailty, and a negative alpha means clusters with many recurrences
die later. Baselines are constant (exponential). Because w^alpha
appears in the terminal hazard the frailty cannot be integrated out
in closed form, so the expectation is taken on a 16-point generalised
Gauss-Laguerre rule matched to the frailty density (Golub-Welsch), for
which E[w] = 1 and E[w^2] = 1 + theta hold to ten digits;
maximisation is a fixed-length coordinate golden-section search with a
guarded comparison, so a flat optimum cannot send the two language arms
into different sub-intervals.
Value
List with estimate (alpha), alpha,
theta, lambda_r, lambda_t, loglik,
n_clusters, n_recurrent, n_terminal,
exposure, naive_lambda_r, naive_lambda_t,
n, method.
References
Liu, Wolfe and Huang (2004), Shared frailty models for recurrent events and a terminal event, Biometrics 60(3):747-756. doi:10.1111/j.0006-341X.2004.00225.x
Examples
set.seed(1)
Jntfr(time = rexp(30), event = rbinom(30, 1, 0.7),
terminal = rbinom(30, 1, 0.3), cluster = rep(1:10, 3))
#> $estimate
#> [1] 0.6523028
#>
#> $alpha
#> [1] 0.6523028
#>
#> $theta
#> [1] 0.001000011
#>
#> $lambda_r
#> [1] 0.6459255
#>
#> $lambda_t
#> [1] 0.3075828
#>
#> $loglik
#> [1] -51.98549
#>
#> $n_clusters
#> [1] 10
#>
#> $n_recurrent
#> [1] 21
#>
#> $n_terminal
#> [1] 10
#>
#> $exposure
#> [1] 32.52618
#>
#> $naive_lambda_r
#> [1] 0.6456337
#>
#> $naive_lambda_t
#> [1] 0.3074446
#>
#> $n
#> [1] 30
#>
#> $method
#> [1] "lambda_R = w lambda_0R, lambda_T = w^alpha lambda_0T, w ~ Gamma(1/theta, 1/theta), Liu, Wolfe & Huang (2004)"
#>
#> attr(,"class")
#> [1] "morie_rich_result" "list"