emfrail_pll function

Profile log-likelihood calculation

emfrail_pll(formula, data, distribution = emfrail_dist(), values)

Arguments

• formula: Same as in emfrail
• data: Same as in emfrail
• distribution: Same as in emfrail
• values: A vector of values on where to calculate the profile likelihood. See details.

Returns

The profile log-likelihood at the specific value of the frailty parameter

Details

This function can be used to calculate the profile log-likelihood for different values of $\theta$. The scale is that of theta as defined in emfrail_dist(). For the gamma and pvf frailty, that is the inverse of the frailty variance.

Note

This function is just a simple wrapper for emfrail() with the control argument a call from emfrail_control with the option opt_fit = FALSE. More flexibility can be obtained by calling emfrail with this option, especially for setting other emfrail_control parameters.

Examples

fr_var <- seq(from = 0.01, to = 1.4, length.out = 20)
pll_gamma <- emfrail_pll(formula = Surv(time, status) ~  rx + sex + cluster(litter),
 data =  rats,
 values = 1/fr_var )
 plot(fr_var, pll_gamma,
     type = "l",
     xlab = "Frailty variance",
     ylab = "Profile log-likelihood")

# check with coxph;
# attention: theta is the the inverse frailty variance in emfrail,
# but theta is the frailty variance in coxph.

pll_cph <- sapply(fr_var, function(th)
  coxph(data =  rats, formula = Surv(time, status) ~ rx + sex + frailty(litter, theta = th),
        method = "breslow")$history[[1]][[3]])

lines(fr_var, pll_cph, col = 2)

# Same for inverse gaussian
pll_if <- emfrail_pll(Surv(time, status) ~  rx + sex + cluster(litter),
                      rats,
                      distribution = emfrail_dist(dist = "pvf"),
                      values = 1/fr_var )

# Same for pvf with a positive pvfm parameter
pll_pvf <- emfrail_pll(Surv(time, status) ~  rx + sex + cluster(litter),
                       rats,
                       distribution = emfrail_dist(dist = "pvf", pvfm = 1.5),
                       values = 1/fr_var )

miny <- min(c(pll_gamma, pll_cph, pll_if, pll_pvf))
maxy <- max(c(pll_gamma, pll_cph, pll_if, pll_pvf))

plot(fr_var, pll_gamma,
     type = "l",
     xlab = "Frailty variance",
     ylab = "Profile log-likelihood",
     ylim = c(miny, maxy))
points(fr_var, pll_cph, col = 2)
lines(fr_var, pll_if, col = 3)
lines(fr_var, pll_pvf, col = 4)

legend(legend = c("gamma (emfrail)", "gamma (coxph)", "inverse gaussian", "pvf, m=1.5"),
       col = 1:4,
       lty = 1,
       x = 0,
       y = (maxy + miny)/2)