Skip to content

Commit

Permalink
Minor fix to predict functions; small changes in the vignette to prev…
Browse files Browse the repository at this point in the history
…ent warnings
  • Loading branch information
mclements committed Oct 7, 2016
1 parent ff69a25 commit 49c857b
Showing 3 changed files with 71 additions and 66 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -6,7 +6,7 @@ Authors@R: c(person("Mark", "Clements", role = c("aut", "cre"),
person("Xing-Rong", "Liu", role = "aut",
email = "xingrong.liu@ki.se"),
person("Paul", "Lambert", role = "ctb", email="pl4@leicester.ac.uk"))
Version: 1.3.3
Version: 1.3.4
Date: 2016-10-07
Depends: R (>= 2.10), methods, survival, splines
Imports: graphics, Rcpp (>= 0.10.2), numDeriv, stats, mgcv, bbmle (>= 1.0.3), fastGHQuad
100 changes: 53 additions & 47 deletions R/pm2-3.R
Original file line number Diff line number Diff line change
@@ -1022,31 +1022,34 @@ setMethod("predict", "stpm2",
##`%c%` <- function(f,g) function(...) g(f(...)) # function composition
setMethod("plot", signature(x="stpm2", y="missing"),
function(x,y,newdata,type="surv",
xlab=NULL,ylab=NULL,line.col=1,ci.col="grey",lty=par("lty"),
add=FALSE,ci=!add,rug=!add,
var=NULL,exposed=incrVar(var),...) {
y <- predict(x,newdata,type=if (type=="fail") "surv" else type,var=var,exposed=exposed,grid=TRUE,se.fit=TRUE)
if (type=="fail") y <- data.frame(Estimate=1-y$Estimate, lower=1-y$upper, upper=1-y$lower)
if (is.null(xlab)) xlab <- deparse(x@timeExpr)
if (is.null(ylab))
ylab <- switch(type,hr="Hazard ratio",hazard="Hazard",surv="Survival",density="Density",
sdiff="Survival difference",hdiff="Hazard difference",cumhaz="Cumulative hazard",
loghazard="log(hazard)",link="Linear predictor",meansurv="Mean survival",
meansurvdiff="Difference in mean survival",odds="Odds",or="Odds ratio",
margsurv="Marginal survival",marghaz="Marginal hazard",marghr="Marginal hazard ratio", haz="Hazard",fail="Failure",
meanhaz="Mean hazard")
xx <- attr(y,"newdata")
xx <- eval(x@timeExpr,xx) # xx[,ncol(xx)]
if (!add) matplot(xx, y, type="n", xlab=xlab, ylab=ylab, ...)
if (ci) polygon(c(xx,rev(xx)), c(y[,2],rev(y[,3])), col=ci.col, border=ci.col)
lines(xx,y[,1],col=line.col,lty=lty)
if (rug) {
Y <- x@y
eventTimes <- Y[Y[,ncol(Y)]==1,ncol(Y)-1]
rug(eventTimes,col=line.col)
}
return(invisible(y))
})
xlab=NULL,ylab=NULL,line.col=1,ci.col="grey",lty=par("lty"),
add=FALSE,ci=!add,rug=!add,
var=NULL,exposed=incrVar(var),...) {
y <- predict(x,newdata,type=if (type=="fail") "surv" else type,var=var,exposed=exposed,
grid=!(x@timeVar %in% names(newdata)), se.fit=ci)
if (type=="fail") y <- if (ci) data.frame(Estimate=1-y$Estimate, lower=1-y$upper, upper=1-y$lower) else data.frame(Estimate=1-y)
if (is.null(xlab)) xlab <- deparse(x@timeExpr)
if (is.null(ylab))
ylab <- switch(type,hr="Hazard ratio",hazard="Hazard",surv="Survival",density="Density",
sdiff="Survival difference",hdiff="Hazard difference",cumhaz="Cumulative hazard",
loghazard="log(hazard)",link="Linear predictor",meansurv="Mean survival",
meansurvdiff="Difference in mean survival",odds="Odds",or="Odds ratio",
margsurv="Marginal survival",marghaz="Marginal hazard",marghr="Marginal hazard ratio", haz="Hazard",fail="Failure",
meanhaz="Mean hazard")
xx <- attr(y,"newdata")
xx <- eval(x@timeExpr,xx) # xx[,ncol(xx)]
if (!add) matplot(xx, y, type="n", xlab=xlab, ylab=ylab, ...)
if (ci) {
polygon(c(xx,rev(xx)), c(y[,2],rev(y[,3])), col=ci.col, border=ci.col)
lines(xx,y[,1],col=line.col,lty=lty)
} else lines(xx,y,col=line.col,lty=lty)
if (rug) {
Y <- x@y
eventTimes <- Y[Y[,ncol(Y)]==1,ncol(Y)-1]
rug(eventTimes,col=line.col)
}
return(invisible(y))
})
if (FALSE) {
lpfunc <- function(delta,fit,dataset,var) {
dataset[[var]] <- dataset[[var]]+delta
@@ -2052,28 +2055,31 @@ setMethod("plot", signature(x="pstpm2", y="missing"),
lwd=par("lwd"),
add=FALSE,ci=!add,rug=!add,exposed=incrVar(var),
var=NULL,...) {
y <- predict(x,newdata,type=if (type=="fail") "surv" else type,var=var,exposed=exposed,grid=TRUE,se.fit=TRUE)
if (type=="fail") y <- data.frame(Estimate=1-y$Estimate, lower=1-y$upper, upper=1-y$lower)
if (is.null(xlab)) xlab <- deparse(x@timeExpr)
if (is.null(ylab))
ylab <- switch(type,hr="Hazard ratio",hazard="Hazard",surv="Survival",density="Density",
sdiff="Survival difference",hdiff="Hazard difference",cumhaz="Cumulative hazard",
loghazard="log(hazard)",link="Linear predictor",meansurv="Mean survival",
meansurvdiff="Difference in mean survival",odds="Odds",or="Odds ratio",
margsurv="Marginal survival",marghaz="Marginal hazard",
marghr="Marginal hazard ratio", haz="Hazard", fail="Failure", meanhaz="Mean hazard")
xx <- attr(y,"newdata")
xx <- eval(x@timeExpr,xx) # xx[,ncol(xx)]
if (!add) matplot(xx, y, type="n", xlab=xlab, ylab=ylab, ...)
if (ci) polygon(c(xx,rev(xx)), c(y[,2],rev(y[,3])), col=ci.col, border=ci.col)
lines(xx,y[,1],col=line.col,lty=lty,lwd=lwd)
if (rug) {
Y <- x@y
eventTimes <- Y[Y[,ncol(Y)]==1,ncol(Y)-1]
rug(eventTimes,col=line.col)
}
return(invisible(y))
})
y <- predict(x,newdata,type=if (type=="fail") "surv" else type,var=var,exposed=exposed,
grid=!(x@timeVar %in% names(newdata)), se.fit=TRUE)
if (type=="fail") y <- if (ci) data.frame(Estimate=1-y$Estimate, lower=1-y$upper, upper=1-y$lower) else data.frame(Estimate=y)
if (is.null(xlab)) xlab <- deparse(x@timeExpr)
if (is.null(ylab))
ylab <- switch(type,hr="Hazard ratio",hazard="Hazard",surv="Survival",density="Density",
sdiff="Survival difference",hdiff="Hazard difference",cumhaz="Cumulative hazard",
loghazard="log(hazard)",link="Linear predictor",meansurv="Mean survival",
meansurvdiff="Difference in mean survival",odds="Odds",or="Odds ratio",
margsurv="Marginal survival",marghaz="Marginal hazard",
marghr="Marginal hazard ratio", haz="Hazard", fail="Failure", meanhaz="Mean hazard")
xx <- attr(y,"newdata")
xx <- eval(x@timeExpr,xx) # xx[,ncol(xx)]
if (!add) matplot(xx, y, type="n", xlab=xlab, ylab=ylab, ...)
if (ci) {
polygon(c(xx,rev(xx)), c(y[,2],rev(y[,3])), col=ci.col, border=ci.col)
lines(xx,y[,1],col=line.col,lty=lty,lwd=lwd)
} else lines(xx,y,col=line.col,lty=lty,lwd=lwd)
if (rug) {
Y <- x@y
eventTimes <- Y[Y[,ncol(Y)]==1,ncol(Y)-1]
rug(eventTimes,col=line.col)
}
return(invisible(y))
})

## sandwich variance estimator (from the sandwich package)

35 changes: 17 additions & 18 deletions vignettes/Introduction.Rnw
Original file line number Diff line number Diff line change
@@ -50,8 +50,6 @@ A more detailed theoretical development is available from the paper by Liu, Pawi

Why would you want to use these models?

\section{}


\section{Mean survival}

@@ -60,8 +58,8 @@ This has a useful interpretation for causal inference.
$E_Z(S(t|Z,X=1))-E_Z(S(t|Z,X=0))$

\begin{verbatim}
fit <- rstpm()
predict(fit,type="meansurvdiff",newdata=data)
fit <- stpm(...)
predict(fit,type="meansurv",newdata=data)
\end{verbatim}

\section{Cure models}
@@ -77,10 +75,8 @@ options(width=80,useFancyQuotes="UTF-8")
require(rstpm2)
@
<<>>=
data(popmort, package="rstpm2")
data(colon,package="rstpm2")
popmort2 <- transform(popmort,exitage=age,exityear=year,age=NULL,year=NULL)
colon2 <- within(colon, {
popmort2 <- transform(rstpm2::popmort,exitage=age,exityear=year,age=NULL,year=NULL)
colon2 <- within(rstpm2::colon, {
status <- ifelse(surv_mm>120.5,1,status)
tm <- pmin(surv_mm,120.5)/12
exit <- dx+tm*365.25
@@ -120,23 +116,26 @@ figures with that in Stata. For the predictions, the Stata model gives:
6. | .86108264 .8542898 .8675839 |
+---------------------------------+
\end{verbatim}
We can estimate the proportion of failures at the end of follow-up using:
We can estimate the proportion of failures prior to the last event time:
<<>>=
newdata.eof <- data.frame(year8594 = unique(colon2$year8594),
tm=max(colon2$tm))
tm=10)
1-predict(fit0, newdata.eof, type="surv", se.fit=TRUE)
1-predict(fit, newdata.eof, type="surv", se.fit=TRUE)
predict(fit, newdata.eof, type="haz", se.fit=TRUE)
@
We can plot the predicted survival estimates:
\begin{center}
<<fig=TRUE>>=
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 85-94"), ylim=0:1,
tms=seq(0,10,length=301)[-1]
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 85-94", tm=tms), ylim=0:1,
xlab="Time since diagnosis (years)", ylab="Relative survival")
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 75-84"),
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 75-84",tm=tms),
add=TRUE,line.col="red",rug=FALSE)
plot(fit,newdata=data.frame(year8594 = "Diagnosed 85-94"),
## warnings: Predicted hazards less than zero for cure
plot(fit,newdata=data.frame(year8594 = "Diagnosed 85-94",tm=tms),
add=TRUE,ci=FALSE,lty=2,rug=FALSE)
plot(fit,newdata=data.frame(year8594="Diagnosed 75-84"),
plot(fit,newdata=data.frame(year8594="Diagnosed 75-84",tm=tms),
add=TRUE,rug=FALSE,line.col="red",ci=FALSE,lty=2)
legend("topright",c("85-94 without cure","75-84 without cure",
"85-94 with cure","75-84 with cure"),
@@ -148,16 +147,16 @@ And the hazard curves:

\begin{center}
<<fig=TRUE>>=
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 85-94"),
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 85-94", tm=tms),
ylim=c(0,0.5), type="hazard",
xlab="Time since diagnosis (years)",ylab="Excess hazard")
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 75-84"),
plot(fit0,newdata=data.frame(year8594 = "Diagnosed 75-84", tm=tms),
type="hazard",
add=TRUE,line.col="red",rug=FALSE)
plot(fit,newdata=data.frame(year8594 = "Diagnosed 85-94"),
plot(fit,newdata=data.frame(year8594 = "Diagnosed 85-94", tm=tms),
type="hazard",
add=TRUE,ci=FALSE,lty=2,rug=FALSE)
plot(fit,newdata=data.frame(year8594="Diagnosed 75-84"),
plot(fit,newdata=data.frame(year8594="Diagnosed 75-84", tm=tms),
type="hazard",
add=TRUE,rug=FALSE,line.col="red",ci=FALSE,lty=2)
legend("topright",c("85-94 without cure","75-84 without cure",

0 comments on commit 49c857b

Please sign in to comment.