General semiparametric shared frailty model

Description

frailtySurv provides a suite of functions for generating clustered survival data, fitting multivariate shared frailty models under a wide range of frailty distributions, and visualizing the output. The semi-parametric estimators have better asymptotic properties than most existing implementations, including consistent and asymptotically-normal estimators. Moreover, this is the first package that implements semi-parametric estimators with inverse Gaussian and PVF frailty models.

The frailtySurv package provides functions

• genfrail for data generation

• fitfrail for model fitting

• simfrail for survival simulation

Details

Author(s)

John V. Monaco, Malka Gorfine, and Li Hsu.

References

Gorfine M, Zucker DM, Hsu L (2006) Prospective survival analysis with a general semiparametric shared frailty model: A pseudo full likelihood approach. Biometrika, 93(3), 735-741.

Zucker DM, Gorfine M, Hsu L (2008) Pseudo-full likelihood estimation for prospective survival analysis with a general semiparametric shared frailty model: Asymptotic theory. Journal of Statistical Planning and Inference, 138(7), 1998-2016.

See Also

genfrail, fitfrail, simfrail

Diabetic Retinopathy Study (DRS)

Description

The Diabetic Retinopathy Study (DRS) was performed to determine whether the onset of blindness in 197 high-risk diabetic patients could be delayed by laser treatment. The treatment was administered to one randomly-selected eye in each patient, leaving the other eye untreated. Thus, there are 394 observations, which are clustered by patient since the level of risk will tend to vary between patients. A failure occurred when visual acuity dropped to below 5/200. All patients had a visual acuity of at least 20/100 at the beginning of the study.

Usage

data("drs")

Format

A data frame with 394 rows and 8 columns. There are two rows for each subject, one row for each eye:

subject_id

unique identifier for each subject

eye

subject's eye, where 1=right and 2=left

time

the observed follow-up time

status

outcome at the end of the observation period, where 1=blindness and 0 indicates censorship

treated

a binary covariate, where 1=treated or 0=untreated

age_at_onset

age (in years) at the onset of diabetes

laser_type

type of laser used for treatment, where 1=xenon, 2=argon

diabetes_type

type of diabetes, where 1=juvenile (age at dx < 20) and 2=adult

Source

https://www.mayo.edu/research/documents/diabeteshtml/doc-10027460

Examples

## Not run: 
data(drs)

# Clustered by subject
fit.drs <- fitfrail(Surv(time, status) ~ treated + cluster(subject_id), 
                    drs, frailty="gamma")

fit.drs

# Variance estimates
vcov(fit.drs)

# Plot the estimated cumulative baseline hazard
plot(fit.drs, type="cumhaz")

## End(Not run)

Fit a shared frailty model

Description

Fit an extended Cox proportional hazards model with unobserved shared frailty variate and unspecified baseline hazard function, using a semiparametric estimation technique. See Gorfine et al.~(2006) and Zucker et al.~(2008) for details.

Usage

fitfrail(formula, dat, control, frailty, weights = NULL, se = FALSE, ...)

Arguments

Value

A fitfrail object representing the shared frailty model.

Convergence

The initial values of the regression coefficients are provided by coxph. Convergence is reached when either the relative reduction or absolute reduction in loglikelihood or score equations (depending on the fitmethod used) are below a threshold. If the maxit iterations are performed before convergence, then the algorithm terminates with a warning.

Author(s)

The estimation method was developed by Malka Gorfine, Li Hsu, and David Zucker; implemented by John V. Monaco.

References

Gorfine M, Zucker DM, Hsu L (2006) Prospective survival analysis with a general semiparametric shared frailty model: A pseudo full likelihood approach. Biometrika, 93(3), 735-741.

Monaco JV, Gorfine M, Hsu L (2018) General Semiparametric Shared Frailty Model: Estimation and Simulation with frailtySurv Journal of Statistical Software, 86(4), 1-42

Zucker DM, Gorfine M, Hsu L (2008) Pseudo-full likelihood estimation for prospective survival analysis with a general semiparametric shared frailty model: Asymptotic theory. Journal of Statistical Planning and Inference, 138(7), 1998-2016.

See Also

vcov.fitfrail, genfrail, simfrail, survfit, coxph

Examples

## Not run: 
# 
# Generate synthetic survival data with regression coefficients 
# beta = c(log(2),log(3)) and theta = 2, where the shared frailty 
# values from a gamma distribution with expectation 1 and variance theta.
# 
dat <- genfrail(N=300, K=2, beta=c(log(2),log(3)), 
                frailty="gamma", theta=2,
                censor.rate=0.35,
                Lambda_0=function(t, tau=4.6, C=0.01) (C*t)^tau)

