[程序汇编]Markov Models using R and Matlab

Authors:	Artur Araújo, Luís Meira-Machado, Javier Roca-Pardiñas
Title:	TPmsm: Estimation of the Transition Probabilities in 3-State Models
Abstract:	One major goal in clinical applications of multi-state models is the estimation of transition probabilities. The usual nonparametric estimator of the transition matrix for non-homogeneous Markov processes is the Aalen-Johansen estimator (Aalen and Johansen 1978). However, two problems may arise from using this estimator: first, its standard error may be large in heavy censored scenarios; second, the estimator may be inconsistent if the process is non-Markovian. The development of the R package TPmsm has been motivated by several recent contributions that account for these estimation problems. Estimation and statistical inference for transition probabilities can be performed using TPmsm. The TPmsm package provides seven different approaches to three-state illness-death modeling. In two of these approaches the transition probabilities are estimated conditionally on current or past covariate measures. Two real data examples are included for illustration of software usage.
	Page views:: 2160. Submitted: 2012-07-05. Published: 2014-12-25.
Paper:	TPmsm: Estimation of the Transition Probabilities in 3-State Models     [size=1em]Download PDF (Downloads: 2382)
Supplements:	TPmsm_1.2.0.tar.gz: R source package [size=1em]Download(Downloads: 127; 90KB) v62i04.R: R example code from the paper [size=1em]Download(Downloads: 131; 2KB)

Authors:	Nicole Ferguson, Somnath Datta, Guy Brock
Title:	msSurv: An R Package for Nonparametric Estimation of Multistate Models
Abstract:	We present an R package, msSurv, to calculate the marginal (that is, not conditional on any covariates) state occupation probabilities, the state entry and exit time distributions, and the marginal integrated transition hazard for a general, possibly non-Markov, multistate system under left-truncation and right censoring. For a Markov model, msSurv also calculates and returns the transition probability matrix between any two states. Dependent censoring is handled via modeling the censoring hazard through observable covariates. Pointwise confidence intervals for the above mentioned quantities are obtained and returned for independent censoring from closed-form variance estimators and for dependent censoring using the bootstrap.
	Page views:: 4452. Submitted: 2011-02-08. Published: 2012-09-22.
Paper:	msSurv: An R Package for Nonparametric Estimation of Multistate Models     [size=1em]Download PDF (Downloads: 4483)
Supplements:	msSurv_1.1-2.tar.gz: R source package [size=1em]Download(Downloads: 491; 484KB) v50i14.R: R example code from the paper [size=1em]Download(Downloads: 464; 6KB)

Authors:	Martin D. King, Fernando Calamente, Chris A. Clark, David G. Gadian
Title:	Markov Chain Monte Carlo Random Effects Modeling in Magnetic Resonance Image Processing Using the BRugs Interface to WinBUGS
Abstract:	A common feature of many magnetic resonance image (MRI) data processing methods is the voxel-by-voxel (a voxel is a volume element) manner in which the processing is performed. In general, however, MRI data are expected to exhibit some level of spatial correlation, rendering an independent-voxels treatment inefficient in its use of the data. Bayesian random effect models are expected to be more efficient owing to their information-borrowing behaviour. To illustrate the Bayesian random effects approach, this paper outlines a Markov chain Monte Carlo (MCMC) analysis of a perfusion MRI dataset, implemented in R using the BRugs package. BRugs provides an interface to WinBUGS and its GeoBUGS add-on. WinBUGS is a widely used programme for performing MCMC analyses, with a focus on Bayesian random effect models. A simultaneous modeling of both voxels (restricted to a region of interest) and multiple subjects is demonstrated. Despite the low signal-to-noise ratio in the magnetic resonance signal intensity data, useful model signal intensity profiles are obtained. The merits of random effects modeling are discussed in comparison with the alternative approaches based on region-of-interest averaging and repeated independent voxels analysis. This paper focuses on perfusion MRI for the purpose of illustration, the main proposition being that random effects modeling is expected to be beneficial in many other MRI applications in which the signal-to-noise ratio is a limiting factor.
	Page views:: 3298. Submitted: 2010-10-01. Published: 2011-10-27.
Paper:	Markov Chain Monte Carlo Random Effects Modeling in Magnetic Resonance Image Processing Using the BRugs Interface to WinBUGS     [size=1em]Download PDF (Downloads: 3337)
Supplements:	v44i02-replication.zip: Replication code/data for examples in the paper [size=1em]Download(Downloads: 761; 12KB)

1. #Section 4
   
2. 
   
3. library("TPmsm")
   
4. 
   
5. setThreadsTP(1)
   
6. 
   
7. seed <- c(2718, 3141, 5436, 6282, 8154, 9423)
   
