R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Compute cumulative sum of estimates.

ci.cum

R Documentation

Compute cumulative sum of estimates.

Description

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.

Usage

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

Arguments

`obj`	A model object (of class `lm`, `glm`, `coxph`, `survreg`, `lme`,`mer`,`nls`,`gnlm`, `MIresult` or `polr`).
`ctr.mat`	Contrast matrix defining the parameter functions from the parameters of the model.
`subset`	Subset of the parameters of the model to which `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 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 log-rates when applied to the coefficients of the model. 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.

The ci.Exp argumen ensures that the confidence intervals for

Value

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.

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 )

# Cut the follow-up 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 log-rates: 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 KM-estimator
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) )

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Epi)

Attaching package: 'Epi'

The following object is masked from 'package:base':

    merge.data.frame

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Epi/ci.cum.Rd_%03d_medium.png", width=480, height=480)
> ### Name: ci.cum
> ### Title: Compute cumulative sum of estimates.
> ### Aliases: ci.cum
> ### Keywords: models regression
> 
> ### ** 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 )
NOTE: entry is assumed to be 0 on the tfd timescale.
> summary( lungL )
    
Transitions:
     To
From  0   1  Records:  Events: Risk time:  Persons:
   0 63 165       228      165      69593       228
> 
> # Cut the follow-up every 10 days
> sL <- splitLexis( lungL, "tfd", breaks=seq(0,1100,10) )
> str( sL )
Classes 'Lexis' and 'data.frame':	7071 obs. of  15 variables:
 $ lex.id   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ tfd      : num  0 10 20 30 40 50 60 70 80 90 ...
 $ lex.dur  : num  10 10 10 10 10 10 10 10 10 10 ...
 $ lex.Cst  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ lex.Xst  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ inst     : num  3 3 3 3 3 3 3 3 3 3 ...
 $ time     : num  306 306 306 306 306 306 306 306 306 306 ...
 $ status   : num  2 2 2 2 2 2 2 2 2 2 ...
 $ age      : num  74 74 74 74 74 74 74 74 74 74 ...
 $ sex      : num  1 1 1 1 1 1 1 1 1 1 ...
 $ ph.ecog  : num  1 1 1 1 1 1 1 1 1 1 ...
 $ ph.karno : num  90 90 90 90 90 90 90 90 90 90 ...
 $ pat.karno: num  100 100 100 100 100 100 100 100 100 100 ...
 $ meal.cal : num  1175 1175 1175 1175 1175 ...
 $ wt.loss  : num  NA NA NA NA NA NA NA NA NA NA ...
 - attr(*, "breaks")=List of 1
  ..$ tfd: num  0 10 20 30 40 50 60 70 80 90 ...
 - attr(*, "time.scales")= chr "tfd"
 - attr(*, "time.since")= chr ""
> summary( sL )
    
Transitions:
     To
From    0   1  Records:  Events: Risk time:  Persons:
   0 6906 165      7071      165      69593       228
> 
> # 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 log-rates: 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 KM-estimator
> 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) )
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>