Scenario II (sim_Markov_group_multinomial.r)
library(nnet)
n.iter<-1
coef.l5 <- matrix(NA,n.iter,30)
z.l5<-matrix(NA,n.iter,30)
p.l5<-matrix(NA,n.iter,30)
sig.l5<-matrix(NA,n.iter,30)
coef.l5.5<-matrix(NA,n.iter,1)
coef.l5.7<-matrix(NA,n.iter,1)
coverage.l5.5<- matrix(NA,n.iter,1)
coverage.l5.7<- matrix(NA,n.iter,1)
coef.l4 <- matrix(NA,n.iter,30)
z.l4<-matrix(NA,n.iter,30)
p.l4<-matrix(NA,n.iter,30)
sig.l4<-matrix(NA,n.iter,30)
coef.l4.5<-matrix(NA,n.iter,1)
coef.l4.7<-matrix(NA,n.iter,1)
coverage.l4.5<- matrix(NA,n.iter,1)
coverage.l4.7<- matrix(NA,n.iter,1)
coef.l3 <- matrix(NA,n.iter,30)
z.l3<-matrix(NA,n.iter,30)
p.l3<-matrix(NA,n.iter,30)
sig.l3<-matrix(NA,n.iter,30)
coef.l3.5<-matrix(NA,n.iter,1)
coef.l3.7<-matrix(NA,n.iter,1)
coverage.l3.5<- matrix(NA,n.iter,1)
coverage.l3.7<- matrix(NA,n.iter,1)
coef.l2 <- matrix(NA,n.iter,30)
z.l2<-matrix(NA,n.iter,30)
p.l2<-matrix(NA,n.iter,30)
sig.l2<-matrix(NA,n.iter,30)
coef.l2.5<-matrix(NA,n.iter,1)
coef.l2.7<-matrix(NA,n.iter,1)
coverage.l2.5<- matrix(NA,n.iter,1)
coverage.l2.7<- matrix(NA,n.iter,1)
coef.tot <- matrix(NA,n.iter,36)
z.tot<-matrix(NA,n.iter,36)
p.tot<-matrix(NA,n.iter,36)
sig.tot<-matrix(NA,n.iter,36)
coef.tot.7<-matrix(NA,n.iter,1)
coverage.tot.7<- matrix(NA,n.iter,1)
for (iter in 1:n.iter){
cat("iter=", iter)
##MARKOV with a group effect######
n<-500
state<-matrix(0,n,10)
ile<-matrix(1,n,1)
time<-matrix(0,n,10)
group<-rbinom(n,1,0.5)
###
time21<-time23<-time24<-time25<-time26<-time27<-matrix(0,n,10)
time31<-time32<-time34<-time35<-time36<-time37<-matrix(0,n,10)
time41<-time42<-time43<-time45<-time46<-time47<-matrix(0,n,10)
time51<-time52<-time53<-time54<-time56<-time57<-matrix(0,n,10)
for ( i in 1: n){
state[i,1]<-5
time[i,1]<-0
for (j in 1:10){
time21[i,j]<-rexp(1,rate=exp(-1))
#time22[i,j]<-rexp(1,rate=exp(-1))
time23[i,j]<-rexp(1,rate=exp(-1))
time24[i,j]<-rexp(1,rate=exp(-1))
time25[i,j]<-rexp(1,rate=exp(-1))
time26[i,j]<-rexp(1,rate=exp(-1))
time27[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
time31[i,j]<-rexp(1,rate=exp(-1))
time32[i,j]<-rexp(1,rate=exp(-1))
#time33[i,j]<-rexp(1,rate=exp(-1))
time34[i,j]<-rexp(1,rate=exp(-1))
time35[i,j]<-rexp(1,rate=exp(-1))
time36[i,j]<-rexp(1,rate=exp(-1))
time37[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
time41[i,j]<-rexp(1,rate=exp(-1))
time42[i,j]<-rexp(1,rate=exp(-1))
time43[i,j]<-rexp(1,rate=exp(-1))
#time44[i,j]<-rexp(1,rate=exp(-1))
time45[i,j]<-rexp(1,rate=exp(-1))
time46[i,j]<-rexp(1,rate=exp(-1))
time47[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
time51[i,j]<-rexp(1,rate=exp(-1))
time52[i,j]<-rexp(1,rate=exp(-1))
time53[i,j]<-rexp(1,rate=exp(-1))
time54[i,j]<-rexp(1,rate=exp(-1))
#time55[i,j]<-rexp(1,rate=exp(-1))
time56[i,j]<-rexp(1,rate=exp(-1))
time57[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
}
}
for ( i in 1:n){
for (j in 2:10){
if (j>2){
prev2.state<-state[i,j-2]
}
prev.state<-state[i,j-1]
if (prev.state==1){
state[i,j]<-1}
if (prev.state==6){
state[i,j]<-6}
if (prev.state==7){
state[i,j]<-7}
if (prev.state==5){
time[i,j]<-min(time51[i,j],time52[i,j],time53[i,j],time54[i,j],9999,time56[i,j],time57[i,j])
times<-c(time51[i,j],time52[i,j],time53[i,j],time54[i,j],9999,time56[i,j],time57[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time51[i,j],time52[i,j],time53[i,j],time54[i,j],9999,time56[i,j],time57[i,j]))
}
if (prev.state==4){
time[i,j]<-min(time41[i,j],time42[i,j],time43[i,j],9999,time45[i,j],time46[i,j],time47[i,j])
times<-c(time41[i,j],time42[i,j],time43[i,j],999,time45[i,j],time46[i,j],time47[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time41[i,j],time42[i,j],time43[i,j],9999,time45[i,j],time46[i,j],time47[i,j]))
}
if (prev.state==3){
time[i,j]<-min(time31[i,j],time32[i,j],9999,time34[i,j],time35[i,j],time36[i,j],time37[i,j])
times<-c(time31[i,j],time32[i,j],9999,time34[i,j],time35[i,j],time36[i,j],time37[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time31[i,j],time32[i,j],9999,time34[i,j],time35[i,j],time36[i,j],time37[i,j]))
}
if (prev.state==2){
time[i,j]<-min(time21[i,j],9999,time23[i,j],time24[i,j],time25[i,j],time26[i,j],time27[i,j])
times<-c(time21[i,j],9999,time23[i,j],time24[i,j],time25[i,j],time26[i,j],time27[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time21[i,j],9999,time23[i,j],time24[i,j],time25[i,j],time26[i,j],time27[i,j]))
}
}#end for i
}# end for j
ile<-matrix(0,n,1)
for ( i in 1:n){
ile[i,1]<-sum(as.numeric(time[i,]>0))
}
##total nr trans##
total.trans<-sum(ile)
###
###DATA WITH GROUP####
id.data<-matrix(0,total.trans,1)
group.data<-matrix(0,total.trans,1)
from<-matrix(0,total.trans,1)
to<-matrix(0,total.trans,1)
tstart<-tstop<-matrix(0,total.trans,1)
ile.juz<-0
for (i in 1: n){
for ( j in 1: ile[i,]){
id.data[ile.juz+j]<-i
group.data[ile.juz+j]<-group[i]
if (j==1){
from[ile.juz+j]<-state[i,j]
to[ile.juz+j]<-state[i,j+1]
tstart[ile.juz+j]<-0
tstop[ile.juz+j]<-time[i,j+1]
}
if(j>1){
from[ile.juz+j]<-to[ile.juz+j-1]
to[ile.juz+j]<-state[i,j+1]
tstart[ile.juz+j]<-tstop[ile.juz+j-1]
tstop[ile.juz+j]<-tstop[ile.juz+j-1]+time[i,j]
}
}
ile.juz<-ile.juz+ile[i,]
}
heart.sym<-as.data.frame(cbind(id.data,group.data,from,to,tstart,tstop))
names(heart.sym)<-c("id","group","from","to","tstart","tstop")
heart.sym$time<-heart.sym$tstop
###
head(heart.sym)
####history dependence ###
###creating history- a list####
heart.sym$hist<-list("N")
num<-1
for (i in 2:(length(heart.sym$id))){
current.id<-heart.sym$id[i]
id.prev<-heart.sym$id[i-1]
if(current.id!=id.prev){num<-1}
if (current.id==id.prev){
if(heart.sym$hist[i-1]=="N"){
num<-2
heart.sym$hist[[i]][num-1]<-as.character(heart.sym$from[i-1])}else{
num<-length(heart.sym$hist[[i-1]])
for( k in 1:num){
heart.sym$hist[[i]][k]<-as.character(heart.sym$hist[[i-1]][k])}
heart.sym$hist[[i]][num+1]<-as.character(heart.sym$from[i-1])
}
}}
head(heart.sym)
rownames(heart.sym)<-seq(1:length(heart.sym$id))
###counting the length of the history###
count<-matrix(0,1,length(heart.sym$id))
for ( i in 1: length(heart.sym$id)){
count[,i]<-as.integer(length(heart.sym$hist[[i]]))
}
heart.sym$count<-t(count)
###
lasteq5<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="5")
lasteq5[i,]<-1}
heart.sym$lasteq5<-lasteq5
###last 2
lasteq2<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="2")
lasteq2[i,]<-1}
heart.sym$lasteq2<-lasteq2
##last 3###
lasteq3<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="3")
lasteq3[i,]<-1}
heart.sym$lasteq3<-lasteq3
##last 4###
lasteq4<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="4")
lasteq4[i,]<-1}
heart.sym$lasteq4<-lasteq4
####history 5 selected####
heart.sym.l5<-heart.sym[heart.sym$lasteq5==1,]
#heart.sym.l5$id<-seq(1:length(heart.sym.l5$id))
id.set.l5<-unique(heart.sym.l5$id)
n.id.l5 <-length(id.set.l5)
jack.se.l5<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l5)
for (i in 1:n.id.l5) {
id.c<-heart.sym.l5$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l5[heart.sym.l5$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l5[r] <- sqrt(((n.id.l5 - 1)/n.id.l5) * sum((u - mean(u))^2))
}
coef.l5[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5))
z.l5[iter,]<-(coef.l5[iter,1:30])/jack.se.l5
p.l5[iter,]<-2-2*pnorm(abs(z.l5[iter,]))
sig.l5[iter,]<-as.numeric(p.l5[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l5.5[iter]<-coef.l5[iter,22]
coef.l5.7[iter]<-coef.l5[iter,24]
coverage.l5.5[iter]<- sig.l5[iter,22]
coverage.l5.7[iter]<- sig.l5[iter,24]
names<-colnames(coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5)))
model.show<-coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5))
colnames(model.show)<-names
jack.se.show<-matrix(jack.se.l5,nrow=6,byrow=TRUE)
colnames(jack.se.show)<-names
model.show
jack.se.show
######we estimate transition times for the subgroup with pooling individuals###
###history 4 selected####
heart.sym.l4<-heart.sym[heart.sym$lasteq4==1,]
id.set.l4<-unique(heart.sym.l4$id)
n.id.l4 <-length(id.set.l4)
jack.se.l4<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l4)
for (i in 1:n.id.l4) {
id.c<-heart.sym.l5$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l4[heart.sym.l4$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l4[r] <- sqrt(((n.id.l4 - 1)/n.id.l4) * sum((u - mean(u))^2))
}
coef.l4[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l4))
z.l4[iter,]<-(coef.l4[iter,1:30])/jack.se.l4
p.l4[iter,]<-2-2*pnorm(abs(z.l4[iter,]))
sig.l4[iter,]<-as.numeric(p.l4[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l4.5[iter]<-coef.l4[iter,22]
coef.l4.7[iter]<-coef.l4[iter,24]
coverage.l4.5[iter]<- sig.l4[iter,22]
coverage.l4.7[iter]<- sig.l4[iter,24]
#####last 3
heart.sym.l3<-heart.sym[heart.sym$lasteq3==1,]
id.set.l3<-unique(heart.sym.l3$id)
n.id.l3 <-length(id.set.l3)
jack.se.l3<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l3)
for (i in 1:n.id.l3) {
id.c<-heart.sym.l5$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l3[heart.sym.l3$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l3[r] <- sqrt(((n.id.l3 - 1)/n.id.l3) * sum((u - mean(u))^2))
}
coef.l3[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l3))
z.l3[iter,]<-(coef.l3[iter,1:30])/jack.se.l3
p.l3[iter,]<-2-2*pnorm(abs(z.l3[iter,]))
sig.l3[iter,]<-as.numeric(p.l3[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l3.5[iter]<-coef.l3[iter,22]
coef.l3.7[iter]<-coef.l3[iter,24]
coverage.l3.5[iter]<- sig.l3[iter,22]
coverage.l3.7[iter]<- sig.l3[iter,24]
###last 2
heart.sym.l2<-heart.sym[heart.sym$lasteq2==1,]
id.set.l2<-unique(heart.sym.l2$id)
n.id.l2 <-length(id.set.l2)
jack.se.l2<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l2)
for (i in 1:n.id.l2) {
id.c<-heart.sym.l2$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l2[heart.sym.l2$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l2[r] <- sqrt(((n.id.l2 - 1)/n.id.l2) * sum((u - mean(u))^2))
}
coef.l2[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l2))
z.l2[iter,]<-(coef.l2[iter,1:30])/jack.se.l2
p.l2[iter,]<-2-2*pnorm(abs(z.l2[iter,]))
sig.l2[iter,]<-as.numeric(p.l2[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l2.5[iter]<-coef.l2[iter,22]
coef.l2.7[iter]<-coef.l2[iter,24]
coverage.l2.5[iter]<- sig.l2[iter,22]
coverage.l2.7[iter]<- sig.l2[iter,24]
###total to asses better 7####
id.set.tot<-unique(heart.sym$id)
n.id.tot <-length(id.set.tot)
jack.se.tot<-rep(0,36)
for ( r in 1:36){
u <- rep(0,n.id.tot)
for (i in 1:n.id.tot) {
id.c<-heart.sym$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym[heart.sym$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.tot[r] <- sqrt(((n.id.tot - 1)/n.id.tot) * sum((u - mean(u))^2))
}
coef.tot[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym))
z.tot[iter,]<-(coef.tot[iter,1:36])/jack.se.tot
p.tot[iter,]<-2-2*pnorm(abs(z.tot[iter,]))
sig.tot[iter,]<-as.numeric(p.tot[iter,]<0.05)
coef.tot.7[iter]<-coef.tot[iter,30]
coverage.tot.7[iter]<- sig.tot[iter,30]
}
iter= 1# weights: 42 (30 variable)
initial value 496.207088
iter 10 value 450.266832
iter 20 value 443.185755
iter 30 value 442.008973
iter 40 value 441.943422
final value 441.943075
converged
# weights: 42 (30 variable)
initial value 496.207088
iter 10 value 450.574806
iter 20 value 442.711864
iter 30 value 441.630182
iter 40 value 441.551080
final value 441.550672
converged
# weights: 42 (30 variable)
initial value 496.207088
iter 10 value 450.191340
iter 20 value 442.363568
iter 30 value 441.278253
iter 40 value 441.199687
final value 441.199294
converged
# weights: 42 (30 variable)
initial value 496.207088
iter 10 value 451.592453
iter 20 value 443.171687
iter 30 value 441.999192
iter 40 value 441.924244
final value 441.923722
converged
....
Showing results for multinomial model with last state equal T:
names<-colnames(coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5)))
model.show<-coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5))
colnames(model.show)<-names
jack.se.show<-matrix(jack.se.l5,nrow=6,byrow=TRUE)
colnames(jack.se.show)<-names
model.show
Intercept) as.factor(from)3 as.factor(from)4 as.factor(group)1 time
2 -14.93685767 15.98771262 14.71544875 0.1546722 -0.35837483
3 0.14035175 -16.48440951 0.39385083 0.1019955 -0.20621497
4 -0.20784156 0.03801484 -17.48415044 0.8131844 0.21955774
5 0.46694934 -0.40290016 -0.59334587 0.1428662 -0.09929488
6 -0.05077117 0.59304129 0.01524984 0.2568383 -0.17309359
7 -0.32232849 0.34966424 -0.06349466 1.0689938 0.12355032
Jackknife SE:
jack.se.show
(Intercept) as.factor(from)3 as.factor(from)4 as.factor(group)1 time
[1,] 5.0690517 0.5656984 0.5542705 0.5026805 0.5053540
[2,] 0.4846059 4.9992243 4.2356059 0.5741753 0.5097811
[3,] 0.5803959 0.5374176 4.9597613 0.5476981 4.9874812
[4,] 0.5467425 0.5323433 0.5253471 0.5546750 0.4854521
[5,] 0.5575211 0.4612572 0.4613411 0.4250448 0.3183462
[6,] 0.2858649 0.3520166 0.2856358 0.3130389 0.2773092
list.l5<-list( coef.l5=coef.l5,z.l5=z.l5,p.l5=p.l5,sig.l5=sig.l5,coef.l5.5=coef.l5.5,coef.l5.7=coef.l5.7,coverage.l5.5=coverage.l5.5,coverage.l5.7=coverage.l5.7)
list.l4<-list( coef.l4=coef.l4,z.l4=z.l4,p.l4=p.l4,sig.l4=sig.l4,coef.l4.5=coef.l4.5,coef.l4.7=coef.l4.7,coverage.l4.5=coverage.l4.5,coverage.l4.7=coverage.l4.7)
list.l3<-list( coef.l3=coef.l3,z.l3=z.l3,p.l3=p.l3,sig.l3=sig.l3,coef.l3.5=coef.l3.5,coef.l3.7=coef.l3.7,coverage.l3.5=coverage.l3.5,coverage.l3.7=coverage.l3.7)
list.l2<-list( coef.l2=coef.l2,z.l2=z.l2,p.l2=p.l2,sig.l2=sig.l2,coef.l2.5=coef.l2.5,coef.l2.7=coef.l2.7,coverage.l2.5=coverage.l2.5,coverage.l2.7=coverage.l2.7)
list.l.tot<-list( coef.tot=coef.tot,z.tot=z.tot,p.tot=p.tot,sig.tot=sig.tot,coef.tot.7=coef.tot.7,coverage.tot.7=coverage.tot.7)
#save(list.l5, file = file.path("list.l5_Markov_group.RData"))
#save(list.l4, file = file.path("list.l4_Markov_group.RData"))
#save(list.l3, file = file.path("list.l3_Markov_group.RData"))
#save(list.l2, file = file.path("list.l2_Markov_group.RData"))
#save(list.l.tot, file = file.path("list.l.tot_Markov_group.RData"))