Recurrent Events Regression Based on Counts and Rate Function

This function fits recurrent event data (event counts) by gamma frailty model with spline rate function. The default model is the gamma frailty model with one piece constant baseline rate function, which is equivalent to negative binomial regression with the same shape and rate parameter in the gamma prior. Spline (including piecewise constant) baseline hazard rate function can be specified for the model fitting.

Usage

rateReg(
  formula,
  data,
  subset,
  na.action,
  start = list(),
  control = list(),
  contrasts = NULL,
  ...
)

rateReg.control(
  df = NULL,
  degree = 0L,
  knots = NULL,
  Boundary.knots = NULL,
  periodic = FALSE,
  verbose = TRUE,
  ...
)

Arguments

formula: Recur object produced by function Recur. The terminal events and risk-free episodes specified in Recur will be ignored since the model does not support them.
data: An optional data frame, list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), usually the environment from which function rateReg is called.
subset: An optional vector specifying a subset of observations to be used in the fitting process.
na.action: A function that indicates what should the procedure do if the data contains NAs. The default is set by the na.action setting of options. The "factory-fresh" default is na.omit. Other possible values inlcude na.fail, na.exclude, and na.pass. See help(na.fail) for details.
start: An optional list of starting values for the parameters to be estimated in the model. See more in Section details.
control: An optional list of parameters to control the maximization process of negative log likelihood function and adjust the baseline rate function. See more in Section details.
contrasts: An optional list, whose entries are values (numeric matrices or character strings naming functions) to be used as replacement values for the contrasts replacement function and whose names are the names of columns of data containing factors. See contrasts.arg of model.matrix.default for details.
...: Other arguments passed to rateReg.control() and stats::constrOptim().
df: A nonnegative integer to specify the degree of freedom of baseline rate function. If argument knots or degree are specified, df will be neglected whether it is specified or not.
degree: A nonnegative integer to specify the degree of spline bases.
knots: A numeric vector that represents all the internal knots of baseline rate function. The default is NULL, representing no any internal knots.
Boundary.knots: A length-two numeric vector to specify the boundary knots for baseline rate funtion. By default, the left boundary knot is the smallest origin time and the right one takes the largest censoring time from data.
periodic: A logical value indicating if periodic splines should be used.
verbose: A logical value with default TRUE. Set it to FALSE to supress messages from this function.

Value

A rateReg object, whose slots include

call: Function call of rateReg.
formula: Formula used in the model fitting.
nObs: Number of observations.
spline: A list contains
- spline: The name of splines used.
- knots: Internal knots specified for the baseline rate function.
- Boundary.knots: Boundary knots specified for the baseline rate function.
- degree: Degree of spline bases specified in baseline rate function.
- df: Degree of freedom of the model specified.
estimates: Estimated coefficients of covariates and baseline rate function, and estimated rate parameter of gamma frailty variable.
control: The control list specified for model fitting.
start: The initial guess specified for the parameters to be estimated.
na.action: The procedure specified to deal with missing values in the covariate.
xlevels: A list that records the levels in each factor variable.
contrasts: Contrasts specified and used for each factor variable.
convergCode: code returned by function optim, which is an integer indicating why the optimization process terminated. help(optim) for details.
logL: Log likelihood of the fitted model.
fisher: Observed Fisher information matrix.

Details

Function Recur in the formula response by default first checks the dataset and will report an error if the dataset does not fall into recurrent event data framework. Subject's ID will be pinpointed if its observation violates any checking rule. See Recur for all the checking rules.

Function rateReg first constructs the design matrix from the specified arguments: formula, data, subset, na.action and constrasts before model fitting. The constructed design matrix will be checked again to fit the recurrent event data framework if any observation with missing covariates is removed.

The model fitting process involves minimization of negative log likelihood function, which calls function constrOptim internally. help(constrOptim) for more details.

The argument start is an optional list that allows users to specify the initial guess for the parameter values for the minimization of negative log likelihood function. The available numeric vector elements in the list include