8. 
   
9. setPackageSeedTP(seed)
   
10. 
    
11. 
    
12. 
    
13. sim_data_exp <- dgpTP(n = 1000, corr = 0, dist = "exponential", dist.par = c(1, 1), model.cens = "uniform", cens.par = 3, state2.prob = 0.5)
    
14. 
    
15. transAJ(object = sim_data_exp, s = 0.5108, t = 0.9163, conf = TRUE, conf.level = 0.95, n.boot = 1000)
    
16. 
    
17. transPAJ(object = sim_data_exp, s = 0.5108, t = 0.9163, conf = TRUE, conf.level = 0.95, n.boot = 1000)
    
18. 
    
19. 
    
20. setPackageSeedTP(seed)
    
21. 
    
22. sim_data_exp2 <- dgpTP(n = 1000, corr = 1, dist = "exponential", dist.par = c(1, 1), model.cens = "uniform", cens.par = 3, state2.prob = 0.5)
    
23. 
    
24. transKMW(object = sim_data_exp2, s = 0.5108, t = 0.9163, conf = TRUE, conf.level = 0.95, n.boot = 1000)
    
25. 
    
26. transKMPW(object = sim_data_exp2, s = 0.5108, t = 0.9163, conf = TRUE, conf.level = 0.95, n.boot = 1000)
    
27. 
    
28. 
    
29. 
    
30. #Section 5
    
31. 
    
32. data("colonTP", package = "TPmsm")
    
33. 
    
34. head(head(colonTP[ , c(1:4, 7)]))
    
35. 
    
36. 
    
37. 
    
38. colon_obj <- with(colonTP, survTP(time1, event1, Stime, event, age))
    
39. 
    
40. colon_obj_TP <- transKMW(object = colon_obj, s = 365, t = 1096, conf = TRUE, conf.level = 0.95)
    
41. 
    
42. colon_obj_TP
    
43. 
    
44. 
    
45. 
    
46. colon_obj2_TP <- transKMPW(object = colon_obj, s = 365, t = 1096, conf = TRUE, conf.level = 0.95)
    
47. 
    
48. colon_obj2_TP
    
49. 
    
50. 
    
51. 
    
52. #Figure 2
    
53. 
    
54. colon_obj_TP <- transKMW(object = colon_obj, s = 365, conf = TRUE, conf.level = 0.95)
    
55. 
    
56. plot(colon_obj_TP, col = seq_len(5), lty = 1, ylab = "p_hj(365,t)")
    
57. 
    
58. 
    
59. 
    
60. #Figure 3
    
61. 
    
62. plot(colon_obj_TP, tr.choice = "1 2", conf.int = TRUE, ylim = c(0, 0.2), legend = FALSE, ylab = "p12(365,t)")
    
63. 
    
64. 
    
65. CTP_obj <- transIPCW(colon_obj, s = 365, t = 1096, x = c(40, 68), conf = TRUE, n.boot = 1000, method.boot = "percentile")
    
66. 
    
67. CTP_obj
    
68. 
    
69. 
    
70. 
    
71. #Figure 4
    
72. 
    
73. plot(CTP_obj, plot.type = "c", tr.choice = "1 1", conf.int = TRUE, xlab = "Age", legend = FALSE, ylab = "p11(365,1096|age)")
    
74. 
    
75. 
    
76. 
    
77. #Figure 5
    
78. 
    
79. plot(CTP_obj, plot.type = "c", tr.choice = "1 2", conf.int = TRUE, xlab = "Age", legend = FALSE, ylab = "p12(365,1096|age)")
    
80. 
    
81. 
    
82. 
    
83. #Figure 6
    
84. 
    
85. plot(CTP_obj, plot.type = "c", col = seq_len(5), lty = 1, xlab = "Age", ylab = "p_hj(365,1096|age)")
    
86. 
    
87. 
    
88. 
    
89. data("bladderTP", package = "TPmsm")
    
90. 
    
91. head(bladderTP)
    
92. 
    
93. 
    
94. 
    
95. bladderTP_obj <- with(bladderTP, survTP(time1, event1, Stime, event))
    
96. 
    
97. LS_obj <- transLS(object = bladderTP_obj, s = 3, t = 8, h = c(0.0001, 1), nh = 100, ncv = 100, conf = TRUE)
    
98. 
    
99. LS_obj
    
100. 
     
101. 
     
102. 
     
103. LS2_obj <- transLS(object = bladderTP_obj, s = 3, t = 60, h = c(0.0001, 1), nh = 100, ncv = 100, conf = TRUE)
     
104. 
     
105. 
     
106. 
     
107. #Figure 7
     
108. 
     
