Nothing
R/AUC.test.R
In markovMSM: Methods for Checking the Markov Condition in Multi-State Survival Data
Defines functions
AUC.test
Documented in
AUC.test
#' Global test and local test for checking the Markov condition on Multi-state
#' Models based on the Area Under the two curves (AUC) given by the non-Markov
#' estimators of the transition probabilities to the Markov Aalen-Johansen
#' estimators.
#' @description This function is used to obtain the local and global tests for
#' checking the Markov condition.
#' @param data A data frame in long format containing the subject \code{id};
#' \code{from} corresponding to the starting state;the receiving state,
#' \code{to}; the transition number, \code{trans}; the starting time of the
#' transition given by \code{Tstart}; the stopping time of the transition,
#' \code{Tstop}, and \code{status} for the status variable, with 1 indicating
#' an event (transition), 0 a censoring.
#' @param from The starting state of the transition probabilities.
#' @param to The last receiving state considered for the estimation of the
#' transition probabilities. All the probabilities among the first and the
#' last states are also computed.
#' @param type Type of test for checking the Markov condition: "local" or
#' "global". By default type='global'.
#' @param times For the local test, times represents the starting times of the
#' transition probabilities. In case of a global test, the argument is given by
#' times between the minimum time and the third quartile times used
#' in the formula of this test. Default to NULL.
#' @param quantiles Quantiles used in the formula of the Global test for the
#' AUC methods. By default takes the percentiles 0.05, 0.10, 0.20, 0.30, 0.40.
#' @param tmat The transition matrix for multi-state model.
#' @param replicas Number of replicas for the Monte Carlo simulation to
#' standardization of the T-statistic given by the difference of the areas of
#' AJ and LM transition probabilities estimates.
#' @param limit Percentile of the event time used as the upper bound for the
#' computation of the AUC-based test.
#' @param positions List of possible transitions; x[[i]] consists of a vector of
#' state numbers reachable from state i.
#' @param namesStates A character vector containing the names of either the
#' competing risks or the states in the multi-state model specified by the
#' competing risks or illness-death model. names should have the same length as
#' the list x (for transMat), or either K or K+1 (for trans.comprisk), or 3
#' (for trans.illdeath).
#' @param timesNames Either 1) a matrix or data frame of dimension n x S
#' (n being the number of individuals and S the number of states in the
#' multi-state model), containing the times at which the states are visited or
#' last follow-up time, or 2) a character vector of length S containing the
#' column names indicating these times. In the latter cases, some elements of
#' time may be NA, see Details
#' @param status Either 1) a matrix or data frame of dimension n x S,
#' containing, for each of the states, event indicators taking the value 1 if
#' the state is visited or 0 if it is not (censored), or 2) a character
#' vector of length S containing the column names indicating these status
#' variables. In the latter cases, some elements of status may be NA,
#' see Details.
#'
#' @return In case of the AUC global test, an object with a list with the
#' following outcomes:
#' \item{globalTest}{p-value of AUC global tests for each transition. These
#' values are obtained through the minimum of the means of each two contiguous
#' quantiles times of the AUC global tests.}
#' \item{localTest}{AUC local tests of the transition probability for each times
#' and transitions.}
#' \item{quantiles}{Quantiles times used for the AUC global tests.}
#' \item{times}{Times used for the AUC global tests.}
#' \item{DIF}{Differences between the AJ and the LMAJ estimates for each
#' transition probabilites from the starting state until the receiving state
#' given by only one replica where 's' represent each of the quantile times.}
#' \item{from}{The starting state considered for the AUC global tests.}
#' \item{to}{The last receiving state considered for the the AUC Local tests.}
#' \item{ET.qiAll}{The lower limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qi2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{ET.qsAll}{The upper limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qs2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{replicas}{Number of replicas for the Monte Carlo simulation.}
#' \item{limit}{Percentil of the times used in the AUC global tests.}
#'
#' In case of the AUC local test, an object with a list with the following
#' outcomes:
#' \item{localTest}{p-value of AUC local tests for each times and transitions.}
#' \item{trans}{The transition matrix describing the states and transitions of
#' the multi-state model.}
#' \item{times}{Times selected for the AUC Local tests.}
#' \item{DIF}{Differences between the AJ and the LMAJ estimates for each
#' transition probabilites from the starting state until the receiving state
#' given by only one replica where 's' represent each of the quantile times.}
#' \item{from}{The starting state considered for the AUC Local tests.}
#' \item{to}{The last receiving state considered for the the AUC Local tests.}
#' \item{ET.qiAll}{The lower limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qi2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{ET.qsAll}{The upper limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qs2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{replicas}{Number of replicas for the Monte Carlo simulation.}
#' \item{limit}{Percentil of the times used in the AUC local tests.}
#' @references Soutinho G, Meira-Machado L (2021). Methods for checking the
#' Markov condition in multi-state survival data. \emph{Computational Statistics}.
#' de Una-alvarez J, Meira-Machado L (2015). Nonparametric estimation of
#' transition probabilities in the non-Markov illness-death model: A comparative
#' study. \emph{Biometrics}, 71(2), 364-375.
#' Putter H, Spitoni C (2018). Non-parametric estimation of transition
#' probabilities in non-Markov multi-state models: The landmark Aalen-Johansen
#' estimator. \emph{Statistical Methods in Medical Research}, 27, 2081-2092.
#' Meira-Machado L, Sestelo M (2019). Estimation in the progressive illness-death
#' model: A nonexhaustive review. \emph{Biometrical Journal}, 61(2), 245-263.
#' @examples
#' \donttest{
#' set.seed(1234)
#' library(markovMSM)
#' data("colonMSM")
#' db_wide<-colonMSM
#' positions<-list(c(2, 3), c(3), c())
#' namesStates = c("Alive", "Rec", "Death")
#' tmat <-transMatMSM(positions, namesStates)
#' timesNames = c(NA, "time1","Stime")
#' status=c(NA, "event1","event")
#' trans = tmat
#' db_long<- prepMSM(data=db_wide, trans, timesNames, status)
#' res<-AUC.test(data=db_long, times=180, from=1, to=3, type='local',
#' replicas = 20, tmat = tmat)
#' res$localTest
#' res2<-AUC.test(data=db_long, times=180, from=2, to=3, type='local',
#' replicas = 20, tmat = tmat)
#' res2$localTest
#'
#' res3<-AUC.test(data=db_long, from=1, to=3, replicas = 20, tmat=tmat)
#' round(res3$globalTest,3)
#'
#' res4<-AUC.test(data=db_long, from=2, to=3, type='global', replicas = 20,
#' tmat=tmat)
#' round(res4$globalTest,3)
#'
#' round(res4$localTest,3)
#' round(res4$localTest,3)
#'
#' #AUC global for the individuals with the treatment "Obs" of covariate `"rx"
#' #for colonMSM data set
#'
#' set.seed(12345)
#'
#' db_wide.obs<-db_wide[db_wide$rx=='Obs',]
#'
#' db_long.obs <- prepMSM(data = db_wide.obs, trans, timesNames, status)
#'
#' res3a<-AUC.test(data=db_long.obs, times=365, from=1, to=3,
#' type='local', replicas= 20, tmat = tmat)
#'
#' res3a$localTest
#'
#' set.seed(12345)
#'
#' res4a<-AUC.test(data=db_long.obs, times=365, from=2, to=3,
#' type='local', replicas= 20, tmat = tmat)
#'
#' res4a$localTest
#'
#' #' data("ebmt4")
#' db_wide <- ebmt4
#' positions <- list(c(2, 3, 5, 6), c(4, 5, 6), c(4, 5, 6),
#' c(5, 6), c(), c())
#' state.names <- c("Tx", "Rec", "AE", "Rec+AE", "Rel", "Death")
#' tmat <-transMatMSM(positions, state.names)
#' timesNames <- c(NA, "rec", "ae","recae", "rel", "srv")
#' status <- c(NA, "rec.s", "ae.s", "recae.s","rel.s", "srv.s")
#' trans <- tmat
#'
#' db_long <- prepMSM(data=db_wide, trans, timesNames, status)
#' db_long[1:10,]
#'
#' res5<-AUC.test(data=db_long, from=1, to=5, type='global',
#' quantiles=c(.05, .10, .20, .30, 0.40),
#' tmat = tmat, replicas = 5,
#' positions=positions, namesStates=state.names,
#' timesNames=timesNames, status=status)
#'
#' round(res5$globalTest, 4)
#' round(res5$localTests,4)
#'
#' res6<-AUC.test(data = prothr, from=2, to=3,
#' type='global', replicas= 5, limit=0.90,
#' quantiles=c(.05, .10, .20, .30, 0.40))
#' round(res6$globalTest,4)
#'
#' round(res6$localTests,4)
#'
#' }
#' @author Gustavo Soutinho and Luis Meira-Machado.
#'
#' @importFrom "survival" coxph Surv survfit strata untangle.specials
#' @importFrom "graphics" legend abline axis legend lines matplot par plot polygon
#' @importFrom "stats" pchisq pnorm quantile sd na.omit terms approxfun as.formula rnorm rpois weighted.mean
#' @importFrom "utils" capture.output flush.console
#' @importFrom "mstate" transMat msprep msfit probtrans LMAJ xsect cutLMms msdata2etm trans2Q events
#' @importFrom "stats" model.matrix model.frame model.response model.offset
#' @importFrom "stats" delete.response delete.response
#' @export AUC.test
#' @export prepMSM
#' @export transMatMSM
#' @export eventsMSM
AUC.test<- function(data, from=1, to=3, type='global',
times=NULL, quantiles=c(.05,.10, .20, .30, 0.40),
tmat=NULL, replicas=10, limit=0.90,
positions=list(c(2, 3), c(3), c()),
namesStates = c("Alive", "Rec", "Death"),
timesNames = c(NA, "time1","Stime"),
status=c(NA, "event1","event")){
db_long<-data
if (class(db_long)[1]!="msdata")
stop("Argument 'data' must be a mstate object")
db_wide<-NULL
if(type=='global'){
i<-from
#i<-2
j<-to
#j<-5
M<-replicas
DIF<-NULL
areaAll<-NULL
ETAll<-NULL
p.valueAll<-NULL
ET_bootAll<-NULL
ET.qiAll<-NULL
ET.qsAll<-NULL
ET.qi2All<-NULL
ET.qs2All<-NULL
times<-sort(times)
if(is.null(quantiles) & is.null(times)){
stop("One of the arguments quantiles or times must be NULL.")
}
if(!is.null(quantiles) & !is.null(times)){
stop("Only one of the arguments quantiles or times must be NULL.")
}
if(length(times)!=5 & !is.null(times)){
stop("The length of times must be 5.")
}
#min e quantile 0.75
if(!is.null(times)){
if(!is.null(db_wide)){
if(length(namesStates)==3){
timeToQuantiles<-timesNames[2]
min<-min(as.numeric(db_wide[, timeToQuantiles]))
val0.75<-as.numeric(quantile(db_wide[, timeToQuantiles], 0.75))
}else{#inicio difere de 3
if(i==1){#from
#head(db_wide)
valToQuantiles<-rep(0,nrow(db_wide))
for(a1 in 1:nrow(db_wide)){
#a1<-1
valToQuantiles[a1]<-min(db_wide[a1, timesNames[-1]])
}
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}else{
valToQuantiles<-db_wide[,timesNames[i]]
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}
min<-min(as.numeric(valToQuantiles))
val0.75<-as.numeric(quantile(valToQuantiles, 0.75))
}#fim de casos em que length(names) difere de 3
}else{#caso db_wide=NULL
tempQ<-unique(db_long$Tstop)
min<-min(as.numeric(tempQ))
val0.75<-as.numeric(quantile(tempQ, 0.75))
}
if(length(which(times<val0.75 & times>min & !is.null(times)))!=5){
stop("All the times must be greater than the minumum time and the third quantile of time.")
}
}
#
if(is.null(times))
times.to.test<-length(quantiles)
if(is.null(quantiles))
times.to.test<-length(times)
#start.time <- Sys.time() #marcacao do tempo inicial de processamento
cat("This method is based on the calculation of 5 percentile times and can lead to a high execution time.
Partial execution times will be provided during execution.", '\n')
for(h in 1:times.to.test){
#h<-1
#print(h) - estava este
#for(i3 in 1:length(colnames(db_wide))){
#i3<-1
# if (colnames(db_wide)[i3]==timeToQuantiles)
# indexS<-i3
#}
#s<-as.numeric(quantile(db_wide[,indexS], quantiles))[h]
if(!is.null(db_wide)){
if(length(namesStates)==3){
if(is.null(times)){
timeToQuantiles<-timesNames[2]
s<-as.numeric(quantile(db_wide[, timeToQuantiles], quantiles))[h]
s_quant<-as.numeric(quantile(db_wide[, timeToQuantiles], quantiles))
}else{
s<-as.numeric(times)[h]
s_quant<-as.numeric(times)
}
}else{ #inicio difere de 3
if(is.null(times)){
if(i==1){#from
#head(db_wide)
valToQuantiles<-rep(0,nrow(db_wide))
for(a1 in 1:nrow(db_wide)){
#a1<-1
valToQuantiles[a1]<-min(db_wide[a1, timesNames[-1]])
}
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}else{
valToQuantiles<-db_wide[,timesNames[i]]
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}
s<-as.numeric(quantile(valToQuantiles, quantiles))[h]
s_quant<-as.numeric(quantile(valToQuantiles, quantiles))
}else{
s<-as.numeric(times)[h]
s_quant<-as.numeric(times)
}
}#fim de casos em que length(names) difere de 3
}else{#......................caso db_wide=NULL
if(is.null(times)){
#head(db_long)
#tempQ<-unique(db_long$Tstop)
#summary(tempQ)
if(length(namesStates)==3){
if(nrow(db_long)==2152){
#head(db_long)
tempQ<-unique(db_long$Tstop)
#summary(tempQ)
s<-as.numeric(quantile(tempQ, quantiles))[h]
tmat <- attr(db_long, "trans")
s_quant<-as.numeric(quantile(tempQ, quantiles))
}else{
values<-rep(0,length(unique(db_long$id)))
for(i3 in 1:length(unique(db_long$id))){
#i3<-1
values[i3]<-min(unique(db_long[db_long$id==i3,'Tstop']))
}
tempQ<-as.numeric(values)
#tempQ==as.numeric(db_wide[, timeToQuantiles])
#tempQ<-unique(db_long$Tstop)
min<-min(as.numeric(tempQ))
val0.75<-as.numeric(quantile(tempQ, 0.75))
s<-as.numeric(quantile(tempQ, quantiles))[h]
tmat <- attr(db_long, "trans")
s_quant<-as.numeric(quantile(tempQ, quantiles))
}
}else{#diferente de illness death model
#from=1
if(i==1){
values<-rep(0,length(unique(db_long$id)))
for(i3 in 1:length(unique(db_long$id))){
#i3<-2
#db_long[db_long$id==i3,]
#db_long[db_long$id==i3 ,'Tstop']
values[i3]<-min(unique(db_long[db_long$id==i3 ,'Tstop']))
}
valToQuantiles<-values
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
s<-as.numeric(quantile(valToQuantiles, quantiles))[h]
s_quant<-as.numeric(quantile(valToQuantiles, quantiles))
}else{#from !=1 ....
db2<-db_long[db_long$from==i,]
ids2<-unique(db2$id)
values2<-rep(0,length(ids2))
for(i4 in 1:length(ids2)){
#i4<-3
#db_long[db_long$id==ids2[i4],]
#db_long[db_long$id==ids2[i4] & db_long$from==i,'Tstop']
values2[i4]<-min(unique(db_long[db_long$id==ids2[i4] & db_long$from==i,'Tstop']))
}
valToQuantiles<-values2
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
s<-as.numeric(quantile(valToQuantiles, quantiles))[h]
s_quant<-as.numeric(quantile(valToQuantiles, quantiles))
}
} #fim diferente de illness death model (3 estados)
}else{#valores em concreto quando nao se conhece WIDE
s<-as.numeric(times)[h]
s_quant<-as.numeric(times)
}
}
#h<-2
#s<-c(10, 30, 100, 200, 300, 650)[h] #10 erro
#cat("Time =", s,"\n")
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans),
data = db_long))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood",
trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
#resAJ[[i]]
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, limit)
tempos <- tempos[tempos < t.limit]
#trans<-as.numeric(na.omit(as.numeric(resAJ$trans)))
#resAJ[[i]][1:10,]
LMpt0 <- suppressWarnings(LMAJ2(db_long, s=s, from=i))
#LMpt0[200:250,]
#length(LMpt0$time)
#tempos[tempos %in% LMpt0$time]
#LMpt0$time[LMpt0$time %in% tempos]
#tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
#head(resAJf[,2:(j+1)])
#head(LMpt0f[,2:(j+1)])
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,] #; pr
#head(pr)
tempos2 <- diff(tComuns)
#length(tempos2)
#tempos2[1]*pr[1,]
#tempos2[2]*pr[2,]
#(pr * tempos2)[1:2,]
area<-colSums((pr * tempos2))
areaAll<-rbind(areaAll, cbind(s, as.data.frame(t(area))))
#areaAll
areasBoot<-matrix(NA, nrow = M, ncol = j) #ncol(resAJ$trans)
prTimes1<-cbind(tComuns, pr.ini)
DIF<-rbind(DIF, cbind(s, prTimes1))
#none<-NULL (para filtro...)
#for(n1 in from:to){
# #n1<-1
# if(n1==from)
# none<-c(none,unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]-1))
# else
# none<-c(none, unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]))
#}
ET_boot<-NULL
start.time <- Sys.time() #marcacao do tempo inicial de processamento
for(k in 1:M){
if(k==6){
ord1<-c('first', 'second', 'third', 'fourth', 'fifth', 'sixth', 'seventh', 'eighth', 'ninth', 'tenth',
'eleventh', 'twelfth')
ord2<-c('one', 'two', 'three', 'four', 'five', 'six', 'seven', 'eight', 'nine', 'ten',
'eleven', 'twelve')
#M<-100
#(as.numeric(difftime(Sys.time(), start.time))/5*M)/60
if(h==1 | h==2){
if(nrow(db_long)==15512){
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*1.5*M)/5)/60, '\r', sep='')
}else{
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*2.5*M)/5)/60, '\r', sep='')
}
}else{
if(h==3){
if(nrow(db_long)==15512){
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*M)/5)/60, '\r', sep='')
}else{
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*2.1*M)/5)/60, '\r', sep='')
}
}else{#nao ser h=3
if(nrow(db_long)==15512){
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*M)/5)/60, '\r', sep='')
}else{
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*1.3*M)/5)/60, '\r', sep='')
}
}#final de nao ser h=3
} #final de ser 3 ou nao ser
}
#k<-1
#cat("Monte Carlo sample =", k,"\n")
#head(db_wide)
#unique(db_wide$id)
if(!is.null(db_wide)){
n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
db_wide2 <- db_wide[pos,]
#head(db_wide2)
#positions<-list(c(2, 3), c(3), c())
#namesStates = names = c("Alive", "Rec", "Death")
tmat <-transMatMSM(positions, namesStates)
#times = c(NA, "time1","Stime")
#status=c(NA, "event1","event")
trans = tmat
db_long2<- prepMSM(data=db_wide2, trans, timesNames, status)
#db_long2<- msprep(data = db_wide2, trans = tmat,
# time = c(NA, "rec", "ae","recae", "rel", "srv"),
# status = c(NA, "rec.s", "ae.s", "recae.s","rel.s", "srv.s"),
# keep = c("match", "proph", "year", "agecl"))
}else{
n1<-length(unique(db_long$id))
#n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
#table(pos)
pos2<-sort(pos)
db_long2<-NULL
for(i4 in 1:n1){
db_long2 <- rbind(db_long2,db_long[db_long$id %in% pos2[i4],])
}
#db_long2[db_long2$id==3,]
}
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans),
data = db_long2))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood", trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, limit)
tempos <- tempos[tempos < t.limit]
#
msdata<-db_long2
xss <- suppressWarnings(xsect(msdata, s))
infrom <- xss$id[xss$state %in% from]
msdatas <- suppressWarnings(cutLMms(msdata, LM = s))
msdatasfrom <- msdatas[msdatas$id %in% infrom, ]
#print(nrow(msdatasfrom))
if(nrow(msdatasfrom)<2){
print('Error')
}else{
LMpt0 <- suppressWarnings(LMAJ2(db_long2, s=s, from=i))
#tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,]; pr
tempos2b <- diff(tComuns)
areasBoot[k,]<-colSums((pr * tempos2b))
#para o grafico
prTimes<-cbind(tComuns, pr.ini)
prTimes1$tComuns %in% prTimes$tComuns
dataframe<-as.data.frame(prTimes1$tComuns)
colnames(dataframe)<-'tComuns'
ET_boot_rep<-merge(x=dataframe, y=prTimes, by='tComuns',all.x = T)
ET_boot<-rbind(ET_boot,cbind(ET_boot_rep,rep(k,nrow(ET_boot_rep)), rep(s,nrow(ET_boot_rep))))
#ET_boot[200:350,]
#cat(k, '\n')
#cat(h, '\n')
}
#
}
sd_areasBoot<-rep(0, j) # ncol(resAJ$trans)
for(l in 1: j){#ncol(resAJ$trans)
#l<-1
sd_areasBoot[l]<-sd(areasBoot[,l], na.rm = T)
}
ET<-area/sd_areasBoot; ET
ETAll<-rbind(ETAll, cbind(s, as.data.frame(t(ET))))
#none<-unique(none)
#if(length(none)==1){
# p.value <- 2*(1-pnorm(abs(ET))); p.value
# for(n2 in from:to){
# p.value[n2]<-'None'
# }
#}else{
# p.value <- 2*(1-pnorm(abs(ET))); p.value
#}
p.value <- 2*(1-pnorm(abs(ET))); p.value
p.valueAll<-rbind(p.valueAll, cbind(s, as.data.frame(t(p.value ))))
p.valueAll
#
colnames(ET_boot)[(j+2)]<-'replica'
colnames(ET_boot)[(j+3)]<-'s'
#head(ET_boot)
len<-unique(ET_boot$tComuns)
ET.qi<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n in 2:(j+1)){
#n<-2
#ET_boot[ET_boot$tComuns==len[m],2]
ET.qi[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = 0.025,na.rm = T))
ET.qs[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = .975,na.rm = T))
}
}
ET.qi[,1]<-len
ET.qs[,1]<-len
ET_boot[1:20,]
ET_boot2<-ET_boot
for(p in 1:M){
#p<-1
#ET_boot[ET_boot$replica==p,]
for(q in 2:(j+1)){
#q<-2
#ET_boot[ET_boot$replica==p,q]
for (r in 1:length(len)){
#r<-14
ET_boot2[ET_boot2$replica==p,q][r]<-ifelse(is.na(ET_boot2[ET_boot2$replica==p,q][r]),
ET_boot2[ET_boot2$replica==p,q][(r-1)],
ET_boot2[ET_boot2$replica==p,q][r])
}
}
}
#ET_boot2[ET_boot2$replica==p,]
ET.qi2<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs2<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m2 in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n2 in 2:(j+1)){
#n<-2
#ET_boot2[ET_boot2$tComuns==len[m],4]
ET.qi2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = 0.025,na.rm = T))
ET.qs2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = .975,na.rm = T))
}
}
ET.qi2[,1]<-len
ET.qs2[,1]<-len
ET.qiAll<-rbind(ET.qiAll,cbind(as.data.frame(ET.qi), rep(s,length(len))))
ET.qi2All<-rbind(ET.qi2All,cbind(as.data.frame(ET.qi2), rep(s,length(len))))
ET.qsAll<-rbind(ET.qsAll,cbind(as.data.frame(ET.qs), rep(s,length(len))))
ET.qs2All<-rbind(ET.qs2All,cbind(as.data.frame(ET.qs2), rep(s,length(len))))
}#fim times
p.value.f<-p.valueAll
#round(p.value.f,4)
minAll<-rep(0, j)
if(is.null(times)){
minInt<-rep(0, (length(quantiles)-1))
}else{
minInt<-rep(0, (length(times)-1))
}
for (i2 in 1:j){
#i2<-1
for(j2 in 1: (nrow(p.valueAll)-1)){
#j2<-1
p.valueAll[,(i2+1)]<-ifelse(p.valueAll[,(i2+1)]== 'NaN', NA, p.valueAll[,(i2+1)])
minInt[j2]<-mean(p.valueAll[j2:(j2+1),(i2+1)],na.rm = T)
}
#is.na(minInt)
minAll[i2]<-min(minInt, na.rm=T)
}
minAll
#head(ET.qiAll)
names(ET.qiAll)[length(names(ET.qiAll))]<-'s'
names(ET.qiAll)[1]<-'tComuns'
names(ET.qiAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qi2All)
names(ET.qi2All)[length(names(ET.qi2All))]<-'s'
names(ET.qi2All)[1]<-'tComuns'
names(ET.qi2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qsAll)
names(ET.qsAll)[length(names(ET.qsAll))]<-'s'
names(ET.qsAll)[1]<-'tComuns'
names(ET.qsAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qs2All)
names(ET.qs2All)[length(names(ET.qs2All))]<-'s'
names(ET.qs2All)[1]<-'tComuns'
names(ET.qs2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
if(!is.null(db_wide)){
p.value.f2<-p.value.f[,c(1, (from+1):ncol(p.value.f))]
minAllf<-minAll[from:length(minAll)]
}else{
p.value.f2<-p.value.f
minAllf<-minAll
}
nomes<-rep("s",length(names(p.value.f2)))
names(p.value.f2)
for(i5 in 2:length(names(p.value.f2))){
#i5<-2
#nchar(names(p.valueAll.f)[i5])
nomes[i5]<-paste(from,'->',substr(names(p.value.f2)[i5], 7, nchar(names(p.value.f2)[i5])),sep='')
}
names(p.value.f2)<-nomes
#p.value.f2
minAllf2<-as.data.frame(t(minAllf))
names(minAllf2)<-nomes[2:length(nomes)]
if(is.null(times)){
res <-list(globalTest=minAllf2, localTests=p.value.f2,
quantiles=as.numeric(s_quant), times=NULL, DIF=DIF, from=from, to=to,
ET.qiAll=ET.qiAll, ET.qi2All=ET.qi2All, ET.qsAll=ET.qsAll, ET.qs2All=ET.qs2All,
replicas=M, limit=limit)
}else{
res <-list(globalTest=minAllf2, localTests=p.value.f2,
quantiles=NULL,times=times, DIF=DIF, from=from, to=to,
ET.qiAll=ET.qiAll, ET.qi2All=ET.qi2All, ET.qsAll=ET.qsAll, ET.qs2All=ET.qs2All,
replicas=M, limit=limit)
}
class(res) <- c("globalTest", "markovMSM")
res$call <- match.call()
#return(res)
options(warn=-1)
suppressWarnings(res)
return(invisible(res))
}else{
#local
i<-from
#i<-1
j<-to
#j<-6
M<-replicas
DIF<-NULL
areaAll<-NULL
ETAll<-NULL
p.valueAll<-NULL
AJall<-NULL
LMAJall<-NULL
ET.qiAll<-NULL
ET.qsAll<-NULL
ET.qi2All<-NULL
ET.qs2All<-NULL
if(is.null(db_wide))
tmat <- attr(db_long, "trans")
cat("This method may require high execution time.
Partial execution times will be provided during execution.", '\n')
for(h in 1:length(times)){
#h<-2
#print(h)
s<-times[h]
#cat("Time =", s,"\n")
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = db_long))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood", trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, 0.90)
tempos <- tempos[tempos < t.limit]
#trans<-as.numeric(na.omit(as.numeric(resAJ$trans)))
resAJ[[1]] #from 1 to 6
#resAJ[[i]][1:10,1:8]
LMpt0 <- suppressWarnings(LMAJ2(db_long, s=s, from=i))
#LMpt0[1:40,1:8]
length(LMpt0$time)
tempos[tempos %in% LMpt0$time]
LMpt0$time[LMpt0$time %in% tempos]
tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,]#; pr
#head(pr)
pr<-pr[-1,]
#head(pr)
tempos2 <- diff(tComuns)
#length(tempos2)
#tempos2[1]*pr[1,]
#tempos2[2]*pr[2,]
#pr[1:10,]
#tempos2[1:10]
#(pr * tempos2)[1:10,]
area<-colSums((pr * tempos2), na.rm = T)
areaAll<-rbind(areaAll, cbind(s, as.data.frame(t(area))))
#areaAll
areasBoot<-matrix(NA, nrow = M, ncol = j) #ncol(resAJ$trans)
prTimes1<-cbind(tComuns, pr.ini)
DIF<-rbind(DIF, cbind(s, prTimes1))
AJall<-rbind(AJall, cbind(s, resAJf[1:(j+1)]))
LMAJall<-rbind(LMAJall, cbind(s, LMpt0f[1:(j+1)]))
none<-NULL
for(n1 in from:to){
#n1<-1
if(n1==from)
none<-c(none,unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]-1))
else
none<-c(none, unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]))
}
ET_boot<-NULL
start.time <- Sys.time() #marcacao do tempo inicial de processamento
for(k in 1:M){
#k<-1
if(k==6){
cat("Expected computation time in minutes for ", s , ": " , ((as.numeric(difftime(Sys.time(), start.time))*2.3*M)/5)/60, '\r', sep='')
}
#cat("Monte Carlo sample =", k,"\n")
#head(db_wide)
#unique(db_wide$id)
if(!is.null(db_wide)){
n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
db_wide2 <- db_wide[pos,]
#head(db_wide2)
#positions<-list(c(2, 3), c(3), c())
#namesStates = names = c("Alive", "Rec", "Death")
tmat <-transMatMSM(positions, namesStates)
#times = c(NA, "time1","Stime")
#status=c(NA, "event1","event")
trans = tmat
db_long2<- prepMSM(data=db_wide2, trans, timesNames, status)
#db_long2<- msprep(data = db_wide2, trans = tmat,
# time = c(NA, "rec", "ae","recae", "rel", "srv"),
# status = c(NA, "rec.s", "ae.s", "recae.s","rel.s", "srv.s"),
# keep = c("match", "proph", "year", "agecl"))
}else{
n1<-length(unique(db_long$id))
#n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
#table(pos)
pos2<-sort(pos)
db_long2<-NULL
for(i4 in 1:n1){
db_long2 <- rbind(db_long2,db_long[db_long$id %in% pos2[i4],])
}
#db_long2[db_long2$id==3,]
}
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = db_long2))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood", trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, 0.90)
tempos <- tempos[tempos < t.limit]
LMpt0 <- suppressWarnings(LMAJ2(db_long2, s=s, from=i))
#tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
#resAJf<-resAJ[[1]][resAJ[[1]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
#o que estava
#pr<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
#pr<-pr[-1,]; pr
#tempos2 <- diff(tComuns)
#areasBoot[k,]<-colSums((pr * tempos2), na.rm = T) #,
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,]; pr
tempos2b <- diff(tComuns)
areasBoot[k,]<-colSums((pr * tempos2b))
#para o grafico
prTimes<-cbind(tComuns, pr.ini)
prTimes1$tComuns %in% prTimes$tComuns
dataframe<-as.data.frame(prTimes1$tComuns)
colnames(dataframe)<-'tComuns'
ET_boot_rep<-merge(x=dataframe, y=prTimes, by='tComuns',all.x = T)
ET_boot<-rbind(ET_boot,cbind(ET_boot_rep,rep(k,nrow(ET_boot_rep)), rep(s,nrow(ET_boot_rep))))
}
sd_areasBoot<-rep(0, j) # ncol(resAJ$trans)
for(l in 1: j){#ncol(resAJ$trans)
#l<-3
sd_areasBoot[l]<-sd(areasBoot[,l], na.rm = T)
}
ET<-area/sd_areasBoot; ET
ETAll<-rbind(ETAll, cbind(s, as.data.frame(t(ET))))
#ETAll
#none<-unique(none)
#if(length(none)==1){
# p.value <- 2*(1-pnorm(abs(ET))); p.value
# for(n2 in from:to){
# p.value[n2]<-'None'
# }
#}else{
# p.value <- 2*(1-pnorm(abs(ET))); p.value
#}
p.value <- 2*(1-pnorm(abs(ET))); p.value
p.valueAll<-rbind(p.valueAll, cbind(s, as.data.frame(t(p.value ))))
p.valueAll
#
colnames(ET_boot)[(j+2)]<-'replica'
colnames(ET_boot)[(j+3)]<-'s'
#head(ET_boot)
len<-unique(ET_boot$tComuns)
ET.qi<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n in 2:(j+1)){
#n<-2
#ET_boot[ET_boot$tComuns==len[m],2]
ET.qi[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = 0.025,na.rm = T))
ET.qs[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = .975,na.rm = T))
}
}
ET.qi[,1]<-len
ET.qs[,1]<-len
ET_boot[1:20,] #with missing values. how to replace with the previous value
ET_boot2<-ET_boot
for(p in 1:M){
#p<-1
#ET_boot[ET_boot$replica==p,]
for(q in 2:(j+1)){
#q<-2
#ET_boot[ET_boot$replica==p,q]
for (r in 1:length(len)){
#r<-14
ET_boot2[ET_boot2$replica==p,q][r]<-ifelse(is.na(ET_boot2[ET_boot2$replica==p,q][r]),
ET_boot2[ET_boot2$replica==p,q][(r-1)],
ET_boot2[ET_boot2$replica==p,q][r])
}
}
}
#ET_boot2[ET_boot2$replica==p,]
ET.qi2<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs2<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m2 in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n2 in 2:(j+1)){
#n<-2
#ET_boot2[ET_boot2$tComuns==len[m],4]
ET.qi2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = 0.025,na.rm = T))
ET.qs2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = .975,na.rm = T))
}
}
ET.qi2[,1]<-len
ET.qs2[,1]<-len
ET.qiAll<-rbind(ET.qiAll,cbind(as.data.frame(ET.qi), rep(s,length(len))))
ET.qi2All<-rbind(ET.qi2All,cbind(as.data.frame(ET.qi2), rep(s,length(len))))
ET.qsAll<-rbind(ET.qsAll,cbind(as.data.frame(ET.qs), rep(s,length(len))))
ET.qs2All<-rbind(ET.qs2All,cbind(as.data.frame(ET.qs2), rep(s,length(len))))
}#fim tempos
if(!is.null(db_wide)){
p.valueAll.f<-p.valueAll[,c(1, (from+1):ncol(p.valueAll))]
}else{
p.valueAll.f<-p.valueAll
}
round(p.valueAll.f,4)
#head(ET.qiAll)
names(ET.qiAll)[length(names(ET.qiAll))]<-'s'
names(ET.qiAll)[1]<-'tComuns'
names(ET.qiAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qi2All)
names(ET.qi2All)[length(names(ET.qi2All))]<-'s'
names(ET.qi2All)[1]<-'tComuns'
names(ET.qi2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qsAll)
names(ET.qsAll)[length(names(ET.qsAll))]<-'s'
names(ET.qsAll)[1]<-'tComuns'
names(ET.qsAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qs2All)
names(ET.qs2All)[length(names(ET.qs2All))]<-'s'
names(ET.qs2All)[1]<-'tComuns'
names(ET.qs2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
nomes<-rep("s",length(names(p.valueAll.f)))
names(p.valueAll.f)
for(i5 in 2:length(names(p.valueAll.f))){
#i5<-2
#nchar(names(p.valueAll.f)[i5])
nomes[i5]<-paste(from,'->',substr(names(p.valueAll.f)[i5], 7, nchar(names(p.valueAll.f)[i5])),sep='')
}
names(p.valueAll.f)<-nomes
#p.valueAll.f
res<-list(localTest=p.valueAll.f, trans = tmat, times=times, DIF=DIF, AJall=AJall, LMAJall=LMAJall, from=from,
to=to, ET.qiAll=ET.qiAll, ET.qi2All=ET.qi2All, ET.qsAll=ET.qsAll, ET.qs2All=ET.qs2All, replicas=M,
limit=limit)
class(res) <- c("localTest", "markovMSM")
res$call <- match.call()
#return(res)
options(warn=-1)
suppressWarnings(res)
return(invisible(res))
}
}
Try the markovMSM package in your browser
Any scripts or data that you put into this service are public.
markovMSM documentation built on Feb. 16, 2023, 8:34 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions AUC.test
Documented in AUC.test
#' Global test and local test for checking the Markov condition on Multi-state
#' Models based on the Area Under the two curves (AUC) given by the non-Markov
#' estimators of the transition probabilities to the Markov Aalen-Johansen
#' estimators.
#' @description This function is used to obtain the local and global tests for
#' checking the Markov condition.
#' @param data A data frame in long format containing the subject \code{id};
#' \code{from} corresponding to the starting state;the receiving state,
#' \code{to}; the transition number, \code{trans}; the starting time of the
#' transition given by \code{Tstart}; the stopping time of the transition,
#' \code{Tstop}, and \code{status} for the status variable, with 1 indicating
#' an event (transition), 0 a censoring.
#' @param from The starting state of the transition probabilities.
#' @param to The last receiving state considered for the estimation of the
#' transition probabilities. All the probabilities among the first and the
#' last states are also computed.
#' @param type Type of test for checking the Markov condition: "local" or
#' "global". By default type='global'.
#' @param times For the local test, times represents the starting times of the
#' transition probabilities. In case of a global test, the argument is given by
#' times between the minimum time and the third quartile times used
#' in the formula of this test. Default to NULL.
#' @param quantiles Quantiles used in the formula of the Global test for the
#' AUC methods. By default takes the percentiles 0.05, 0.10, 0.20, 0.30, 0.40.
#' @param tmat The transition matrix for multi-state model.
#' @param replicas Number of replicas for the Monte Carlo simulation to
#' standardization of the T-statistic given by the difference of the areas of
#' AJ and LM transition probabilities estimates.
#' @param limit Percentile of the event time used as the upper bound for the
#' computation of the AUC-based test.
#' @param positions List of possible transitions; x[[i]] consists of a vector of
#' state numbers reachable from state i.
#' @param namesStates A character vector containing the names of either the
#' competing risks or the states in the multi-state model specified by the
#' competing risks or illness-death model. names should have the same length as
#' the list x (for transMat), or either K or K+1 (for trans.comprisk), or 3
#' (for trans.illdeath).
#' @param timesNames Either 1) a matrix or data frame of dimension n x S
#' (n being the number of individuals and S the number of states in the
#' multi-state model), containing the times at which the states are visited or
#' last follow-up time, or 2) a character vector of length S containing the
#' column names indicating these times. In the latter cases, some elements of
#' time may be NA, see Details
#' @param status Either 1) a matrix or data frame of dimension n x S,
#' containing, for each of the states, event indicators taking the value 1 if
#' the state is visited or 0 if it is not (censored), or 2) a character
#' vector of length S containing the column names indicating these status
#' variables. In the latter cases, some elements of status may be NA,
#' see Details.
#'
#' @return In case of the AUC global test, an object with a list with the
#' following outcomes:
#' \item{globalTest}{p-value of AUC global tests for each transition. These
#' values are obtained through the minimum of the means of each two contiguous
#' quantiles times of the AUC global tests.}
#' \item{localTest}{AUC local tests of the transition probability for each times
#' and transitions.}
#' \item{quantiles}{Quantiles times used for the AUC global tests.}
#' \item{times}{Times used for the AUC global tests.}
#' \item{DIF}{Differences between the AJ and the LMAJ estimates for each
#' transition probabilites from the starting state until the receiving state
#' given by only one replica where 's' represent each of the quantile times.}
#' \item{from}{The starting state considered for the AUC global tests.}
#' \item{to}{The last receiving state considered for the the AUC Local tests.}
#' \item{ET.qiAll}{The lower limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qi2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{ET.qsAll}{The upper limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qs2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{replicas}{Number of replicas for the Monte Carlo simulation.}
#' \item{limit}{Percentil of the times used in the AUC global tests.}
#'
#' In case of the AUC local test, an object with a list with the following
#' outcomes:
#' \item{localTest}{p-value of AUC local tests for each times and transitions.}
#' \item{trans}{The transition matrix describing the states and transitions of
#' the multi-state model.}
#' \item{times}{Times selected for the AUC Local tests.}
#' \item{DIF}{Differences between the AJ and the LMAJ estimates for each
#' transition probabilites from the starting state until the receiving state
#' given by only one replica where 's' represent each of the quantile times.}
#' \item{from}{The starting state considered for the AUC Local tests.}
#' \item{to}{The last receiving state considered for the the AUC Local tests.}
#' \item{ET.qiAll}{The lower limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qi2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{ET.qsAll}{The upper limit of the diferences between the AJ and the LMAJ
#' estimates given by the Monte Carlo simulation in each transition for each "s"
#' quantile times. \code{ET.qs2All} means the same but missing values were
#' replaced by the previous diferences of estimators.}
#' \item{replicas}{Number of replicas for the Monte Carlo simulation.}
#' \item{limit}{Percentil of the times used in the AUC local tests.}
#' @references Soutinho G, Meira-Machado L (2021). Methods for checking the
#' Markov condition in multi-state survival data. \emph{Computational Statistics}.
#' de Una-alvarez J, Meira-Machado L (2015). Nonparametric estimation of
#' transition probabilities in the non-Markov illness-death model: A comparative
#' study. \emph{Biometrics}, 71(2), 364-375.
#' Putter H, Spitoni C (2018). Non-parametric estimation of transition
#' probabilities in non-Markov multi-state models: The landmark Aalen-Johansen
#' estimator. \emph{Statistical Methods in Medical Research}, 27, 2081-2092.
#' Meira-Machado L, Sestelo M (2019). Estimation in the progressive illness-death
#' model: A nonexhaustive review. \emph{Biometrical Journal}, 61(2), 245-263.
#' @examples
#' \donttest{
#' set.seed(1234)
#' library(markovMSM)
#' data("colonMSM")
#' db_wide<-colonMSM
#' positions<-list(c(2, 3), c(3), c())
#' namesStates = c("Alive", "Rec", "Death")
#' tmat <-transMatMSM(positions, namesStates)
#' timesNames = c(NA, "time1","Stime")
#' status=c(NA, "event1","event")
#' trans = tmat
#' db_long<- prepMSM(data=db_wide, trans, timesNames, status)
#' res<-AUC.test(data=db_long, times=180, from=1, to=3, type='local',
#' replicas = 20, tmat = tmat)
#' res$localTest
#' res2<-AUC.test(data=db_long, times=180, from=2, to=3, type='local',
#' replicas = 20, tmat = tmat)
#' res2$localTest
#'
#' res3<-AUC.test(data=db_long, from=1, to=3, replicas = 20, tmat=tmat)
#' round(res3$globalTest,3)
#'
#' res4<-AUC.test(data=db_long, from=2, to=3, type='global', replicas = 20,
#' tmat=tmat)
#' round(res4$globalTest,3)
#'
#' round(res4$localTest,3)
#' round(res4$localTest,3)
#'
#' #AUC global for the individuals with the treatment "Obs" of covariate `"rx"
#' #for colonMSM data set
#'
#' set.seed(12345)
#'
#' db_wide.obs<-db_wide[db_wide$rx=='Obs',]
#'
#' db_long.obs <- prepMSM(data = db_wide.obs, trans, timesNames, status)
#'
#' res3a<-AUC.test(data=db_long.obs, times=365, from=1, to=3,
#' type='local', replicas= 20, tmat = tmat)
#'
#' res3a$localTest
#'
#' set.seed(12345)
#'
#' res4a<-AUC.test(data=db_long.obs, times=365, from=2, to=3,
#' type='local', replicas= 20, tmat = tmat)
#'
#' res4a$localTest
#'
#' #' data("ebmt4")
#' db_wide <- ebmt4
#' positions <- list(c(2, 3, 5, 6), c(4, 5, 6), c(4, 5, 6),
#' c(5, 6), c(), c())
#' state.names <- c("Tx", "Rec", "AE", "Rec+AE", "Rel", "Death")
#' tmat <-transMatMSM(positions, state.names)
#' timesNames <- c(NA, "rec", "ae","recae", "rel", "srv")
#' status <- c(NA, "rec.s", "ae.s", "recae.s","rel.s", "srv.s")
#' trans <- tmat
#'
#' db_long <- prepMSM(data=db_wide, trans, timesNames, status)
#' db_long[1:10,]
#'
#' res5<-AUC.test(data=db_long, from=1, to=5, type='global',
#' quantiles=c(.05, .10, .20, .30, 0.40),
#' tmat = tmat, replicas = 5,
#' positions=positions, namesStates=state.names,
#' timesNames=timesNames, status=status)
#'
#' round(res5$globalTest, 4)
#' round(res5$localTests,4)
#'
#' res6<-AUC.test(data = prothr, from=2, to=3,
#' type='global', replicas= 5, limit=0.90,
#' quantiles=c(.05, .10, .20, .30, 0.40))
#' round(res6$globalTest,4)
#'
#' round(res6$localTests,4)
#'
#' }
#' @author Gustavo Soutinho and Luis Meira-Machado.
#'
#' @importFrom "survival" coxph Surv survfit strata untangle.specials
#' @importFrom "graphics" legend abline axis legend lines matplot par plot polygon
#' @importFrom "stats" pchisq pnorm quantile sd na.omit terms approxfun as.formula rnorm rpois weighted.mean
#' @importFrom "utils" capture.output flush.console
#' @importFrom "mstate" transMat msprep msfit probtrans LMAJ xsect cutLMms msdata2etm trans2Q events
#' @importFrom "stats" model.matrix model.frame model.response model.offset
#' @importFrom "stats" delete.response delete.response
#' @export AUC.test
#' @export prepMSM
#' @export transMatMSM
#' @export eventsMSM
AUC.test<- function(data, from=1, to=3, type='global',
times=NULL, quantiles=c(.05,.10, .20, .30, 0.40),
tmat=NULL, replicas=10, limit=0.90,
positions=list(c(2, 3), c(3), c()),
namesStates = c("Alive", "Rec", "Death"),
timesNames = c(NA, "time1","Stime"),
status=c(NA, "event1","event")){
db_long<-data
if (class(db_long)[1]!="msdata")
stop("Argument 'data' must be a mstate object")
db_wide<-NULL
if(type=='global'){
i<-from
#i<-2
j<-to
#j<-5
M<-replicas
DIF<-NULL
areaAll<-NULL
ETAll<-NULL
p.valueAll<-NULL
ET_bootAll<-NULL
ET.qiAll<-NULL
ET.qsAll<-NULL
ET.qi2All<-NULL
ET.qs2All<-NULL
times<-sort(times)
if(is.null(quantiles) & is.null(times)){
stop("One of the arguments quantiles or times must be NULL.")
}
if(!is.null(quantiles) & !is.null(times)){
stop("Only one of the arguments quantiles or times must be NULL.")
}
if(length(times)!=5 & !is.null(times)){
stop("The length of times must be 5.")
}
#min e quantile 0.75
if(!is.null(times)){
if(!is.null(db_wide)){
if(length(namesStates)==3){
timeToQuantiles<-timesNames[2]
min<-min(as.numeric(db_wide[, timeToQuantiles]))
val0.75<-as.numeric(quantile(db_wide[, timeToQuantiles], 0.75))
}else{#inicio difere de 3
if(i==1){#from
#head(db_wide)
valToQuantiles<-rep(0,nrow(db_wide))
for(a1 in 1:nrow(db_wide)){
#a1<-1
valToQuantiles[a1]<-min(db_wide[a1, timesNames[-1]])
}
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}else{
valToQuantiles<-db_wide[,timesNames[i]]
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}
min<-min(as.numeric(valToQuantiles))
val0.75<-as.numeric(quantile(valToQuantiles, 0.75))
}#fim de casos em que length(names) difere de 3
}else{#caso db_wide=NULL
tempQ<-unique(db_long$Tstop)
min<-min(as.numeric(tempQ))
val0.75<-as.numeric(quantile(tempQ, 0.75))
}
if(length(which(times<val0.75 & times>min & !is.null(times)))!=5){
stop("All the times must be greater than the minumum time and the third quantile of time.")
}
}
#
if(is.null(times))
times.to.test<-length(quantiles)
if(is.null(quantiles))
times.to.test<-length(times)
#start.time <- Sys.time() #marcacao do tempo inicial de processamento
cat("This method is based on the calculation of 5 percentile times and can lead to a high execution time.
Partial execution times will be provided during execution.", '\n')
for(h in 1:times.to.test){
#h<-1
#print(h) - estava este
#for(i3 in 1:length(colnames(db_wide))){
#i3<-1
# if (colnames(db_wide)[i3]==timeToQuantiles)
# indexS<-i3
#}
#s<-as.numeric(quantile(db_wide[,indexS], quantiles))[h]
if(!is.null(db_wide)){
if(length(namesStates)==3){
if(is.null(times)){
timeToQuantiles<-timesNames[2]
s<-as.numeric(quantile(db_wide[, timeToQuantiles], quantiles))[h]
s_quant<-as.numeric(quantile(db_wide[, timeToQuantiles], quantiles))
}else{
s<-as.numeric(times)[h]
s_quant<-as.numeric(times)
}
}else{ #inicio difere de 3
if(is.null(times)){
if(i==1){#from
#head(db_wide)
valToQuantiles<-rep(0,nrow(db_wide))
for(a1 in 1:nrow(db_wide)){
#a1<-1
valToQuantiles[a1]<-min(db_wide[a1, timesNames[-1]])
}
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}else{
valToQuantiles<-db_wide[,timesNames[i]]
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
}
s<-as.numeric(quantile(valToQuantiles, quantiles))[h]
s_quant<-as.numeric(quantile(valToQuantiles, quantiles))
}else{
s<-as.numeric(times)[h]
s_quant<-as.numeric(times)
}
}#fim de casos em que length(names) difere de 3
}else{#......................caso db_wide=NULL
if(is.null(times)){
#head(db_long)
#tempQ<-unique(db_long$Tstop)
#summary(tempQ)
if(length(namesStates)==3){
if(nrow(db_long)==2152){
#head(db_long)
tempQ<-unique(db_long$Tstop)
#summary(tempQ)
s<-as.numeric(quantile(tempQ, quantiles))[h]
tmat <- attr(db_long, "trans")
s_quant<-as.numeric(quantile(tempQ, quantiles))
}else{
values<-rep(0,length(unique(db_long$id)))
for(i3 in 1:length(unique(db_long$id))){
#i3<-1
values[i3]<-min(unique(db_long[db_long$id==i3,'Tstop']))
}
tempQ<-as.numeric(values)
#tempQ==as.numeric(db_wide[, timeToQuantiles])
#tempQ<-unique(db_long$Tstop)
min<-min(as.numeric(tempQ))
val0.75<-as.numeric(quantile(tempQ, 0.75))
s<-as.numeric(quantile(tempQ, quantiles))[h]
tmat <- attr(db_long, "trans")
s_quant<-as.numeric(quantile(tempQ, quantiles))
}
}else{#diferente de illness death model
#from=1
if(i==1){
values<-rep(0,length(unique(db_long$id)))
for(i3 in 1:length(unique(db_long$id))){
#i3<-2
#db_long[db_long$id==i3,]
#db_long[db_long$id==i3 ,'Tstop']
values[i3]<-min(unique(db_long[db_long$id==i3 ,'Tstop']))
}
valToQuantiles<-values
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
s<-as.numeric(quantile(valToQuantiles, quantiles))[h]
s_quant<-as.numeric(quantile(valToQuantiles, quantiles))
}else{#from !=1 ....
db2<-db_long[db_long$from==i,]
ids2<-unique(db2$id)
values2<-rep(0,length(ids2))
for(i4 in 1:length(ids2)){
#i4<-3
#db_long[db_long$id==ids2[i4],]
#db_long[db_long$id==ids2[i4] & db_long$from==i,'Tstop']
values2[i4]<-min(unique(db_long[db_long$id==ids2[i4] & db_long$from==i,'Tstop']))
}
valToQuantiles<-values2
valToQuantiles<-valToQuantiles[valToQuantiles!=0]
s<-as.numeric(quantile(valToQuantiles, quantiles))[h]
s_quant<-as.numeric(quantile(valToQuantiles, quantiles))
}
} #fim diferente de illness death model (3 estados)
}else{#valores em concreto quando nao se conhece WIDE
s<-as.numeric(times)[h]
s_quant<-as.numeric(times)
}
}
#h<-2
#s<-c(10, 30, 100, 200, 300, 650)[h] #10 erro
#cat("Time =", s,"\n")
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans),
data = db_long))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood",
trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
#resAJ[[i]]
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, limit)
tempos <- tempos[tempos < t.limit]
#trans<-as.numeric(na.omit(as.numeric(resAJ$trans)))
#resAJ[[i]][1:10,]
LMpt0 <- suppressWarnings(LMAJ2(db_long, s=s, from=i))
#LMpt0[200:250,]
#length(LMpt0$time)
#tempos[tempos %in% LMpt0$time]
#LMpt0$time[LMpt0$time %in% tempos]
#tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
#head(resAJf[,2:(j+1)])
#head(LMpt0f[,2:(j+1)])
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,] #; pr
#head(pr)
tempos2 <- diff(tComuns)
#length(tempos2)
#tempos2[1]*pr[1,]
#tempos2[2]*pr[2,]
#(pr * tempos2)[1:2,]
area<-colSums((pr * tempos2))
areaAll<-rbind(areaAll, cbind(s, as.data.frame(t(area))))
#areaAll
areasBoot<-matrix(NA, nrow = M, ncol = j) #ncol(resAJ$trans)
prTimes1<-cbind(tComuns, pr.ini)
DIF<-rbind(DIF, cbind(s, prTimes1))
#none<-NULL (para filtro...)
#for(n1 in from:to){
# #n1<-1
# if(n1==from)
# none<-c(none,unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]-1))
# else
# none<-c(none, unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]))
#}
ET_boot<-NULL
start.time <- Sys.time() #marcacao do tempo inicial de processamento
for(k in 1:M){
if(k==6){
ord1<-c('first', 'second', 'third', 'fourth', 'fifth', 'sixth', 'seventh', 'eighth', 'ninth', 'tenth',
'eleventh', 'twelfth')
ord2<-c('one', 'two', 'three', 'four', 'five', 'six', 'seven', 'eight', 'nine', 'ten',
'eleven', 'twelve')
#M<-100
#(as.numeric(difftime(Sys.time(), start.time))/5*M)/60
if(h==1 | h==2){
if(nrow(db_long)==15512){
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*1.5*M)/5)/60, '\r', sep='')
}else{
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*2.5*M)/5)/60, '\r', sep='')
}
}else{
if(h==3){
if(nrow(db_long)==15512){
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*M)/5)/60, '\r', sep='')
}else{
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*2.1*M)/5)/60, '\r', sep='')
}
}else{#nao ser h=3
if(nrow(db_long)==15512){
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*M)/5)/60, '\r', sep='')
}else{
cat("Expected computation time in minutes for ", ord1[h], ' (of ', ord2[times.to.test],') percentiles: ',
((as.numeric(difftime(Sys.time(), start.time))*1.3*M)/5)/60, '\r', sep='')
}
}#final de nao ser h=3
} #final de ser 3 ou nao ser
}
#k<-1
#cat("Monte Carlo sample =", k,"\n")
#head(db_wide)
#unique(db_wide$id)
if(!is.null(db_wide)){
n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
db_wide2 <- db_wide[pos,]
#head(db_wide2)
#positions<-list(c(2, 3), c(3), c())
#namesStates = names = c("Alive", "Rec", "Death")
tmat <-transMatMSM(positions, namesStates)
#times = c(NA, "time1","Stime")
#status=c(NA, "event1","event")
trans = tmat
db_long2<- prepMSM(data=db_wide2, trans, timesNames, status)
#db_long2<- msprep(data = db_wide2, trans = tmat,
# time = c(NA, "rec", "ae","recae", "rel", "srv"),
# status = c(NA, "rec.s", "ae.s", "recae.s","rel.s", "srv.s"),
# keep = c("match", "proph", "year", "agecl"))
}else{
n1<-length(unique(db_long$id))
#n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
#table(pos)
pos2<-sort(pos)
db_long2<-NULL
for(i4 in 1:n1){
db_long2 <- rbind(db_long2,db_long[db_long$id %in% pos2[i4],])
}
#db_long2[db_long2$id==3,]
}
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans),
data = db_long2))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood", trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, limit)
tempos <- tempos[tempos < t.limit]
#
msdata<-db_long2
xss <- suppressWarnings(xsect(msdata, s))
infrom <- xss$id[xss$state %in% from]
msdatas <- suppressWarnings(cutLMms(msdata, LM = s))
msdatasfrom <- msdatas[msdatas$id %in% infrom, ]
#print(nrow(msdatasfrom))
if(nrow(msdatasfrom)<2){
print('Error')
}else{
LMpt0 <- suppressWarnings(LMAJ2(db_long2, s=s, from=i))
#tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,]; pr
tempos2b <- diff(tComuns)
areasBoot[k,]<-colSums((pr * tempos2b))
#para o grafico
prTimes<-cbind(tComuns, pr.ini)
prTimes1$tComuns %in% prTimes$tComuns
dataframe<-as.data.frame(prTimes1$tComuns)
colnames(dataframe)<-'tComuns'
ET_boot_rep<-merge(x=dataframe, y=prTimes, by='tComuns',all.x = T)
ET_boot<-rbind(ET_boot,cbind(ET_boot_rep,rep(k,nrow(ET_boot_rep)), rep(s,nrow(ET_boot_rep))))
#ET_boot[200:350,]
#cat(k, '\n')
#cat(h, '\n')
}
#
}
sd_areasBoot<-rep(0, j) # ncol(resAJ$trans)
for(l in 1: j){#ncol(resAJ$trans)
#l<-1
sd_areasBoot[l]<-sd(areasBoot[,l], na.rm = T)
}
ET<-area/sd_areasBoot; ET
ETAll<-rbind(ETAll, cbind(s, as.data.frame(t(ET))))
#none<-unique(none)
#if(length(none)==1){
# p.value <- 2*(1-pnorm(abs(ET))); p.value
# for(n2 in from:to){
# p.value[n2]<-'None'
# }
#}else{
# p.value <- 2*(1-pnorm(abs(ET))); p.value
#}
p.value <- 2*(1-pnorm(abs(ET))); p.value
p.valueAll<-rbind(p.valueAll, cbind(s, as.data.frame(t(p.value ))))
p.valueAll
#
colnames(ET_boot)[(j+2)]<-'replica'
colnames(ET_boot)[(j+3)]<-'s'
#head(ET_boot)
len<-unique(ET_boot$tComuns)
ET.qi<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n in 2:(j+1)){
#n<-2
#ET_boot[ET_boot$tComuns==len[m],2]
ET.qi[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = 0.025,na.rm = T))
ET.qs[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = .975,na.rm = T))
}
}
ET.qi[,1]<-len
ET.qs[,1]<-len
ET_boot[1:20,]
ET_boot2<-ET_boot
for(p in 1:M){
#p<-1
#ET_boot[ET_boot$replica==p,]
for(q in 2:(j+1)){
#q<-2
#ET_boot[ET_boot$replica==p,q]
for (r in 1:length(len)){
#r<-14
ET_boot2[ET_boot2$replica==p,q][r]<-ifelse(is.na(ET_boot2[ET_boot2$replica==p,q][r]),
ET_boot2[ET_boot2$replica==p,q][(r-1)],
ET_boot2[ET_boot2$replica==p,q][r])
}
}
}
#ET_boot2[ET_boot2$replica==p,]
ET.qi2<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs2<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m2 in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n2 in 2:(j+1)){
#n<-2
#ET_boot2[ET_boot2$tComuns==len[m],4]
ET.qi2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = 0.025,na.rm = T))
ET.qs2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = .975,na.rm = T))
}
}
ET.qi2[,1]<-len
ET.qs2[,1]<-len
ET.qiAll<-rbind(ET.qiAll,cbind(as.data.frame(ET.qi), rep(s,length(len))))
ET.qi2All<-rbind(ET.qi2All,cbind(as.data.frame(ET.qi2), rep(s,length(len))))
ET.qsAll<-rbind(ET.qsAll,cbind(as.data.frame(ET.qs), rep(s,length(len))))
ET.qs2All<-rbind(ET.qs2All,cbind(as.data.frame(ET.qs2), rep(s,length(len))))
}#fim times
p.value.f<-p.valueAll
#round(p.value.f,4)
minAll<-rep(0, j)
if(is.null(times)){
minInt<-rep(0, (length(quantiles)-1))
}else{
minInt<-rep(0, (length(times)-1))
}
for (i2 in 1:j){
#i2<-1
for(j2 in 1: (nrow(p.valueAll)-1)){
#j2<-1
p.valueAll[,(i2+1)]<-ifelse(p.valueAll[,(i2+1)]== 'NaN', NA, p.valueAll[,(i2+1)])
minInt[j2]<-mean(p.valueAll[j2:(j2+1),(i2+1)],na.rm = T)
}
#is.na(minInt)
minAll[i2]<-min(minInt, na.rm=T)
}
minAll
#head(ET.qiAll)
names(ET.qiAll)[length(names(ET.qiAll))]<-'s'
names(ET.qiAll)[1]<-'tComuns'
names(ET.qiAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qi2All)
names(ET.qi2All)[length(names(ET.qi2All))]<-'s'
names(ET.qi2All)[1]<-'tComuns'
names(ET.qi2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qsAll)
names(ET.qsAll)[length(names(ET.qsAll))]<-'s'
names(ET.qsAll)[1]<-'tComuns'
names(ET.qsAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qs2All)
names(ET.qs2All)[length(names(ET.qs2All))]<-'s'
names(ET.qs2All)[1]<-'tComuns'
names(ET.qs2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
if(!is.null(db_wide)){
p.value.f2<-p.value.f[,c(1, (from+1):ncol(p.value.f))]
minAllf<-minAll[from:length(minAll)]
}else{
p.value.f2<-p.value.f
minAllf<-minAll
}
nomes<-rep("s",length(names(p.value.f2)))
names(p.value.f2)
for(i5 in 2:length(names(p.value.f2))){
#i5<-2
#nchar(names(p.valueAll.f)[i5])
nomes[i5]<-paste(from,'->',substr(names(p.value.f2)[i5], 7, nchar(names(p.value.f2)[i5])),sep='')
}
names(p.value.f2)<-nomes
#p.value.f2
minAllf2<-as.data.frame(t(minAllf))
names(minAllf2)<-nomes[2:length(nomes)]
if(is.null(times)){
res <-list(globalTest=minAllf2, localTests=p.value.f2,
quantiles=as.numeric(s_quant), times=NULL, DIF=DIF, from=from, to=to,
ET.qiAll=ET.qiAll, ET.qi2All=ET.qi2All, ET.qsAll=ET.qsAll, ET.qs2All=ET.qs2All,
replicas=M, limit=limit)
}else{
res <-list(globalTest=minAllf2, localTests=p.value.f2,
quantiles=NULL,times=times, DIF=DIF, from=from, to=to,
ET.qiAll=ET.qiAll, ET.qi2All=ET.qi2All, ET.qsAll=ET.qsAll, ET.qs2All=ET.qs2All,
replicas=M, limit=limit)
}
class(res) <- c("globalTest", "markovMSM")
res$call <- match.call()
#return(res)
options(warn=-1)
suppressWarnings(res)
return(invisible(res))
}else{
#local
i<-from
#i<-1
j<-to
#j<-6
M<-replicas
DIF<-NULL
areaAll<-NULL
ETAll<-NULL
p.valueAll<-NULL
AJall<-NULL
LMAJall<-NULL
ET.qiAll<-NULL
ET.qsAll<-NULL
ET.qi2All<-NULL
ET.qs2All<-NULL
if(is.null(db_wide))
tmat <- attr(db_long, "trans")
cat("This method may require high execution time.
Partial execution times will be provided during execution.", '\n')
for(h in 1:length(times)){
#h<-2
#print(h)
s<-times[h]
#cat("Time =", s,"\n")
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = db_long))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood", trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, 0.90)
tempos <- tempos[tempos < t.limit]
#trans<-as.numeric(na.omit(as.numeric(resAJ$trans)))
resAJ[[1]] #from 1 to 6
#resAJ[[i]][1:10,1:8]
LMpt0 <- suppressWarnings(LMAJ2(db_long, s=s, from=i))
#LMpt0[1:40,1:8]
length(LMpt0$time)
tempos[tempos %in% LMpt0$time]
LMpt0$time[LMpt0$time %in% tempos]
tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,]#; pr
#head(pr)
pr<-pr[-1,]
#head(pr)
tempos2 <- diff(tComuns)
#length(tempos2)
#tempos2[1]*pr[1,]
#tempos2[2]*pr[2,]
#pr[1:10,]
#tempos2[1:10]
#(pr * tempos2)[1:10,]
area<-colSums((pr * tempos2), na.rm = T)
areaAll<-rbind(areaAll, cbind(s, as.data.frame(t(area))))
#areaAll
areasBoot<-matrix(NA, nrow = M, ncol = j) #ncol(resAJ$trans)
prTimes1<-cbind(tComuns, pr.ini)
DIF<-rbind(DIF, cbind(s, prTimes1))
AJall<-rbind(AJall, cbind(s, resAJf[1:(j+1)]))
LMAJall<-rbind(LMAJall, cbind(s, LMpt0f[1:(j+1)]))
none<-NULL
for(n1 in from:to){
#n1<-1
if(n1==from)
none<-c(none,unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]-1))
else
none<-c(none, unique(prTimes1[,(1+n1)]==resAJf[,(1+n1)]))
}
ET_boot<-NULL
start.time <- Sys.time() #marcacao do tempo inicial de processamento
for(k in 1:M){
#k<-1
if(k==6){
cat("Expected computation time in minutes for ", s , ": " , ((as.numeric(difftime(Sys.time(), start.time))*2.3*M)/5)/60, '\r', sep='')
}
#cat("Monte Carlo sample =", k,"\n")
#head(db_wide)
#unique(db_wide$id)
if(!is.null(db_wide)){
n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
db_wide2 <- db_wide[pos,]
#head(db_wide2)
#positions<-list(c(2, 3), c(3), c())
#namesStates = names = c("Alive", "Rec", "Death")
tmat <-transMatMSM(positions, namesStates)
#times = c(NA, "time1","Stime")
#status=c(NA, "event1","event")
trans = tmat
db_long2<- prepMSM(data=db_wide2, trans, timesNames, status)
#db_long2<- msprep(data = db_wide2, trans = tmat,
# time = c(NA, "rec", "ae","recae", "rel", "srv"),
# status = c(NA, "rec.s", "ae.s", "recae.s","rel.s", "srv.s"),
# keep = c("match", "proph", "year", "agecl"))
}else{
n1<-length(unique(db_long$id))
#n1 <- dim(db_wide)[1]
pos <- sample.int(n1, n1, replace=TRUE)
#table(pos)
pos2<-sort(pos)
db_long2<-NULL
for(i4 in 1:n1){
db_long2 <- rbind(db_long2,db_long[db_long$id %in% pos2[i4],])
}
#db_long2[db_long2$id==3,]
}
c0 <- suppressWarnings(coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = db_long2))
msf0 <- suppressWarnings(msfit(object = c0, vartype = "greenwood", trans = tmat))
resAJ<- suppressWarnings(probtrans(msf0, predt=s))
tempos <- resAJ[[1]][,1]
t.limit <- quantile(tempos, 0.90)
tempos <- tempos[tempos < t.limit]
LMpt0 <- suppressWarnings(LMAJ2(db_long2, s=s, from=i))
#tempos[tempos %in% LMpt0$time] == LMpt0$time[LMpt0$time %in% tempos]
tComuns<-tempos[tempos %in% LMpt0$time]
resAJf<-resAJ[[i]][resAJ[[i]][,1] %in% tComuns,]
#resAJf<-resAJ[[1]][resAJ[[1]][,1] %in% tComuns,]
LMpt0f<-LMpt0[LMpt0$time %in% tComuns,]
#o que estava
#pr<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
#pr<-pr[-1,]; pr
#tempos2 <- diff(tComuns)
#areasBoot[k,]<-colSums((pr * tempos2), na.rm = T) #,
pr.ini<-resAJf[,2:(j+1)]-LMpt0f[,2:(j+1)]
pr<-pr.ini[-1,]; pr
tempos2b <- diff(tComuns)
areasBoot[k,]<-colSums((pr * tempos2b))
#para o grafico
prTimes<-cbind(tComuns, pr.ini)
prTimes1$tComuns %in% prTimes$tComuns
dataframe<-as.data.frame(prTimes1$tComuns)
colnames(dataframe)<-'tComuns'
ET_boot_rep<-merge(x=dataframe, y=prTimes, by='tComuns',all.x = T)
ET_boot<-rbind(ET_boot,cbind(ET_boot_rep,rep(k,nrow(ET_boot_rep)), rep(s,nrow(ET_boot_rep))))
}
sd_areasBoot<-rep(0, j) # ncol(resAJ$trans)
for(l in 1: j){#ncol(resAJ$trans)
#l<-3
sd_areasBoot[l]<-sd(areasBoot[,l], na.rm = T)
}
ET<-area/sd_areasBoot; ET
ETAll<-rbind(ETAll, cbind(s, as.data.frame(t(ET))))
#ETAll
#none<-unique(none)
#if(length(none)==1){
# p.value <- 2*(1-pnorm(abs(ET))); p.value
# for(n2 in from:to){
# p.value[n2]<-'None'
# }
#}else{
# p.value <- 2*(1-pnorm(abs(ET))); p.value
#}
p.value <- 2*(1-pnorm(abs(ET))); p.value
p.valueAll<-rbind(p.valueAll, cbind(s, as.data.frame(t(p.value ))))
p.valueAll
#
colnames(ET_boot)[(j+2)]<-'replica'
colnames(ET_boot)[(j+3)]<-'s'
#head(ET_boot)
len<-unique(ET_boot$tComuns)
ET.qi<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n in 2:(j+1)){
#n<-2
#ET_boot[ET_boot$tComuns==len[m],2]
ET.qi[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = 0.025,na.rm = T))
ET.qs[m,n]<-as.numeric(quantile(ET_boot[ET_boot$tComuns==len[m],n],probs = .975,na.rm = T))
}
}
ET.qi[,1]<-len
ET.qs[,1]<-len
ET_boot[1:20,] #with missing values. how to replace with the previous value
ET_boot2<-ET_boot
for(p in 1:M){
#p<-1
#ET_boot[ET_boot$replica==p,]
for(q in 2:(j+1)){
#q<-2
#ET_boot[ET_boot$replica==p,q]
for (r in 1:length(len)){
#r<-14
ET_boot2[ET_boot2$replica==p,q][r]<-ifelse(is.na(ET_boot2[ET_boot2$replica==p,q][r]),
ET_boot2[ET_boot2$replica==p,q][(r-1)],
ET_boot2[ET_boot2$replica==p,q][r])
}
}
}
#ET_boot2[ET_boot2$replica==p,]
ET.qi2<-matrix(NA, nrow = length(len), ncol = (j+1))
ET.qs2<-matrix(NA, nrow = length(len), ncol = (j+1))
for(m2 in 1:length(len)){
#m<-100
#ET_boot[ET_boot$tComuns==len[m],]
for(n2 in 2:(j+1)){
#n<-2
#ET_boot2[ET_boot2$tComuns==len[m],4]
ET.qi2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = 0.025,na.rm = T))
ET.qs2[m2,n2]<-as.numeric(quantile(ET_boot2[ET_boot2$tComuns==len[m2],n2],probs = .975,na.rm = T))
}
}
ET.qi2[,1]<-len
ET.qs2[,1]<-len
ET.qiAll<-rbind(ET.qiAll,cbind(as.data.frame(ET.qi), rep(s,length(len))))
ET.qi2All<-rbind(ET.qi2All,cbind(as.data.frame(ET.qi2), rep(s,length(len))))
ET.qsAll<-rbind(ET.qsAll,cbind(as.data.frame(ET.qs), rep(s,length(len))))
ET.qs2All<-rbind(ET.qs2All,cbind(as.data.frame(ET.qs2), rep(s,length(len))))
}#fim tempos
if(!is.null(db_wide)){
p.valueAll.f<-p.valueAll[,c(1, (from+1):ncol(p.valueAll))]
}else{
p.valueAll.f<-p.valueAll
}
round(p.valueAll.f,4)
#head(ET.qiAll)
names(ET.qiAll)[length(names(ET.qiAll))]<-'s'
names(ET.qiAll)[1]<-'tComuns'
names(ET.qiAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qi2All)
names(ET.qi2All)[length(names(ET.qi2All))]<-'s'
names(ET.qi2All)[1]<-'tComuns'
names(ET.qi2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qsAll)
names(ET.qsAll)[length(names(ET.qsAll))]<-'s'
names(ET.qsAll)[1]<-'tComuns'
names(ET.qsAll)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
#
#head(ET.qs2All)
names(ET.qs2All)[length(names(ET.qs2All))]<-'s'
names(ET.qs2All)[1]<-'tComuns'
names(ET.qs2All)[2:(length(3:ncol(DIF))+1)]<-names(DIF)[3:ncol(DIF)]
nomes<-rep("s",length(names(p.valueAll.f)))
names(p.valueAll.f)
for(i5 in 2:length(names(p.valueAll.f))){
#i5<-2
#nchar(names(p.valueAll.f)[i5])
nomes[i5]<-paste(from,'->',substr(names(p.valueAll.f)[i5], 7, nchar(names(p.valueAll.f)[i5])),sep='')
}
names(p.valueAll.f)<-nomes
#p.valueAll.f
res<-list(localTest=p.valueAll.f, trans = tmat, times=times, DIF=DIF, AJall=AJall, LMAJall=LMAJall, from=from,
to=to, ET.qiAll=ET.qiAll, ET.qi2All=ET.qi2All, ET.qsAll=ET.qsAll, ET.qs2All=ET.qs2All, replicas=M,
limit=limit)
class(res) <- c("localTest", "markovMSM")
res$call <- match.call()
#return(res)
options(warn=-1)
suppressWarnings(res)
return(invisible(res))
}
}
Try the markovMSM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.