Scenario IV (sim_non_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]


if (prev2.state==5){


time25[i,j]<-time25[i,j]/(group[i]+1)
time35[i,j]<-time35[i,j]/(group[i]+1)
time45[i,j]<-time45[i,j]/(group[i]+1)

        }

if (prev2.state==4){

time24[i,j]<-time24[i,j]/(group[i]+1)
time34[i,j]<-time34[i,j]/(group[i]+1)
time54[i,j]<-time54[i,j]/(group[i]+1)

        }

if (prev2.state==3){

time23[i,j]<-time23[i,j]/(group[i]+1)
time43[i,j]<-time43[i,j]/(group[i]+1)
time53[i,j]<-time53[i,j]/(group[i]+1)

        }
    }




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
###
####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])
}
}}


#####creating history


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]

######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 (7)#

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 476.747987 
iter  10 value 431.933646
iter  20 value 424.189870
iter  30 value 423.130775
iter  40 value 423.060801
final  value 423.059918 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 433.356309
iter  20 value 424.703050
iter  30 value 423.666997
iter  40 value 423.570330
final  value 423.569537 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 432.273103
iter  20 value 424.536869
iter  30 value 423.513515
iter  40 value 423.418895
final  value 423.417887 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 431.257315
iter  20 value 424.137881
iter  30 value 423.154080
iter  40 value 423.066380
final  value 423.065316 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 430.855419
iter  20 value 424.640749
iter  30 value 423.600830
iter  40 value 423.524824
final  value 423.524152 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 432.375229
iter  20 value 424.217145
iter  30 value 422.993218
iter  40 value 422.895711
final  value 422.894611 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 432.428738
iter  20 value 424.115228
iter  30 value 423.134155
iter  40 value 423.004464
final  value 423.002782 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 434.204872
iter  20 value 424.584257
iter  30 value 423.343592
iter  40 value 423.240819
final  value 423.240024 
converged
# weights:  42 (30 variable)
initial  value 476.747987 
iter  10 value 431.199322
iter  20 value 424.696606
iter  30 value 423.700171
iter  40 value 423.612941
final  value 423.612031 
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 -15.2893928       14.7929774       15.7293827       -0.24395013 0.3325791
3   0.2919006      -16.6952682        0.4077411       -0.40012116 0.1487177
4  -0.5179829        0.1592948      -16.5820002        0.01348601 0.4403449
5   0.3346142       -0.3547416        0.7267854       -0.09079121 0.1717782
6   0.3856212       -0.6485112        0.1653082       -0.98617957 0.1711206
7  -0.2223938       -0.2254897       -0.1327517        0.70863071 0.2872732

Jackknife SE:
jack.se.show

(Intercept) as.factor(from)3 as.factor(from)4 as.factor(group)1      time
[1,]   2.9278207        0.5421481        0.7927232         0.5576469 0.5590972
[2,]   0.5506199        2.9393823        5.4270839         0.6940920 0.5796293
[3,]   0.6154395        0.5647047        2.9785828         0.6222915 6.9418369
[4,]   0.6510179        0.6447650        0.6376941         0.5658553 0.5417642
[5,]   0.6438745        0.4832308        0.5345415         0.4929255 0.3881445
[6,]   0.3926737        0.4770983        0.3703444         0.3865607 0.3625086


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_nonMarkov_group.RData"))
         save(list.l4,  file = file.path("list.l4_nonMarkov_group.RData"))
         save(list.l3,  file = file.path("list.l3_nonMarkov_group.RData"))
         save(list.l2,  file = file.path("list.l2_nonMarkov_group.RData"))
         save(list.l.tot,  file = file.path("list.l.tot_nonMarkov_group.RData"))