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"))