R：如何制作自己的 msdata 对象以用于 mstate 包

问：

我正在尝试制作一个生存模型，但我可能会有周期。

mstate 包说我不能将 msprep 函数用于循环数据，因为该函数是递归的。msprep 函数返回一个 msdata 对象。

因此，我正在编写自己的函数来从周期性数据构建 msdata 对象。下面列出了该功能。

我可以让对象看起来就像 mstate 引用中的 msdata 对象一样，但是当我运行 msfit 函数时，出现错误：

> msf0 <- msfit(object = c0, vartype = "greenwood", trans = tmat)
Warning messages:
1: In min(diff(time)) : no non-missing arguments to min; returning Inf
2: In min(diff(time)) : no non-missing arguments to min; returning Inf
3: In min(diff(time)) : no non-missing arguments to min; returning Inf

绘制结果看起来像是无稽之谈。

我在这里扫描了 CRAN 中的源代码：https://github.com/cran/mstate/blob/a6567548c1a79bd381e02a474168b8e9cb9de27f/R/msfit.R

我能找到的唯一关于“min”的参考资料与一些数据副本有关......但我不确定他们在做什么。我怀疑问题是否归结为我的集合中缺少数据，但 df 看起来确实像 mstate 样本。

library(mstate)
ids <- c(1551, 1555, 1551, 1555, 1560, 1551, 1555, 1560, 1561, 1562, 1551, 1555, 1560, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1561, 1562, 1551, 1555, 1562, 1551, 1555, 1562, 1551, 1555, 1562, 1551, 1555, 1562, 1551, 1555, 1562)
fails <- c(0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
succeeds <- c(1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
status <- c("a", "a", "a", "a", "a", "a", "a", "a", "b", "a", "a", "a", "c", "b", "a", "a", "e", "b", "a", "a", "c", "b", "a", "a", "c", "d", "a", "a", "c", "d", "e", "a", "c", "d", "c", "a", "c", "f", "c", "a", "b", "f", "c", "a", "b", "f", "b", "a", "b", "f", "b", "c", "g", "b", "g", "e", "b", "g", "f", "b", "g", "f", "b", "g", "f", "g")
df <- data.frame(ids, fails, succeeds, status)


splits <- split(df,df$ids)

# Build the transMat object
states <- as.character(unique(df$status))
states <- c("#EntersSystem#",states,"#Succeeds#","#Fails#")
transitioncounts <- matrix(0,nrow=length(states),ncol=length(states))
rownames(transitioncounts) <- states
colnames(transitioncounts) <- states
for(s in splits){
    s$status <- as.character(s$status)
    for(i in 1:nrow(s)){
        if(i == 1){
            # first observation moves from #EntersSystem#
            previousState <- "#EntersSystem#"
        } else {
            previousState <- s[i-1,]$status
        }
        nextState <- s[i,]$status
        transitioncounts[previousState,nextState] <- transitioncounts[previousState,nextState] + 1
    }

    # last observation moves to #Succeeds# or #Fails#
    if(s[1,]$succeeds == 1){
        nextState <- "#Succeeds#"
    } else {
        nextState <- "#Fails#"
    }
    previousState <- s[nrow(s),]$status
    transitioncounts[previousState,nextState] <- transitioncounts[previousState,nextState] + 1
}



# Build the transition list
transitionlist <- list()
for(i in 1:nrow(transitioncounts)){
    v <- as.integer(vector())
    for(j in 1:ncol(transitioncounts)){
        if((i != j)&(transitioncounts[i,j]>0)){
            v[length(v)+1] <- j
        }
    }
    transitionlist[[length(transitionlist)+1]] <- v
}

tmat <- transMat(x = transitionlist, names = states)


# Build the multistate data frame object
msdf <- data.frame()
for(s in splits){
    #s <- plits[[1]]
    s$status <- as.character(s$status)

    for(i in 1:nrow(s)){
        #i <- 2
        if(i == 1){
            # first observation moves from #EntersSystem#
            fromState <- "#EntersSystem#"
            Tstart <- 0
        } else {
            fromState <- s[i-1,]$status
        }
        nextState <- which(states == s[i,]$status)

        if(fromState != rownames(tmat)[nextState]){
            # there is a transition, so add a transition row
            to <- which(!is.na(tmat[fromState,]))
            status <- ifelse(to==nextState,1,0)
            from <- which(rownames(tmat)==fromState)
            id <- s[i,]$ids
            Tstop <- i
            Time <- Tstop - Tstart
            trans <- tmat[from,to]

            msdf <- msdf <- rbind(msdf,data.frame(id = as.numeric(id),
                                                  from = as.numeric(from),
                                                  to = as.numeric(to),
                                                  trans = as.numeric(trans),
                                                  Tstart = as.numeric(Tstart),
                                                  Tstop = as.numeric(Tstop),
                                                  time = as.numeric(Time),
                                                  status = as.numeric(status),
                                                  row.names = NULL))

            Tstart <- Tstop # prime the loop
        }
    }

    # last observation moves to #Succeeds# or #Fails#
    if(s[1,]$succeeds == 1){
        nextState <- which(states == "#Succeeds#")
    } else {
        nextState <- which(states == "#Fails#")
    }
    fromState <- s[nrow(s),]$status
    to <- which(!is.na(tmat[fromState,]))
    status <- ifelse(to==nextState,1,0)
    from <- which(rownames(tmat)==fromState)
    id <- s[1,]$ids
    Tstop <- nrow(s)+1
    Time <- Tstop - Tstart
    trans <- tmat[from,to]

    msdf <- rbind(msdf,data.frame(id = as.numeric(id),
                                  from = as.numeric(from),
                                  to = as.numeric(to),
                                  trans = as.numeric(trans),
                                  Tstart = as.numeric(Tstart),
                                  Tstop = as.numeric(Tstop),
                                  time = as.numeric(Time),
                                  status = as.numeric(status),
                                  row.names = NULL))
}

head(msdf)

attr(msdf, "trans") <- tmat
class(msdf) <- c("msdata","data.frame")

events(msdf) # looks good

## Cox Model with No Covariates
c0 <- coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = msdf, method = "breslow")
msf0 <- msfit(object = c0, vartype = "greenwood", trans = tmat) # borks


#### Code for comparison...

library(mstate)
data("ebmt4")
ebmt <- ebmt4

tmat2 <- transMat(x = list(c(2, 3, 5, 6), c(4, 5, 6), c(4, 5, 6), c(5, 6), c(), c()), names = c("Tx", "Rec", "AE", "Rec+AE", "Rel", "Death"))

msebmt <- msprep(data = ebmt, trans = tmat2, 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"))

events(msebmt)

c0 <- coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = msebmt, method = "breslow")
msff0 <- msfit(object = c0, vartype = "greenwood", trans = tmat2) # doesn't bork