来自不同患者就诊序列的马尔可夫转移矩阵

Question

我正在尝试根据不同患者的就诊顺序创建马尔可夫转换矩阵。在我的马尔可夫模型中，状态是不同的医生，而联系是患者的访问。患者可以留在同一个提供者处，也可以过渡到另一个提供者进行下一次访问。使用该信息，我需要创建一个转换矩阵。

这是excel中的一部分数据。数据包括对近 100 家不同提供商的 3 万多次访问。

这是excel中的部分数据。 data

如何使用这个 excel 数据（或 csv）并创建一个马尔可夫转移矩阵作为访问次数，例如： ....

我需要的矩阵如下所示：

enter image description here

如何使用 R 将我的数据转换为转换矩阵？

我是 R 的新手，真的需要帮助。

谢谢

Answer 1

这是一种适用于您的示例数据的方法。

我将使用readxl 获取数据并使用data.table 来操作它。

读取数据：

library(readxl)
library(data.table)

data <- setDT(read_excel("~/Desktop/Book2.xlsx"))[!is.na(PatId)]

#read_excel doesn't have the option to specify integers... silly...
data[ , (names(data)) := lapply(.SD, as.integer)]

预分配转移矩阵：

provs <- data[ , sort(unique(SeenByProv))]
nprov <- length(provs)

markov <- matrix(nrow = nprov, ncol = nprov,
                 dimnames = list(provs, provs))

逐行赋值

for (pr in provs){
  markov[as.character(pr), ] <-
    data[ , {nxt <- SeenByProv[which(SeenByProv == pr) + 1L]
    .(prov = provs, count = 
        sapply(provs, function(pr2) sum(nxt == pr2, na.rm = TRUE)))}, by = PatId
    ][, sum(count), by = prov]$V1
}

这可能会在一些地方加速，但它确实有效。

Answer 2

我想在不使用 data.table 的情况下比较我的方法，发现它快了 45 倍（而且可能更容易理解）。

首先，我根据接受的答案对 data.table 解决方案进行计时：

rm(list=ls())
library(readxl)
library(data.table)

############## Using data.table method() ######################
data <- setDT(read_excel("Book2.xlsx"))[!is.na(PatId)]
data[ , (names(data)) := lapply(.SD, as.integer)]
provs <- data[ , sort(unique(SeenByProv))]
nprov <- length(provs)
markov <- matrix(nrow = nprov, ncol = nprov, dimnames = list(provs, provs))

system.time(      ## Timing the main loop
  for (pr in provs){
    markov[as.character(pr), ] <-
      data[ , {nxt <- SeenByProv[which(SeenByProv == pr) + 1L]
      .(prov = provs, count =
          sapply(provs, function(pr2) sum(nxt == pr2, na.rm = TRUE)))}, by = PatId
      ][, sum(count), by = prov]$V1
  }
)
#   user  system elapsed 
#  3.128   0.000   3.135 
table(markov)
#markov
#   0    1    2    3    4    5    6    7    8    9   10   11   13   22  140 
#3003  308   89   34   14   11    6    4    1    3    4    1    1    1    1 

接下来只使用基本 R 调用：

############## Using all base R calls method() ###################
tm_matrix<-matrix(0, nrow = nprov, ncol = nprov, dimnames = list(provs, provs))
d<-read_excel("Book2.xlsx")
d<-d[!is.na(d$PatId),] # Note: Data is already ordered by PatId, DaysOfStudy

baseR<-function(tm_matrix){
  d1<-cbind(d[-nrow(d),-3],d[-1,-3]); # Form the transitions and drop the DaysofStudy
  colnames(d1)<-c("SeenByProv","PatId","NextProv","PatId2");
  d1<-d1[d1$PatId==d1$PatId2,];       # Drop those transition between different patients
  d1$SeenByProv<-as.character(d1$SeenByProv); # transform to strings to use as rownames
  d1$NextProv  <-as.character(d1$NextProv);   # and column names
  for (i in 1:nrow(d1)){                      # Fill in the transition matrix
    tm_matrix[d1$SeenByProv[i],d1$NextProv[i]]<-tm_matrix[d1$SeenByProv[i],d1$NextProv[i]]+1
  };
  return(tm_matrix)
}
system.time(tm_matrix<-baseR(tm_matrix))
#   user  system elapsed 
#  0.072   0.000   0.072 

table(tm_matrix)
#tm_matrix
#   0    1    2    3    4    5    6    7    8    9   10   11   13   22  140 
#3003  308   89   34   14   11    6    4    1    3    4    1    1    1    1 

all.equal(markov,tm_matrix)
#[1] TRUE

我的 base-R 方法快 3.135/0.072 = 43.54