drizopoulos / JMbayes

Joint Models for Longitudinal and Survival Data using MCMC
58 stars 23 forks source link

Bugs in the function aucJM #91

Closed ChenXinaha closed 2 years ago

ChenXinaha commented 2 years ago

Hello, I think there may be some problems in the function aucJM, in the follow code,because " if (any(dupl <- duplicated(Time))) { Time[dupl] <- Time[dupl] + runif(length(Time[dupl]), 1e-07, 1e-06) }" only keep the unique event time, which will result in removing subjects whose raw event time less than "Thoriz" using "Ti <= Thoriz" and errors in the estimation of auc. Further, I don't understand the effect of "Time[dupl] <- Time[dupl] + runif(length(Time[dupl]), 1e-07, 1e-06)".

I hope you can answer that for me, thanks very much !

names(Time) <- names(event) <- as.character(unique(id))
if (any(dupl <- duplicated(Time))) {
    Time[dupl] <- Time[dupl] + runif(length(Time[dupl]), 1e-07, 1e-06)
}
if (!all(names(pi.u.t) == names(Time)))
    stop("mismatch between 'Time' variable names and survival probabilities names.")
auc <- if (length(Time) > 1) {
    pairs <- combn(as.character(unique(id)), 2)
    Ti <- Time[pairs[1, ]]
    Tj <- Time[pairs[2, ]]
    di <- event[pairs[1, ]]
    dj <- event[pairs[2, ]]
    pi.u.t.i <- pi.u.t[pairs[1, ]]
    pi.u.t.j <- pi.u.t[pairs[2, ]]
    ind1 <- (Ti <= Thoriz & (di == 1 | di == 3)) & Tj > Thoriz
    ind2 <- (Ti <= Thoriz & (di == 0 | di == 2)) & Tj > Thoriz
    ind3 <- (Ti <= Thoriz & (di == 1 | di == 3)) & (Tj <= Thoriz & (dj == 0 | dj == 2))
    ind4 <- (Ti <= Thoriz & (di == 0 | di == 2)) & (Tj <= Thoriz & (dj == 0 | dj == 2))
    names(ind1) <- names(ind2) <- names(ind3) <- names(ind4) <- paste(names(Ti), names(Tj), sep = "_")
    ind <- ind1 | ind2 | ind3 | ind4
    if (any(ind2)) {
        nams <- strsplit(names(ind2[ind2]), "_")
        nams_i <- sapply(nams, "[", 1)
        unq_nams_i <- unique(nams_i)
        pi2 <- if (is_counting) {
            survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ], idVar = idVar, 
                         last.time = Time[unq_nams_i], survTimes = Thoriz, 
                         M = M, LeftTrunc_var = all.vars(TermsT)[1L])
        } else {
            survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ], idVar = idVar, 
                      last.time = Time[unq_nams_i], survTimes = Thoriz, M = M)
        }
        pi2 <- 1 - sapply(pi2$summaries, "[", 1, 2)
        ind[ind2] <- ind[ind2] * pi2[nams_i]
    }
    if (any(ind3)) {
        nams <- strsplit(names(ind3[ind3]), "_")
        nams_j <- sapply(nams, "[", 2)
        unq_nams_j <- unique(nams_j)
        pi3 <- if (is_counting) {
            survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ], idVar = idVar, 
                      last.time = Time[unq_nams_j], survTimes = Thoriz, 
                      M = M, LeftTrunc_var = all.vars(TermsT)[1L])
        } else {
            survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ], idVar = idVar, 
                      last.time = Time[unq_nams_j], survTimes = Thoriz, M = M)
        }
        pi3 <- sapply(pi3$summaries, "[", 1, 2)
        ind[ind3] <- ind[ind3] * pi3[nams_j]
    }
    if (any(ind4)) {
        nams <- strsplit(names(ind4[ind4]), "_")
        nams_i <- sapply(nams, "[", 1)
        nams_j <- sapply(nams, "[", 2)
        unq_nams_i <- unique(nams_i)
        unq_nams_j <- unique(nams_j)
        if (is_counting) {
            pi4_i <- survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ], 
                               idVar = idVar, last.time = Time[unq_nams_i], survTimes = Thoriz, 
                               M = M, 
                               LeftTrunc_var = all.vars(TermsT)[1L])
            pi4_j <- survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ], 
                               idVar = idVar, last.time = Time[unq_nams_j], survTimes = Thoriz, 
                               M = M, 
                               LeftTrunc_var = all.vars(TermsT)[1L])
        } else {
            pi4_i <- survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ], 
                               idVar = idVar, last.time = Time[unq_nams_i], survTimes = Thoriz, 
                               M = M)
            pi4_j <- survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ], 
                               idVar = idVar, last.time = Time[unq_nams_j], survTimes = Thoriz, 
                               M = M)
        }
        pi4_i <- 1 - sapply(pi4_i$summaries, "[", 1, 2)
        pi4_j <- sapply(pi4_j$summaries, "[", 1, 2)
        ind[ind4] <- ind[ind4] * pi4_i[nams_i] * pi4_j[nams_j]
    }
    sum((pi.u.t.i < pi.u.t.j) * c(ind), na.rm = TRUE) / sum(ind, na.rm = TRUE)
} else {
    NA
}
drizopoulos commented 2 years ago

