Epi: ci.cum – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

ci.cum

Compute cumulative sum of estimates.

Description

Computes the cumulative sum of parameter functions and the standard error of it. Used for computing the cumulative rate, or the survival function based on a glm with parametric baseline.

Usage

ci.cum( obj,
    ctr.mat = NULL,
     subset = NULL,
       intl = 1,
      alpha = 0.05,
        Exp = TRUE,
     ci.Exp = FALSE,
     sample = FALSE )
ci.surv( obj,
    ctr.mat = NULL,
     subset = NULL,
       intl = 1,
      alpha = 0.05,
        Exp = TRUE,
     sample = FALSE )

Arguments

`obj`	A model object (of class `lm`, `glm`.
`ctr.mat`	Matrix or data frame. If `ctr.mat` is a matrix, it should be a contrast matrix to be multiplied to the parameter vector, i.e. the desired linear function of the parameters. If it is a data frame it should have columns corresponding to a prediction data frame for `obj`, see details for `ci.lin`.
`subset`	Subset of the parameters of the model to which a matrix `ctr.mat` should be applied.
`intl`	Interval length for the cumulation. Either a constant or a numerical vector of length `nrow(ctr.mat)`.
`alpha`	Significance level used when computing confidence limits.
`Exp`	Should the parameter function be exponentiated before it is cumulated?
`ci.Exp`	Should the confidence limits for the cumulative rate be computed on the log-scale, thus ensuring that exp(-cum.rate) is always in [0,1]?
`sample`	Should a sample of the original parameters be used to compute a cumulative rate?

Details

The purpose of this function is to the compute cumulative rate (integrated intensity) at a set of points based on a model for the rates. ctr.mat is a matrix which, when premultiplied to the parameters of the model returns the (log)rates at a set of increasing time points. If log-rates are returned from the model, the they should be exponentiated before cumulated, and the variances computed accordingly. Since the primary use is for log-linear Poisson models the Exp parameter defaults to TRUE.

Each row in the object supplied via ctr.mat is assumed to represent a midpoint in an interval. ci.cum will then return the cumulative rates at the end of these intervals. ci.surv will return the survival probability at the start of each of these intervals, assuming the the first interval starts at 0 - the first row of the result is c(1,1,1).

The ci.Exp argument ensures that the confidence intervals for the cumulative rates are always positive, so that exp(-cum.rate) is always in [0,1].

Value

A matrix with 3 columns: Estimate, lower and upper c.i. and standard error, unless se=TRUE, in which case the standard error is returned too. If sample is TRUE, a sampled vector is returned, if sample is numeric a matrix with sample columns is returned, each column a cumulative rate based on a random sample from the distribution of the parameter estimates.

ci.surv returns a 3 column matrix with estimate, lower and upper confidence interval.

Author(s)

Bendix Carstensen, http://bendixcarstensen.com

Examples

# Packages required for this example
library( splines )
library( survival )
data( lung )
par( mfrow=c(1,2) )

# Plot the Kaplan-meier-estimator
plot( survfit( Surv( time, status==2 ) ~ 1, data=lung ) )

# Declare data as Lexis
lungL <- Lexis( exit=list(tfd=time),
                exit.status=(status==2)*1,
                data=lung )
summary( lungL )

# Split the follow-up every 10 days
sL <- splitLexis( lungL, "tfd", breaks=seq(0,1100,10) )
summary( sL )

# Fit a Poisson model with a natural spline for the effect of time (left
# end points of intervals are used as covariate)
mp <- glm( cbind(lex.Xst==1,lex.dur) ~ Ns(tfd,knots=c(0,50,100,200,400,700)),
           family=poisreg, data=sL )

# mp is now a model for the rates along the time scale tfd
# prediction data frame for select time points on this time scale
nd <- data.frame( tfd = seq(5,995,10) ) # *midpoints* of intervals
Lambda <- ci.cum ( mp, nd, intl=10 )
surv   <- ci.surv( mp, nd, intl=10 )

# Put the estimated survival function on top of the KM-estimator
# recall the ci.surv return the survival at *start* of intervals
matshade( nd$tfd-5, surv, col="Red", alpha=0.15 )

# Extract and plot the fitted intensity function
lambda <- ci.pred( mp, nd )*365.25 # mortality 
matshade( nd$tfd, lambda, log="y", ylim=c(0.2,5), plot=TRUE,
          xlab="Time since diagnosis", ylab="Mortality per year" )

# same thing works with gam from mgcv
library( mgcv )
mg <- gam( cbind(lex.Xst==1,lex.dur) ~ s(tfd), family=poisreg, data=sL )
matshade( nd$tfd-5, ci.surv( mg, nd, intl=10 ), plot=TRUE,
          xlab="Days since diagnosis", ylab="P(survival)" )
matshade( nd$tfd  , ci.pred( mg, nd )*365.25  , plot=TRUE, log="y",
          xlab="Days since diagnosis", ylab="Mortality per 1 py" )

Epi

Statistical Analysis in Epidemiology

v2.44

GPL-2

Authors

Bendix Carstensen [aut, cre], Martyn Plummer [aut], Esa Laara [ctb], Michael Hills [ctb]

Initial release

2021-02-28