File: analyticp.R

package info (click to toggle)
r-cran-msm 1.1-1
  • links: PTS, VCS
  • area: main
  • in suites: wheezy
  • size: 1,740 kB
  • sloc: ansic: 2,700; makefile: 14
file content (248 lines) | stat: -rw-r--r-- 26,109 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
source("local.R")
library(msm)

### SIMULATIONS TO TEST ANALYTIC P MATRICES

if (developer.local) { 

nsubj <- 50; nobspt <- 6
sim.df <- data.frame(subject = rep(1:nsubj, each=nobspt), time = seq(0, 20, length=nobspt), x = rnorm(nsubj*nobspt), y = rnorm(nsubj*nobspt)* 5 + 2 )
set.seed(22061976)

## 2 STATES
## 1
(two.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-2)), c(0, 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=two.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-3)), c(0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-3)), c(0, 0)), analyticp=FALSE))

## 1,2
(two.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-2)), c(exp(-2), 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=two.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-3)), c(exp(-1), 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-3)), c(exp(-1), 0)), analyticp=FALSE))

## 3 STATES
## 1,2 
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), exp(-3)), c(0, 0, 0), c(0, 0, 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-1)), c(0, 0, 0), c(0, 0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-1)), c(0, 0, 0), c(0, 0, 0)), analyticp=FALSE))

## 1,4
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), 0), c(0, 0, exp(-3)), c(0, 0, 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0), c(0, 0, exp(-1)), c(0, 0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0), c(0, 0, exp(-1)), c(0, 0, 0)), analyticp=FALSE))

# 4,5 (== 1,4)
nsubj <- 500; nobspt <- 6
sim.df <- data.frame(subject = rep(1:nsubj, each=nobspt), time = seq(0, 20, length=nobspt), 
                     x = rnorm(nsubj*nobspt), y = rnorm(nsubj*nobspt)* 5 + 2 )
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, 0, 0), c(0, 0, exp(-2)), c(exp(-3), 0, 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q, start=rep(2,500))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, 0, 0), c(0, 0, exp(-3)), c(exp(-2), 0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, 0, 0), c(0, 0, exp(-3)), c(exp(-2), 0, 0)), analyticp=FALSE))

## 1,6
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), 0), c(0, 0, 0), c(0, exp(-3), 0))))
nsubj <- 50; nobspt <- 6
sim.df <- data.frame(subject = rep(1:nsubj, each=nobspt), time = seq(0, 20, length=nobspt), x = rnorm(nsubj*nobspt), y = rnorm(nsubj*nobspt)* 5 + 2 )
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q, start=c(rep(3,25),rep(1,25)))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0), c(0, 0, 0), c(0, exp(-1), 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0), c(0, 0, 0), c(0, exp(-1), 0)), analyticp=FALSE))

#1,2,4 (= 3,4,5)
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, 0, 0), c(exp(-3), 0, exp(-3)), c(exp(-3), 0, 0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q, start=rep(2, 50))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, 0, 0), c(exp(-2), 0, exp(-2)), c(2*exp(-2), 0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, 0, 0), c(exp(-2), 0, exp(-2)), c(2*exp(-2), 0, 0)), analyticp=FALSE))

# 1,2,4
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), exp(-6)), c(0, 0, exp(-3)), c(0, 0, 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-2)), c(0, 0, exp(-1)), c(0, 0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-2)), c(0, 0, exp(-1)), c(0, 0, 0)), analyticp=FALSE))

# 1,2,4 (=1,2,6)
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), exp(-6)), c(0, 0, 0), c(0, exp(-3), 0))))
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-2)), c(0, 0, 0), c(0, exp(-2), 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-2)), c(0, 0, 0), c(0, exp(-2), 0)), analyticp=FALSE))

#1,3,5  (= 2,4,5)
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, 0, exp(-3)), c(0, 0, exp(-3)), c(exp(-3), 0, 0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q, start=c(rep(1, 25), rep(2,25)))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, 0, exp(-1)), c(0, 0, exp(-1)), c(exp(-2), 0, 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, 0, exp(-1)), c(0, 0, exp(-1)), c(exp(-2), 0, 0)), analyticp=FALSE))

#1,2,4,6
(three.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), exp(-3)), c(0, 0, exp(-3)), c(0, exp(-3), 0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=three.q, start=c(rep(1, 25), rep(2,25)))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-1)), c(0, 0, exp(-1)), c(0, exp(-2), 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-1)), c(0, 0, exp(-1)), c(0, exp(-2), 0)), analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-1)), c(0, 0, exp(-1)), c(0, exp(-1), 0)), analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), exp(-1)), c(0, 0, exp(-1)), c(0, exp(-1), 0)), analyticp=FALSE))

## FOUR STATES
#1,5,9
(four.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), 0, 0), c(0, 0, exp(-3), 0), c(0, 0, 0, exp(-3)), c(0,0,0,0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=four.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-1)), c(0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-1)), c(0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-2)), c(0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-2)), c(0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0, 0), c(0, 0, exp(-2), 0), c(0, 0, 0, exp(-1)), c(0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-1), 0, 0), c(0, 0, exp(-2), 0), c(0, 0, 0, exp(-1)), c(0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-1)), c(0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-1)), c(0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-3)), c(0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, 0), c(0, 0, exp(-1), 0), c(0, 0, 0, exp(-3)), c(0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

#13569
(four.q <- msm.fixdiag.qmatrix(rbind(c(0, exp(-3), 0, exp(-3)), c(0, 0, exp(-3), exp(-3)), c(0, 0, 0, exp(-3)), c(0,0,0,0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=four.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-2), exp(-2)), c(0, 0, 0, 2*exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-2), exp(-2)), c(0, 0, 0, 2*exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-2), exp(-2)), c(0, 0, 0, exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-2), exp(-2)), c(0, 0, 0, exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-3), exp(-3)), c(0, 0, 0, 2*exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-3), exp(-3)), c(0, 0, 0, 2*exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-3), 0, exp(-3)), c(0, 0, exp(-2), exp(-2)), c(0, 0, 0, 2*exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-3), 0, exp(-3)), c(0, 0, exp(-2), exp(-2)), c(0, 0, 0, 2*exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-3), exp(-3)), c(0, 0, 0, exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0, exp(-2), 0, exp(-2)), c(0, 0, exp(-3), exp(-3)), c(0, 0, 0, exp(-2)), c(0,0,0,0)), fixedpars=FALSE, analyticp=FALSE))