109. plot(LS2_obj, col = seq_len(5), lty = 1, ylab = "p_hj(3,t)")
     
110. 
     
111. 
     
112. 
     
113. #Figure 8
     
114. 
     
115. plot(LS2_obj, tr.choice = "1 2", conf.int = TRUE, ylab = "p12(3,t)", ylim = c(0, 0.35), legend = FALSE)

复制代码

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. # Create survTP object
   
5. data(heartTP)
   
6. heartTP_obj <- with(heartTP, survTP(time1, event1, Stime, event))
   
7. # Compute transition probabilities
   
8. transAJ(object=heartTP_obj, s=33, t=412)
   
9. # Compute transition probabilities with confidence band
   
10. transAJ(object=heartTP_obj, s=33, t=412, conf=TRUE, conf.level=0.9,
    
11. method.boot="percentile")
    
12. # Restore the number of threads
    
13. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. # Create survTP object with age as covariate
   
5. data(heartTP)
   
6. heartTP_obj <- with(heartTP, survTP(time1, event1, Stime, event, age=age))
   
7. # Compute unconditioned transition probabilities
   
8. transIPCW(object=heartTP_obj, s=33, t=412)
   
9. # Compute unconditioned transition probabilities with confidence band
   
10. transIPCW(object=heartTP_obj, s=33, t=412, conf=TRUE, conf.level=0.9,
    
11. method.boot="basic", method.est=2)
    
12. # Compute conditional transition probabilities
    
13. transIPCW(object=heartTP_obj, s=33, t=412, x=0)
    
14. # Compute conditional transition probabilities with confidence band
    
15. transIPCW(object=heartTP_obj, s=33, t=412, x=0, conf=TRUE, conf.level=0.95,
    
16. n.boot=100, method.boot="percentile", method.est=2)
    
17. # Restore the number of threads
    
18. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. # Create survTP object
   
5. data(heartTP)
   
6. 38 transKMW
   
7. heartTP_obj <- with(heartTP, survTP(time1, event1, Stime, event))
   
8. # Compute transition probabilities
   
9. transKMPW(object=heartTP_obj, s=33, t=412)
   
10. # Compute transition probabilities with confidence band
    
11. transKMPW(object=heartTP_obj, s=33, t=412, conf=TRUE, conf.level=0.9,
    
12. method.boot="percentile", method.est=4)
    
13. # Restore the number of threads
    
14. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. # Create survTP object
   
5. data(heartTP)
   
6. heartTP_obj <- with(heartTP, survTP(time1, event1, Stime, event))
   
7. # Compute transition probabilities
   
8. transKMW(object=heartTP_obj, s=33, t=412)
   
9. transLIN 41
   
10. # Compute transition probabilities with confidence band
    
11. transKMW(object=heartTP_obj, s=33, t=412, conf=TRUE, conf.level=0.9,
    
12. method.boot="basic", method.est=2)
    
13. # Restore the number of threads
    
14. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. 
   
5. # Create survTP object with age as covariate
   
6. data(heartTP)
   
7. heartTP_obj <- with(heartTP, survTP(time1, event1, Stime, event, age=age))
   
8. 
   
9. # Compute unconditioned transition probabilities
   
10. transLIN(object=heartTP_obj, s=33, t=412)
    
11. 
    
12. # Compute unconditioned transition probabilities with confidence band
    
13. transLIN(object=heartTP_obj, s=33, t=412, conf=TRUE, conf.level=0.9,
    
14. method.boot="basic")
    
15. 
    
16. # Compute conditional transition probabilities
    
17. transLIN(object=heartTP_obj, s=33, t=412, x=0)
    
18. # Compute conditional transition probabilities with confidence band
    
19. transLIN(object=heartTP_obj, s=33, t=412, x=0, conf=TRUE, conf.level=0.95,
    
20. n.boot=100, method.boot="percentile")
    
21. 
    
22. # Restore the number of threads
    
23. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. # Create survTP object
   
5. data(bladderTP)
   
6. bladderTP_obj <- with(bladderTP, survTP(time1, event1, Stime, event))
   
7. # Compute transition probabilities
   
8. LS0 <- transLS(object=bladderTP_obj, s=5, t=59, h=c(0.25, 2.5), nh=25, ncv=50, conf=FALSE)
   
9. print(LS0)
   
10. # Compute transition probabilities with confidence band
    
11. h <- with( LS0, c( rep(h[1], 2), rep(h[2], 2) ) )
    
12. transLS(object=bladderTP_obj, s=5, t=59, h=h, conf=TRUE,
    
13. conf.level=0.95, method.boot="percentile", boot.cv=FALSE)
    
14. # Restore the number of threads
    
15. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf