swh:1:snp:6a0ac420dcf26f4ce7b9c8b5d3c5816f5620bd13
Tip revision: 40c3033311c24f521975a0b6e8685b3400562362 authored by Mark Clements on 10 May 2022, 11:30:05 UTC
version 1.5.6
version 1.5.6
Tip revision: 40c3033
multistate.R
markov_msm <-
function (x, trans, t = c(0,1), newdata = NULL, init=NULL,
tmvar = NULL,
sing.inf = 1e+10, method="adams", rtol=1e-10, atol=1e-10, slow=FALSE,
min.tm=1e-8,
utility=function(t) rep(1, nrow(trans)),
utility.sd=rep(0,nrow(trans)),
use.costs=FALSE,
transition.costs=function(t) rep(0,sum(!is.na(trans))), # per transition
transition.costs.sd=rep(0,sum(!is.na(trans))),
state.costs=function(t) rep(0,nrow(trans)), # per unit time
state.costs.sd=rep(0,nrow(trans)),
discount.rate = 0,
block.size=500,
spline.interpolation=FALSE,
debug = FALSE,
...)
{
call <- match.call()
inherits <- function(x, ...)
base::inherits(x, ...) ||
(base::inherits(x, c("hazFun","zeroModel","hrModel","stratifiedModel"))
&& base::inherits(x$base, ...))
base.classes <- c("stpm2","pstpm2","glm","survPen","gam","aft","flexsurvreg","aftreg")
stopifnot(all(sapply(x, function(xi) inherits(xi,base.classes) | is.function(xi))))
stopifnot(!is.null(newdata))
stopifnot(sum(!is.na(trans)) == length(x))
if (is.null(init)) init <- "[<-"(rep(0,nrow(trans)),1,1)
stopifnot(length(init) == nrow(trans))
## if newdata are many, then separate into blocks
if (nrow(newdata)>block.size) {
n.blocks <- ceiling(nrow(newdata)/block.size)
indices <- parallel::splitIndices(nrow(newdata), n.blocks)
lst <-
parallel::mclapply(indices,
function(index)
markov_msm(x = x, trans = trans, t = t,
newdata = newdata[index,],
init = init, tmvar = tmvar,
sing.inf = sing.inf,
method = method, rtol = rtol,
atol = atol, slow = slow,
min.tm = min.tm, utility = utility,
use.costs = use.costs,
transition.costs = transition.costs,
transition.costs.sd = transition.costs.sd,
state.costs.sd = state.costs.sd,
state.costs = state.costs,
discount.rate = discount.rate,
block.size = block.size,
spline.interpolation = spline.interpolation,
...))
return(do.call(rbind.markov_msm,lst))
}
if (use.costs)
slow <- TRUE
nt <- length(t)
if (nt < 2)
stop("number of times should be at least two")
stopifnot(length(utility(t[2])) %in% c(1,nrow(trans)))
if (is.null(tmvar) && all(sapply(x,inherits,c("stpm2","pstpm2","aft","survPen","flexsurvreg",
"aftreg"))))
tmvar <- sapply(x,function(object)
if(inherits(object,c("stpm2","pstpm2"))) object@timeVar
else if (inherits(object,"aft")) object@args$timeVar
else if (inherits(object,"flexsurvreg")) "t"
else if (inherits(object,"aftreg"))
local({lhs <- object$call$formula[[2]]
deparse(if(length(lhs)==4) lhs[[4]] else lhs[[3]])})
else object$t1.name)
stopifnot(!is.null(tmvar))
stopifnot(length(tmvar) %in% c(1,length(x)))
if (length(tmvar)==1) tmvar <- rep(tmvar, length(x))
for(i in 1:length(x)) {
if(inherits(x[[i]],"survPen") && !inherits(x[[i]],"survPenWrap"))
x[[i]] <- survPenWrap(x[[i]])
if(inherits(x[[i]],"gam") && !inherits(x[[i]],"gamWrap"))
x[[i]] <- gamWrap(x[[i]])
else if(is.function(x[[i]]))
x[[i]] <- hazFun(x[[i]], tmvar[i])
}
if (spline.interpolation) {
replace <- function(object,tmvar) {
logModels <- c("aftreg","flexsurvreg","stpm2","pstpm2","aft")
logModel <- inherits(object, logModels)
t <- if (logModel) exp(seq(log(max(c(min.tm,min(t)))), log(max(t)), length.out=1000))
else seq(max(c(min.tm,min(t))), max(t), length.out=1000)
if (inherits(object,"hazFun")) object
else if (inherits(object, "addModel"))
structure(lapply(object,replace,tmvar=tmvar), class="addModel")
else makeSplineFun(object, newdata=newdata, tmvar=tmvar, min.tm=min.tm, tm=t,
log=logModel)
}
x <- mapply(replace, x, tmvar, SIMPLIFY=FALSE)
}
ntr <- sum(!is.na(trans))
nobs <- nrow(newdata)
nstates <- nrow(trans)
n <- nstates*nobs
cs <- sapply(x,function(item) length(coef(item)))
cumcs <- cumsum(c(0,cs))
ncoef <- sum(cs)
lambda.discount <- log(1+discount.rate)
matwhich <- function(ind) {
stopifnot(sum(ind,na.rm=TRUE) == 1)
nr <- nrow(ind)
nc <- ncol(ind)
wind <- which(ind)-1
c(wind %% nr + 1, wind %/% nr + 1)
}
coef.surv_list <- function(objects)
unlist(lapply(objects,function(x) coef(x))) # sapply(objects, coef) FAILS
vcov.surv_list <- function(objects) {
vcovs <- lapply(objects, function(x) vcov(x))
lengths <- sapply(vcovs,nrow)
index <- cumsum(c(0,lengths))
m <- matrix(0,sum(lengths),sum(lengths))
for (i in 1:length(objects))
m[(index[i]+1):index[i+1],(index[i]+1):index[i+1]] <- vcovs[[i]]
rownames(m) <- colnames(m) <- names(coef(objects))
m
}
dp <- function(t, y, parms, ...) {
P <- matrix(y[1:(nstates*nobs)],nstates,nobs) # only one state vector per individual
Pu <- array(y[nstates*nobs+1:(nstates*ncoef*nobs)],c(nstates,ncoef,nobs))
Q <- array(0, c(nstates, nstates, nobs))
Qu <- array(0, c(nstates, nstates, ncoef, nobs))
QC <- array(0, c(nstates, nobs))
QuC <- array(0, c(nstates, ncoef, nobs))
PQCu <- array(0, c(nstates, ntr, nobs))
froms <- tos <- rep(0,ntr)
for (i in 1:ntr) {
newdata[[tmvar[i]]] <- pmax(min.tm,t) # special case: t=0
ij <- matwhich(trans == i)
froms[i] <- from <- ij[1]
tos[i] <- to <- ij[2]
Q[from,to,] <- pmax(0,predict(x[[i]], newdata=newdata, type="haz", tmvar=tmvar[i]))
Qu[from,to,1:cs[i]+cumcs[i],] <- t(predict(x[[i]], newdata=newdata, type="gradh",
tmvar=tmvar[i])) # transposed
if (use.costs) {
## collapse across "to" states
QC[from,] <- QC[from,] + Q[from,to,]*transition.costs(t)[[i]]
QuC[from,1:cs[i]+cumcs[i],] <- QuC[from,1:cs[i]+cumcs[i],] + Qu[from,to,1:cs[i]+cumcs[i],]*transition.costs(t)[i]
PQCu[from,i,] <- P[from,]*Q[from,to,]
}
}
Q[is.na(Q)] <- Qu[is.na(Qu)] <- QC[is.na(QC)] <- QuC[is.na(QuC)] <- 0
Q[is.infinite(Q) & Q > 0] <- Qu[is.infinite(Qu) & Qu > 0] <- sing.inf
QC[is.infinite(QC) & QC > 0] <- QuC[is.infinite(QuC) & QuC > 0] <- sing.inf
if (slow) {
dTCdt <- dPdt <- array(0,c(nstates,nobs))
dTCudt <- dPudt <- array(0,c(nstates,ncoef,nobs))
for (i in 1:nobs) {
Qi <- Q[,,i] # nstates*nstates
Pi <- P[,i] # nstates
diag(Qi) <- -rowSums(Qi)
dPdt[,i] <- t(Pi)%*%Qi # nstates
Qui <- Qu[,,,i] # nstates*nstates*ncoef
Pui <- Pu[,,i] # nstates*ncoef
for (j in 1:ncoef)
diag(Qui[,,j]) <- -rowSums(Qui[,,j])
dPudt[,,i] <- t(t(Pui)%*%Qi) + apply(Qui,3,function(x) t(Pi)%*%x) # nstates*ncoef
if (use.costs) {
QCi <- QC[,i] # nstates
QuCi <- QuC[,,i] # nstates*ncoef
dTCdt[,i] <- Pi*(state.costs(t) + QCi) # nstates
dTCudt[,,i] <- Pui*(state.costs(t) + QCi) +
Pi*QuCi # nstates*ncoef
}
}
} else {
Qu2 <- array(Qu,c(nstates*nstates,ncoef,nobs))
QuC2 <- array(QuC,c(nstates*nstates,ncoef,nobs))
lst <- .Call("multistate_ddt",P,Pu,Q,Qu2,PACKAGE="rstpm2")
dPdt <- lst[[1]]
dPudt <- lst[[2]]
dTCdt <- array(0,c(nstates,nobs)) # not implemented!
dTCudt <- array(0,c(nstates,ncoef,nobs)) # not implemented!
}
if (use.costs)
list(c(dPdt,
dPudt,
P*utility(t)*exp(-lambda.discount*t),
Pu*utility(t)*exp(-lambda.discount*t),
P*exp(-lambda.discount*t), # d/dt for gradient of U wrt utilities *and* gradient of C wrt state costs
dTCdt*exp(-lambda.discount*t),
dTCudt*exp(-lambda.discount*t),
PQCu*exp(-lambda.discount*t) # gradient of C wrt transition costs
))
else
list(c(dPdt,
dPudt,
P*utility(t)*exp(-lambda.discount*t),
Pu*utility(t)*exp(-lambda.discount*t),
P*exp(-lambda.discount*t) # d/dt for gradient of U wrt utilities (assumes u'_m(t)=1)
))
}
class(x) <- "surv_list"
Sigma <- vcov(x)
if (any(is.na(Sigma)))
stop("NAs in the covariance matrix")
if (any(is.na(coef(x))))
stop("NAs in the coefficient vector")
init <- if (use.costs) c(rep(init,nobs),
rep(0,3*nobs*nstates+3*nobs*nstates*ncoef+nobs*nstates*ntr))
else
c(rep(init,nobs),rep(0,2*nobs*nstates+2*nobs*nstates*ncoef))
res <- ode(y = init, times = t, func = dp, method=method, atol=atol, rtol=rtol,
...)
P <- array(res[,1+1:(nstates*nobs)],c(nt,nstates,nobs))
Pu <- array(res[,1+nobs*nstates+1:(nobs*nstates*ncoef)],c(nt,nstates,ncoef,nobs))
L <- array(res[,1+nstates*nobs+nstates*nobs*ncoef+1:(nstates*nobs)],c(nt,nstates,nobs))
Lu <- array(res[,1+2*nstates*nobs+nstates*nobs*ncoef+1:(nobs*nstates*ncoef)],c(nt,nstates,ncoef,nobs))
Pdisc <- array(res[,1+2*nstates*nobs+2*nstates*nobs*ncoef+1:(nstates*nobs)],c(nt,nstates,1,nobs))
state.names <- if(is.null(rownames(trans))) 1:nrow(trans) else rownames(trans)
dimnames(P) <- dimnames(L) <-
list(time=t, state= state.names, obs=rownames(newdata))
dimnames(Pdisc) <-
list(time=t, state= state.names, coef="Pdisc", obs=rownames(newdata))
coef.names <- unlist(lapply(x, function(xi) names(coef(xi))))
dimnames(Pu) <- dimnames(Lu) <-
list(time=t,
state=state.names,
coef=coef.names,
obs=rownames(newdata))
if (use.costs) {
costs <- array(res[,1+3*nstates*nobs+2*nstates*nobs*ncoef+1:(nstates*nobs)],c(nt,nstates,nobs))
costsu <- array(res[,1+4*nstates*nobs+2*nstates*nobs*ncoef+1:(nobs*nstates*ncoef)],c(nt,nstates,ncoef,nobs))
costsu.trans <- array(res[,1+4*nstates*nobs+3*nstates*nobs*ncoef+1:(nobs*nstates*ntr)],c(nt,nstates,ntr,nobs))
dimnames(costs) <- dimnames(P)
dimnames(costsu) <- dimnames(Pu)
dimnames(costsu.trans) <- list(time=t,
state=state.names,
trans=rownames(trans),
obs=rownames(newdata))
}
structure(list(time = t, P=P, Pu=Pu, L=L, Lu=Lu, Pdisc=Pdisc,
res=if (debug) res else NULL,
use.costs=use.costs,
costs=if (use.costs) costs else NULL,
costsu=if (use.costs) costsu else NULL,
costsu.trans=if (use.costs) costsu.trans else NULL,
vcov=Sigma, newdata=newdata, trans=trans,
utility.sd=utility.sd,
state.costs.sd=state.costs.sd,
transition.costs.sd=transition.costs.sd,
coefficients=coef(x),
call=call),
class="markov_msm")
}
## Pu, vcov
## Lu || Pdisc, bdiag(vcov,diag(utility.sd^2))
## costsu || Pdisc || costsu.trans, bdiag(vcov, diag(state.costs.sd^2), diag(transition.costs.sd^2))
bdiag <- function(...,initial=0) {
l <- list(...)
for (elt in l)
if (is.array(elt) && length(dim(elt))>2)
stop("not defined for arrays with more than 2 dimensions")
l <- lapply(l, as.matrix)
nr <- sapply(l,nrow)
nc <- sapply(l,ncol)
cumnr <- c(0,cumsum(nr))
cumnc <- c(0,cumsum(nc))
out <- matrix(initial,sum(nr),sum(nc))
for (i in 1:length(l))
out[(cumnr[i]+1):cumnr[i+1], (cumnc[i]+1):cumnc[i+1]] <- l[[i]]
## if rownames/colnames are all missing, then NULL else convert the NULLs to blank strings
replace <- function(l,f,len) {
lf <- lapply(l,f)
nulls <- sapply(lf, is.null)
if (all(nulls)) NULL
else unlist(lapply(1:length(l), function(i) if (nulls[i]) rep("",len(l[[i]])) else lf[[i]]))
}
rownames(out) <- replace(l,rownames, nrow)
colnames(out) <- replace(l,colnames, ncol)
out
}
## all(rownames(bdiag(matrix(1:4,2),c(a=5,b=6),7L,t(8:9),initial=0L))==c("","","a","b","",""))
## is.integer(bdiag(matrix(1:4,2),a=5:6,7L,t(8:9),initial=0L))
## bdiag(list(1:2),list(3:4+0),initial=list())
hazFun <- function(f, tmvar="t", ...) {
if(all(c("newdata","t") %in% names(formals(f))))
newf <- function(newdata) f(t=newdata[[tmvar]],newdata=newdata,...)
else if("t" %in% names(formals(f)))
newf <- function(newdata) f(t=newdata[[tmvar]], ...)
else if("newdata" %in% names(formals(f)))
newf <- function(newdata) f(newdata=newdata, ...)
else
newf <- function(newdata) f(...)
out <- list(haz=newf)
class(out) <- "hazFun"
out
}
predict.hazFun <- function(object, newdata, type=c("haz","gradh"), ...) {
type <- match.arg(type)
val <- if (type=="haz") object$haz(newdata) else 0
if (length(val)==1 && length(val)<nrow(newdata))
val <- rep(val,nrow(newdata))
val
}
coef.hazFun <- function(object, ...) c(hazFun=0)
vcov.hazFun <- function(object, ...) matrix(0,1,1,FALSE,list("hazFun","hazFun"))
splinefunx <- function(x, y, method="natural", constant.left=FALSE, constant.right=FALSE, ...) {
xstar <- x
ystar <- y
if (constant.right) {
xstar <- c(xstar,tail(x,1)+(1:10)*1e-7)
ystar <- c(ystar,rep(tail(y,1),10))
}
if (constant.left) {
xstar <- c(x[1]-(10:1)*1e-7,xstar)
ystar <- c(rep(y[1],10),ystar)
}
stats::splinefun(xstar, ystar, ..., method=method)
}
smoothpwc <- function(midts, rates, tmvar="t", offsetvar="", ...) {
log.smoother <- splinefunx(midts, log(rates), constant.right=TRUE)
haz <- function(newdata) {
t <- newdata[[tmvar]] + (if (offsetvar!="") newdata[[offsetvar]] else 0)
exp(log.smoother(t))
}
structure(list(haz=haz), class="smoothpwc")
}
coef.smoothpwc <- function(object, ...) c(smoothpwc=0)
vcov.smoothpwc <- function(object, ...) matrix(0,1,1,FALSE,list("pwc","pwc"))
predict.smoothpwc <- function(object, newdata, type=c("haz","gradh"), ...) {
type <- match.arg(type)
val <- if (type=="haz") object$haz(newdata) else 0
if (length(val)==1 && length(val)<nrow(newdata))
val <- rep(val,nrow(newdata))
val
}
as.data.frame.markov_msm <- function(x, row.names=NULL, optional=FALSE,
ci=TRUE,
P.conf.type="logit", L.conf.type="log", C.conf.type="log",
P.range=c(0,1), L.range=c(0,Inf),
C.range=c(0,Inf),
state.weights=NULL,
obs.weights=NULL,
...) {
if(!is.null(obs.weights))
x <- standardise(x, weights=obs.weights, ...)
if (!is.null(state.weights))
return(state_weights(x=x, row.names=row.names,
optional=optional, ci=ci,
P.conf.type=P.conf.type, L.conf.type=L.conf.type,
C.conf.type=C.conf.type,
P.range=P.range, L.range=L.range, C.range=C.range,
state.weights=state.weights, ...))
state.names <- rownames(x$trans)
perm <- function(x) aperm(x, c(2,3,1))
## probabilities
seP <- apply(x$Pu, c(1,2), function(m) sqrt(colSums(m * (vcov(x) %*% m)))) # c(nobs,nt,nstates)
seP <- array(seP, dim(x$Pu)[c(4,1,2)]) # c(nobs,nt,nstates))
seP <- aperm(seP,c(2,3,1))
## utilities
keep.ith <- function(a,i) {
slice <- a[,i,,]
a <- a*0
a[,i,,] <- slice
a
}
Lu <- do.call(abind, c(list(3, x$Lu),
lapply(1:nrow(x$trans), function(i) keep.ith(x$Pdisc,i))))
vcovL <- bdiag(vcov(x), diag(x$utility.sd^2))
seL <- apply(Lu, c(1,2), function(m) sqrt(colSums(m * (vcovL %*% m)))) # c(nobs,nt,nstates)
seL <- array(seL, dim(Lu)[c(4,1,2)]) # c(nobs,nt,nstates))
seL <- aperm(seL,c(2,3,1))
out <- data.frame(P=as.vector(perm(x$P)),
seP=as.vector(perm(seP)),
L=as.vector(perm(x$L)),
seL=as.vector(perm(seL)))
if (!is.null(x$costsu)) {
## error in the following
costsu <- do.call(abind, c(list(3,x$costsu),
lapply(1:nrow(x$trans), function(i) keep.ith(x$Pdisc,i)),
list(x$costsu.trans)))
vcovC <- bdiag(vcov(x), diag(x$state.costs.sd^2), diag(x$transition.costs.sd^2))
seC <- apply(costsu, c(1,2), function(m) sqrt(colSums(m * (vcovC %*% m)))) # c(nobs,nt,nstates)
seC <- array(seC, dim(costsu)[c(4,1,2)]) # c(nobs,nt,nstates))
seC <- aperm(seC,c(2,3,1))
out$C <- as.vector(perm(x$costs))
out$seC <- as.vector(perm(seC))
}
dimnames. <- dimnames(perm(x$P))
dimnames.$obs <- 1:length(dimnames.$obs)
y <- expand.grid(dimnames.)
y$time <- as.numeric(attr(y$time,"levels"))[y$time]
out <- cbind(y,out)
if (ci) {
tmp <- surv.confint(out$P,out$seP, conf.type=P.conf.type, min.value=P.range[1], max.value=P.range[2])
out$P.lower=tmp$lower
out$P.upper=tmp$upper
tmp <- surv.confint(out$L,out$seL, conf.type=L.conf.type, min.value=L.range[1], max.value=L.range[2])
out$L.lower=tmp$lower
out$L.upper=tmp$upper
if (!is.null(x$costs)) {
tmp <- surv.confint(out$C,out$seC, conf.type=C.conf.type, min.value=C.range[1], max.value=C.range[2])
out$C.lower=tmp$lower
out$C.upper=tmp$upper
}
}
newdata <- as.data.frame(x$newdata[as.numeric(out$obs), , drop=FALSE])
names(newdata) <- names(x$newdata)
out <- "rownames<-"(cbind(newdata, out),rownames(out))
if(!is.null(row.names)) rownames(out) <- row.names
out
}
# if var(P1)=g1*Sigma*g1' and var(P2)=g2*Sigma*g2' then var(P1+P2)=(g1+g2)*Sigma*(g1+g2)'
state_weights <- function(x, row.names=NULL, optional=FALSE,
ci=TRUE,
P.conf.type="logit", L.conf.type="log", C.conf.type="log",
P.range=c(0,1), L.range=c(0,Inf),
C.range=c(0,Inf),
state.weights,
...) {
stopifnot(!missing(state.weights))
state.names <- rownames(x$trans)
## stopifnot(length(state.names) == length(state.weights))
w <- function(m) colSums(m * state.weights)
perm <- function(x) aperm(x, c(2,3,1))
## probabilities
seP <- apply(x$Pu, c(1,4), function(m) sqrt(colSums(w(m) * (vcov(x) %*% w(m))))) # c(nobs,nt)
## utilities
keep.ith <- function(a,i) {
slice <- a[,i,,]
a <- a*0
a[,i,,] <- slice
a
}
Lu <- do.call(abind, c(list(3, x$Lu),
lapply(1:nrow(x$trans), function(i) keep.ith(x$Pdisc,i))))
vcovL <- bdiag(vcov(x), diag(x$utility.sd^2))
seL <- apply(Lu, c(1,4), function(m) sqrt(colSums(w(m) * (vcovL %*% w(m))))) # c(nobs,nt)
out <- data.frame(P=as.vector(apply(x$P,1,w)),
seP=as.vector(t(seP)),
L=as.vector(apply(x$L,1,w)),
seL=as.vector(t(seL)))
if (!is.null(x$costsu)) {
## error in the following
costsu <- do.call(abind, c(list(3,x$costsu),
lapply(1:nrow(x$trans), function(i) keep.ith(x$Pdisc,i)),
list(x$costsu.trans)))
vcovC <- bdiag(vcov(x), diag(x$state.costs.sd^2), diag(x$transition.costs.sd^2))
seC <- apply(costsu, c(1,4), function(m) sqrt(colSums(w(m) * (vcovC %*% w(m))))) # c(nobs,nt)
out$C <- as.vector(apply(x$costs,1,w))
out$seC <- as.vector(t(seC))
}
dimnames. <- dimnames(x$P)[c(3,1)]
dimnames.$obs <- 1:length(dimnames.$obs)
y <- expand.grid(dimnames.)
y$time <- as.numeric(attr(y$time,"levels"))[y$time]
out <- cbind(y,out)
if (ci) {
tmp <- surv.confint(out$P,out$seP, conf.type=P.conf.type, min.value=P.range[1], max.value=P.range[2])
out$P.lower=tmp$lower
out$P.upper=tmp$upper
tmp <- surv.confint(out$L,out$seL, conf.type=L.conf.type, min.value=L.range[1], max.value=L.range[2])
out$L.lower=tmp$lower
out$L.upper=tmp$upper
if (!is.null(x$costs)) {
tmp <- surv.confint(out$C,out$seC, conf.type=C.conf.type, min.value=C.range[1], max.value=C.range[2])
out$C.lower=tmp$lower
out$C.upper=tmp$upper
}
}
newdata <- as.data.frame(x$newdata[as.numeric(out$obs), , drop=FALSE])
names(newdata) <- names(x$newdata)
out <- "rownames<-"(cbind(newdata, out),rownames(out))
if(!is.null(row.names)) rownames(out) <- row.names
out
}
standardise <- function(x, ...)
UseMethod("standardise")
standardise.markov_msm <- function(x,
weights = rep(1,nrow(x$newdata)),
normalise = TRUE, ...) {
call <- match.call
state.names <- rownames(x$trans)
add.dim <- function(x) array(x,dim=c(dim(x),1))
if (length(weights)==0) weights <- 1
stopifnot(length(weights) %in% c(1,nrow(x$newdata)))
if (normalise) weights <- weights/sum(weights)
wtd.mean <- function(x) sum(x*weights)
Sigma <- vcov(x)
P <- add.dim(apply(x$P,c(1,2),wtd.mean))
Pu <- add.dim(apply(x$Pu,1:3,wtd.mean))
L <- add.dim(apply(x$L,c(1,2),wtd.mean))
Lu <- add.dim(apply(x$Lu,1:3,wtd.mean))
Pdisc <- add.dim(apply(x$Pdisc,1:3,wtd.mean))
dimnames(P) <- dimnames(L) <- list(time=x$time, state=state.names, obs=1)
oldnames <- dimnames(x$Pu)
dimnames(Pdisc) <- list(time=x$time, state=state.names, coef="Pdisc", obs=1)
dimnames(Pu) <- dimnames(Lu) <-
list(time=x$time, state=state.names, coef=oldnames[[3]], obs=1)
## update x
x$newdata <- x$newdata[0,]
x$P <- P
x$Pu <- Pu
x$L <- L
x$Lu <- Lu
x$Pdisc <- Pdisc
x$call <- call
x
}
vcov.survPen <- function(object) object$Vp
"coef<-.survPen" <- function(this,value) {
this$coefficients <- value
this
}
survPenWrap <- function(object) {
stopifnot(inherits(object, "survPen"))
if (!inherits(object, "survPenWrap"))
class(object) <- c("survPenWrap", class(object))
object
}
predict.survPenWrap <- function(object, newdata, type=NULL, min.eps=1e-6, ...) {
newdata[[object$t1.name]] <- pmax(min.eps, newdata[[object$t1.name]])
if (is.null(type)) NextMethod("predict", object)
else if (type=="haz") predict(object, newdata=newdata, do.surv=FALSE)$haz
else if (type=="gradh") predict(object, newdata=newdata, type="lpmatrix") *
predict(object, newdata=newdata, do.surv=FALSE)$haz
else NextMethod("predict", object)
}
predict.glm <- function (object, newdata=NULL, type=NULL, ...) {
if (is.null(type)) NextMethod("predict", object)
if(type=="gradh" && !(object$family$link %in% c("identity","log")))
stop("Currently only implemented for log and identity links")
## stopifnot() # Poisson family with log link?
## assumes response is a rate
if(!is.na(pmatch(type,"hazard")))
stats::predict.glm(object, newdata=newdata, type="response", ...)
else if (type=="gradh" && object$family$link == "log")
stats::predict.glm(object, newdata=newdata, type="response", ...) *
lpmatrix.lm(object, newdata=newdata)
else if (type=="gradh" && object$family$link == "identity")
lpmatrix.lm(object, newdata=newdata)
else NextMethod("predict", object)
}
gamWrap <- function(object) {
stopifnot(inherits(object, "gam"))
if (!inherits(object, "gamWrap"))
class(object) <- c("gamWrap", class(object))
object
}
predict.gamWrap <- function (object, newdata=NULL, type=NULL, ...) {
if (is.null(type)) NextMethod("predict", object)
if(type=="gradh" && object$family$link != "log")
stop("Currently only implemented for a log link")
## stopifnot() # Poisson family with log link?
## assumes response is a rate
if (type=="haz") predict(object, newdata=newdata, type="response", ...)
else if (type=="gradh") {
as.vector(predict(object, newdata=newdata, type="response", ...)) *
predict(object, newdata=newdata, type="lpmatrix", ...)
} else NextMethod("predict", object)
}
predict.flexsurvreg <- function(object, newdata, type=NULL, t=NULL, se.fit=FALSE, tmvar="t", ...) {
if (is.null(t)) t <- unique(newdata[[tmvar]])
if (is.null(type)) return(summary(object, newdata=newdata, t=t, ci=se.fit, ...))
if (type=="haz")
sapply(summary(object, newdata=newdata, t=t, ci=se.fit, type="haz"),
"[[", "est")
else if (type=="gradh") {
g <- function(coef) {
object$res[,"est"] <- object$res.t[,"est"] <- coef
predict(object, newdata, type="haz", t=t, ci=FALSE, ...)
}
## numDeriv::jacobian(g, x=coef(object), method="simple")
grad(g, coef(object))
}
else summary(object, newdata=newdata, type=type, t=t, ci=se.fit, ...)
}
vcov.aftreg <- function(object, ...)
object$var
predict.aftreg <- function (object, type = c("haz", "cumhaz", "density", "surv", "gradh"),
newdata, t=NULL, tmvar=NULL, na.action=na.pass, fd=FALSE, ...)
{
## utility function
mref <- function(m,i,j) i+(j-1)*nrow(m)
if (!inherits(object, "aftreg"))
stop("Works only with 'aftreg' objects.")
type <- match.arg(type)
if (type=="gradh" && fd) {
coef <- coef(object)
g1 <- function(coef) {
object$coefficients <- coef
predict.aftreg(object, newdata=newdata, type="haz", t=t,
tmvar=tmvar, na.action=na.action, fd=FALSE, ...)
}
## return(numDeriv::jacobian(g, coef))
## return(numDeriv::jacobian(g, coef, method="simple"))
return(grad(g1, coef))
}
if (is.null(t)) {
if (is.null(tmvar)) {
lhs <- object$call$formula[[2]]
expr <- if (length(lhs)==4) lhs[[3]] else lhs[[2]]
t <- eval(expr, newdata, parent.frame())
rm(lhs,expr)
} else
t <- newdata[[tmvar]]
}
if (length(t)<nrow(newdata)) t <- rep(t,length=nrow(newdata))
if (type=="haz" && fd) {
S <- predict.aftreg(object, newdata=newdata, type="surv", t=t,
tmvar=tmvar, na.action=na.action, ...)
g2 <- function(t)
predict.aftreg(object, newdata=newdata, type="surv", t=t,
tmvar=tmvar, na.action=na.action, fd=FALSE, ...)
## return(-numDeriv::jacobian(g, t)/S)
return(-grad(g2, t)/S)
}
## BUG: add names to object$levels
if (any(object$isF))
names(object$levels) <- object$covars[object$isF]
Terms <- object$terms
if (!inherits(Terms, "terms"))
stop("invalid terms component of object")
strats <- attr(Terms, "specials")$strata
intercept <- attr(Terms, "intercept")
Terms <- delete.response(Terms)
newframe <- stats::model.frame(Terms, data = newdata,
na.action = na.action,
xlev = object$levels)
if (length(strats)) {
if (!("strata" %in% names(object)))
stop("predictions with strata requires 'aftreg(..., x=TRUE)'")
temp <- untangle.specials(Terms, "strata", 1)
if (length(temp$vars) == 1)
strata.keep <- newframe[[temp$vars]]
else strata.keep <- survival::strata(newframe[, temp$vars], shortlabel = TRUE)
strata.keep <- factor(as.character(strata.keep),levels=levels(object$strata))
strata <- as.numeric(strata.keep)
## remove strats from object$terms
newTerms <- Terms[-temp$terms]
attr(newTerms, "intercept") <- attr(Terms, "intercept")
Terms <- newTerms
} else {
stopifnot(object$n.stata == 1)
strata <- rep(1, nrow(newframe))
}
"%mv%" <- function(a,b) as.vector(a %*% b)
newframe <- stats::model.frame(Terms, data = newdata,
na.action = na.action,
xlev = object$levels)
object$terms <- Terms
ncov <- length(object$means)
if (object$pfixed) {
p <- object$shape[strata]
lambda <- exp(object$coefficients[ncov + strata])
} else {
p <- exp(object$coefficients[ncov + strata * 2])
lambda <- exp(object$coefficients[ncov + strata * 2 - 1])
}
if (ncov) {
x <- model.matrix(object, newframe)[,-1,drop=FALSE] # is the intercept always first?
x <- sweep(x, 2, object$means)
param.scale <- if (object$param=="lifeAcc") -1 else 1
lambda <- lambda *
exp(param.scale*(x %mv% object$coefficients[1:ncov]))
} else {
x <- NULL
}
xx <- t
if (object$dist == "weibull") {
dist <- "Weibull"
haza <- eha::hweibull
Haza <- eha::Hweibull
Surviv <- stats::pweibull
Dens <- stats::dweibull
}
else if (object$dist == "loglogistic") {
dist <- "Loglogistic"
haza <- eha::hllogis
Haza <- eha::Hllogis
Surviv <- eha::pllogis
Dens <- eha::dllogis
}
else if (object$dist == "lognormal") {
dist = "Lognormal"
haza <- eha::hlnorm
Haza <- eha::Hlnorm
Surviv <- stats::plnorm
Dens <- stats::dlnorm
}
else if (object$dist == "ev") {
dist = "Extreme value"
haza <- eha::hEV
Haza <- eha::HEV
Surviv <- eha::pEV
Dens <- eha::dEV
}
else if (object$dist == "gompertz") {
dist = "Gompertz"
haza <- eha::hgompertz
Haza <- eha::Hgompertz
Surviv <- eha::pgompertz
Dens <- eha::dgompertz
## canonical parameterisation
p <- p/lambda
}
if (type=="haz")
return(haza(xx, scale = lambda, shape = p))
if (type=="cumhaz")
return(Haza(xx, scale = lambda, shape = p))
if (type=="density") {
if (object$dist == "lognormal") {
sdlog <- 1/p
meanlog <- log(lambda)
den <- stats::dlnorm(xx, meanlog, sdlog)
}
else
den <- Dens(xx, scale = lambda, shape = p)
return(den)
}
if (type=="surv") {
if (object$dist == "lognormal") {
sdlog <- 1/p
meanlog <- log(lambda)
sur <- stats::plnorm(xx, meanlog, sdlog, lower.tail = FALSE)
}
else {
sur <- Surviv(xx, scale = lambda, shape = p,
lower.tail = FALSE)
}
return(sur)
}
if (type=="gradh" && !fd) {
mapper <- list(ev="extreme",
weibull="weibull",
loglogistic="loglogistic",
Lognormal="lognormal")
if (object$dist %in% names(mapper)) {
distname <- mapper[[object$dist]]
dd <- survival::survreg.distributions[[distname]]
if (is.null(dd$itrans)) {
trans <- function(x) x
itrans <- function(x) x
dtrans <- function(x) 1
}
else {
trans <- dd$trans
itrans <- dd$itrans
dtrans <- dd$dtrans
}
if (!is.null(dd$dist))
dd <- survival::survreg.distributions[[dd$dist]]
pred <- log(lambda)
scale <- 1/p
u <- (trans(t)-pred) / scale # check dimensions
density <- dd$density(u, list()) # F, 1-F, f, f'/f, f''/f
S <- density[,2]
f <- density[,3]
f.prime <- density[,4]*f
## NOTE: x may be NULL
grad.beta <- cbind(x,-1)*(f.prime*S+f^2)*dtrans(t)/(S*scale)^2
out <- matrix(0,nrow(newdata),length(coef(object)))
nc <- ncol(grad.beta)
if (ncov)
out[,1:(nc-1)] <- grad.beta[,1:(nc-1)]
out[mref(out,1:nrow(newdata),nc+strata*2-2)] <- grad.beta[,nc]
if (!object$pfixed) {
grad.logscale <- (-f.prime*u*S/dtrans(t) -
f^2*u/dtrans(t) -
f*S/dtrans(t)) / (S*scale/dtrans(t))^2*scale
out[mref(out,1:nrow(newdata),nc+strata*2-1)] <- -grad.logscale
}
return(out)
} else {
## gompertz
h <- haza(xx, scale = lambda, shape = p)
grad.beta <- cbind(x,-1)*h*xx/p
grad.logscale <- h
out <- matrix(0,nrow(newdata),length(coef(object)))
nc <- ncol(grad.beta)
if (!is.null(x))
out[,1:(nc-1)] <- grad.beta*x
out[mref(out,1:nrow(newdata),nc+strata*2-2)] <- grad.beta[,nc]
out[mref(out,1:nrow(newdata),nc+strata*2-1)] <- -grad.logscale
return(out)
}
}
}
surv.confint <- function(p, se, conf.type=c("log","log-log","plain","logit","arcsin"),
conf.int=0.95, min.value=0, max.value=1) {
stopifnot(length(p) == length(se))
conf.type <- match.arg(conf.type)
zval <- qnorm(1 - (1 - conf.int)/2)
selog <- ifelse(p==0, NaN, se/p) # transform to selog
out <- if (conf.type == "plain") {
se2 <- se * zval
list(lower = pmax(p - se2, min.value), upper = pmin(p + se2, max.value))
}
else if (conf.type == "log") {
xx <- ifelse(p == 0, NaN, p)
selog2 <- zval * selog
temp1 <- exp(log(xx) - selog2)
temp2 <- exp(log(xx) + selog2)
list(lower = temp1, upper = pmin(temp2, max.value))
}
else if (conf.type == "log-log") {
xx <- ifelse(p == 0, NaN, p)
selog2 <- zval * selog/log(xx)
temp1 <- exp(-exp(log(-log(xx)) - selog2))
temp2 <- exp(-exp(log(-log(xx)) + selog2))
list(lower = temp1, upper = temp2)
}
else if (conf.type == "logit") {
xx <- ifelse(p == 0 | p == 1, NaN, p)
selog2 <- zval * selog * (1 + xx/(1 - xx))
temp1 <- 1 - 1/(1 + exp(log(xx/(1 - xx)) - selog2))
temp2 <- 1 - 1/(1 + exp(log(xx/(1 - xx)) + selog2))
list(lower = temp1, upper = temp2)
}
else if (conf.type == "arcsin") {
xx <- ifelse(p == 1, NaN, p)
selog2 <- 0.5 * zval * selog * sqrt(xx/(1 - xx))
list(lower = (sin(pmax(0, asin(sqrt(xx)) - selog2)))^2,
upper = (sin(pmin(pi/2, asin(sqrt(xx)) + selog2)))^2)
}
if (any(index <- is.na(p) | is.na(se)))
out$lower[index] <- out$upper[index] <- NA
out
}
print.markov_msm <- function(x,
digits=5,
se=FALSE, ci=TRUE,
P.conf.type="logit", L.conf.type="log",
C.conf.type="log",
P.range=c(0,1), L.range=c(0,Inf),
C.range=c(0,Inf),
...) {
df <- as.data.frame(x, ci=ci,
P.conf.type=P.conf.type, L.conf.type=L.conf.type,
C.conf.type=C.conf.type,
P.range=P.range, L.range=L.range, C.range=C.range)
if (!se) df$seP <- df$seL <- NULL
print(df, digits=digits, ...)
invisible(x)
}
zeroModel <- function(object) {
if (inherits(object, "zeroModel")) object
else structure(list(base=object), class="zeroModel")
}
predict.zeroModel <- function(object, type=c("haz","gradh"), ...) {
type <- match.arg(type)
0.0*predict(object$base, type=type, ...)
}
## we keep the same length/dim to support comparisons between different interventions
coef.zeroModel <- function(object, ...) 0.0*coef(object$base, ...)
vcov.zeroModel <- function(object, ...) 0.0*vcov(object$base, ...)
hrModel <- function(object, hr=1, ci=NULL, seloghr=NULL) {
stopifnot(is.null(ci) || (is.numeric(ci) && length(ci)==2))
stopifnot(is.numeric(hr) && length(hr)==1 && hr>0)
newobject <- if (inherits(object, "hrModel")) object else structure(list(base=object),
class="hrModel")
if (!is.null(ci) && is.null(seloghr))
seloghr <- log(ci[2]/ci[1])/2/qnorm(0.975)
attr(newobject, "loghr") <- log(hr)
attr(newobject, "seloghr") <- if (is.null(seloghr)) 0 else seloghr
newobject
}
predict.hrModel <- function(object, type=c("haz","gradh"), ...) {
type <- match.arg(type)
hr <- exp(attr(object,"loghr"))
pred1 <- predict(object$base, type="haz", ...)
if (type=="haz") hr*pred1
else
cbind(predict(object$base, type="gradh", ...)*hr, pred1*hr)
}
## This has different lengths/dimensions to the base model:
## wrap base intervention in hrModel(..., hr=1)
coef.hrModel <- function(object, ...) c(coef(object$base), attr(object,"loghr"))
vcov.hrModel <- function(object, ...) {
out <- rbind(cbind(vcov(object$base),0),0)
out[nrow(out),ncol(out)] <- attr(object, "seloghr")^2
out
}
## aftModel requires dh/dt - performed using finite differences
## Slow? Cool for completeness - but unlikely to be used
## axiom: h := operator 'h
## axiom: hstar := h(t*exp(eta),beta)*exp(eta)
## axiom: D(hstar,beta)
## axiom: D(hstar,eta)
## h(t)=h(exp(u)) where t=exp(u)
## => dh/dt = dh/du * du/dt where du/dt=1/t
## and dh/du ~= (h(t*exp(eps)) - h(t/exp(eps)))/2/eps
aftModel <- function(object, af=1, ci=NULL, selogaf=NULL) {
stopifnot(is.null(ci) || (is.numeric(ci) && length(ci)==2))
stopifnot(is.numeric(af) && length(af)==1 && af>0)
if (!is.null(ci) && is.null(selogaf))
selogaf <- log(ci[2]/ci[1])/2/qnorm(0.975)
newobject <- if (inherits(object, "aftModel")) object else structure(list(base=object),
class="aftModel")
attr(newobject, "logaf") <- log(af)
attr(newobject, "selogaf") <- if (is.null(selogaf)) 0 else selogaf
newobject
}
predict.aftModel <- function(object, newdata, type=c("haz","gradh"), tmvar=NULL, eps=1e-5, ...) {
type <- match.arg(type)
af <- exp(attr(object,"logaf"))
time <- newdata[[tmvar]]
timeStar <- newdata[[tmvar]] <- af*time
pred1 <- predict(object$base, type="haz", newdata=newdata, ...)
if (type=="haz") af*pred1
else {
newdata2 <- newdata
newdata2[[tmvar]] <- timeStar*exp(eps)
predu <- predict(object$base, type="haz", newdata=newdata2, ...)
newdata2[[tmvar]] <- timeStar/exp(eps)
predl <- predict(object$base, type="haz", newdata=newdata2, ...)
dhdt <- (predu-predl)/2/eps/timeStar
cbind(predict(object$base, type="gradh", newdata=newdata, ...)*af,
timeStar*af*dhdt + af*pred1)
}
}
## This has different lengths/dimensions to the base model:
## wrap base intervention in aftModel(..., af=1)
coef.aftModel <- function(object, ...) c(coef(object$base), attr(object,"logaf"))
vcov.aftModel <- function(object, ...) {
out <- rbind(cbind(vcov(object$base),0),0)
out[nrow(out),ncol(out)] <- attr(object, "selogaf")^2
out
}
## h'(t):
## gsm: S(t|x)=G(eta(t,x)) => H(t|x)=-log(S(t|x))
## => h(t|x)= -G'(eta(t,x))/G(eta(t,x))*eta_{,1}(t,x)
## => h'(t|x) = mess
## export
splineFun <- function(time, rate, method="natural", scale=1, ...) {
fun <- stats::splinefun(time, log(rate), method=method, ...)
hazFun(function(t) scale*exp(fun(t)))
}
## export
addModel <- function (...) {
structure (as.list (...), class="addModel")
}
predict.addModel <- function (object, ...) {
Reduce("+", lapply(object, predict, ...))
}
coef.addModel <- function (object)
do.call(c,lapply(object, coef))
vcov.addModel <- function (object)
do.call(bdiag,lapply(object, vcov))
subset.markov_msm <- function(x, subset, ...) {
e <- substitute(subset)
r <- eval(e, x$newdata, parent.frame())
if (!is.logical(r))
stop("'subset' must be logical")
r <- r & !is.na(r)
# update x
x$P <- x$P[,,r,drop=FALSE]
x$Pu <- x$Pu[,,,r,drop=FALSE]
x$L <- x$L[,,r,drop=FALSE]
x$Lu <- x$Lu[,,,r,drop=FALSE]
x$Pdisc <- x$Pdisc[,,,r,drop=FALSE]
x$newdata <- x$newdata[r,,drop=FALSE]
if (x$use.costs) {
x$costs <- x$costs[,,r,drop=FALSE]
x$costsu <- x$costsu[,,,r,drop=FALSE]
}
x$res <- NULL # not strictly needed
x$call <- match.call()
x
}
diff.markov_msm <- function(x, y, ...) {
stopifnot(inherits(x,"markov_msm"))
stopifnot(inherits(y,"markov_msm"))
## vcov, time, trans should all be the same
stopifnot(all(x$vcov == y$vcov))
stopifnot(all(x$time == y$time))
stopifnot(all(x$trans == y$trans, na.rm=TRUE))
z <- x # copy of x
z$P <- x$P-y$P
z$Pu <- x$Pu-y$Pu
z$L <- x$L-y$L
z$Lu <- x$Lu-y$Lu
z$Pdisc <- x$Pdisc-y$Pdisc
z$call <- match.call()
z$newdata[x$newdata != y$newdata] <- NA
if (z$use.costs) {
z$costs <- x$costs-y$costs
z$costsu <- x$costsu-y$costsu
}
class(z) <- c("markov_msm_diff","markov_msm") # not strictly "markov_msm"...
z
}
as.data.frame.markov_msm_diff <- function(x, row.names=NULL, optional=FALSE,
P.conf.type="plain", L.conf.type="plain",
C.conf.type="plain",
P.range=c(-Inf, Inf), L.range=c(-Inf, Inf),
C.range=c(-Inf,Inf), ...)
as.data.frame.markov_msm(x,
P.conf.type=P.conf.type, L.conf.type=L.conf.type,
C.conf.type=C.conf.type,
P.range=P.range, L.range=L.range,
C.range=C.range, ...)
plot.markov_msm_diff <- function(x, y, ...)
stop("Plots for markov_msm_diff have not yet been implemented")
plot.markov_msm <- function(x, y, stacked=TRUE, which=c("P","L"),
xlab="Time", ylab=NULL, col=2:6, border=col,
ggplot2=FALSE, lattice=FALSE, alpha=0.2,
strata=NULL,
...) {
stopifnot(inherits(x,"markov_msm"))
which <- match.arg(which)
## ylab defaults
if(is.null(ylab)) {
ylab <- if(which=='P') "Probability" else "Length of stay"
if(inherits(x,"markov_msm_diff"))
ylab <- if(which=='P') "Difference in probabilities"
else "Difference in lengths of stay"
if(inherits(x,"markov_msm_ratio"))
ylab <- if(which=='P') "Ratio of probabilities"
else "Ratio of lengths of stay"
}
if (nrow(x$newdata)>1) {
warning("More than one set of covariates; covariates have been standardised")
x <- standardise(x)
}
if (ggplot2)
ggplot.markov_msm(x, which=which, stacked=stacked, xlab=xlab, ylab=ylab,
alpha=alpha, ...)
else if (lattice)
xyplot.markov_msm(x, which=which, stacked=stacked, xlab=xlab, ylab=ylab,
col=col, border=border, strata=strata, ...)
else {
if (!missing(y)) warning("y argument is ignored")
df <- as.data.frame(x)
states <- unique(df$state)
if (stacked) {
if (which == "L")
stop("Stacked plot is not sensible for length of stay")
out <- graphics::plot(range(x$time),0:1, type="n", xlab=xlab, ylab=ylab, ...)
lower <- 0
for (i in length(states):1) { # put the last state at the bottom
df2 <- df[df$state==states[i],]
if (length(lower)==1) lower <- rep(0,nrow(df2))
upper <- lower+df2[[which]]
graphics::polygon(c(df2$time,rev(df2$time)), c(lower,rev(upper)),
border=border[i], col=col[i])
lower <- upper
}
graphics::box()
invisible(out)
}
else stop('Unstacked plot not implemented in base graphics; use ggplot2=TRUE or lattice=TRUE')
}
}
ggplot.markov_msm <- function(data, mapping=NULL,
which=c("P","L"),
stacked = TRUE, alpha=0.2,
xlab=NULL, ylab=NULL,
..., environment=parent.frame()) {
if (requireNamespace("ggplot2", quietly=TRUE)) {
which <- match.arg(which)
df <- as.data.frame(data, ci=!stacked)
if (stacked)
ggplot2::ggplot(df, if(is.null(mapping))
ggplot2::aes_string(x='time', y=which, fill='state')
else mapping) +
ggplot2::geom_area() +
ggplot2::xlab(xlab) + ggplot2::ylab(ylab)
else {
lower <- paste0(which,'.lower')
upper <- paste0(which,'.upper')
ggplot2::ggplot(df, if (is.null(mapping))
ggplot2::aes_string(x='time', y=which, ymin=lower,
ymax=upper)
else mapping) +
ggplot2::geom_line() + ggplot2::geom_ribbon(alpha=alpha) +
ggplot2::facet_grid(stats::reformulate(".","state")) +
ggplot2::xlab(xlab) + ggplot2::ylab(ylab)
}
}
else stop("Suggested package 'ggplot2' is not available")
}
xyplot.markov_msm <- function(x, data, strata=NULL,
which=c("P","L"),
stacked = TRUE,
col=2:6, border=col,
..., environment=parent.frame()) {
if (requireNamespace("lattice", quietly=TRUE)) {
which <- match.arg(which)
## if (!missing(data) && is.null(strata)) strata <- data # copy from data to strata
df <- as.data.frame(x, ci=!stacked)
if (!is.factor(df$state)) df$state <- as.factor(df$state)
states <- levels(df$state)
if (stacked) {
rhs_string <- if (!is.null(strata)) paste('time',deparse(rhs(strata)),sep='|')
else 'time'
lattice::xyplot(stats::reformulate(rhs_string,which), data=df,
panel = function(x, y, subscripts, group, ...) {
df2 <- df[subscripts,]
lattice::panel.xyplot(x,y, type="n")
lower <- 0
for (i in length(states):1) {
df3 <- df2[df2$state==states[i], , drop=FALSE]
if (length(lower)==1) lower <- rep(0,nrow(df3))
upper <- lower+df3[[which]]
lattice::panel.polygon(c(df3$time,rev(df3$time)),
c(lower,rev(upper)),
border=border[i], col=col[i])
lower <- upper
}
},
subscripts=TRUE, ...)
}
else {
lower <- paste0(which,'.lower')
upper <- paste0(which,'.upper')
rhs_string <- if (!is.null(strata)) paste('time|state',deparse(rhs(strata)),sep='+')
else 'time|state'
lattice::xyplot(stats::reformulate(rhs_string,which), data=df,
panel = function(x, y, subscripts, group, ...) {
lattice::panel.xyplot(x,y, type="n")
lattice::panel.polygon(c(x,rev(x)),
c(df[[lower]][subscripts],
rev(df[[upper]][subscripts])),
border=NULL, col="grey")
lattice::panel.lines(x,y)
},
subscripts=TRUE, ...)
}
}
else stop("Suggested package 'lattice' is not available")
}
## g(beta) = log(losA) - log(losB)
## g'(beta) = losA'/losA - losB'/losB # which is not defined when los=0
## NB: data are on a log scale!
ratio_markov_msm <- function(x, y, ...) {
stopifnot(inherits(x,"markov_msm"))
stopifnot(inherits(y,"markov_msm"))
## vcov, time and trans should all be the same
stopifnot(all(x$vcov == y$vcov))
stopifnot(all(x$time == y$time))
stopifnot(all(x$trans == y$trans, na.rm=TRUE))
z <- x # copy of x
z$P <- log(x$P/y$P)
z$L <- log(x$L/y$L) # logs of the ratios
z$Pu <- apply(x$Pu, 3, function(slice) slice/x$P) - apply(y$Pu, 3, function(slice) slice/y$P)
dim(z$Pu) <- dim(x$Pu)
dimnames(z$Pu) <- dimnames(x$Pu)
z$Lu <- apply(x$Lu, 3, function(slice) slice/x$L) - apply(y$Lu, 3, function(slice) slice/y$L)
dim(z$Lu) <- dim(x$Lu)
dimnames(z$Lu) <- dimnames(x$Lu)
z$Pdisc <- log(x$Pdisc/y$Pdisc)
z$call <- match.call()
z$newdata[x$newdata != y$newdata] <- NA
class(z) <- c("markov_msm_ratio","markov_msm") # not strictly "markov_msm"...
z
}
as.data.frame.markov_msm_ratio <- function(x, row.names=NULL, optional=FALSE, ...) {
## data are on a log scale!!
z <- as.data.frame.markov_msm_diff(x, ...)
names <- c("P","L","P.lower","P.upper","L.lower","L.upper")
if (!is.null(x$costs))
names <- c(names, c("C","C.lower","C.upper"))
for (name in names)
z[[name]] <- exp(z[[name]])
z
}
## prev_markov_msm <- function(x, w, ...) {
## stopifnot(inherits(x,"markov_msm"))
## stopifnot(nrow(x$trans) == length(w))
## stopifnot(all(w %in% 0:1))
## x <- msm2
## w <- c(1,1,0)
## pi <- piL <- P <- L <- array(0,dim(x$P))
## sumP <- sumL <- array(0,dim(x$P)[-2])
## for (i in 1:length(w)) {
## P[,i,] <- x$P[,i,]*w[i]
## L[,i,] <- x$L[,i,]*w[i]
## sumP <- sumP + P[,i,]
## sumL <- sumL + L[,i,]
## }
## for (i in 1:length(w)) {
## pi[,i,] <- P[,i,]/sumP
## piL[,i,] <- P[,i,]/sumL
## }
## dw <- diag(w)
## logit <- stats::poisson()$linkinv
## z <- list(P=logit(P), L=log(L)) # logit and log
## z$Pu <- apply(x$Pu, 3, function(slice) slice/x$P) - apply(y$Pu, 3, function(slice) slice/y$P)
## dim(z$Pu) <- dim(x$Pu)
## dimnames(z$Pu) <- dimnames(x$Pu)
## z$Lu <- apply(x$Lu, 3, function(slice) slice/x$L) - apply(y$Lu, 3, function(slice) slice/y$L)
## dim(z$Lu) <- dim(x$Lu)
## dimnames(z$Lu) <- dimnames(x$Lu)
## z <- c(list(time=x$time,
## vcov=x$vcov,
## trans=x$trans,
## res=NULL),
## z)
## z$call <- match.call()
## z$newdata <- x$newdata
## z$newdata[x$newdata != y$newdata] <- NA
## class(z) <- c("markov_msm_prev","markov_msm") # not strictly "markov_msm"...
## z
## }
## as.data.frame.markov_msm_prev <- function(x, ...) {
## ## data are on a log scale!!
## z <- as.data.frame.markov_msm_diff(x, ...)
## for (name in c("P","L","P.lower","P.upper","L.lower","L.upper"))
## z[[name]] <- exp(z[[name]])
## z
## }
abind <- function(index, ...) { # bind on nth slice for a bag of arrays
x <- list(...)
stopifnot(all(sapply(x,is.array)))
dim1 <- dim(x[[1]])
## check index
stopifnot(index==floor(index) && index>=1 || index<=length(dim))
## check all arrays have the same number of dimensions
stopifnot(all(sapply(x[-1],function(xi) length(dim(xi))==length(dim1))))
## case: arrays of dimension 1
if (length(dim1)==1) return(array(unlist(x)))
## check all arrays for the same dimensions (excluding the index dimension)
for (i in ((1:length(dim1))[-index]))
stopifnot(all(sapply(x[-1],function(xi) dim(xi)[i]==dim1[i])))
## re-order to put the index dimension last
ord <- c((1:length(dim1))[-index],index)
x <- lapply(x, function(xi) aperm(xi,ord)) # REPLACEMENT
y <- do.call("c",lapply(x,as.vector)) # magic...
dims <- dim1[ord]
dimnames. <- dimnames(x[[1]])
dims[length(dims)] <- sum(sapply(x,function(xi) dim(xi)[length(dims)]))
if (!is.null(dimnames.[[length(dims)]]))
dimnames.[[length(dims)]] <- unlist(lapply(x,function(xi) dimnames(xi)[length(dims)]))
dim(y) <- dims
dimnames(y) <- dimnames.
aperm(y,order(ord))
}
## abind(1,array(1),array(2),array(3:4))
## abind(1,array(1,c(1,1)),array(2,c(1,1)),array(3:4,c(2,1)))
## abind(2,array(1,c(1,1)),array(2,c(1,1)),array(3:4,c(1,2)))
## abind(1,array(1,c(1,1,1)),array(2,c(1,1,1)),array(3:4,c(2,1,1)))
add_dim <- function(a,index,dimname="") {
stopifnot(is.array(a))
dims <- dim(a)
ndim <- length(dims)
dimnames. <- dimnames(a)
stopifnot(index==floor(index) && index>=1 && index<=ndim+1)
newdim <- if (index==1) c(1,dims)
else if (index==ndim+1) c(dims,1)
else c(dims[1:(index-1)], 1, dims[index:ndim])
dimnames. <- if (index==1) c(dimname,dimnames.)
else if (index==ndim+1) c(dimnames.,dimname)
else c(dimnames.[1:(index-1)], dimname, dimnames.[index:ndim])
array(a,newdim,dimnames.)
}
## add_dim(array(c(1,2,3,4),dimnames=list(11:14)),2,"blank")
## add_dim(array(1:4),2)
## add_dim(array(1:4,c(2,2)),1)[1,,]
## add_dim(array(1:4,c(2,2)),2)[,1,]
## add_dim(array(1:4,c(2,2)),3)[,,1]
bindlast <- function(...) { # bind on last slice for a bag of arrays
x <- list(...)
stopifnot(all(sapply(x,is.array)))
dim1 <- dim(x[[1]])
stopifnot(all(sapply(x[-1],function(xi) length(dim(xi))==length(dim1))))
for (i in 1:(length(dim1)-1))
stopifnot(all(sapply(x[-1],function(xi) dim(xi)[i]==dim1[i])))
y <- do.call("c",lapply(x,as.vector))
dims <- dim1
dimnames. <- dimnames(x[[1]])
dims[length(dims)] <- sum(sapply(x,function(xi) dim(xi)[length(dims)]))
if (!is.null(dimnames.[[length(dims)]]))
dimnames.[[length(dims)]] <- unlist(lapply(x,function(xi) dimnames(xi)[length(dims)]))
dim(y) <- dims
dimnames(y) <- dimnames.
y
}
rbind.markov_msm <- function(..., deparse.level = 1) {
x <- list(...)
stopifnot(all(sapply(x,inherits,"markov_msm")))
## class, vcov, time and trans should all be the same
stopifnot(all(sapply(x[-1],function(xi) all(class(xi) == class(x[[1]])))))
stopifnot(all(sapply(x[-1],function(xi) all(xi$vcov==x[[1]]$vcov))))
stopifnot(all(sapply(x[-1],function(xi) all(xi$time==x[[1]]$time))))
stopifnot(all(sapply(x[-1],function(xi) all(xi$trans==x[[1]]$trans, na.rm=TRUE))))
bind <- function(name) do.call(bindlast, lapply(x,"[[",name))
z <- x[[1]]
z$P <- bind("P")
z$Pu <- bind("Pu")
z$L <- bind("L")
z$Lu <- bind("Lu")
z$Pdisc <- bind("Pdisc")
z$call <- match.call()
z$newdata <- do.call(rbind,lapply(x,"[[", "newdata"))
## is this a good idea? It will be lost if done more than once...
z$newdata$.index <- unlist(lapply(1:length(x),
function(i) rep(i,nrow(x[[i]]$newdata))))
## state.names <- rownames(z$trans)
z
}
transform.markov_msm <- function(`_data`, ...) {
`_data`$newdata <- transform(`_data`$newdata, ...)
`_data`
}
reorder <- function(x,order) {
l <- attr(x,"levels")
factor(l[x],l[order])
}
nrow.markov_msm <- function(x, ...) nrow(as.data.frame(x, ...))
vcov.markov_msm <- function(object, ...) object$vcov
## Collapse across states
## issue with variance calculations...
collapse_markov_msm <- function(object, which=NULL, sep="; ") {
## example: which=c(1,2) => combine states 1 and 2
stopifnot(inherits(object, "markov_msm"))
stopifnot(!is.null(which))
if(is.character(which))
which <- pmatch(which, rownames(object$trans))
stopifnot(all(which %in% 1:nrow(object$trans)))
stopifnot(nrow(object$trans)>=2)
call <- match.call()
which <- sort(unique(which)) # in case of duplicates
if(length(which)==1) return(object) # no change
## algorithm:
n.newstates <- nrow(object$trans) - length(which) + 1
## initalise
newstates <- vector("list", n.newstates)
j <- 0
base <- NULL
for (i in 1:nrow(object$trans)) {
if (i %in% which) {
if (is.null(base)) {
j <- j + 1
base <- j
}
newstates[[base]] <- c(newstates[[base]],i)
}
else {
j <- j + 1
newstates[[j]] <- i
}
}
index <- vector('numeric',nrow(object$trans))
for (j in 1:length(newstates))
for (k in newstates[[j]])
index[k] <- j
state.names <- rownames(object$trans)
if (is.null(state.names)) state.names <- 1:nrow(object$trans)
new.state.names <- tapply(state.names, index, paste0, collapse=sep)
trans <- object$trans
for (i in rev(which[-1])) {
trans <- trans[-i,,drop=FALSE]
trans <- trans[,-i,drop=FALSE]
}
## sum <- function(x) tapply(x, index, base::sum)
sum3 <- function(x,index,DIM=2) {
dims <- dim(x)
dims[DIM] <- length(unique(index))
y <- array(0,dims)
dimnames <- dimnames(x)
dimnames[[DIM]] <- new.state.names
for (i in 1:length(index)) {
if (DIM==1)
y[index[i],,] <- y[index[i],,] + x[i,,]
else if (DIM==2)
y[,index[i],] <- y[,index[i],] + x[,i,]
}
dimnames(y) <- dimnames
y
}
sum4 <- function(x,index) {
## assumes DIM=2
dims <- dim(x)
dims[2] <- length(unique(index))
dimnames <- dimnames(x)
dimnames[[2]] <- new.state.names
y <- array(0,dims)
for (i in 1:length(index))
y[,index[i],,] <- y[,index[i],,] + x[,i,,]
dimnames(y) <- dimnames
y
}
P <- sum3(object$P,index)
L <- sum3(object$L,index)
Pu <- sum4(object$Pu,index)
Lu <- sum4(object$Lu,index)
Pdisc <- sum4(object$Pdisc,index)
if ("C" %in% names(object)) {
C <- sum3(object$C,index)
Cu <- sum4(object$Cu,index)
}
index <- 1
for (i in 1:nrow(object$trans)) {
if (any(i==trans,na.rm=TRUE)) {
trans[which(i==trans)] <- index
index <- index + 1
}
}
colnames(trans) <- rownames(trans) <- new.state.names
z <- object
z$P <- P
z$L <- L
z$Pu <- Pu
z$Lu <- Lu
z$Pdisc <- Pdisc
z$trans <- trans
## costs??
z$call <- call
z
}
stratifiedModel <- function (object,strata) {
## if(inherits (object, "stratifiedModel"))
## return(object)
if(! inherits (object, "stratifiedModel"))
class (object) <- c("stratifiedModel", class (object))
object$strata.name <- substitute (strata)
stopifnot(is.name(object$strata.name))
object$strata.index <- NULL
object
}
"[[.stratifiedModel" <- function (object, index) {
object$strata.index <- index
object
}
## does the following only work for index 1..K?
predict.stratifiedModel <- function (object, type, newdata, ...) {
if (! is.null(object$strata.index))
newdata[[object$strata.name]] <- object$strata.index
## otherwise assume that the strata is specified in newdata
NextMethod ("predict", object)
}
## predict.test <- function(object, type, newdata) return(newdata$strata)
## m <- stratifiedModel("class<-"(list(),"test"), strata)
## predict(m[[2]],"",data.frame(x=1))
coef.stratifiedModel <- function(object, ...) NextMethod("coef", object)
vcov.stratifiedModel <- function(object, ...) NextMethod("vcov", object)
vsplinefun <- function(x,y,...) {
stopifnot(is.matrix(y))
splinefuns <- lapply(1:ncol(y), function(j) stats::splinefun(x,y[,j],...))
function(x) sapply(splinefuns, function(obj) obj(x))
}
makeSplineFun <- function(object,tm,newdata=data.frame(one=1),tmvar="", min.tm=1e-8, log=FALSE) {
tm <- pmax(tm,min.tm)
trans <- if (log) base::log else base::identity
## itrans <- if (log) base::exp else base::identity
newdata <- as.data.frame(newdata)
times <- data.frame(t=tm); if (tmvar != "") names(times) <- tmvar
hazard <- lapply(1:nrow(newdata), function(i)
splinefun(trans(tm),
predict(object,type="hazard",newdata=merge(newdata[i,,drop=FALSE], times))))
gradh <- lapply(1:nrow(newdata), function(i)
vsplinefun(trans(tm),
predict(object,type="gradh",newdata=merge(newdata[i,,drop=FALSE], times))))
out <- list(
object=object,
hazard=function(t) sapply(hazard, function(obj) obj(trans(t))),
gradh=function(t) t(sapply(gradh, function(obj) obj(trans(t)))))
class(out) <- "SplineFun"
out
}
predict.SplineFun <- function(object, newdata, type=c("hazard","gradh"), tmvar="", ...) {
type <- match.arg(type)
time <- newdata[[tmvar]]
stopifnot(all(time[1] == time[-1]))
if (type=="hazard") object$hazard(time[1]) else object$gradh(time[1])
}
vcov.SplineFun <- function(object) vcov(object$object)
coef.SplineFun <- function(object) coef(object$object)
## Non-parametric baseline: SDE approach due to Ryalen and colleagues
markov_sde <- function(models, trans, newdata, init=NULL, nLebesgue=1e4+1, los=FALSE, nOut=300,
weights=1) {
transfun <- function(tmat) {
indices <- sort(as.vector(tmat)); indices <- setdiff(indices,NA)
nStates <- nrow(tmat)
out <- do.call(rbind,
lapply(indices, function(i) {
index2 <- which(tmat == i)
from <- (index2-1) %% nStates +1
to <- (index2-1) %/% nStates + 1
data.frame(from=from,to=to)
}))
matrix(as.integer(as.matrix(out)),nrow(out))-1L
}
## TODO check parameters
nStates <- nrow(trans)
nTrans <- sum(!is.na(trans))
if (is.null(init)) {
init <- c(1,rep(0,nStates-1))
}
stopifnot(length(init)==nStates)
## init <- rep(init,nrow(newdata))
n <- sum(sapply(models,attr,"orig.max.clust"))
cumHazList <- lapply(models, function(object)
-t(log(predict(object, newdata=newdata, se=FALSE)$S0)))
timesList <- lapply(models, function(model) model$cum[,1])
eventTimes <- times <- sort(unique(unlist(timesList)))
hazList <- lapply(1:nrow(newdata), function(i) {
hazMatrix <- matrix(0,nrow=nTrans,ncol=length(times))
for (j in 1:nTrans) {
hazMatrix[j,match(timesList[[j]],times)] <- diff(c(0,cumHazList[[j]][,i]))
}
hazMatrix
})
hazMatrix <- do.call(rbind,hazList)
if (length(weights)==1 && weights != 0)
weights <- rep(weights,nrow(newdata))/(weights*nrow(newdata))
vcov <- matrix(1,1,1)
if (!los) {
out <- .Call("plugin_P_by", n, nrow(newdata), hazMatrix, init, transfun(trans), weights,
vcov,
PACKAGE="rstpm2")
out$P <- out$X
out$P.se <- sqrt(out$variance)
} else {
out <- .Call("plugin_P_L_by",
n, nrow(newdata), hazMatrix, init, transfun(trans), times, weights, nOut,
vcov, nLebesgue,
PACKAGE="rstpm2")
PIndex <- 1:(nrow(out$X)/2)
tr <- function(x) array(as.vector(x), dim=c(nStates,nrow(newdata),ncol(x)))
out$P <- tr(out$X[PIndex,])
out$L <- tr(out$X[-PIndex,])
out$P.se <- sqrt(tr(out$variance[PIndex,]))
out$L.se <- sqrt(tr(out$variance[-PIndex,]))
}
out$n <- n
out$times <- if (los) out$time else times
out$newdata <- newdata
out$trans <- trans
out$los <- los
out$init <- init
out$weights <- weights
class(out) <- "markov_sde"
if (!all(weights==0)) {
stand <- out # copy! This may be a bad idea...
stand$newdata <- newdata
if (!los) {
stand$P <- stand$Y
stand$P.se <- sqrt(stand$varY)
} else {
PIndex <- 1:(nrow(out$Y)/2)
stand$P <- stand$Y[PIndex,]
stand$L <- stand$Y[-PIndex,]
stand$P.se <- sqrt(stand$varY[PIndex,])
stand$L.se <- sqrt(stand$varY[-PIndex,])
}
## tidy up
stand$X <- stand$Y
stand$variance <- stand$varY
out$X <- out$variance <- out$Y <- out$varY <- stand$Y <- stand$varY <- NULL
stand$newdata <- stand$newdata[1,,drop=FALSE]
class(stand) <- "markov_sde"
out$stand <- stand
}
out
}
standardise.markov_sde <- function(x, ...) {
x$stand
}
plot.markov_sde <- function(x, y, stacked=TRUE, which=c("P","L"), index=NULL,
xlab="Time", ylab=NULL, col=2:6, border=col,
ggplot2=FALSE, lattice=FALSE, alpha=0.2,
strata=NULL,
...) {
stopifnot(inherits(x,"markov_sde"))
which <- match.arg(which)
if (!missing(y)) warning("y argument is ignored")
## ylab defaults
if (is.null(ylab))
ylab <- if(which=='P') "Probability" else "Length of stay"
if (ggplot2)
ggplot.markov_msm(x, which=which, stacked=stacked, xlab=xlab, ylab=ylab,
alpha=alpha, ...)
else if (lattice)
xyplot.markov_msm(x, which=which, stacked=stacked, xlab=xlab, ylab=ylab,
col=col, border=border, strata=strata, ...)
else {
if (is.null(index) && nrow(x$newdata)>1) {
warning("More than one set of covariates; defaults to weighted estimator")
x <- x$stand # Warning: replacement
index <- 1
}
if (is.null(index)) index <- 1
df <- merge(x$newdata[index,,drop=FALSE], as.data.frame(x))
states <- unique(df$state)
if (stacked) {
out <- graphics::plot(range(x$times, na.rm=TRUE),0:1, type="n", xlab=xlab, ylab=ylab, ...)
lower <- 0
for (i in length(states):1) { # put the last state at the bottom
df2 <- df[df$state==states[i],]
if (length(lower)==1) lower <- rep(0,nrow(df2))
upper <- lower+df2[[which]]
graphics::polygon(c(df2$time,rev(df2$time)), c(lower,rev(upper)),
border=border[i], col=col[i])
lower <- upper
}
graphics::box()
invisible(out)
}
else stop('Unstacked plot not implemented in base graphics; use ggplot2=TRUE or lattice=TRUE')
}
}
as.data.frame.markov_sde <- function(x, row.names=NULL, optional=NULL, ci=TRUE,
P.conf.type="logit", L.conf.type="log",
P.range=c(0,1), L.range=c(0,Inf),
...) {
if (any(x$weights<0)) {
P.conf.type <- L.conf.type <- "plain"
P.range <- L.range <- c(-Inf,Inf)
}
.id. <- 1:nrow(x$newdata)
nStates <- nrow(x$trans)
state.names <- rownames(x$trans)
stateNames <- if (!is.null(rownames(x$trans))) rownames(x$trans) else 1:nrow(x$trans)
out <- expand.grid(state=stateNames, .id.=.id., time=x$times)
out <- cbind(x$newdata[out$.id.,],out)
names(out)[1:ncol(x$newdata)] <- colnames(x$newdata)
out$P <- as.vector(x$P)
out$P.se <- as.vector(x$P.se)
if (ci) {
tmp <- surv.confint(out$P,out$P.se, conf.type=P.conf.type, min.value=P.range[1], max.value=P.range[2])
out$P.lower <- tmp$lower
out$P.upper <- tmp$upper
}
if (x$los) {
out$L <- as.vector(x$L)
out$L.se <- as.vector(x$L.se)
if (ci) {
tmp <- surv.confint(out$L,out$L.se, conf.type=L.conf.type, min.value=L.range[1], max.value=L.range[2])
out$L.lower <- tmp$lower
out$L.upper <- tmp$upper
}
}
out <- out[order(out$.id.,out$state,out$time),]
out$.id. <- NULL
if(!is.null(row.names)) rownames(out) <- row.names
out
}