## FIVE STATES
#1_6_11_16
(five.q <- msm.fixdiag.qmatrix(rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=five.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-4),0), c(0,0,0,0,exp(-5)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-4),0), c(0,0,0,0,exp(-5)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-3),0), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-3),0), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-4),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-4),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-4),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-4),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-4),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-4),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-3)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-3)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-3),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-3),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),0,0), c(0,0,0,exp(-2),0), c(0,0,0,0,exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

#1_4_6_8_11_12_16
(five.q <- msm:::msm.fixdiag.qmatrix(rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-2)), c(0,0,0,0,0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=five.q)
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,exp(-5)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,exp(-5)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-3),exp(-3)), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-3),exp(-3)), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-4)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-4),0,exp(-4)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-4),0,exp(-4)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-3)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,exp(-3)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-3),exp(-3)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-3),exp(-3)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-3),0,exp(-3)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,exp(-3)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,exp(-2)), c(0,0,exp(-2),0,exp(-2)), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,2*exp(-2)), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

#1_6_7_11_12
(five.q <- msm:::msm.fixdiag.qmatrix(rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0))))
set.seed(22061976)
sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=five.q)

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),exp(-3),0), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-3),exp(-3),0), c(0,0,0,exp(-4),exp(-4)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,2*exp(-2),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-3),exp(-3)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,2*exp(-2),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-3),exp(-3)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,2*exp(-2),0,0,0), c(0,0,exp(-3),exp(-3),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,2*exp(-2),0,0,0), c(0,0,exp(-3),exp(-3),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-3),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))

(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=TRUE))
(sim.mod <- msm(state ~ time, subject=subject, data=sim2.df, qmatrix = rbind(c(0,exp(-2),0,0,0), c(0,0,exp(-2),exp(-2),0), c(0,0,0,exp(-2),exp(-2)), c(0,0,0,0,0), c(0,0,0,0,0)), fixedpars=TRUE, analyticp=FALSE))


### APPLICATIONS 
## one-way cav misc
## Takes 34% of the time using new method. 
oneway4.q <- rbind(c(0, 0.148, 0, 0.0171), c(0, 0, 0.202, 0.081), c(0, 0, 0, 0.126), c(0, 0, 0, 0))
rownames(oneway4.q) <- colnames(oneway4.q) <- c("Well","Mild","Severe","Death")
ematrix <- rbind(c(0, 0.1, 0, 0),c(0.1, 0, 0.1, 0),c(0, 0.1, 0, 0),c(0, 0, 0, 0))
system.time(misc.msm <- msm(state ~ years, subject = PTNUM, data = cav,
                qmatrix = oneway4.q, ematrix=ematrix, death = 4, fixedpars=FALSE, analyticp=FALSE,
                control = list(trace=1, REPORT=1), method="BFGS"))
misc.msm
system.time(misc.msm <- msm(state ~ years, subject = PTNUM, data = cav,
                qmatrix = oneway4.q, ematrix=ematrix, death = 4, fixedpars=FALSE, analyticp=TRUE,
                control = list(trace=1, REPORT=1), method="BFGS"))
misc.msm

## psoriatic
## half of the time using new method
psor <- read.table("../data/psor.txt", header=TRUE)
psor.q <- rbind(c(0,0.1,0,0),c(0,0,0.1,0),c(0,0,0,0.1),c(0,0,0,0))
system.time(psor.msm <- msm(state ~ months, subject=ptnum, data=psor,
                            qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn,
                            constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)),
                            fixedpars=FALSE, control = list(REPORT=1,trace=2), method="BFGS",analyticp=FALSE))
system.time(psor.msm <- msm(state ~ months, subject=ptnum, data=psor,
                            qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn,
                            constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)),
                            fixedpars=FALSE, control = list(REPORT=1,trace=2), method="BFGS",analyticp=TRUE))


## Jean-Luc
c2.df <- read.table("~/msm/tests/jeanluc/donneesaveccancerPT.txt", header=TRUE)
qx <- rbind( c(0, 0.005, 0, 0, 0), c(0, 0, 0.01, 0.02,0), c(0, 0, 0, 0.04, 0.03), c(0, 0, 0, 0, 0), c(0, 0, 0, 0, 0))
system.time(c2.msm <- msm(state~years, subject=PTNUM, data=c2.df,
              qmatrix=qx, death=c(4, 5), method="BFGS", # fixedpars = 1:5, 
              control=list(trace=2, REPORT=1, fnscale=100000)))
c2.msm
system.time(c2.msm <- msm(state~years, subject=PTNUM, data=c2.df, analyticp=FALSE,
              qmatrix=qx, death=c(4, 5), method="BFGS", # fixedpars = 1:5, 
              control=list(trace=2, REPORT=1, fnscale=100000)))
c2.msm

}