# Fit a shared frailty model
fit <- fitfrail(Surv(time, status) ~ Z1 + Z2 + cluster(family), 
                dat, frailty="gamma")
fit

# The Lambda.fun function can give the estimated cumulative baseline hazard at
# any time
fit$Lambda.fun(seq(0, 100, by=10))

# Fit the DRS data, clustered on patient
data(drs)
fit.drs <- fitfrail(Surv(time, status) ~ treated + cluster(subject_id), 
                    drs, frailty="gamma")
fit.drs

## End(Not run)

# 
# A small example with c(log(2),log(3)) coefficients, Gamma(2) frailty, and
# 0.10 censorship.
# 
dat <- genfrail(N=30, K=2, beta=c(log(2),log(3)), 
                frailty="gamma", theta=2,
                censor.rate=0.10,
                Lambda_0=function(t, tau=4.6, C=0.01) (C*t)^tau)

# Fit a shared frailty model
fit <- fitfrail(Surv(time, status) ~ Z1 + Z2 + cluster(family), 
                dat, frailty="gamma", se=TRUE)
fit

# Summarize the survival curve
head(summary(fit))

Control parameters for fitfrail

Description

This function creates a list of control parameters needed by fitfrail.

Usage

fitfrail.control(fitmethod = "loglik", 
                 abstol = 0, reltol = 1e-6, maxit = 100,
                 int.abstol = 0, int.reltol = 1, int.maxit = 1000, 
                 init.beta="coxph", init.theta=NULL, verbose = FALSE)

Arguments

Value

A list of control parameters.

Author(s)

John V. Monaco, Malka Gorfine, Li Hsu

See Also

fitfrail

Generate survival data

Description

Generate clustered survival data from a shared frailty model, with hazard function given by

S(t)=\exp [-\Lambda_0(t) \omega_i \exp (\beta Z_{ij})]

where \Lambda_0 is the cumulative baseline hazard, \omega_i is the frailty value of cluster i, \beta is the regression coefficient vector, and Z_ij is the covariate vector for individual i in cluster j.

The baseline hazard can be specified by the inverse cumualative baseline hazard, cumulative baseline hazard, or simply the baseline hazard. Frailty values can be sampled from gamma, power variance function (PVF), log-normal, inverse Gaussian, and positive stable distributions.

Usage

genfrail(N = 300, K = 2, K.param = c(2, 0), beta = c(log(2)),
         frailty = "gamma", theta = c(2), 
         covar.distr = "normal", covar.param = c(0, 1), covar.matrix = NULL,
         censor.distr = "normal", censor.param = c(130, 15), 
         censor.rate = NULL, censor.time = NULL,
         lambda_0 = NULL, Lambda_0 = NULL, Lambda_0_inv = NULL, 
         round.base = NULL, control, ...)

Arguments

Value

A data.frame with row-observations is returned.

Author(s)

John V. Monaco, Malka Gorfine, and Li Hsu.

See Also

fitfrail

Examples

# Generate the same dataset 3 different ways

# Using the baseline hazard (least efficient)
set.seed(1234)
dat.1 <- genfrail(N = 300, K = 2, 
                  beta = c(log(2),log(3)),
                  frailty = "gamma", theta = 2,
                  lambda_0=function(t, tau=4.6, C=0.01) (tau*(C*t)^tau)/t)

# Using the cumulative baseline hazard
set.seed(1234)
dat.2 <- genfrail(N = 300, K = 2, 
                  beta = c(log(2),log(3)),
                  frailty = "gamma", theta = 2, 
                  Lambda_0 = function(t, tau=4.6, C=0.01) (C*t)^tau)

# Using the inverse cumulative baseline hazard (most efficient)
set.seed(1234)
dat.3 <- genfrail(N = 300, K = 2, 
                  beta = c(log(2),log(3)),
                  frailty = "gamma", theta = 2, 
                  Lambda_0_inv=function(t, tau=4.6, C=0.01) (t^(1/tau))/C)

# Generate data with PVF frailty, truncated Poisson cluster sizes, normal
# covariates, and 0.35 censorship from a lognormal distribution
set.seed(1234)
dat.4 <- genfrail(N = 100, K = "poisson", K.param=c(5, 1), 
                  beta = c(log(2),log(3)),
                  frailty = "pvf", theta = 0.3, 
                  covar.distr = "lognormal", 
                  censor.rate = 0.35) # Use the default baseline hazard

# Cluster sizes have size >= 2, summarized by
summary(dat.4)