Check JMbayes2

—— Dimitris Rizopoulos Professor of Biostatistics Erasmus University Medical Center The Netherlands


From: ChenXinaha @.> Sent: Sunday, April 24, 2022 11:11:15 AM To: drizopoulos/JMbayes @.> Cc: Subscribed @.***> Subject: [drizopoulos/JMbayes] Bugs in the function aucJM (Issue #91)

Waarschuwing: Deze e-mail is afkomstig van buiten de organisatie. Klik niet op links en open geen bijlagen, tenzij u de afzender herkent en weet dat de inhoud veilig is. Caution: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.

Hello, I think there may be some problems in the function aucJM, in the follow code,because " if (any(dupl <- duplicated(Time))) { Time[dupl] <- Time[dupl] + runif(length(Time[dupl]), 1e-07, 1e-06) }" only keep the unique event time, which will result in removing subjects whose raw event time less than "Thoriz" using "Ti <= Thoriz" and errors in the estimation of auc. Further, I don't understand the effect of "Time[dupl] <- Time[dupl] + runif(length(Time[dupl]), 1e-07, 1e-06)".

I hope you can answer that for me, thanks very much !

names(Time) <- names(event) <- as.character(unique(id))

if (any(dupl <- duplicated(Time))) {

Time[dupl] <- Time[dupl] + runif(length(Time[dupl]), 1e-07, 1e-06)

}

if (!all(names(pi.u.t) == names(Time)))

stop("mismatch between 'Time' variable names and survival probabilities names.")

auc <- if (length(Time) > 1) {

pairs <- combn(as.character(unique(id)), 2)

Ti <- Time[pairs[1, ]]

Tj <- Time[pairs[2, ]]

di <- event[pairs[1, ]]

dj <- event[pairs[2, ]]

pi.u.t.i <- pi.u.t[pairs[1, ]]

pi.u.t.j <- pi.u.t[pairs[2, ]]

ind1 <- (Ti <= Thoriz & (di == 1 | di == 3)) & Tj > Thoriz

ind2 <- (Ti <= Thoriz & (di == 0 | di == 2)) & Tj > Thoriz

ind3 <- (Ti <= Thoriz & (di == 1 | di == 3)) & (Tj <= Thoriz & (dj == 0 | dj == 2))

ind4 <- (Ti <= Thoriz & (di == 0 | di == 2)) & (Tj <= Thoriz & (dj == 0 | dj == 2))

names(ind1) <- names(ind2) <- names(ind3) <- names(ind4) <- paste(names(Ti), names(Tj), sep = "_")

ind <- ind1 | ind2 | ind3 | ind4

if (any(ind2)) {

    nams <- strsplit(names(ind2[ind2]), "_")

    nams_i <- sapply(nams, "[", 1)

    unq_nams_i <- unique(nams_i)

    pi2 <- if (is_counting) {

        survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ], idVar = idVar,

                     last.time = Time[unq_nams_i], survTimes = Thoriz,

                     M = M, LeftTrunc_var = all.vars(TermsT)[1L])

    } else {

        survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ], idVar = idVar,

                  last.time = Time[unq_nams_i], survTimes = Thoriz, M = M)

    }

    pi2 <- 1 - sapply(pi2$summaries, "[", 1, 2)

    ind[ind2] <- ind[ind2] * pi2[nams_i]

}

