Nothing
########################################################################
## This program is Open Source Software: you can redistribute it ##
## and/or modify it under the terms of the GNU General Public License ##
## as published by the Free Software Foundation, either version 3 of ##
## the License, or (at your option) any later version. ##
## ##
## This program is distributed in the hope that it will be useful, ##
## but WITHOUT ANY WARRANTY; without even the implied warranty of ##
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ##
## General Public License for more details. ##
## ##
## You should have received a copy of the GNU General Public License ##
## along with this program. If not, see http://www.gnu.org/licenses/. ##
########################################################################
#require("methods")
# classes
#' Class \code{"ArtCohort"}
#'
#' TestIs a S4 class for the artificial cohort generated by \code{simulateCohort}.
#' @name ArtCohort
#' @aliases ArtCohort-class [,ArtCohort-method summary,ArtCohort-method
#' update,ArtCohort-method head,ArtCohort-method tail,ArtCohort-method
#' @docType class
#' @section Objects from the Class: Objects are created by calls to the
#' function \code{simulateCohort}.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{simulateCohort}},
#' \code{\link{transition.structure-class}},
#' \code{\link{transitionProbabilities}}, {\link{cumulativeIncidence}}
#' @keywords classes
#' @param
#' \code{x} object an ArtCohort
#' \code{i}, \code{j}, \code{drop} same as for data.frame
#' \code{...} passed on to data.frame method
#' \code{newsize} size of the updated cohort
#' \code{addbaseline} baseline for new part of cohort
#' \code{newInitialStates} initial states for new part of cohort
#' @slot
#' states.number Object of class "numeric": number of states
#' size Object of class "numeric": cohort size
#' baseline Object of class "matrix": baseline matrix
#' follow.up Object of class "numeric": maximum follow-up time
#' parameters Object of class "transition.structure": input parameters
#' parametersCovariances
#' Object of class "transition.structure": input covariance matrices
#' timeToTransition Object of class "matrix": input timeToTransition matrix. logical components
#' transitionFunctions
#' Object of class "transition.structure": input hazard functions
#' time.to.state Object of class "data.frame": entry times for each patient into each of the states
#' @examples
#'
#' showClass("ArtCohort")
#'
NULL
# #' GUIgems
# #' @name GUIgems
# #' package \code{GUIgems}
# #' Graphical User Interface for R package gems
# #' @author Zofia Baranczuk, Janne Estill, Olivia Keiser, Luisa Salazar Vizcaya, Nello Blaser
# #' @note
# #'
# NULL
#' Class \code{"PosteriorProbabilities"}
#'
#' This S4 class summarizes the posterior probabilities over time for objects
#' of class \code{"ArtCohort"}
#'
#' @name PosteriorProbabilities
#' @aliases PosteriorProbabilities [,PosteriorProbabilities-method head,PosteriorProbabilities-method
#' tail,PosteriorProbabilities-method
#' @docType class
#' @slot states Object of class "character": names of states
#' @slot times Object of class "numeric": time points at which probabilities are evaluated
#' @slot probabilities Object of class "matrix": posterior Probabilities to be in each state at each time
#' @slot lower Object of class "matrix": lower prediction bound to be in each state at each time
#' @slot upper Object of class "matrix": upper prediction bound to be in each state at each time
#' @slot type Object of class "character": describes type of probability
#' @section Objects from the Class: Objects are created by calls to the
#' function \code{simulateCohort}.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{transitionProbabilities}},
#' \code{\link{cumulativeIncidence}}, \code{\link{ArtCohort}}
#' @keywords classes
#' @examples
#'
#' showClass("PosteriorProbabilities")
#'
NULL
#' Class \code{"transition.structure"}
#'
#' This S4 class provides a structure to specify different characteristics of
#' transitions, such as transition functions functions, parameters or parameter
#' covariances.
#'
#' @import methods
#' @name transition.structure
#' @aliases transition.structure-class [[,transition.structure-method
#' [[<-,transition.structure-method print,transition.structure-method
#' @docType class
#' @section Objects from the Class: Objects are created by calls to the
#' functions \code{generateHazardMatrix}, \code{generateParameterMatrix},
#' \code{generateParameterCovarianceMatrix}.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{generateHazardMatrix}},
#' \code{\link{generateParameterMatrix}},
#' \code{\link{generateParameterCovarianceMatrix}}
#' @keywords classes
#' @examples
#'
#' showClass("transition.structure")
#'
NULL
setClass("transition.structure", representation(states.number = "numeric", list.matrix = "matrix") )
setClass("FuncInteractions")
setClass("FuncParametersDist", representation(func = "function", mu = "list"), contains = "FuncInteractions")
setClass("ArtCohort",
representation(
states.number = "numeric",
size = "numeric",
baseline = "matrix",
baselineFunction = "function",
follow.up = "numeric",
parameters= "transition.structure",
parametersCovariances = "transition.structure",
timeToTransition = "matrix",
transitionFunctions = "transition.structure",
time.to.state = "data.frame"
))
setClass("PosteriorProbabilities", representation(states = "character", times = "numeric", probabilities = "matrix",
lower= "matrix", upper = "matrix", type="character"))
setMethod("[[", "transition.structure", function(x, i, j){
x@list.matrix[[i,j]]
})
setMethod("[[<-", "transition.structure", function(x, i, j, value){
x@list.matrix[[i,j]] <- value
return(x)
})
setMethod("[", "PosteriorProbabilities", function(x, i, j){
x@probabilities[i,j]
})
setMethod("[", "ArtCohort", function(x, i, j){
x@time.to.state[i,j]
})
possibleTransitions = function(object) 0
implicitGeneric("possibleTransitions")
setGeneric("possibleTransitions")
setMethod( "possibleTransitions", "transition.structure", function(object){
aux1 <- auxcounter(object@states.number)
aux1[aux1 >0] <- 1
aux2 <- matrix(0,object@states.number,object@states.number)
aux2[object@list.matrix == "impossible"] = 1
aux2[is.null(object@list.matrix)] = 1
aux3 <- aux1-aux2
aux4 <- which(aux3 ==1, arr.ind = TRUE)
colnames(aux4) <- c("From", "To")
rownames(aux4) <- character()
aux4
})
#' plot
#' @description
#' Plots S4 objects
#' @param ci a binary variable, cumulative incidence
#' @param states names of states in the plot
#' @param lwd vector of lentht 2, with first component
#' for the point estimate and second component
#' for the confidence interval
#' @param col vector of lentht 2, with first component
#' for the point estimate and second component
#' for the confidence interval
#' @param lty vector of lentht 2, with first component
#' for the point estimate and second component
#' for the confidence interval
#' @param main same as in any plot
#' @param xlab same as in any plot
#' @param ylab same as in any plot
#' @param x a cohort simulated by \code{simulateCohort}.
#' @param ... arguments passed on to main method
# #' @usage plot(x, ci = FALSE, main = paste(x@type,
# #' "after starting in State", x@states[1], "at time 0"),
# #' states = 1:dim(x@probabilities)[2], lwd = c(2, 2),
# #' col = c("blue", "green3"), lty = c(1, 2), xlab = "Time",
# #' ylab = "Probability", ...)
#' @export
#' @import graphics
#' @import stats
#' @import grDevices
setMethod("plot", "PosteriorProbabilities", function(x, ci=FALSE, main = paste(x@type, "after starting in State", x@states[1], "at time 0"), states=1:dim(x@probabilities)[2],
lwd=c(2,2), col=c('blue','green3'), lty=c(1,2), xlab="Time", ylab="Probability", ...){
if (ci){
plotPrevalence(x@times, x@probabilities, x@states, lower=x@lower, upper=x@upper, main=main, states=states,
lwd=lwd, col=col, lty=lty, xlab=xlab, ylab=ylab, ...)
if (sum(complete.cases(x@lower))==0) {
warning("Too few simulations for prediction intervals")
}
par(mfrow=c(1,1))
}
else {
plotPrevalence(x@times, x@probabilities, x@states, lower=NA, upper=NA, main=main, states=states,
lwd=lwd, col=col, lty=lty, xlab=xlab, ylab=ylab, ...)
par(mfrow=c(1,1))
}
})
setMethod( "update", "ArtCohort", function(object, newsize, addbaseline=matrix(NA, nrow = newsize - object@size), newInitialStates=rep(1, newsize - object@size)){
# View((object@parametersCovariances)@list.matrix)
if (dim((object@parametersCovariances)@list.matrix) == c(0,0)) parameterCovariance0 = FALSE
else parameterCovariance0 = object@parametersCovariances
print("is.na.data.frame")
if( all(is.na.data.frame(addbaseline)) &
!(all(is.na.data.frame(object@baseline))))
print("fun(size)")
print(newsize - object@size)
print(object@baselineFunction(newsize - object@size))
addbaseline <- object@baselineFunction(newsize - object@size)
print("start simulating the cohort")
aux1 =
simulateCohort(
transitionFunctions = object@transitionFunctions,
parameters = object@parameters,
cohortSize = newsize - object@size ,
parameterCovariances = parameterCovariance0,
timeToTransition= object@timeToTransition,
baseline = addbaseline,
baselineFunction =object@baselineFunction,
initialState=newInitialStates,
to=object@follow.up#5
)
newcohort <- object
print("binding baselines")
View(as.data.frame(object@baseline))
View(as.matrix((addbaseline)))
newcohort@baseline <- as.matrix(rbind(as.data.frame(object@baseline), as.data.frame(as.matrix(addbaseline))))
newcohort@size <- newsize
newcohort@follow.up <-object@follow.up
print("binding time.to.state")
aux2 <- rbind(object@time.to.state, aux1@time.to.state)
print("aux2")
rownames(aux2) <- paste("Patient", 1:newsize)
newcohort@time.to.state <- aux2
newcohort
})
setMethod("head", "ArtCohort", function(x, ...){
head(x@time.to.state, ...)
})
setMethod("tail", "ArtCohort", function(x, ...){
tail(x@time.to.state, ...)
})
setMethod("head", "PosteriorProbabilities", function(x, ...){
head(x@probabilities, ...)
})
setMethod("tail", "PosteriorProbabilities", function(x, ...){
tail(x@probabilities, ...)
})
setMethod("print", "transition.structure", function(x){
mat <- x@list.matrix
mchar <- lapply(mat, function(x) {
y <- x
if (is.function(x)) {
y <- paste("function(", paste(names(formals(x)), collapse=", "), ")", sep="")
}
return(y)
})
dim(mchar) <- dim(mat)
dimnames(mchar) <- dimnames(mat)
print(mchar)
})
setMethod( "summary", "ArtCohort", function(object){
summ <- matrix(list(
(object@states.number),
(object@size),
(object@baseline),
# (object@baselineFunction),
(object@follow.up),
(object@parameters),
(object@parametersCovariances),
(object@timeToTransition),
(object@transitionFunctions),
(object@time.to.state)),9,2)
aux2 <- getSlots("ArtCohort")
summ[,2] <- noquote(aux2)
colnames(summ) <- c("Mode", "Class")
rownames(summ) <- names(aux2)
print(noquote(matrix( c(class(object), mode(object)),ncol = 2, nrow = 1,dimnames = list("",c("Class", "Mode")))))
summ.all <- list()
summ.all$general <- summ
summ.all$events <- noquote(summary(object@time.to.state)[,2:object@states.number])
aux3 <- unlist(apply(object@time.to.state, 2, function(x) length(which(is.na(x) == 0))))
message("
Number of patients entering each state: ")
print(aux3)
message("
object and events: ")
summ.all
})
###TODO: methods for ArtCohort
#
# setMethod("addState", "ArtCohort", function(object){
#
# })
#
# setMethod("mergeStates", "ArtCohort", function(object){
#
# })
#
#
# setMethod("chooseSubcohort", "ArtCohort", function(object){
#
# })
#
#
# setMethod("plotFun", "ArtCohort", function(object){
#
# })
# setMethod("plotPaths", "ArtCohort", function(object){
#
# })
#
# setMethod("nextPath", "ArtCohort", function(object){
#
# })
#
# setMethod("save", "ArtCohort", function(object){
# is it already there? is it loadable?
# })
#
# setMethod("addState", "ArtCohort", function(x, ...){
#
# })
#' Counts a cost of the interventions in the cohort
#' @name countCost
#' @description
#' Counts a cost of the intervention based on the cost data and the times spent in each state.
#' @param cohort a cohort simulated in simulateCohort()
#' @param cost a list of lists:stateCost, units and qalys defining costs of basic units of cost,
#' cost per state expressed in basic units and QALY for each state.
#' @return \code{c("qaly","cost")}
#' \code{qaly} an average QALY for an individual in the cohort
#' \code{cost} an average cost for an individual in the cohort
#' @export
#' @import utils
countCost <- function(cohort, cost){
stis <- sumTimeInStates(cohort@time.to.state, cohort@follow.up)
View(stis)
# edit(cost[["units"]])
# edit(cost[["stateCost"]])
eval_cost <- c(1:length(cost[["stateCost"]]))
for (f in seq(1:length(cost[["stateCost"]]))){
eval_cost[f] <- with(cost[["units"]], eval(parse(text=(cost[["stateCost"]][f]))))
}
View(eval_cost)
g_qaly <- unlist(stis)*unlist(cost[["qualy"]])
n_patients <- nrow(cohort@time.to.state)
# g_cost <- data.frame(mapply(`*`,df,v))
View(list(cohort@follow.up, sum(stis)))
g_cost <- unlist(stis)*unlist(eval_cost)
return(c("qaly"=sum(g_qaly)/(nrow(cohort@time.to.state)*cohort@follow.up),
"cost"=sum(g_cost)/nrow(cohort@time.to.state)) )
}
#' Builds hazard function matrix based on the hazard function names matrix.
#' @param mydata - cohort@@time.to.state created in simulateCohort()
#' @param end_time - time of the simulation
#' @return time spent in each of the states by each simulated unit
#' @export
timeInStates <- function(mydata, end_time=100){
#mydata - cohort@time.to.state
#res -times in each state for every patient
#TODO:
nr <- nrow(mydata)
nc <- ncol(mydata)
tmp <- mydata*0
tmp[is.na(tmp)] <- 0
tmp <- tmp+end_time+1
tmp1 <- tmp
mydata[,nc+1]<- end_time
#Checking the differences between subsequent stages, with distance 2, 3, etc.
for (f in c(1:(nc))){
tmp1[, 1:(nc+1-f)]<- (-1)*(mydata[,1:(nc+1-f)]-mydata[,(f+1):(nc+1)])
tmp1[is.na(tmp1)]<-end_time+1
tmp1[tmp1<0]<-end_time +1
tmp <- pmin((tmp),(tmp1))
}
res <- tmp
res[res==(end_time+1)]<-0
print(dim(res))
View(res)
return(res)
}
## #' Builds hazard function matrix based on the hazard function names matrix.
## #' @param mypanel, hfNames
## #' @return mypanel
## #' @export
sumTimeInStates <- function(mydata, end_time = 100){
colSums(timeInStates(mydata, end_time))
}
# functions
auxcounter <-
function (statesNumber, statesNamesFrom =paste("State", 1:statesNumber) ,statesNamesTo = statesNamesFrom)
{
### "auxcounter" associates an index to each position in the transition matrix:
# Input: number of states
# Output: matrix M with transitions numbers, where element M_ij is the
# transition number of the transition from State i to state j.
m = matrix(0, ncol = statesNumber, nrow = statesNumber)
for (i in 1:(statesNumber - 1)) {
m[i, (i + 1):statesNumber] <- 1:(statesNumber - i) +
(i - 1) * (statesNumber - i/2)
}
dimnames(m) <- list(from=statesNamesFrom, to=statesNamesTo)#(from = paste("State", 1:statesNumber),
# to = paste("State", 1:statesNumber))
return(m)
}
auxposition <-
function (matrix, value)
{
### "auxposition" finds the position of a value in a matrix:
# Input: a matrix and a value
# Output: position of (1st occurence of) value in matrix
which(matrix == value, arr.ind = TRUE)[1, ]
}
### The function "createCohorts" is the main internal function.
createCohorts <-
function (hazardf, statesNumber,
statesNamesFrom = strsplit((dimnames(hazardf)["from"])[[1]]," "),statesNamesTo = statesNamesFrom,
cohortSize, mu, sigma = matrix(0, nrow = length(unlist(mu)), ncol = length(unlist(mu))), historyl = FALSE,
startingStates = rep(1, cohortSize), absorbing=cohortSize, impossible = NULL, fixpar = NULL,
direct = NULL, bl0 = matrix(0, nrow = cohortSize), to = 100, report.every, sampler.steps)
{
statesNames <- list(from=statesNamesFrom, to=statesNamesTo)
hazardf <- t(hazardf)[t(auxcounter(statesNumber)) > 0]
if (length(hazardf) < (statesNumber * (statesNumber - 1)/2)) {
length(hazardf) <- statesNumber * (statesNumber - 1)/2
}
for (tr in 1:length(hazardf)) {
if (length(hazardf[[tr]]) > 0) {
if (!is.function(hazardf[[tr]])) {
if (!hazardf[[tr]] %in% c("Weibull", "multWeibull", "Exponential", "impossible") &&
!is.na(hazardf[[tr]]))
try(hazardf[[tr]] <- as.function(hazardf[[tr]]))
}
}
}
mu <- as.list(unlist(t(mu), recursive = FALSE))
if (length(mu) > 0) {
mu <- lapply(mu, as.numeric)
}
if(class(sigma) == "transition.structure"){
if (dim(sigma)[1] != length(unlist(mu)))
stop("size of parameters and parameterCovariances are inconsistent")
}
if (length(startingStates) == 1) rep(startingStates, cohortSize)
if (cohortSize>1) {try(bl0 <- as.matrix(bl0))}
else if (is.null(dim(bl0))) {try(bl0 <- t(as.matrix(bl0)))}
if (nrow(bl0) != cohortSize)
warning("baseline is not of the right dimension")
for (i in 1:length(hazardf)) {
if (is.element(i, impossible)) {
if (historyl == FALSE)
hazardf[[i]] = function(t, bl) 99 * to
else hazardf[[i]] = function(t, history, bl) 99 *
to
}
}
parametric = numeric()
k = 0
for (i in 1:max(auxcounter(statesNumber))) {
if (is.character(hazardf[[i]])) {
k = k + 1
parametric[k] = i
if (hazardf[[i]] %in% c("Weibull", "multWeibull", "Exponential")) {
hazardf[[i]] <- switch(hazardf[[i]],
Weibull=function(t, shape, scale, history, bl) rweibull(t, shape, scale),
multWeibull=function(t, w, shapes, scales, history, bl) multPar(t, w, shapes, scales),
Exponential=function(t, rate, history, bl) rexp(t, rate))
if (!historyl) formals(hazardf[[i]])$history <- NULL
}
else {
stop(paste("Transition function for transition", i, "not recognized", sep=" "))
}
}
}
allFunctions = mainFunctions(statesNumber = statesNumber,
Mu = mu, sigma = sigma, cohortSize = cohortSize, history = historyl,
hazardf, impossible = c(impossible, fixpar, direct))
parametric = sort(c(parametric, impossible, direct))
cohorts <- sapply(1:cohortSize, function(i){
if ((i%%report.every ==0)| (i ==1)) message("Simulating patient:" ,i) ;
print(bl0[i,])
historical(gf = allFunctions[[i]],
statesNumber = statesNumber, parametric = parametric,
historyl = historyl, startingState = startingStates[i],
absorbing=absorbing, bl = bl0[i, ], to = to, sampler.steps = sampler.steps)[[1]]})
print("dimnames")
dimnames(cohorts) <- list(statesNamesFrom, paste("Patient", 1:cohortSize))#(paste("State", 1:statesNumber),
#paste("Patient", 1:cohortSize))
print("return(cohorts)")
return(cohorts)
}
data_prep <-
function (progressionObject, maxtime)
{
states_number <- dim(progressionObject)[1]
cohortSize <- dim(progressionObject)[2]
times <- t(progressionObject)
status <- array(TRUE, dim = dim(times))
status[is.na(times)] <- FALSE
times[times>=maxtime]<-maxtime
status[times >= maxtime] <- FALSE
times[status == FALSE] <- maxtime + 1
tmat <- auxcounter(states_number)
tmat[tmat==0] <- NA
prepData <- msprep(times, status, trans = tmat,
start = list(state = 1, time = 0))
return(prepData)
}
fold <-
function (vector, object)
{
kk <- 0
lc <- 0
newobject <- list()
for (i in 1:length(object)) {
k <- 0
newobject[[i]] <- list()
for (j in 1:(length(object[[i]]@mu))) {
l <- length((object[[i]]@mu)[[j]])
newobject[[i]][[j]] = c(vector[(1 + k + kk):(k +
kk + l)])
k <- k + l
lc <- lc + l
}
kk <- lc
}
return(newobject)
}
# Builds hazard function matrix based on the hazard function names matrix.
# @param mypanel, hfNames
# @return mypanel
# @export
addStatesNames<- function(list.matrix, statesNamesFrom, statesNamesTo = statesNamesFrom){
statesNames <- list(from=statesNamesFrom,to=statesNamesTo)
#if (exists(hf@list.matrix)){
dimnames(list.matrix) <- statesNames
#}
return(list.matrix)
}
#' generate template for transition functions
#'
#' This function simplifies generating the matrix of transition functions.
#'
#'
#' @param statesNumber the number of states to be considered.
#' @param statesNamesFrom names of the states in the model.
#' @param statesNamesTo names of the states in the model.
#' @return a \code{transition.structure} of dimension \eqn{N \times N}{N x N},
#' where \eqn{N} is the number of states and with value "impossible" for all
#' potential transitions.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer, Zofia Baranczuk
#' @seealso \code{\link{transition.structure-class}},
#' \code{\link{simulateCohort}}
#' @keywords utilities
#' @export generateHazardMatrix
generateHazardMatrix <- function (statesNumber, statesNamesFrom = paste("State", 1:statesNumber),
statesNamesTo = statesNamesFrom){
# to = paste("State", 1:statesNumber))
statesNames <- list(from=statesNamesFrom,to=statesNamesTo)
hf <- list()
length(hf) <- statesNumber * (statesNumber - 1)/2
for (i in 1:length(hf)) {
hf[[i]] <- "impossible"
}
tmat <- auxcounter(statesNumber, statesNamesFrom, statesNamesTo)
for (trans in 1:(statesNumber * (statesNumber - 1)/2)) {
attributes(hf)$names[trans] <- paste("Transition", trans,
"- from state", auxposition(tmat, trans)[1], "to state",
auxposition(tmat, trans)[2])
}
hfmat <- matrix(list(), nrow = statesNumber, ncol = statesNumber)
for (i in 1:(statesNumber * (statesNumber - 1)/2)) {
hfmat[tmat == i][[1]] <- hf[[i]]
}
dimnames(hfmat) <- statesNames#list(from = paste("State", 1:statesNumber),
#to = paste("State", 1:statesNumber))
hfclass <- new("transition.structure")
hfclass@states.number <- dim(hfmat)[1]
hfclass@list.matrix <- hfmat
return(hfclass)
}
#' generate a template for parameter covariances
#'
#' This function simplifies generating the matrix of parameter covariances from
#' a matrix of mean parameters.
#'
#'
#' @param muM a \code{transition.structure} of dimension \eqn{N \times N}{N x
#' N}, whose components \code{list} the mean values for the parameters in the
#' \code{transitionFunctions}.
#' @param statesNamesFrom a list of names of states in the model. By default list(State 1, State 2, ..., State N)
#' @param statesNamesTo a list of names of states in the model. By default list(State 1, State 2, ..., State N)
#' @return a \code{transition.structure} of dimension \eqn{N \times N}{N x N}
#' of covariance matrices for the \code{parameters}.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso
#' \code{\link{transition.structure}},
#' \code{\link{generateParameterMatrix}}, \code{\link{simulateCohort}}
#' @keywords utilities
#' @export generateParameterCovarianceMatrix
generateParameterCovarianceMatrix <- function (muM, statesNamesFrom = (dimnames(muM@list.matrix)["from"])[[1]],
statesNamesTo = statesNamesFrom )
{
print(muM)
mu = muM@list.matrix
SigmaMat <- matrix(list(), nrow = dim(mu)[1], ncol = dim(mu)[2])
dimnames(SigmaMat) <- list(from=statesNamesFrom, to=statesNamesTo)# dimnames(mu)#list(from = paste("State", 1:dim(mu)[1]),
# to = paste("State", 1:dim(mu)[1]))
lengthMu <- numeric(length(mu))
for (trans in 1:length(mu)) {
lengthMu[trans] <- length(unlist(mu[[trans]]))
if (lengthMu[trans] > 0) {
SigmaMat[[trans]] <- matrix(0, nrow = lengthMu[trans],
ncol = lengthMu[trans])
}
}
Sigmaclass <- new("transition.structure")
Sigmaclass@states.number <- dim(SigmaMat)[1]
Sigmaclass@list.matrix <- SigmaMat
return(Sigmaclass)
}
#' generate a template for mean parameters
#'
#' This function simplifies generating the matrix of mean parameters from a
#' matrix of transition functions.
#'
#'
#' @param hf a \code{transition.structure} of dimension \eqn{N \times N}{N x
#' N}, where \eqn{N} is the number of states.
#' @param statesNamesFrom - a list of names of states in the model. By default list(State 1, State 2, ..., State N)
#' @param statesNamesTo - a list of names of states in the model. By default the same as statesNamesFrom
#' @return a \code{transition.structure} of dimension \eqn{N \times N}{N x N},
#' whose components are \code{lists} of the right length for the parameters in
#' the corresponding hazard function \code{hf}.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{transition.structure-class}},
#' \code{\link{simulateCohort}}
#' @keywords utilities
#' @export generateParameterMatrix
generateParameterMatrix <- function (hf,statesNamesFrom = (dimnames(hf@list.matrix)["from"])[[1]], # strsplit((dimnames(hf@list.matrix)["from"])[[1]]," "),
statesNamesTo = statesNamesFrom)
{
hfMat <- hf@list.matrix
tmat <- auxcounter(dim(hfMat)[1],statesNamesFrom, statesNamesTo)
hf <- t(hfMat)[t(tmat > 0)]
MuMat <- matrix(list(), nrow = dim(hfMat)[1], ncol = dim(hfMat)[2])
dimnames(MuMat) <- list(from=statesNamesFrom, to=statesNamesTo)
lengthMu <- numeric(length(hf))
for (trans in 1:length(hf)) {
if (class(hf[[trans]]) == "character") {
# print("hf[trans]")
# print(hf[trans])
lengthMu[trans] <- length(my_formals(hf[[trans]]))#[[1]])#switch(hf[[trans]], Weibull = 2,
# multWeibull = 3, Exponential = 1, impossible = 0)
}
else if (class(hf[[trans]]) == "function") {
if (length(formals(hf[[trans]])) > 0) {
n1 <- attributes(formals(hf[[trans]]))$names
n2 <- n1[!(n1 %in% c("t", "bl", "history", "process"))]
lengthMu[trans] <- length(n2)
}
}
}
# for (trans in 1:length(hf)) {
# if (lengthMu[trans] > 0) {
# MuMat[tmat == trans][[1]] <- as.list(numeric(lengthMu[trans]))
# }
# }
for (f in 2:nrow(hfMat)){
for (g in 1:(f-1)){
MuMat[[g,f]] <- as.list(rep(0, length(my_formals(hfMat[[g,f]]))))#[[1]])))
names(MuMat[[g,f]])<-names(my_formals(hfMat[[g,f]])) #[[1]])
# }
# else{
# MuMat[f,g] <- list()
}
}
#View(MuMat)
Muclass <- new("transition.structure")
Muclass@states.number <- dim(MuMat)[1]
Muclass@list.matrix <- MuMat
return(Muclass)
}
historical <- function(gf, statesNumber, parametric, historyl, startingState, absorbing, bl, to, sampler.steps){
if (historyl) histH(gf, statesNumber, parametric, startingState, absorbing, bl, to, sampler.steps = sampler.steps)
else histNoH(gf, statesNumber, parametric, startingState, absorbing, bl, to, sampler.steps = sampler.steps)
}
histNoH <- function(gf, statesNumber, parametric, startingState, absorbing, bl, to, sampler.steps){
possible <- auxcounter(statesNumber)
# all possible transition times
npar <- (1:length(gf))[!1:length(gf)%in%parametric]
time0 <- rep(NA,length(gf))
time0[parametric] <- unlist(lapply(gf[parametric], function(ff) samplerP(function(t) ff(t, bl=bl), n=1)))
time0[npar] <- unlist(lapply(gf[npar], function(ff) sampler(n=1, function(t) ff(t, bl=bl), to=to, length = sampler.steps)))
#times in the transition matrix
time = matrix(ncol = statesNumber, nrow = statesNumber)
time[unlist(lapply(1:max(possible), function(i) which(possible == i)))] <- time0
time[possible==0] <- 99*to
minTimeState <- apply(time, 1, function(x) c(min(x), which.min(x)))
kk = startingState
path = rep(NA,statesNumber)
path[startingState]<-0
aux = startingState
while(aux %in% setdiff(startingState:statesNumber, absorbing) && sum(path, na.rm=TRUE)<=to){
aux = minTimeState[2,aux]
path[aux] = minTimeState[1,kk] + sum(path, na.rm=TRUE)
kk = aux
}
return(list(path,NULL))
}
histH <- function(gf, statesNumber, parametric, startingState, absorbing, bl, to, sampler.steps){
possible <- auxcounter(statesNumber)
formerhistory <- array(0,dim = c(max(possible), 4))
eventTimes <- matrix(NA, ncol = statesNumber ,nrow= 1)
aux = possible[startingState,]; aux = aux[aux !=0]
par = match(intersect(aux,parametric),aux)
f = gf[aux]; rm(aux)
npar <- (1:length(f))[!1:length(f)%in%par]
time = numeric(length(f))
time[par] <- unlist(lapply(f[par], function(ff) samplerP(function(t) ff(t, history=0, bl=bl), n=1)))
time[npar] <- unlist(lapply(f[npar], function(ff) sampler(n=1, function(t) ff(t, history=0, bl=bl),to=to, length = sampler.steps)))
#########################################
mt = min(time)
wmin = which.min(time)
first <- possible[startingState, startingState+wmin]
k = auxposition(possible, max(first))[2] #current state of the patient
formerhistory[possible[startingState,startingState+1]+wmin-1,1:4] = c(first, mt, mt, k)
lim = mt
######## loop over states...
while(k %in% setdiff(startingState:statesNumber, absorbing) && lim <= to){
aux = possible[auxposition(possible, max(na.rm = TRUE,formerhistory[,1]))[2],]; aux = aux[aux !=0]
f = gf[aux]
par = match(intersect(aux,parametric),aux) #those index of f (possible transition hazards corespond to parametric)
npar <- (1:length(f))[!1:length(f)%in%par]
timesp = numeric(length(f))
timesp[par] <- unlist(lapply(f[par], function(ff) samplerP(function(t) ff(t, history=formerhistory[,2], bl=bl), n=1)))
timesp[npar] <- unlist(lapply(f[npar], function(ff) sampler(n=1, function(t) ff(t, history=formerhistory[,2], bl=bl),to=to, length = sampler.steps)))
mt = min(timesp)
wmin = which.min(timesp)
auxmin <- aux[wmin]
lim = mt + formerhistory[max(formerhistory[,1]) ,3]
k = auxposition(possible, max(auxmin))[2]
formerhistory[auxmin,1:4] = c(auxmin, mt, lim, k)
}
for(st in startingState:statesNumber){
if(length(which(formerhistory[,4]==st)) == 0)
eventTimes[1,st] = ifelse(st==startingState, 0, NA)
else eventTimes[1,st] = formerhistory[which(formerhistory[,4]==st)[1],3]
}
return(list(eventTimes, NULL))
}
mainFunctions <-
function (statesNumber, Mu, sigma, cohortSize, history, functions,
impossible)
{
objects = setFunctions(statesNumber = statesNumber, Mu = Mu,
history = history, functions, impossible = impossible)
if (class(sigma) == "logical") return(simFunctions.NoUnc(objects, covariances = sigma, history = history,
statesNumber = statesNumber, impossible = impossible,
cohortSize = cohortSize))
else return(simFunctions(objects, covariances = sigma, history = history,
statesNumber = statesNumber, impossible = impossible,
cohortSize = cohortSize))
}
msmDataPrep <-
function (mstateDataPrep)
{
a <- mstateDataPrep
cohortSize <- max(a$id)
initState <- sapply(1:cohortSize, function(ID) min(a$from[a$id==ID]))
unsortedId <- c(a$id[a$status == 1], seq(1, cohortSize))
unsortedTime <- c(a$Tstop[a$status == 1], rep(0, cohortSize))
unsortedState <- c(a$to[a$status == 1], initState)
unsortedData <- data.frame(unsortedId, unsortedTime, unsortedState)
sortedData <- unsortedData[order(unsortedId, unsortedTime),
]
msmData <- data.frame(id = sortedData$unsortedId, time = sortedData$unsortedTime,
state = sortedData$unsortedState)
}
multPar <-
function (t, w, shape, scale)
{
ii <- apply(rmultinom(t, 1, c(w, 1 - sum(w))), 2, which.max)
mW <- rweibull(t, shape[ii], scale[ii])
attr(mW, "weights") <- w
attr(mW, "index") <- ii
return(mW)
}
plotPrevalence <-
function (times, prevalence, stateNames, lower, upper, typ = "separate", main = "State-wise probability over time", states,
lwd=c(2,2), col=c('blue','green3'), lty=c(1,2), xlab="Time", ylab="Probability" , ...)
{
if (typ == "separate") {
if (length(col)==1) col <- rep(col, 2)
if (length(lty)==1) lty <- rep(col, 2)
if (length(lwd)==1) lwd <- rep(lwd, 2)
par(mfrow = c(2, ceiling(length(states)/2)))
if (length(states)==1) par(mfrow = c(1,1))
par(oma = c(0, 0, 2, 0))
for (i in states) {
plot(times, prevalence[, i], type = "l", col = col[1],
lwd = lwd, ylim = c(0, 1), ylab = ylab,
xlab = xlab, main = stateNames[[i]], ...)
try(lines(times, lower[,i], col=col[2], lty=lty[2]), silent=TRUE)
try(lines(times, upper[,i], col=col[2], lty=lty[2]), silent=TRUE)
}
mtext(main, outer = TRUE, cex = 1.5)
}
else if (typ == "stacked") {
par(mfrow = c(1, 1))
totalPrev <- prevalence[, 1]
plot(times, totalPrev, col = 1, type = "l", ylim = c(0,
1), ylab = "Probability", xlab = "Time")
for (i in 2:dim(prevalence)[2]) {
totalPrev <- totalPrev + prevalence[, i]
lines(times, totalPrev, col = i)
}
}
}
posteriorProbabilities = function(object, times, M=100) {
if (class(object)=="ArtCohort") cohorts <- t(object@time.to.state)
else cohorts <- t(object)
statesNumber <- dim(cohorts)[1]
stateNames = rownames(cohorts)
print(toString(statesNumber))
print(rownames(cohorts))
# Prediction interval:
if (dim(cohorts)[2]>=1000) {
dd <- data.frame(t(cohorts))
dd$cc <- c(rep(1:M, times=floor(dim(cohorts)[2])/M) , 1:(dim(cohorts)[2]%%M+1))[1:dim(cohorts)[2]]
dd <- split(dd,dd$cc)
prev <- list()
ppp <- NULL
for (ciind in 1:M){
cohorts <- t(dd[[ciind]][,1:statesNumber])
if (ncol(cohorts)==0) {
prev[[ciind]] <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
}
if (ncol(cohorts)==1) {
prev[[ciind]] <- matrix(0, ncol = nrow(cohorts), nrow = length(times))
for (tInd in 1:(nrow(cohorts)-1)) {
prev[[ciind]][times < min(c(max(times)+1,cohorts[(tInd+1):nrow(cohorts)]), `na.rm`=TRUE) & times >= cohorts[tInd] , tInd] <- 1
}
prev[[ciind]][times >= cohorts[nrow(cohorts)] , nrow(cohorts)] <- 1
}
if (ncol(cohorts)>=2) {
prep1 <- data_prep(cohorts, max(times))
prep1 <- prep1[!(prep1$Tstart == prep1$Tstop), ]
prepData <- msmDataPrep(prep1)
prev[[ciind]] <- prevalence(prepData, times, dim(cohorts)[1])
}
ppp <- rbind(ppp,prev[[ciind]])
}
ddd <- `data.frame`(ppp)
ddd$times <- times
prev <- `as.matrix`(plyr::ddply(ddd, plyr::.(times), function(x) colMeans(x)))[, 1:statesNumber]
#prev <- as.matrix(plyr::ddply(ddd, .(times), function(x) apply(x, 2, function(y){quantile(y, .5)})))[, 1:statesNumber]
lower <- `as.matrix`(plyr::ddply(ddd, plyr::.(times), function(x) apply(x, 2, function(y){quantile(y, .025)})))[, 1:statesNumber]
upper <- as.matrix(plyr::ddply(ddd, plyr::.(times), function(x) apply(x, 2, function(y){quantile(y, .975)})))[, 1:statesNumber]
}
else {
if (ncol(cohorts)==0) {
prev <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
}
if (ncol(cohorts)==1) {
prev <- matrix(0, ncol = nrow(cohorts), nrow = length(times))
for (tInd in 1:(nrow(cohorts)-1)) {
prev[times < min(c(max(times)+1,cohorts[(tInd+1):nrow(cohorts)]), `na.rm`=TRUE) & times >= cohorts[tInd] , tInd] <- 1
}
prev[times >= cohorts[nrow(cohorts)] , nrow(cohorts)] <- 1
}
if (ncol(cohorts)>=2) {
prep1 <- data_prep(cohorts, max(times))
prep1 <- prep1[!(prep1$Tstart == prep1$Tstop), ]
prepData <- msmDataPrep(prep1)
prev<- prevalence(prepData, times, dim(cohorts)[1])
}
lower <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
upper <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
}
print(dim(prev))
length(stateNames)
print(stateNames)
print(length(list(paste("Time", times), stateNames)))
dimnames(prev) <- list(paste("Time", times), stateNames)
dimnames(lower) <- list(paste("Time", times), stateNames)
dimnames(upper) <- list(paste("Time", times), stateNames)
pp = new("PosteriorProbabilities")
pp@states = stateNames
pp@times = times
pp@probabilities = prev
pp@lower = lower
pp@upper = upper
pp@type = "Transition probabilities"
return(pp)
}
#' Cumulative Incidence
#'
#' Calculates the cumulative incidence and prediction intervals after first
#' state
#'
#'
#' @param object either the output of \code{\link{simulateCohort}} or the
#' \code{matrix} with the \code{probabilities} slot of that output.
#' @param times a time vector.
#' @param M number of groups for calculating confidence intervals.
#' @return an object of class \code{"PosteriorProbabilities"}, containing the
#' statenames, timepoints and the cumulative incidence with pointwise
#' prediction intervals over time.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{PosteriorProbabilities}},
#' \code{\link{ArtCohort}}, \code{\link{simulateCohort}}
#' @keywords utilities
#' @export cumulativeIncidence
cumulativeIncidence <- function (object, times, M = 100) {
if (class(object) == "ArtCohort")
cohorts <- t(object@time.to.state)
else cohorts <- t(object)
statesNumber <- dim(cohorts)[1]
stateNames <- rownames(cohorts)
if (dim(cohorts)[2] >= 1000) {
dd <- data.frame(t(cohorts))
dd$cc <- c(rep(1:M, times = floor(dim(cohorts)[2])/M),
1:(dim(cohorts)[2]%%M + 1))[1:dim(cohorts)[2]]
dd <- split(dd, dd$cc)
inc <- list()
ppp <- NULL
for (ciind in 1:M) {
cohorts <- t(dd[[ciind]][, 1:statesNumber])
if (ncol(cohorts) == 0) {
inc[[ciind]] <- matrix(NA, ncol = nrow(cohorts),
nrow = length(times))
}
if (ncol(cohorts) == 1) {
inc[[ciind]] <- matrix(0, ncol = nrow(cohorts),
nrow = length(times))
for (tInd in 1:(nrow(cohorts))) {
inc[[ciind]][times >= cohorts[tInd], tInd] <- 1
}
}
if (ncol(cohorts) >= 2) {
prep1 <- data_prep(cohorts, max(times))
prep1 <- prep1[!(prep1$Tstart == prep1$Tstop),
]
prepData <- msmDataPrep(prep1)
inc[[ciind]] <- incidence(prepData, times, dim(cohorts)[1])
}
ppp <- rbind(ppp, inc[[ciind]])
}
ddd <- `data.frame`(ppp)
ddd$times <- times
inc <- `as.matrix`(plyr::ddply(ddd, plyr::.(times), function(x) colMeans(x)))[,
1:statesNumber]
lower <- `as.matrix`(plyr::ddply(ddd, plyr::.(times), function(x) apply(x,
2, function(y) {
quantile(y, 0.025)
})))[, 1:statesNumber]
upper <- `as.matrix`(plyr::ddply(ddd, plyr::.(times), function(x) apply(x,
2, function(y) {
quantile(y, 0.975)
})))[, 1:statesNumber]
}
else {
if (ncol(cohorts) == 0) {
inc <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
}
if (ncol(cohorts) == 1) {
inc <- matrix(0, ncol = nrow(cohorts), nrow = length(times))
for (tInd in 1:(nrow(cohorts))) {
inc[times >= cohorts[tInd], tInd] <- 1
}
}
if (ncol(cohorts) >= 2) {
prep1 <- data_prep(cohorts, max(times))
prep1 <- prep1[!(prep1$Tstart == prep1$Tstop), ]
prepData <- msmDataPrep(prep1)
inc <- incidence(prepData, times, dim(cohorts)[1])
}
lower <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
upper <- matrix(NA, ncol = nrow(cohorts), nrow = length(times))
}
dimnames(inc) <- list(paste("Time", times), stateNames)
dimnames(lower) <- list(paste("Time", times), stateNames)
dimnames(upper) <- list(paste("Time", times), stateNames)
pp = new("PosteriorProbabilities")
pp@states = stateNames
pp@times = times
pp@probabilities = inc
pp@lower = lower
pp@upper = upper
pp@type = "Cumulative incidence"
return(pp)
}
#' transition probabilities
#'
#' Calculates the probabilities and prediction intervals after first state
#'
#'
#' @param object either the output of \code{\link{simulateCohort}} or the
#' \code{matrix} with the \code{probabilities} slot of that output.
#' @param times a time vector.
#' @param M number of groups for calculating confidence intervals.
#' @return an object of class \code{"PosteriorProbabilities"}, containing the
#' statenames, timepoints and the transition probabilities with pointwise
#' prediction intervals over time.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{PosteriorProbabilities}},
#' \code{\link{ArtCohort-class}}, \code{\link{simulateCohort}}
#' @keywords utilities
#' @export transitionProbabilities
transitionProbabilities <- posteriorProbabilities
prepareF <-
function (func, known, historyl, historyp = numeric(2), blp = numeric(2))
{
lsim = length(known)
invisible(func)
if (historyl == TRUE) {
ofnh = function(t, sim = known, history = historyp, bl = blp) {
hisaux = matrix(0, ncol = 1, nrow = length(history))
hisaux[1:length(history), 1] = history[1:length(history)]
l = list()
l[[1]] = t
l[2:(length(sim) + 1)] = sim[1:length(sim)]
l[[length(sim) + 2]] = hisaux
l[[length(sim) + 3]] = bl
return(`do.call`(func, l))
}
return(ofnh)
}
else {
ofnl = function(t, sim = known, bl = blp) {
l = list(lsim + 2)
l[[1]] = t
l[2:(length(sim) + 1)] = sim[1:length(sim)]
l[[length(sim) + 2]] = bl
return(`do.call`(func, l))
}
return(ofnl)
}
}
## #' @export
prevalence <-
function (data, times, states_number)
{
states <- array(0, dim = c(length(times), max(data$id)))
for (i in 1:max(data$id)) {
ti <- data$time[data$id == i]
st <- data$state[data$id == i]
for (j in 1:length(ti)) {
states[times >= ti[j], i] <- st[j]
}
}
prev <- matrix(0, ncol = states_number, nrow = length(times))
for (i in 1:states_number) {
prev[, i] <- apply((states == i), 1, sum)
}
prev <- prev/length(unique(data$id))
return(prev)
}
incidence <-
function (data, times, states_number)
{
spl <- split(data, factor(data$state, levels=1:states_number))
inc <- matrix(0, ncol = states_number, nrow = length(times))
for (i in 1:states_number){
if (dim(spl[[i]])[1]>0){
sti <- sapply(times, function(x) x>=spl[[i]]$time)
if (dim(spl[[i]])[1]>1) inc[,i] <- colSums(sti)
else inc[,i] <- sti
}
else {
inc[,i] <- 0
}
}
inc <- inc/length(unique(data$id))
return(inc)
}
## #' @export
sampler <-
function (n, f, to = 100, length = 1000)
{
time = seq(0, to, length = length)
rate = f(time)
rate[rate == Inf] <- 0
rate[is.nan(rate)] <- 0
rate[is.na(rate)] <- 0
rate[rate < .Machine$double.eps] <- .Machine$double.eps
samples = msm::rpexp(n, rate, time)
if (is.nan(samples)) {
print(rate)
}
samples[is.nan(samples)] = 99*to
return(samples)
}
samplerP <-
function (f, n)
{
f = match.fun(f)
return(f(t = n))
}
setFunctions <-
function (statesNumber, Mu, history, functions, impossible)
{
trans = max(auxcounter(statesNumber))
sample = list(trans)
narg = numeric()
k = 0
for (i in 1:trans) {
sample[[i]] = new("FuncParametersDist")
sample[[i]]@func = functions[[i]]
if (!is.element(i, impossible)) {
if (history == TRUE)
narg[i] = length(formals(functions[[i]])) - 2 -
1
else narg[i] = length(formals(functions[[i]])) -
1 - 1
if (i == 1)
sample[[i]]@mu = Mu[1:narg[i]]
else sample[[i]]@mu = Mu[(k + 1):(k + narg[i])]
k = k + narg[i]
}
}
return(sample)
}
simFunctions <-
function (so, covariances, history, statesNumber, impossible,
cohortSize)
{
possible = max(auxcounter(statesNumber))
lso = length(so)
prov = setdiff(1:possible, impossible)
for (i in prov) {
if (i == prov[1]) {
aux1 = unfold(so[[i]]@mu)
aux2 = rep(1, length(aux1))
}
else {
uf = unfold(so[[i]]@mu)
aux1 = c(aux1, uf)
aux2 = c(aux2, rep(i, length(uf)))
}
}
if (sum(covariances !=0)){
simpar = MASS::mvrnorm(cohortSize, aux1, covariances)
}
else{
simpar <- matrix(aux1, nrow=cohortSize, ncol=length(aux1), byrow=TRUE)
}
simpar = rbind(aux2, simpar)
auxso = so[prov]
simpar1 <- lapply(1:cohortSize, function(k) fold(simpar[k+1,], auxso))
lp = list(lso)
ltp = list(cohortSize)
for (j in 1:cohortSize) {
k = 0
for (i in 1:lso) {
if (!is.element(i, impossible)) {
k = k + 1
lp[[i]] = prepareF(func = so[[i]]@func, known = simpar1[[j]][[k]],
historyl = history)
}
else {
lp[[i]] = so[[i]]@func
}
}
ltp[[j]] = lp
}
return(ltp)
}
simFunctions.NoUnc <- function (so, covariances, history, statesNumber, impossible, cohortSize)
{
possible = max(auxcounter(statesNumber))
lso = length(so)
prov = setdiff(1:possible, impossible)
for (i in prov) {
if (i == prov[1]) {
aux1 = unfold(so[[i]]@mu)
aux2 = rep(1, length(aux1))
}
else {
uf = unfold(so[[i]]@mu)
aux1 = c(aux1, uf)
aux2 = c(aux2, rep(i, length(uf)))
}
}
simpar = rbind(aux2, aux1)
auxso = so[prov]
simpar1 = list()
simpar1[[1]] = fold(simpar[2, ], auxso)
lp = list()
ltp = list()
k = 0
for (i in 1:lso) {
if (!is.element(i, impossible)) {
k = k + 1
lp[[i]] = prepareF(func = so[[i]]@func, known = simpar1[[1]][[k]],
historyl = history)
}
else {
lp[[i]] = so[[i]]@func
}
}
ltp[1:cohortSize] = list(lp)
return(ltp)
}
#' Simulate cohort
#'
#' Simulates a cohort of patients from a set of functions associated to each
#' possible transition in a multistate model. The multistate model is not
#' required to be a Markov model and may take the history of previous events
#' into account. In the basic version, it allows to simulate from
#' transition-specific hazard function, whose parameters are multivariable
#' normally distributed. For each state, all transition-specific hazard
#' functions and their parameters need to be specified. For simulating one
#' transition, all possible event times are simulated and the minimum is
#' chosen. Then simulation continues from the corresponding state until an
#' absorbing state of time \code{to} is reached.
#'
#' The \code{transitionFunctions} contains hazard functions or time to event
#' function associated to each possible transition. The elements of this
#' \code{list} can be either expressed as an explicit R \code{function} or as a
#' \code{character} ("impossible", "Weibull", "multWeibull", "exponential") in
#' order to express impossible transitions or parametric forms for the
#' distributions of time to event. If the functions should depend on time,
#' baseline characteristics or be \emph{history-dependent}, the function
#' arguments \emph{t}, \emph{bl} or \emph{history} can be used. Time \emph{t}
#' refers to the time since entry into the current state. For the time since
#' the initial state, use \code{t+sum(history)}.
#'
#' The components of the \code{parameters} argument \code{list} the mean values
#' for the parameters in the \code{transitionFunction}. If the corresponding
#' \code{transitionFunction} is a \code{function}, the parameters should appear
#' in the same order as in the \code{function}, leaving out \emph{t}, \emph{bl}
#' and \emph{history}. If the corresponding \code{transitionFunction} is the
#' \code{character} "Weibull", the first argument is the shape and the second
#' one the scale. If the corresponding \code{transitionFunction} is the
#' \code{character} "multWeibull", specify weights, shapes, scales in this
#' order.
#'
#' Note that when using the \code{parameterCovariances} argument it is the
#' users responsibility to ensure that the functions are parametrized such that
#' \code{parameters} for each transition are multivariate normally distributed
#' and mutually independent.
#'
#' @param transitionFunctions a \code{transition.structure} of dimension
#' \eqn{N \times N}{N x N} that contains the hazard functions %history: vector, each
#' element corresponds to a transition number, 0 if transition does not occur
#' @param parameters a \code{transition.structure} of dimension
#' \eqn{N \times N}{N x N} that contains the parameters
#' @param cohortSize a \code{numeric} indicating the number of patients to be
#' simulated.
#' @param parameterCovariances a \code{transition.structure} of dimension
#' \eqn{N \times N}{N x N} of covariance matrices for the \code{parameters}.
#' @param timeToTransition a \code{logical} \code{matrix}; TRUE for all
#' transitions whose \code{transitionFunction} is specified as the time until
#' transition instead of as a hazard function or as a \code{character}.
#' @param baseline a \code{matrix} or \code{data.frame} of dimension
#' \eqn{cohortSize \times M}{cohortSize x M} with \eqn{M} baseline
#' characteristics of subjects to be simulated.
#' @param baselineFunction a \code{function} generating a \code{data.frame} of dimension
#' \eqn{cohortSize \times M}{cohortSize x M} with \eqn{M} baseline
#' characteristics of subjects to be simulated.
#' @param initialState a \code{numeric} of length \code{cohortSize} with the
#' initial state for each subject simulated.
#' @param absorbing a \code{numeric} containing all absorbing states.
#' @param to final time of the simulation.
#' @param report.every a \code{numeric} to check progress of simulation.
#' @param sampler.steps a \code{numeric} indicating number of steps for
#' discretization of hazard functions
#' @return an object of class \code{"ArtCohort"} with \code{time.to.state} slot
#' of dimension \eqn{cohortSize \times N}{cohortSize x N} with entry times for
#' each patient into each of the states.
#' @author Luisa Salazar Vizcaya, Nello Blaser, Thomas Gsponer
#' @seealso \code{\link{generateHazardMatrix}},
#' @seealso \code{\link{generateParameterMatrix}},
#' @seealso \code{\link{generateParameterCovarianceMatrix}},
#' @seealso \code{\link{ArtCohort}}, \code{\link{transitionProbabilities}},
#' @seealso \code{\link{cumulativeIncidence}}
#' @keywords main function
#' @examples
#'
#' # Here is an example model with 3 states and 2 possible transitions.
#'
#' # number of states in the model
#' statesNumber <- 3
#'
#' # cohort size
#' cohortSize <- 100
#'
#' # specification of hazard functions
#' hazardf <- generateHazardMatrix(statesNumber)
#' hazardf[[1,2]] <- function(t, r1, r2)
#' {
#' ifelse(t<=2, r1 , r2)
#' }
#' hazardf[[2,3]] <- "Weibull"
#'
#' # list of parameters for the hazard functions
#' mu <- generateParameterMatrix(hazardf)
#' mu[[1,2]] <- list(r1=0.33, r2=0.03) # r1, r2
#' mu[[2,3]] <- list(shape=1,scale=0.84) # shape, scale
#'
#' # time
#' maxTime <- 10
#'
#' # simulate the cohort
#' cohort <- simulateCohort(
#' transitionFunctions = hazardf,
#' parameters = mu,
#' cohortSize = cohortSize,
#' to=maxTime)
#'
#' # output
#' head(cohort)
#'
#' # transition probability
#' tr <- transitionProbabilities(cohort, times=seq(0,4,.1))
#' plot(tr, ci=FALSE)
#'
#' # cumulative incidence
#' inc <- cumulativeIncidence(cohort, times=seq(0,4,.1))
#' plot(inc, ci=FALSE, states=c(2,3))
#' @import graphics
#' @export simulateCohort
simulateCohort <-
function(transitionFunctions,
parameters,
cohortSize=1000,
parameterCovariances=generateParameterCovarianceMatrix(parameters),
timeToTransition=array(FALSE, dim = dim(transitionFunctions@list.matrix)),
baseline=matrix(NA, nrow = cohortSize),
baselineFunction = baselineFunction_empty,
initialState=rep(1, cohortSize),
absorbing=transitionFunctions@states.number,
to=100, report.every = 100, sampler.steps = 1000){
try(cohortSize <- as.integer(cohortSize))
if (length(cohortSize) > 1)
warning("cohortSize has length > 1 and only the first element will be used")
cohortSize <- cohortSize[1]
try(to <- as.numeric(to))
if (length(to) > 1)
warning("to has length > 1 and only the first element will be used")
to <- to[1]
statesNumber <- dim(transitionFunctions@list.matrix)[1]
print("class(parameters)")
print(class(parameters))
if (class(parameterCovariances) == "transition.structure")
stopifnot(cohortSize>0,
is.list(transitionFunctions@list.matrix),
is.list(parameters@list.matrix),
is.list(parameterCovariances@list.matrix),
dim(transitionFunctions@list.matrix)[1]==dim(transitionFunctions@list.matrix)[2],
identical(dim(transitionFunctions@list.matrix), dim(parameters@list.matrix)),
identical(dim(parameters@list.matrix), dim(parameterCovariances@list.matrix)),
identical(dim(transitionFunctions@list.matrix), dim(timeToTransition)),
initialState %in% 1:statesNumber,
absorbing %in% 1:statesNumber,
length(initialState) %in% c(1, cohortSize)
)
else stopifnot(cohortSize>0,
is.list(transitionFunctions@list.matrix),
is.list(parameters@list.matrix),
is.logical(parameterCovariances),
dim(transitionFunctions@list.matrix)[1]==dim(transitionFunctions@list.matrix)[2],
identical(dim(transitionFunctions@list.matrix), dim(parameters@list.matrix)),
identical(dim(transitionFunctions@list.matrix), dim(timeToTransition)),
initialState %in% 1:statesNumber,
absorbing %in% 1:statesNumber,
length(initialState) %in% c(1, cohortSize)
)
transitionFunction = transitionFunctions@list.matrix
if (class(parameterCovariances) != "transition.structure"){
if (parameterCovariances == FALSE) parameterCovariance = parameterCovariances
else message("parameterCovariances needs to be of class transition.structure or the logical FALSE")
}
else parameterCovariance = parameterCovariances@list.matrix
parameters.in = parameters@list.matrix
historyDep <- length(grep(c("history"), transitionFunction))>0
impossible=matrix(FALSE, nrow=dim(transitionFunction)[1], ncol=dim(transitionFunction)[1])
fixpar=matrix(FALSE, nrow=dim(transitionFunction)[1], ncol=dim(transitionFunction)[1])
hf <- transitionFunction
for (trans in 1:length(transitionFunction)){
if (length(hf[[trans]]) > 0 & is.character(hf[[trans]])) {
if (hf[[trans]]=="impossible"){
impossible[[trans]] <- TRUE
}
}
if (is.function(hf[[trans]])) {
params <- names(formals(hf[[trans]]))
historyIndex <- which(params=="history")
blIndex <- which(params=="bl")
timeIndex <- which(params=="t")
if(length(params[-c(historyIndex, blIndex, timeIndex)])==0) fixpar[[trans]] <- TRUE
if (historyDep){
newParams <- c("t", params[-c(historyIndex, blIndex, timeIndex)], "history", "bl")
}
if (length(c(historyIndex, blIndex, timeIndex))>0 & !historyDep){
newParams <- c("t", params[-c(historyIndex, blIndex, timeIndex)], "bl")
}
if (length(c(historyIndex, blIndex, timeIndex))==0 & !historyDep){
newParams <- c("t", params, "bl")
}
fct <- function(x) {}
forms <- vector("pairlist", length=length(newParams))
for (par in 1:length(forms)){
forms[[par]] <- formals(fct)$x
}
names(forms) <- newParams
formals(hf[[trans]]) <- forms
rm(fct, forms)
}
}
if (class(parameterCovariances) == "logical") sigma = parameterCovariance
else {
sigma <- matrix(0, ncol=length(unlist(parameters.in)), nrow=length(unlist(parameters.in)))
start <- 0
for (matIndex in 1:dim(parameterCovariance)[1]) {
leng <- dim(parameterCovariance[[matIndex]])[1]
if (is.null(leng)) leng <- 0
if (leng>0) {
sigma[(start+1):(start+leng), (start+1):(start+leng)] <- parameterCovariance[[matIndex]]
}
start <- start+leng
}
}
impossible <- auxcounter(statesNumber)[impossible]
fixpar <- auxcounter(statesNumber)[fixpar]
direct <- auxcounter(statesNumber)[timeToTransition]
cohort <- createCohorts(hazardf=hf,
statesNumber=statesNumber,
cohortSize=cohortSize,
mu=parameters.in,
sigma = sigma,
historyl = historyDep,
impossible = impossible,
fixpar = fixpar,
direct = direct,
bl0=baseline,
startingStates=initialState,
absorbing=absorbing,
to=to, report.every = report.every, sampler.steps = sampler.steps)
cohort[cohort>to] <- NA
simulatedcohort <- new("ArtCohort")
simulatedcohort@baseline <- as.matrix(baseline)
simulatedcohort@follow.up <- to
simulatedcohort@parameters <- parameters
simulatedcohort@timeToTransition <- timeToTransition
simulatedcohort@transitionFunctions <- transitionFunctions
simulatedcohort@states.number <- dim(cohort)[1]
simulatedcohort@size <- dim(cohort)[2]
simulatedcohort@time.to.state <- as.data.frame(t(cohort))
if (class(parameterCovariances) == "transition.structure")
simulatedcohort@parametersCovariances <- parameterCovariances
return(simulatedcohort)
}
unfold <-
function (object)
{
vector <- numeric()
k <- 0
for (i in 1:length(object)) {
vector[(1 + k):(length(object[[i]]) + k)] <- c(object[[i]])
k <- length(object[[i]]) + k
}
return(vector)
}
msprep <- function (time, status, data, trans, start, id, keep)
{
if (!(is.matrix(time) | (is.data.frame(time)))) {
if (!is.character(time))
stop("argument \"time\" should be a character vector")
if (missing(data))
stop("missing \"data\" argument not allowed when time argument is character vector")
startings <- which(apply(!is.na(trans), 2, sum) == 0)
wh <- which(is.na(time))
if (!all(wh %in% startings))
stop("no NA's allowed in the \"time\" argument for non-starting states")
tcols <- match(time[!is.na(time)], names(data))
if (any(is.na(tcols)))
stop("at least one of elements of \"time\" not in data")
time <- matrix(NA, nrow(data), length(time))
whcols <- (1:ncol(time))[!(1:ncol(time)) %in% wh]
time[, whcols] <- as.matrix(data[, tcols])
}
if (!(is.matrix(status) | (is.data.frame(status)))) {
if (!is.character(status))
stop("argument \"status\" should be a character vector")
if (missing(data))
stop("missing \"data\" argument not allowed when status argument is character vector")
startings <- which(apply(!is.na(trans), 2, sum) == 0)
wh <- which(is.na(status))
if (!all(wh %in% startings))
stop("no NA's allowed in the \"status\" argument for non-starting states")
dcols <- match(status[!is.na(status)], names(data))
if (any(is.na(dcols)))
stop("at least one of elements of \"status\" not in data")
status <- matrix(NA, nrow(data), length(status))
whcols <- (1:ncol(status))[!(1:ncol(status)) %in% wh]
status[, whcols] <- as.matrix(data[, dcols])
}
time <- as.matrix(time)
status <- as.matrix(status)
if (!all(dim(time) == dim(status)))
stop("unequal dimensions of \"time\" and \"status\" data")
n <- nrow(time)
K <- dim(trans)[1]
if ((dim(trans)[2] != K) | (dim(time)[2] != K))
stop("dimension of \"trans\" does not match with length of \"time\" and \"status\"")
idname <- "id"
if (missing(id))
id <- 1:n
else {
if (!is.vector(id))
stop("argument \"id\" is not a vector")
else {
if (!is.character(id)) {
if (length(id) != n)
stop("argument \"id\" not of correct length")
}
else {
if (length(id) == 1) {
if (n == 1)
stop("cannot determine whether \"id\" argument indicates ")
else {
idname <- id
id <- data[[id]]
}
}
else {
if (length(id) != n)
stop("argument \"id\" not of correct length")
id <- factor(id)
}
}
}
}
idlevels <- NULL
if (is.factor(id))
idlevels <- levels(id)
if (!missing(start)) {
startstate <- start$state
starttime <- start$time
if (length(startstate) != length(starttime))
stop("starting states and times not of equal length")
if (length(startstate) > 1) {
if (length(startstate) != n)
stop("length of starting states and times different from no of subjects in data")
}
else {
startstate <- rep(startstate, n)
starttime <- rep(starttime, n)
}
}
else {
startstate <- rep(1, n)
starttime <- rep(0, n)
}
msres <- msprepEngine(time = time, status = status, id = id,
starttime = starttime, startstate = startstate, trans = trans,
originalStates = (1:nrow(trans)), longmat = NULL)
msres <- as.data.frame(msres)
names(msres) <- c(idname, "from", "to", "trans", "Tstart",
"Tstop", "status")
msres$time <- msres$Tstop - msres$Tstart
msres <- msres[, c(1:6, 8, 7)]
msres <- msres[order(msres[, 1], msres[, 5], msres[, 2],
msres[, 3]), ]
row.names(msres) <- 1:nrow(msres)
if (!is.null(idlevels))
msres[, 1] <- factor(msres[, 1], 1:length(idlevels),
labels = idlevels)
if (!missing(keep)) {
if (!(is.matrix(keep) | (is.data.frame(keep)))) {
if (is.character(keep)) {
if (missing(data))
stop("argument \"data\" is missing, with no default")
nkeep <- length(keep)
kcols <- match(keep, names(data))
if (any(is.na(kcols)))
stop("at least one of elements of \"keep\" not in data")
keepname <- keep
keep <- data[, kcols]
}
else {
nkeep <- 1
keepname <- names(keep)
if (length(keep) != n)
stop("argument \"keep\" has incorrect dimension")
}
}
else {
nkeep <- ncol(keep)
keepname <- names(keep)
if (nrow(keep) != n)
stop("argument \"keep\" has incorrect dimension")
if (nkeep == 1)
keep <- keep[, 1]
}
if (is.null(keepname))
keepname <- paste("keep", as.character(1:nkeep),
sep = "")
if (nkeep > 0) {
if (is.factor(msres[, 1]))
msres[, 1] <- factor(msres[, 1])
tbl <- table(msres[, 1])
if (nkeep == 1) {
ddcovs <- rep(keep, tbl)
ddcovs <- as.data.frame(ddcovs)
names(ddcovs) <- keepname
}
else {
ddcovs <- lapply(1:nkeep, function(i) rep(keep[,
i], tbl))
ddcovs <- as.data.frame(ddcovs)
names(ddcovs) <- keepname
}
msres <- cbind(msres, ddcovs)
}
}
attr(msres, "trans") <- trans
class(msres) <- c("msdata", "data.frame")
return(msres)
}
msprepEngine <- function (time, status, id, starttime, startstate, trans, originalStates, longmat)
{
#if (is.null(nrow(time)))
# return(longmat)
#if (nrow(time) == 0)
# return(longmat)
if (length(time) == 0)
return(longmat)
states.to <- apply(!is.na(trans), 1, sum)
absorbing <- which(states.to == 0)
states.from <- apply(!is.na(trans), 2, sum)
startings <- which(states.from == 0)
newstate <- startstate
newtime <- starttime
to.remove <- NULL
for (starting in startings) {
subjs <- which(startstate == starting)
nstart <- length(subjs)
tostates <- which(!is.na(trans[starting, ]))
transs <- trans[starting, tostates]
nreach <- length(tostates)
if ((nstart > 0) & (nreach > 0)) {
Tstart <- starttime[subjs]
###
if (!is.null(dim(time))) {
Tstop <- time[subjs, tostates, drop = FALSE]
Tstop[Tstop < Tstart] <- Inf
stat <- status[subjs, tostates, drop = FALSE]
smallesttime <- apply(Tstop, 1, min)
hlp <- Tstop * 1/stat
hlp[Tstop == 0 & stat == 0] <- Inf
nexttime <- apply(hlp, 1, min)
censored <- which(is.infinite(apply(hlp, 1, min)))
}
else {
Tstop <- time[tostates, drop = FALSE]
Tstop[Tstop < Tstart] <- Inf
stat <- status[tostates, drop = FALSE]
smallesttime <- min(Tstop)
hlp <- Tstop * 1/stat
hlp[Tstop == 0 & stat == 0] <- Inf
nexttime <- min(hlp)
censored <- which(is.infinite(min(hlp)))
}
wh <- which(smallesttime < nexttime)
whminc <- setdiff(wh, censored)
if (length(whminc) > 0) {
whsubjs <- id[subjs[whminc]]
whsubjs <- paste(whsubjs, collapse = " ")
warning("From starting state ", originalStates[starting],
", subject ", whsubjs, " has smallest transition time with status=0, larger transition time with status=1")
}
nexttime[censored] <- smallesttime[censored]
if (!is.null(dim(time))){
if (ncol(hlp) > 1) {
hlpsrt <- t(apply(hlp, 1, sort))
warn1 <- which(hlpsrt[, 1] - hlpsrt[, 2] == 0)
if (length(warn1) > 0) {
isw <- id[subjs[warn1]]
isw <- paste(isw, collapse = " ")
hsw <- hlpsrt[warn1, 1]
hsw <- paste(hsw, collapse = " ")
warning("Starting from state ", originalStates[starting],
", simultaneous transitions possible for subjects ",
isw, " at times ", hsw, "; smallest receiving state chosen")
}
}
}
if (length(censored) > 0) {
if (!is.null(dim(hlp))) nextstate <- apply(hlp[-censored, , drop = FALSE], 1, which.min)
else nextstate <- which.min(hlp[-censored , drop = FALSE])
reachAbsorb <- (1:nstart)[-censored][which(tostates[nextstate] %in%
absorbing)]
}
else {
if (!is.null(dim(time))) nextstate <- apply(hlp, 1, which.min)
else nextstate <- which.min(hlp)
reachAbsorb <- (1:nstart)[which(tostates[nextstate] %in%
absorbing)]
}
statmat <- matrix(0, nstart, nreach)
if (length(censored) > 0)
statmatmin <- statmat[-censored, , drop = FALSE]
else statmatmin <- statmat
if (nrow(statmatmin) > 0)
statmatmin <- t(sapply(1:nrow(statmatmin), function(i) {
x <- statmatmin[i, ]
x[nextstate[i]] <- 1
return(x)
}))
if (length(censored) > 0)
statmat[-censored, ] <- statmatmin
else statmat <- statmatmin
mm <- matrix(c(rep(id[subjs], rep(nreach, nstart)),
rep(originalStates[starting], nreach * nstart),
rep(originalStates[tostates], nstart), rep(transs,
nstart), rep(Tstart, rep(nreach, nstart)),
rep(nexttime, rep(nreach, nstart)), as.vector(t(statmat))),
nreach * nstart, 7)
longmat <- rbind(longmat, mm)
to.remove <- c(to.remove, subjs[c(censored, reachAbsorb)])
if (length(censored) > 0)
newstate[subjs[-censored]] <- tostates[nextstate]
else newstate[subjs] <- tostates[nextstate]
if (length(censored) > 0)
newtime[subjs[-censored]] <- nexttime[-censored]
else newtime[subjs] <- nexttime
}
}
if (length(to.remove) > 0) {
if (!is.null(dim(time))){
time <- time[-to.remove, ]
status <- status[-to.remove, ]
}
else {
time <- time[-to.remove]
status <- status[-to.remove]
}
newtime <- newtime[-to.remove]
newstate <- newstate[-to.remove]
id <- id[-to.remove]
}
K <- nrow(trans)
idx <- rep(1, K)
idx[startings] <- 0
idx <- cumsum(idx)
newstate <- idx[newstate]
if (!is.null(dim(time))) {
Recall(time = time[, -startings], status = status[, -startings],
id = id, starttime = newtime, startstate = newstate,
trans = trans[-startings, -startings], originalStates = originalStates[-startings],
longmat = longmat)
}
else {
Recall(time = time[-startings], status = status[ -startings],
id = id, starttime = newtime, startstate = newstate,
trans = trans[-startings, -startings], originalStates = originalStates[-startings],
longmat = longmat)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.