beta: Coefficient(s) of covariates, set to be all 0.1 by default.
theta: Parameter in Gamma(theta, 1 / theta) for frailty random effect, set to be 0.5 by default.
alpha: Coefficient(s) of baseline rate function, set to be all 0.05 by default.

The argument control allows users to control the process of minimization of negative log likelihood function passed to constrOptim and specify the boundary knots of baseline rate function.

References

Fu, H., Luo, J., & Qu, Y. (2016). Hypoglycemic events analysis via recurrent time-to-event (HEART) models. Journal Of Biopharmaceutical Statistics, 26(2), 280--298.

Examples

library(reda)

## constant rate function
(constFit <- rateReg(Recur(time, ID, event) ~ group + x1, data = simuDat))
#> Call: 
#> rateReg(formula = Recur(time, ID, event) ~ group + x1, data = simuDat)
#> 
#> Coefficients of covariates: 
#> groupTreat         x1 
#> -0.6133832  0.3248600 
#> 
#> Frailty parameter:  0.5878151 
#> 
#> Boundary knots: 
#> 0, 168
#> 
#> Coefficients of pieces:
#> M-spline1 
#>  5.124268 

## six pieces' piecewise constant rate function
(piecesFit <- rateReg(Recur(time, ID, event) ~ group + x1,
                      data = simuDat, subset = ID %in% 1:50,
                      knots = seq.int(28, 140, by = 28)))
#> Call: 
#> rateReg(formula = Recur(time, ID, event) ~ group + x1, data = simuDat, 
#>     subset = ID %in% 1:50, knots = seq.int(28, 140, by = 28))
#> 
#> Coefficients of covariates: 
#> groupTreat         x1 
#> -0.8030745  0.3361279 
#> 
#> Frailty parameter:  0.687028 
#> 
#> Internal knots: 
#> 28, 56, 84, 112, 140
#> 
#> Boundary knots: 
#> 0, 168
#> 
#> Coefficients of pieces:
#> M-spline1 M-spline2 M-spline3 M-spline4 M-spline5 M-spline6 
#> 1.0354803 1.0354803 0.7060105 1.1647606 1.1907062 1.7539445 

## fit rate function with cubic spline
(splineFit <- rateReg(Recur(time, ID, event) ~ group + x1, data = simuDat,
                      knots = c(56, 84, 112), degree = 3))
#> Call: 
#> rateReg(formula = Recur(time, ID, event) ~ group + x1, data = simuDat, 
#>     knots = c(56, 84, 112), degree = 3)
#> 
#> Coefficients of covariates: 
#> groupTreat         x1 
#> -0.6357172  0.3061369 
#> 
#> Frailty parameter:  0.5885601 
#> 
#> Internal knots: 
#> 56, 84, 112
#> 
#> Boundary knots: 
#> 0, 168
#> 
#> Coefficients of spline bases:
#> M-spline1 M-spline2 M-spline3 M-spline4 M-spline5 M-spline6 M-spline7 
#> 0.2620605 0.8194386 0.7333654 0.7877011 1.0425639 1.2607529 0.4621612 

## more specific summary
summary(constFit)
#> Call: 
#> rateReg(formula = Recur(time, ID, event) ~ group + x1, data = simuDat)
#> 
#> Coefficients of covariates: 
#>                coef exp(coef) se(coef)       z Pr(>|z|)  
#> groupTreat -0.61338   0.54152  0.28527 -2.1502  0.03154 *
#> x1          0.32486   1.38384  0.16641  1.9522  0.05092 .
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Parameter of frailty: 
#>         parameter        se
#> Frailty 0.5878151 0.1102496
#> 
#> Boundary knots:
#> 0, 168
#> 
#> Degree of spline bases: 0 
#> 
#> Coefficients of spline bases:
#>             coef se(coef)
#> M-spline1 5.1243   0.9632
#> 
#> Loglikelihood:  -1676.421 
summary(piecesFit)
#> Call: 
#> rateReg(formula = Recur(time, ID, event) ~ group + x1, data = simuDat, 
#>     subset = ID %in% 1:50, knots = seq.int(28, 140, by = 28))
#> 
#> Coefficients of covariates: 
#>                coef exp(coef) se(coef)       z Pr(>|z|)  
#> groupTreat -0.80307   0.44795  0.38513 -2.0852  0.03705 *
#> x1          0.33613   1.39952  0.23242  1.4462  0.14812  
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Parameter of frailty: 
#>         parameter        se
#> Frailty  0.687028 0.1736583
#> 
#> Internal knots: 
#> 28, 56, 84, 112, 140
#> 
#> Boundary knots:
#> 0, 168
#> 
#> Degree of spline bases: 0 
#> 
#> Coefficients of spline bases:
#>              coef se(coef)
#> M-spline1 1.03548   0.2784
#> M-spline2 1.03548   0.2784
#> M-spline3 0.70601   0.2033
#> M-spline4 1.16476   0.3127
#> M-spline5 1.19071   0.3251
#> M-spline6 1.75394   0.4926
#> 
#> Loglikelihood:  -989.3347 
summary(splineFit)
#> Call: 
#> rateReg(formula = Recur(time, ID, event) ~ group + x1, data = simuDat, 
#>     knots = c(56, 84, 112), degree = 3)
#> 
#> Coefficients of covariates: 
#>                coef exp(coef) se(coef)       z Pr(>|z|)  
#> groupTreat -0.63572   0.52956  0.28550 -2.2267  0.02597 *
#> x1          0.30614   1.35817  0.16650  1.8387  0.06596 .
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Parameter of frailty: 
#>         parameter        se
#> Frailty 0.5885601 0.1103332
#> 
#> Internal knots: 
#> 56, 84, 112
#> 
#> Boundary knots:
#> 0, 168
#> 
#> Degree of spline bases: 3 
#> 
#> Coefficients of spline bases:
#>              coef se(coef)
#> M-spline1 0.26206   0.1106
#> M-spline2 0.81944   0.2732
#> M-spline3 0.73337   0.3056
#> M-spline4 0.78770   0.3259
#> M-spline5 1.04256   0.3862
#> M-spline6 1.26075   0.3965
#> M-spline7 0.46216   0.1938
#> 
#> Loglikelihood:  -1663.231 

## model selection based on AIC or BIC
AIC(constFit, piecesFit, splineFit)
#> Warning: Models are not all fitted to the same number of observations. Consider
#> BIC instead?
#>           df      AIC
#> constFit   4 3360.843
#> piecesFit  9 1996.669
#> splineFit 10 3346.462
BIC(constFit, piecesFit, splineFit)
#>           df      BIC
#> constFit   4 3377.701
#> piecesFit  9 2030.003
#> splineFit 10 3388.608

## estimated covariate coefficients
coef(piecesFit)
#> groupTreat         x1 
#> -0.8030745  0.3361279 
coef(splineFit)
#> groupTreat         x1 
#> -0.6357172  0.3061369 

## confidence intervals for covariate coefficients
confint(piecesFit)
#>                  2.5%       97.5%
#> groupTreat -1.5579228 -0.04822615
#> x1         -0.1194082  0.79166391
confint(splineFit, "x1", 0.9)
#>            5%       95%
#> x1 0.03227031 0.5800034
confint(splineFit, 1, 0.975)
#>                1.25%      98.75%
#> groupTreat -1.275639 0.004204251

## estimated baseline rate function
splinesBase <- baseRate(splineFit)
plot(splinesBase, conf.int = TRUE)


## estimated baseline mean cumulative function (MCF) from a fitted model
piecesMcf <- mcf(piecesFit)
plot(piecesMcf, conf.int = TRUE, col = "blueviolet")


## estimated MCF for given new data
newDat <- data.frame(x1 = rep(0, 2), group = c("Treat", "Contr"))
splineMcf <- mcf(splineFit, newdata = newDat, groupName = "Group",
                 groupLevels = c("Treatment", "Control"))
plot(splineMcf, conf.int = TRUE, lty = c(1, 5))


## example of further customization by ggplot2
library(ggplot2)
plot(splineMcf) +
    geom_ribbon(aes(x = time, ymin = lower,
                    ymax = upper, fill = Group),
                data = splineMcf@MCF, alpha = 0.2) +
    xlab("Days")