if (any(ind3)) {

    nams <- strsplit(names(ind3[ind3]), "_")

    nams_j <- sapply(nams, "[", 2)

    unq_nams_j <- unique(nams_j)

    pi3 <- if (is_counting) {

        survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ], idVar = idVar,

                  last.time = Time[unq_nams_j], survTimes = Thoriz,

                  M = M, LeftTrunc_var = all.vars(TermsT)[1L])

    } else {

        survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ], idVar = idVar,

                  last.time = Time[unq_nams_j], survTimes = Thoriz, M = M)

    }

    pi3 <- sapply(pi3$summaries, "[", 1, 2)

    ind[ind3] <- ind[ind3] * pi3[nams_j]

}

if (any(ind4)) {

    nams <- strsplit(names(ind4[ind4]), "_")

    nams_i <- sapply(nams, "[", 1)

    nams_j <- sapply(nams, "[", 2)

    unq_nams_i <- unique(nams_i)

    unq_nams_j <- unique(nams_j)

    if (is_counting) {

        pi4_i <- survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ],

                           idVar = idVar, last.time = Time[unq_nams_i], survTimes = Thoriz,

                           M = M,

                           LeftTrunc_var = all.vars(TermsT)[1L])

        pi4_j <- survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ],

                           idVar = idVar, last.time = Time[unq_nams_j], survTimes = Thoriz,

                           M = M,

                           LeftTrunc_var = all.vars(TermsT)[1L])

    } else {

        pi4_i <- survfitJM(object, newdata = newdata2[id %in% unq_nams_i, ],

                           idVar = idVar, last.time = Time[unq_nams_i], survTimes = Thoriz,

                           M = M)

        pi4_j <- survfitJM(object, newdata = newdata2[id %in% unq_nams_j, ],

                           idVar = idVar, last.time = Time[unq_nams_j], survTimes = Thoriz,

                           M = M)

    }

    pi4_i <- 1 - sapply(pi4_i$summaries, "[", 1, 2)

    pi4_j <- sapply(pi4_j$summaries, "[", 1, 2)

    ind[ind4] <- ind[ind4] * pi4_i[nams_i] * pi4_j[nams_j]

}

sum((pi.u.t.i < pi.u.t.j) * c(ind), na.rm = TRUE) / sum(ind, na.rm = TRUE)

} else {

NA

}

— Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdrizopoulos%2FJMbayes%2Fissues%2F91&data=05%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cadbdcdaea3a5404eb39208da25d26356%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637863882798256059%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=A2WzC3ZgwEybdy%2BBcl3mRsI4gogNVerp1EG%2FvGiCQmY%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FADE7TTYSANNHSTI2DOLK7GTVGUF3HANCNFSM5UF62AIA&data=05%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cadbdcdaea3a5404eb39208da25d26356%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637863882798256059%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=iNU2jw9MQtrzxYchhlVFQOZQ%2FcM7ZIgExthzUMDhEtA%3D&reserved=0. You are receiving this because you are subscribed to this thread.Message ID: @.***>