Computes the cumulative sum of parameter functions and the standard error of it. Optionally the exponential is applied to the parameter functions before it is cumulated.
1 2 3 4 5 6 7 8 
obj 
A model object (of class

ctr.mat 
Contrast matrix defining the parameter functions from the parameters of the model. 
subset 
Subset of the parameters of the model to which

intl 
Interval length for the cumulation. Either a constant or
a numerical vector of length 
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 logscale, 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? 
The purpose of this function is to compute cumulative rate based on a
model for the rates. If the model is a multiplicative model for the
rates, the purpose of ctr.mat
is to return a vector of rates or
logrates when applied to the coefficients of the model. If logrates
are returned from the model, the they should be exponentiated before
cumulated, and the variances computed accordingly. Since the primary
use is for loglinear Poisson models the Exp
parameter defaults
to TRUE.
The ci.Exp
argumen ensures that the confidence intervals for
A matrix with 4 columns: Estimate, lower and upper c.i. and standard
error. If sample
is TRUE, a sampled vector is reurned, 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.
Bendix Carstensen, http://BendixCarstensen.com
See also ci.lin
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41  # Packages required for this example
library( splines )
library( survival )
data( lung )
par( mfrow=c(1,2) )
# Plot the Kaplanmeierestimator
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 )
# Cut the followup every 10 days
sL < splitLexis( lungL, "tfd", breaks=seq(0,1100,10) )
str( sL )
summary( sL )
# Fit a Poisson model with a natural spline for the effect of time.
# Extract the variables needed
D < status(sL, "exit")
Y < dur(sL)
tB < timeBand( sL, "tfd", "left" )
MM < ns( tB, knots=c(50,100,200,400,700), intercept=TRUE )
mp < glm( D ~ MM  1 + offset(log(Y)),
family=poisson, eps=10^8, maxit=25 )
# Contrast matrix to extract effects, i.e. matrix to multiply with the
# coefficients to produce the logrates: unique rows of MM, in time order.
T.pt < sort( unique( tB ) )
T.wh < match( T.pt, tB )
Lambda < ci.cum( mp, ctr.mat=MM[T.wh,], intl=diff(c(0,T.pt)) )
# Put the estimated survival function on top of the KMestimator
matlines( c(0,T.pt[1]), exp(Lambda[,1:3]), lwd=c(3,1,1), lty=1, col="Red" )
# Extract and plot the fitted intensity function
lambda < ci.lin( mp, ctr.mat=MM[T.wh,], Exp=TRUE )
matplot( T.pt, lambda[,5:7]*10^3, type="l", lwd=c(3,1,1), col="black", lty=1,
log="y", ylim=c(0.2,20) )

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.