# An oscillating baseline hazard
set.seed(1234)
dat.5 <- genfrail(lambda_0=function(t, tau=4.6, C=0.01, A=2, f=0.1) 
                              A^sin(f*pi*t) * (tau*(C*t)^tau)/t)

# Uniform censorship with 0.25 censoring rate
set.seed(1234)
dat.6 <- genfrail(N = 300, K = 2, 
                  beta = c(log(2),log(3)),
                  frailty = "gamma", theta = 2, 
                  censor.distr = "uniform", 
                  censor.param = c(50, 150), 
                  censor.rate = 0.25,
                  Lambda_0_inv=function(t, tau=4.6, C=0.01) (t^(1/tau))/C)

Control parameters for genfrail

Description

This function creates a list of control parameters needed by genfrail.

Usage

genfrail.control(censor.reltol = 1e-4,
                 censor.subdivisions = 1000L,
                 crowther.reltol = 1e-4,
                 crowther.subdivisions = 1000L)

Arguments

Value

A list of control parameters.

Author(s)

John V. Monaco, Malka Gorfine, Li Hsu

See Also

genfrail

Hard drive failure dataset

Description

This dataset contains the observed follow-up times and SMART statistics of 52k unique hard drives.

Daily snapshots of a large backup storage provider over 2 years were made publicly available. On each day, the Self-Monitoring, Analysis, and Reporting Technology (SMART) statistics of operational drives are recorded. When a hard drive is no longer operational, it is marked as a failure and removed from the subsequent daily snapshots. New hard drives are also continuously added to the population. In total, there are over 52k unique hard drives over approximately 2 years and 2885 (5.5%) failures.

Usage

data("hdfail")

Format

A data frame with 52422 observations on the following 8 variables.

serial

unique serial number of the hard drive

model

hard drive model

time

the observed followup time

status

failure indicator

temp

temperature in Celsius

rsc

binary covariate, where 1 indicates sectors that encountered read, write, or verification errors

rer

binary covariate, where 1 indicates a non-zero rate of errors that occur in hardware when reading from data from disk.

psc

binary covariate, where 1 indicates there were sectors waiting to be remapped due to an unrecoverable error.

Source

https://www.backblaze.com/cloud-storage/resources/hard-drive-test-data

Examples

## Not run: 
data(hdfail)

# Select only Western Digital hard drives
dat <- subset(hdfail, grepl("WDC", model))

fit.hd <- fitfrail(Surv(time, status) ~ temp + rer + rsc 
                                      + psc + cluster(model), 
                   dat, frailty="gamma", fitmethod="score")

fit.hd

## End(Not run)

Plot method for fitfrail objects

Description

Plot the cumulative baseline hazard estimates or the parameter trace from model estimation.

Usage

## S3 method for class 'fitfrail'
plot(x, type = c("cumhaz", "trace"), ...)

Arguments

Value

The plot object.

Author(s)

John. V Monaco, Malka Gorfine, Li Hsu

See Also

fitfrail

Examples

## Not run: 
data(drs)
fit.drs <- fitfrail(Surv(time, status) ~ treated + cluster(subject_id), 
                    drs, frailty="gamma")

# Plot the parameter and log-likelihood trace
plot(fit.drs, type="trace") 

# This may take a while to run. 
# Use parameter B to specify the number of repetitions in the weighted bootstrap
plot(fit.drs, type="cumhaz", CI=0.95)

## End(Not run)

Plot method for simfrail objects

Description

Plot the estimated parameter residuals or the mean estimated cumulative baseline hazard.

Usage

## S3 method for class 'simfrail'
plot(x, type = c("residuals", "cumhaz"), ...)

Arguments

Value

The plot object.

Author(s)

John. V Monaco, Malka Gorfine, Li Hsu

See Also

simfrail

Examples

## Not run: 
set.seed(2015)
sim <- simfrail(1000,
    genfrail.args=alist(beta=c(log(2),log(3)), frailty="gamma", 
                        censor.rate=0.30, N=300, K=2, theta=2,
                        covar.distr="uniform", covar.param=c(0, 1),
                        Lambda_0=function(t, tau=4.6, C=0.01) (C*t)^tau), 
    fitfrail.args=alist(formula=Surv(time, status) ~ Z1 + Z2 
                                                   + cluster(family), 
                        frailty="gamma"),
    Lambda.times=1:120)

# Make a boxplot of residuals
plot(sim, type="residuals")

# Plot the mean estimated cumulative baseline hazard and empirical 0.95 CI
plot(sim, type="cumhaz")

## End(Not run)

Simulate survival data and fit models

Description

Generates simulated clustered survival data by repeatedly generating data, using a shared frailty model, and fitting the models. Respective arguments are passed to genfrail and coxph, and the resulting parameter estimates are aggregated and summarized.

