※ 引述《yummy7922 (crucify)》之铭言：
: ※ 引述《yummy7922 (crucify)》之铭言：
: : [问题类型]:
: : 程式咨询(我想用R 做某件事情，但是我不知道要怎么用R 写出来)
: : [软件熟悉度]:
: : 入门(写过其他程式，只是对语法不熟悉)
: : [问题叙述]:
: : 我的资料是一个重复测量的资料，资料中有当期是否使用药物的资料(若有使用，设
为1
: : 没有为0)、还有记录该笔资料为该位病人的第几笔观察值，我希望能算出，药物转换
的
: : 率，例如有多少人第一期使用a药物，但在第二期时转换成使用b药物。
: : 资料中共有6种药物，一位病人最多有48笔观察值，
: : (6x6)种转换可能 x 47(个时间隔) = 1692 个机率值。
: : 想请教各位高手们，该怎么做比较有效率。
: 不好意思，我没有说明清楚，
: 我想算的机率其实是很直观的条件机率，例如：
: P(第二期使用b药物|第一期使用a药物)
: = (第一期使用a药且第二期使用b药的人数)/(第一期使用a药的人数)
: 但是我不知道该怎么在分组之后，还能够给定条件，
: 计算出第一期使用a药第二期使用b药的人数。
library(data.table)
library(dplyr)
library(magrittr)
library(Matrix)
# data generation
N_patient = 1000
dat = sample(48, N_patient, replace = TRUE) %>% {
cbind(rep(1:N_patient, times=.), sapply(., seq, from = 1) %>% unlist())
} %>% cbind(sample(6, nrow(.), replace = TRUE)) %>%
data.table() %>% setnames(c("id", "obs_time", "dose"))
dat_transform = dat$dose %>%
spMatrix(length(.), 6, 1:length(.), ., rep(1,length(.))) %>%
as.matrix() %>% data.table() %>% setnames(paste0("dose_", 1:6)) %>%
cbind(select(dat, 1:2), .)
# dat_transform为你的资料格式，透过下面可以转成dat
dat = dat_transform %>%
mutate(dose = dose_1 + dose_2 * 2 + dose_3 * 3 + dose_4 * 4 +
dose_5 * 5 + dose_6 * 6) %>% select(c(1,2, 9))
# 你的资料column顺序跟我的不一样，select要改不同数字
st = proc.time()
dat_count = dat %>% group_by(id) %>%
mutate(previous_dose = c(0, dose[-length(dose)])) %>%
filter(obs_time > 1) %>%
group_by(obs_time, dose, previous_dose) %>%
summarise(count = n())
transitMatrix = dat_count %>% group_by(obs_time, previous_dose) %>%
mutate(transitP = count / sum(count)) %>% ungroup() %>%
split(.$obs_time) %>%
lapply(function(x) spMatrix(6, 6, x$previous_dose, x$dose,
x$transitP))
proc.time() - st
# user system elapsed
# 0.72 0.05 0.78
随机过程里面应该有定义这个矩阵：transition probability matrix
翻译应该是转移机率矩阵，下面简称转移矩阵
> transitMatrix[[1]] # 第二期的转移矩阵
[1,] 0.1595092 0.1288344 0.1411043 0.2085890 0.1901840 0.1717791
[2,] 0.1506024 0.2108434 0.1506024 0.1626506 0.2108434 0.1144578
[3,] 0.1676647 0.1137725 0.1616766 0.1497006 0.2035928 0.2035928
[4,] 0.2000000 0.1161290 0.1806452 0.1612903 0.1870968 0.1548387
[5,] 0.1707317 0.2256098 0.1890244 0.1341463 0.1463415 0.1341463
[6,] 0.1296296 0.1851852 0.2160494 0.1543210 0.1728395 0.1419753
(1,1)位置就是上一期服dose 1，下一期也是dose 1的条件机率
(1,2)位置就是上一期服dose 1，下一期服dose 2的条件机率
剩下依此类推
PS: 因此，转移矩阵的每一个列，其和为1
备注：
此处用spMatrix不是因为output是sparse matrix (如果人少，则是)
而是程式比较好写，直接透过矩阵位置做initial，省功夫去排列而已
特此说明。