This function is similar to simfrail, except models are fitted using the coxph.

Usage

simcoxph(reps, genfrail.args, coxph.args, Lambda.times, cores = 0)

Arguments

Value

A simcoxph object that is essentially a data.frame of the resulting parameter estimates. Each row is a single run, and columns are as follows.

Author(s)

John. V Monaco, Malka Gorfine, Li Hsu

See Also

coxph, genfrail, simfrail

Examples

## Not run: 
sim <- simcoxph(reps=100,
                genfrail.args=alist(
                  N=50, K=2,
                  beta=c(log(2),log(3)), 
                  frailty="gamma", theta=2,
                  Lambda_0 = function(t, tau=4.6, C=0.01) (C*t)^tau), 
                coxph.args=alist(
                  formula=Surv(time, status) ~ Z1 + Z2 + cluster(family), 
                  frailty="gamma"),
                Lambda.times=1:120, cores = 0)

# Summarize the results
summary(sim)

# Plot the residuals
plot(sim, "residuals")

## End(Not run)

Simulate survival data and fit models

Description

Generates simulated clustered survival data by repeatedly generating data, using a shared frailty model, and fitting the models. Respective arguments are passed to genfrail and fitfrail, and the resulting parameter estimates are aggregated and summarized.

Usage

simfrail(reps, genfrail.args, fitfrail.args, Lambda.times, 
         vcov.args = list(), cores = 0, skip.SE = FALSE)

Arguments

Value

A simfrail object that is essentially a data.frame of the resulting parameter estimates. Each row is a single run, and columns are as follows.

Author(s)

John. V Monaco, Malka Gorfine, Li Hsu

See Also

genfrail, fitfrail

Examples

## Not run: 
sim <- simfrail(reps=100,
                genfrail.args=alist(
                  N=50, K=2,
                  beta=c(log(2),log(3)), 
                  frailty="gamma", theta=2,
                  Lambda_0 = function(t, tau=4.6, C=0.01) (C*t)^tau), 
                fitfrail.args=alist(
                  formula=Surv(time, status) ~ Z1 + Z2 + cluster(family), 
                  frailty="gamma"),
                Lambda.times=1:120, cores = 0)

# Summarize the results
summary(sim)

# Plot the residuals
plot(sim, "residuals")

## End(Not run)

Summary of the survival curve

Description

Returns a data.frame summarizing the survival curve of the fitted model. If specified, this function uses a weighted bootstrap procedure to calculate SE of the survival curve estimates. Subsequente calls with the same arguments will use the cached SE and avoid performing the weighted bootstrap again.

Usage

## S3 method for class 'fitfrail'
summary(object, type = "survival", Lambda.times = NULL, 
                           censored = FALSE, se = FALSE, CI = 0.95, ...)

Arguments

Value

A data.frame summarizing the survival curve with the following columns.

Note

Similar to summary.survfit function in the survival package.

See Also

fitfrail, vcov.fitfrail

Examples

## Not run: 
dat <- genfrail(N=200, K=2, beta=c(log(2),log(3)), 
                frailty="gamma", theta=2,
                censor.rate=0.35,
                Lambda_0=function(t, tau=4.6, C=0.01) (C*t)^tau)

fit <- fitfrail(Surv(time, status) ~ Z1 + Z2 + cluster(family), 
                dat, frailty="gamma")

surv <- summary(fitfrail, B=50, se=TRUE, CI=0.95)
head(surv)

## End(Not run)

Compute variance/covariance matrix for fitfrail model

Description

Compute the variance/covariance matrix for fitfrail estimated parameters. This can be performed by a an asymptotically-normal and consistent variance estimator or a weighted bootstrap. The resulting covariance matrix is cached in the fitted object and later retrieved if the same arguments to vcov.fitfrail are supplied.

Usage

## S3 method for class 'fitfrail'
vcov(object, boot=FALSE, B=100, Lambda.times=NULL, cores=0, ...)

Arguments

Value

variance/covariance matrix for the fitfrail model parameters

See Also

fitfrail

Examples

## Not run: 
dat <- genfrail(N=200, K=2, beta=c(log(2),log(3)), 
                frailty="gamma", theta=2,
                censor.rate=0.35,
                Lambda_0=function(t, tau=4.6, C=0.01) (C*t)^tau)

fit <- fitfrail(Surv(time, status) ~ Z1 + Z2 + cluster(family), 
                dat, frailty="gamma")

# boot=TRUE will give the weighted bootstrap variance estimates
COV <- vcov(fit, boot=FALSE)
COV

## End(Not run)