Scenario IV (sim_non_Markov_group_pseudo_multivariate.r)

library(etm)
Loading required package: lattice
Loading required package: survival
Loading required package: splines
library(geepack)

n.iter<-1



bias.1<-matrix(NA,1,n.iter)
stand.se.1<-matrix(NA,1,n.iter)
jack.se.1<-matrix(NA,1,n.iter)
stand.cov.1<-matrix(NA,1,n.iter)


bias.1.mult<--matrix(NA,1,n.iter)
stand.se.1.mult<-matrix(NA,1,n.iter)
jack.se.1.mult<-matrix(NA,1,n.iter)
stand.cov.1.mult<-matrix(NA,1,n.iter)


bias.2<-matrix(NA,1,n.iter)
stand.se.2<-matrix(NA,1,n.iter)
jack.se.2<-matrix(NA,1,n.iter)
stand.cov.2<-matrix(NA,1,n.iter)


bias.2.mult<--matrix(NA,1,n.iter)
stand.se.2.mult<-matrix(NA,1,n.iter)
jack.se.2.mult<-matrix(NA,1,n.iter)
stand.cov.2.mult<-matrix(NA,1,n.iter)



bias.3<-matrix(NA,1,n.iter)
stand.se.3<-matrix(NA,1,n.iter)
jack.se.3<-matrix(NA,1,n.iter)
stand.cov.3<-matrix(NA,1,n.iter)


bias.3.mult<--matrix(NA,1,n.iter)
stand.se.3.mult<-matrix(NA,1,n.iter)
jack.se.3.mult<-matrix(NA,1,n.iter)
stand.cov.3.mult<-matrix(NA,1,n.iter)


bias.4<-matrix(NA,1,n.iter)
stand.se.4<-matrix(NA,1,n.iter)
jack.se.4<-matrix(NA,1,n.iter)
stand.cov.4<-matrix(NA,1,n.iter)

bias.4.mult<--matrix(NA,1,n.iter)
stand.se.4.mult<-matrix(NA,1,n.iter)
jack.se.4.mult<-matrix(NA,1,n.iter)
stand.cov.4.mult<-matrix(NA,1,n.iter)

bias.5<-matrix(NA,1,n.iter)
stand.se.5<-matrix(NA,1,n.iter)
jack.se.5<-matrix(NA,1,n.iter)
stand.cov.5<-matrix(NA,1,n.iter)


bias.5.mult<--matrix(NA,1,n.iter)
stand.se.5.mult<-matrix(NA,1,n.iter)
jack.se.5.mult<-matrix(NA,1,n.iter)
stand.cov.5.mult<-matrix(NA,1,n.iter)

bias.6<-matrix(NA,1,n.iter)
stand.se.6<-matrix(NA,1,n.iter)
jack.se.6<-matrix(NA,1,n.iter)
stand.cov.6<-matrix(NA,1,n.iter)


bias.6.mult<--matrix(NA,1,n.iter)
stand.se.6.mult<-matrix(NA,1,n.iter)
jack.se.6.mult<-matrix(NA,1,n.iter)
stand.cov.6.mult<-matrix(NA,1,n.iter)

bias.7<-matrix(NA,1,n.iter)
stand.se.7<-matrix(NA,1,n.iter)
jack.se.7<-matrix(NA,1,n.iter)
stand.cov.7<-matrix(NA,1,n.iter)

bias.7.mult<--matrix(NA,1,n.iter)
stand.se.7.mult<-matrix(NA,1,n.iter)
jack.se.7.mult<-matrix(NA,1,n.iter)
stand.cov.7.mult<-matrix(NA,1,n.iter)




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

##ETM########

library(etm)
# Possible transitions.
tra <- matrix(ncol=7,nrow=7,FALSE)
tra[1:7, 1:7] <- TRUE
tra[1,1]<-FALSE
tra[2,2]<-FALSE
tra[3,3]<-FALSE
tra[4,4]<-FALSE
tra[5,5]<-FALSE
tra[1,]<-FALSE
tra[6,]<-FALSE
tra[7,]<-FALSE



etm.sym <- etm(heart.sym, c("1", "2", "3","4","5","6","7"), tra, NULL,0)

library(lattice)

# plotting
#xyplot(etm.sym, tr.choice=c("2 1", "2 3", "2 4", "2 5","2 6","2 7","3 1","3 2","3 4","3 5","3 6", "3 7", "4 1", "4 2", "4 3", "4 5", "4 6", "4 7","5 1", "5 2","5 3", "5 4","5 6","5 7"),layout=c(6, 7),
#strip=strip.custom(bg="white",factor.levels=c("2 to 1", "2 to 3", "2 to 4", "2 to 5"," 2 to 6", "2 to 7","3 to 1","3 to 2","3 to 4","3 to 5","3 to 6", "3 to 7",
#"4 to 1", "4 to 2", "4 to 3", "4 to 5", "4 to 6", "4 to 7","5 to 1", "5 to 2", "5 to 3", "5 to 4"," 5 to 6", "5 to 7")))
####

#######for GROUPS ########

etm.sym.0<- etm(heart.sym[heart.sym$group==0,], c("1", "2", "3","4","5","6","7"), tra, NULL,0)
etm.sym.1<- etm(heart.sym[heart.sym$group==1,], c("1", "2", "3","4","5","6","7"), tra, NULL,0)

#####calculating Occupation Probabilities##################
##for state 5####
new.model.4<-as.data.frame(summary(etm.sym)[[4]])
new.model.10<-as.data.frame(summary(etm.sym)[[10]])
new.model.16<-as.data.frame(summary(etm.sym)[[16]])
new.model.28<-as.data.frame(summary(etm.sym)[[28]])

####
n.all<-new.model.4$n.risk[1]+new.model.10$n.risk[1]+new.model.16$n.risk[1]+new.model.28$n.risk[1]
tall.i<-(n.all-1)
###

p2.i<-new.model.4$n.risk[1]/tall.i
p3.i<-new.model.10$n.risk[1]/tall.i
p4.i<-new.model.16$n.risk[1]/tall.i
p5.i<-new.model.28$n.risk[1]/tall.i

p5<-p2.i*new.model.4[,1]+p3.i*new.model.10[,1]+p4.i*new.model.16[,1]+p5.i*new.model.28[,1]

#lines(new.model.4[,2],p5,xlab=list("Days",cex=1.5),ylab=list("Occupation Probability 5",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1))


####for group 0 state 5###
new.model.4.0<-as.data.frame(summary(etm.sym.0)[[4]])
new.model.10.0<-as.data.frame(summary(etm.sym.0)[[10]])
new.model.16.0<-as.data.frame(summary(etm.sym.0)[[16]])
new.model.28.0<-as.data.frame(summary(etm.sym.0)[[28]])
n.all.0<-new.model.4.0$n.risk[1]+new.model.10.0$n.risk[1]+new.model.16.0$n.risk[1]+new.model.28.0$n.risk[1]
tall.i.0<-(n.all.0-1)
p2.i.0<-new.model.4.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.10.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.16.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.28.0$n.risk[1]/tall.i.0

p5.0<-p2.i.0*new.model.4.0[,1]+p3.i.0*new.model.10.0[,1]+p4.i.0*new.model.16.0[,1]+p5.i.0*new.model.28.0[,1]

### group 1 state 5
new.model.4.1<-as.data.frame(summary(etm.sym.1)[[4]])
new.model.10.1<-as.data.frame(summary(etm.sym.1)[[10]])
new.model.16.1<-as.data.frame(summary(etm.sym.1)[[16]])
new.model.28.1<-as.data.frame(summary(etm.sym.1)[[28]])
n.all.1<-new.model.4.1$n.risk[1]+new.model.10.1$n.risk[1]+new.model.16.1$n.risk[1]+new.model.28.1$n.risk[1]
tall.i.1<-(n.all.1-1)
p2.i.1<-new.model.4.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.10.1$n.risk[1]/tall.i.1
p4.i.1<-new.model.16.1$n.risk[1]/tall.i.1
p5.i.1<-new.model.28.1$n.risk[1]/tall.i.1

p5.1<-p2.i.1*new.model.4.1[,1]+p3.i.1*new.model.10.1[,1]+p4.i.1*new.model.16.1[,1]+p5.i.1*new.model.28.1[,1]


plot(new.model.4.0[,2],p5.0,xlab=list("Days",cex=1.5),ylab=list("Occupation Probability ",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),lty=1)
lines(new.model.4.1[,2],p5.1,xlab=list("Days",cex=1.5),ylab=list("Occupation Probability",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),col=2,lty=2)

legend(1, 0.50, c("Group 0","Group 1"), col = c(1,2), lty = c(1,2),lwd=3,  merge =TRUE, bg = 'gray90')
title ("Aalen-Johansen estimator for OP for state Transplantable")



####


#state 7###
new.model.6<-as.data.frame(summary(etm.sym)[[6]])
new.model.12<-as.data.frame(summary(etm.sym)[[12]])
new.model.18<-as.data.frame(summary(etm.sym)[[18]])
new.model.24<-as.data.frame(summary(etm.sym)[[24]])
n.all.i<-new.model.6$n.risk[1]+new.model.12$n.risk[1]+new.model.18$n.risk[1]+new.model.24$n.risk[1]
tall.i<-(n.all.1-1)

p2.i<-new.model.6$n.risk[1]/tall.i
p3.i<-new.model.12$n.risk[1]/tall.i
p4.i<-new.model.18$n.risk[1]/tall.i
p5.i<-new.model.24$n.risk[1]/tall.i

p7<-p2.i*new.model.6[,1]+p3.i*new.model.12[,1]+p4.i*new.model.18[,1]+p5.i*new.model.24[,1]

#plot(new.model.6[,2],p7,xlab=list("Days",cex=1.5),ylab=list("Occupation Probability",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1))
#lines(new.model.6[,2],p7,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=2)


#######for GROUPS# 7#############

new.model.6.0<-as.data.frame(summary(etm.sym.0)[[6]])
new.model.12.0<-as.data.frame(summary(etm.sym.0)[[12]])
new.model.18.0<-as.data.frame(summary(etm.sym.0)[[18]])
new.model.24.0<-as.data.frame(summary(etm.sym.0)[[24]])
n.all.0<-new.model.6.0$n.risk[1]+new.model.12.0$n.risk[1]+new.model.18.0$n.risk[1]+new.model.24.0$n.risk[1]
tall.i.0<-(n.all.0-1)
p2.i.0<-new.model.6.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.12.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.18.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.24.0$n.risk[1]/tall.i.0

p7.0<-p2.i.0*new.model.6.0[,1]+p3.i.0*new.model.12.0[,1]+p4.i.0*new.model.18.0[,1]+p5.i.0*new.model.24.0[,1]

new.model.6.1<-as.data.frame(summary(etm.sym.1)[[6]])
new.model.12.1<-as.data.frame(summary(etm.sym.1)[[12]])
new.model.18.1<-as.data.frame(summary(etm.sym.1)[[18]])
new.model.24.1<-as.data.frame(summary(etm.sym.1)[[24]])
n.all.1<-new.model.6.1$n.risk[1]+new.model.12.1$n.risk[1]+new.model.18.1$n.risk[1]+new.model.24.1$n.risk[1]
tall.i.1<-(n.all.1-1)
p2.i.1<-new.model.6.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.12.0$n.risk[1]/tall.i.1
p4.i.1<-new.model.18.0$n.risk[1]/tall.i.1
p5.i.1<-new.model.24.0$n.risk[1]/tall.i.1

p7.1<-p2.i.1*new.model.6.1[,1]+p3.i.1*new.model.12.1[,1]+p4.i.1*new.model.18.1[,1]+p5.i.1*new.model.24.1[,1]

plot(new.model.6.0[,2],p7.0,xlab=list("Time",cex=1.5),ylab=list("Occupation Probability",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),lty=1)
lines(new.model.6.1[,2],p7.1,xlab=list("Time",cex=1.5),type="l",lwd=3,xlim=c(0,7),ylim=c(0,1),col=2,lty=1)

legend(2, 0.80, c("Group 0","Group 1"), col = c(1,2), lty = c(1,2),lwd=3,  merge =TRUE, bg = 'gray90')
title ("Aalen-Johansen estimator for OP for state Transplanted")
#####
##state 6####

new.model.5<-as.data.frame(summary(etm.sym)[[5]])
new.model.11<-as.data.frame(summary(etm.sym)[[11]])
new.model.17<-as.data.frame(summary(etm.sym)[[17]])
new.model.23<-as.data.frame(summary(etm.sym)[[23]])

n.all<-new.model.5$n.risk[1]+new.model.11$n.risk[1]+new.model.17$n.risk[1]+new.model.23$n.risk[1]
tall.i<-(n.all-1)
p2.i<-new.model.5$n.risk[1]/tall.i
p3.i<-new.model.11$n.risk[1]/tall.i
p4.i<-new.model.17$n.risk[1]/tall.i
p5.i<-new.model.23$n.risk[1]/tall.i

p6<-p2.i*new.model.5[,1]+p3.i*new.model.11[,1]+p4.i*new.model.17[,1]+p5.i*new.model.23[,1]

#plot(new.model.5[,2],p6,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=3)

###for groups state 6####

etm.sym.0<- etm(heart.sym[heart.sym$group==0,], c("1", "2", "3","4","5","6","7"), tra, NULL,0)
etm.sym.1<- etm(heart.sym[heart.sym$group==1,], c("1", "2", "3","4","5","6","7"), tra, NULL,0)

new.model.5.0<-as.data.frame(summary(etm.sym.0)[[5]])
new.model.11.0<-as.data.frame(summary(etm.sym.0)[[11]])
new.model.17.0<-as.data.frame(summary(etm.sym.0)[[17]])
new.model.23.0<-as.data.frame(summary(etm.sym.0)[[23]])

n.all.0<-new.model.5.0$n.risk[1]+new.model.11.0$n.risk[1]+new.model.17.0$n.risk[1]+new.model.23.0$n.risk[1]
tall.i.0<-(n.all.0-1)
p2.i.0<-new.model.5.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.11.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.17.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.23.0$n.risk[1]/tall.i.0

p6.0<-p2.i.0*new.model.5.0[,1]+p3.i.0*new.model.11.0[,1]+p4.i.0*new.model.17.0[,1]+p5.i.0*new.model.23.0[,1]


new.model.5.1<-as.data.frame(summary(etm.sym.1)[[5]])
new.model.11.1<-as.data.frame(summary(etm.sym.1)[[11]])
new.model.17.1<-as.data.frame(summary(etm.sym.1)[[17]])
new.model.23.1<-as.data.frame(summary(etm.sym.1)[[23]])

n.all.1<-new.model.5.1$n.risk[1]+new.model.11.1$n.risk[1]+new.model.17.1$n.risk[1]+new.model.23.1$n.risk[1]
tall.i.1<-(n.all.1-1)
p2.i.1<-new.model.5.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.11.1$n.risk[1]/tall.i.1
p4.i.1<-new.model.17.1$n.risk[1]/tall.i.1
p5.i.1<-new.model.23.1$n.risk[1]/tall.i.1

p6.1<-p2.i.1*new.model.5.1[,1]+p3.i.1*new.model.11.1[,1]+p4.i.1*new.model.17.1[,1]+p5.i.1*new.model.23.1[,1]

#plot(new.model.5.0[,2],p6.0,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=1)
#lines(new.model.5.1[,2],p6.1,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=2)

#plot(new.model.5.0[,2],p6.0,xlab=list("Time",cex=1.5),ylab=list("Occupation Probability ",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),lty=1)
#lines(new.model.5.1[,2],p6.1,xlab=list("Time",cex=1.5),ylab=list("Occupation Probability",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),col=2,lty=2)
#legend(1, 0.80, c("Group 0","Group 1"), col = c(1,2), lty = c(1,2),lwd=3,  merge =TRUE, bg = 'gray90')
######

###state 3####

new.model.2<-as.data.frame(summary(etm.sym)[[2]])
new.model.26<-as.data.frame(summary(etm.sym)[[26]])
new.model.15<-as.data.frame(summary(etm.sym)[[15]])
new.model.21<-as.data.frame(summary(etm.sym)[[21]])
n.all<-new.model.2$n.risk[1]+new.model.26$n.risk[1]+new.model.15$n.risk[1]+new.model.21$n.risk[1]
tall.i<-(n.all-1)
p2.i<-new.model.2$n.risk[1]/tall.i
p3.i<-new.model.26$n.risk[1]/tall.i
p4.i<-new.model.15$n.risk[1]/tall.i
p5.i<-new.model.21$n.risk[1]/tall.i

p3<-p2.i*new.model.2[,1]+p3.i*new.model.26[,1]+p4.i*new.model.15[,1]+p5.i*new.model.21[,1]


#lines(new.model.2[,2],p3,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=4)

###groups state 3###

new.model.2.0<-as.data.frame(summary(etm.sym.0)[[2]])
new.model.26.0<-as.data.frame(summary(etm.sym.0)[[26]])
new.model.15.0<-as.data.frame(summary(etm.sym.0)[[15]])
new.model.21.0<-as.data.frame(summary(etm.sym.0)[[21]])
n.all.0<-new.model.2.0$n.risk[1]+new.model.26.0$n.risk[1]+new.model.15.0$n.risk[1]+new.model.21.0$n.risk[1]
tall.i.0<-(n.all.0-1)
p2.i.0<-new.model.2.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.26.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.15.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.21.0$n.risk[1]/tall.i.0

p3.0<-p2.i.0*new.model.2.0[,1]+p3.i.0*new.model.26.0[,1]+p4.i.0*new.model.15.0[,1]+p5.i.0*new.model.21.0[,1]

new.model.2.1<-as.data.frame(summary(etm.sym.1)[[2]])
new.model.26.1<-as.data.frame(summary(etm.sym.1)[[26]])
new.model.15.1<-as.data.frame(summary(etm.sym.1)[[15]])
new.model.21.1<-as.data.frame(summary(etm.sym.1)[[21]])
n.all.1<-new.model.2.1$n.risk[1]+new.model.26.1$n.risk[1]+new.model.15.1$n.risk[1]+new.model.21.1$n.risk[1]
tall.i.1<-(n.all.1-1)
p2.i.1<-new.model.2.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.26.1$n.risk[1]/tall.i.1
p4.i.1<-new.model.15.1$n.risk[1]/tall.i.1
p5.i.1<-new.model.21.1$n.risk[1]/tall.i.1

p3.1<-p2.i.1*new.model.2.1[,1]+p3.i.1*new.model.26.1[,1]+p4.i.1*new.model.15.1[,1]+p5.i.1*new.model.21.1[,1]

#plot(new.model.2.0[,2],p3.0,xlab=list("Time",cex=1.5),ylab=list("Occupation Probability",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),lty=1)

#lines(new.model.2.1[,2],p3.1,xlab=list("Time",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),col=2,lty=2)

#legend(2, 0.80, c("Group 0","Group 1"), col = c(1,2), lty = c(1,2),lwd=3,  merge =TRUE, bg = 'gray90')

####

##state 4#####

new.model.3<-as.data.frame(summary(etm.sym)[[3]])
new.model.9<-as.data.frame(summary(etm.sym)[[9]])
new.model.27<-as.data.frame(summary(etm.sym)[[27]])
new.model.22<-as.data.frame(summary(etm.sym)[[22]])
tall.i<-(n.all.i-1)
p2.i<-new.model.3$n.risk[1]/tall.i
p3.i<-new.model.9$n.risk[1]/tall.i
p4.i<-new.model.27$n.risk[1]/tall.i
p5.i<-new.model.22$n.risk[1]/tall.i


p4<-p2.i*new.model.3[,1]+p3.i*new.model.9[,1]+p4.i*new.model.27[,1]+p5.i*new.model.22[,1]

#lines(new.model.3[,2],p4,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=5)

###for groups state 4####

new.model.3.0<-as.data.frame(summary(etm.sym.0)[[3]])
new.model.9.0<-as.data.frame(summary(etm.sym.0)[[9]])
new.model.27.0<-as.data.frame(summary(etm.sym.0)[[27]])
new.model.22.0<-as.data.frame(summary(etm.sym.0)[[22]])
tall.i.0<-(n.all.0-1)
p2.i.0<-new.model.3.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.9.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.27.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.22.0$n.risk[1]/tall.i.0


p4.0<-p2.i.0*new.model.3.0[,1]+p3.i.0*new.model.9.0[,1]+p4.i.0*new.model.27.0[,1]+p5.i.1*new.model.22.0[,1]



new.model.3.1<-as.data.frame(summary(etm.sym.1)[[3]])
new.model.9.1<-as.data.frame(summary(etm.sym.1)[[9]])
new.model.27.1<-as.data.frame(summary(etm.sym.1)[[27]])
new.model.22.1<-as.data.frame(summary(etm.sym.1)[[22]])
tall.i.1<-(n.all.1-1)
p2.i.1<-new.model.3.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.9.1$n.risk[1]/tall.i.1
p4.i.1<-new.model.27.1$n.risk[1]/tall.i.1
p5.i.1<-new.model.22.1$n.risk[1]/tall.i.1


p4.1<-p2.i.1*new.model.3.1[,1]+p3.i.1*new.model.9.1[,1]+p4.i.1*new.model.27.1[,1]+p5.i.1*new.model.22.1[,1]

#plot(new.model.3.0[,2],p4.0,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=1)

#lines(new.model.3.1[,2],p4.1,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=2)

####



##state 2#####

new.model.25<-as.data.frame(summary(etm.sym)[[25]])
new.model.8<-as.data.frame(summary(etm.sym)[[8]])
new.model.14<-as.data.frame(summary(etm.sym)[[14]])
new.model.20<-as.data.frame(summary(etm.sym)[[20]])
n.all<-new.model.25$n.risk[1]+new.model.8$n.risk[1]+new.model.14$n.risk[1]+new.model.20$n.risk[1]
tall.i<-(n.all-1)
p2.i<-new.model.25$n.risk[1]/tall.i
p3.i<-new.model.8$n.risk[1]/tall.i
p4.i<-new.model.14$n.risk[1]/tall.i
p5.i<-new.model.20$n.risk[1]/tall.i

p2<-p2.i*new.model.25[,1]+p3.i*new.model.8[,1]+p4.i*new.model.14[,1]+p5.i*new.model.20[,1]

#lines(new.model.25[,2],p2,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=6)

###state 2 groups####

new.model.25.0<-as.data.frame(summary(etm.sym.0)[[25]])
new.model.8.0<-as.data.frame(summary(etm.sym.0)[[8]])
new.model.14.0<-as.data.frame(summary(etm.sym.0)[[14]])
new.model.20.0<-as.data.frame(summary(etm.sym.0)[[20]])
n.all.0<-new.model.25.0$n.risk[1]+new.model.8.0$n.risk[1]+new.model.14.0$n.risk[1]+new.model.20.0$n.risk[1]
tall.i.0<-(n.all.0-1)
p2.i.0<-new.model.25.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.8.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.14.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.20.0$n.risk[1]/tall.i.0

p2.0<-p2.i.0*new.model.25.0[,1]+p3.i.0*new.model.8.0[,1]+p4.i.0*new.model.14.0[,1]+p5.i.0*new.model.20.0[,1]


new.model.25.1<-as.data.frame(summary(etm.sym.1)[[25]])
new.model.8.1<-as.data.frame(summary(etm.sym.1)[[8]])
new.model.14.1<-as.data.frame(summary(etm.sym.1)[[14]])
new.model.20.1<-as.data.frame(summary(etm.sym.1)[[20]])
n.all.1<-new.model.25.1$n.risk[1]+new.model.8.1$n.risk[1]+new.model.14.1$n.risk[1]+new.model.20.1$n.risk[1]
tall.i.1<-(n.all.1-1)
p2.i.1<-new.model.25.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.8.1$n.risk[1]/tall.i.1
p4.i.1<-new.model.14.1$n.risk[1]/tall.i.1
p5.i.1<-new.model.20.1$n.risk[1]/tall.i.1

p2.1<-p2.i.1*new.model.25.1[,1]+p3.i.1*new.model.8.1[,1]+p4.i.1*new.model.14.1[,1]+p5.i.1*new.model.20.1[,1]

#plot(new.model.25.0[,2],p2.0,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=1)
#lines(new.model.25.1[,2],p2.1,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=6)

######state 1#####

new.model.1<-as.data.frame(summary(etm.sym)[[1]])
new.model.7<-as.data.frame(summary(etm.sym)[[7]])
new.model.13<-as.data.frame(summary(etm.sym)[[13]])
new.model.19<-as.data.frame(summary(etm.sym)[[19]])

n.all<-new.model.1$n.risk[1]+new.model.7$n.risk[1]+new.model.13$n.risk[1]+new.model.19$n.risk[1]
tall.i<-(n.all-1)
p2.i<-new.model.1$n.risk[1]/tall.i
p3.i<-new.model.7$n.risk[1]/tall.i
p4.i<-new.model.13$n.risk[1]/tall.i
p5.i<-new.model.19$n.risk[1]/tall.i

p1<-p2.i*new.model.1[,1]+p3.i*new.model.7[,1]+p4.i*new.model.13[,1]+p5.i*new.model.19[,1]

#lines(new.model.19[,2],p1,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=7)

###state 1 groups###

new.model.1.0<-as.data.frame(summary(etm.sym.0)[[1]])
new.model.7.0<-as.data.frame(summary(etm.sym.0)[[7]])
new.model.13.0<-as.data.frame(summary(etm.sym.0)[[13]])
new.model.19.0<-as.data.frame(summary(etm.sym.0)[[19]])

n.all.0<-new.model.1.0$n.risk[1]+new.model.7.0$n.risk[1]+new.model.13.0$n.risk[1]+new.model.19.0$n.risk[1]
tall.i.0<-n.all.0
p2.i.0<-new.model.1.0$n.risk[1]/tall.i.0
p3.i.0<-new.model.7.0$n.risk[1]/tall.i.0
p4.i.0<-new.model.13.0$n.risk[1]/tall.i.0
p5.i.0<-new.model.19.0$n.risk[1]/tall.i.0

p1.0<-p2.i.0*new.model.1.0[,1]+p3.i.0*new.model.7.0[,1]+p4.i.0*new.model.13.0[,1]+p5.i.0*new.model.19.0[,1]


new.model.1.1<-as.data.frame(summary(etm.sym.1)[[1]])
new.model.7.1<-as.data.frame(summary(etm.sym.1)[[7]])
new.model.13.1<-as.data.frame(summary(etm.sym.1)[[13]])
new.model.19.1<-as.data.frame(summary(etm.sym.1)[[19]])

n.all.1<-new.model.1.1$n.risk[1]+new.model.7.1$n.risk[1]+new.model.13.1$n.risk[1]+new.model.19.1$n.risk[1]
tall.i.1<-n.all.1
p2.i.1<-new.model.1.1$n.risk[1]/tall.i.1
p3.i.1<-new.model.7.1$n.risk[1]/tall.i.1
p4.i.1<-new.model.13.1$n.risk[1]/tall.i.1
p5.i.1<-new.model.19.1$n.risk[1]/tall.i.1

p1.1<-p2.i.1*new.model.1.1[,1]+p3.i.1*new.model.7.1[,1]+p4.i.1*new.model.13.1[,1]+p5.i.1*new.model.19.1[,1]

#plot(new.model.19.0[,2],p1.0,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=1)
#lines(new.model.19.1[,2],p1.1,xlab=list("Days",cex=1.5),type="l",lwd=3,xlim=c(0,10),ylim=c(0,1),col=7)


#plot(new.model.19.0[,2],p1.0,xlab=list("Time",cex=1.5),ylab=list("Occupation Probability ",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),lty=1)
#lines(new.model.19.1[,2],p1.1,xlab=list("Time",cex=1.5),ylab=list("Occupation Probability",cex=1.5),type="l",lwd=3,xlim=c(0,5),ylim=c(0,1),col=2,lty=2)
#legend(1, 0.80, c("Group 0","Group 1"), col = c(1,2), lty = c(1,2),lwd=3,  merge =TRUE, bg = 'gray90')

#######################
###chosing times for calculating pseudo-values#####
###state 5
new.model.fr.4<-as.data.frame(summary(etm.sym)[[4]])
new.model.fr.10<-as.data.frame(summary(etm.sym)[[10]])
new.model.fr.16<-as.data.frame(summary(etm.sym)[[16]])
new.model.fr.28<-as.data.frame(summary(etm.sym)[[28]])
###state7
new.model.fr.6<-as.data.frame(summary(etm.sym)[[6]])
new.model.fr.12<-as.data.frame(summary(etm.sym)[[12]])
new.model.fr.18<-as.data.frame(summary(etm.sym)[[18]])
new.model.fr.24<-as.data.frame(summary(etm.sym)[[24]])
###state 6
new.model.fr.5<-as.data.frame(summary(etm.sym)[[5]])
new.model.fr.11<-as.data.frame(summary(etm.sym)[[11]])
new.model.fr.17<-as.data.frame(summary(etm.sym)[[17]])
new.model.fr.23<-as.data.frame(summary(etm.sym)[[23]])
####state 1
new.model.fr.1<-as.data.frame(summary(etm.sym)[[1]])
new.model.fr.7<-as.data.frame(summary(etm.sym)[[7]])
new.model.fr.13<-as.data.frame(summary(etm.sym)[[13]])
new.model.fr.19<-as.data.frame(summary(etm.sym)[[19]])
####state 4###
new.model.fr.3<-as.data.frame(summary(etm.sym)[[3]])
new.model.fr.9<-as.data.frame(summary(etm.sym)[[9]])
new.model.fr.22<-as.data.frame(summary(etm.sym)[[22]])
new.model.fr.27<-as.data.frame(summary(etm.sym)[[27]])
###state 3####
new.model.fr.2<-as.data.frame(summary(etm.sym)[[2]])
new.model.fr.15<-as.data.frame(summary(etm.sym)[[15]])
new.model.fr.21<-as.data.frame(summary(etm.sym)[[21]])
new.model.fr.26<-as.data.frame(summary(etm.sym)[[26]])
###state 2###
new.model.fr.8<-as.data.frame(summary(etm.sym)[[8]])
new.model.fr.14<-as.data.frame(summary(etm.sym)[[14]])
new.model.fr.20<-as.data.frame(summary(etm.sym)[[20]])
new.model.fr.25<-as.data.frame(summary(etm.sym)[[25]])
####time quantiles###
time.quantiles<-quantile(new.model.fr.8$time,probs = seq(0, 1, 0.1))

quantiles.states<-matrix(0,1,11)


for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.states[i]<-which(abs(new.model.fr.7$time-your.number)==min(abs(new.model.fr.7$time-your.number)))[1]

}


##occupation prob 5 all individuals#######
n.all<-new.model.fr.4$n.risk[1]+new.model.fr.10$n.risk[1]+new.model.fr.16$n.risk[1]+new.model.fr.28$n.risk[1]
tall.i<-n.all
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.5<-p2.i*new.model.fr.4[quantiles.states,1]+p3.i*new.model.fr.10[quantiles.states,1]+p4.i*new.model.fr.16[quantiles.states,1]+p5.i*new.model.fr.28[quantiles.states,1]
t.together.5<-rep(t.5,n)
#### state 7####
##occupation prob 7#######
n.all.i<-new.model.fr.6$n.risk[1]+new.model.fr.12$n.risk[1]+new.model.fr.18$n.risk[1]+new.model.fr.24$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.7<-p2.i*new.model.fr.6[quantiles.states,1]+p3.i*new.model.fr.12[quantiles.states,1]+p4.i*new.model.fr.18[quantiles.states,1]+p5.i*new.model.fr.24[quantiles.states,1]
t.together.7<-rep(t.7,n)
#####state 6###############
n.all.i<-new.model.fr.5$n.risk[1]+new.model.fr.11$n.risk[1]+new.model.fr.17$n.risk[1]+new.model.fr.23$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.6<-p2.i*new.model.fr.5[quantiles.states,1]+p3.i*new.model.fr.11[quantiles.states,1]+p4.i*new.model.fr.17[quantiles.states,1]+p5.i*new.model.fr.23[quantiles.states,1]

t.together.6<-rep(t.6,n)
###state 1####
n.all.i<-new.model.fr.1$n.risk[1]+new.model.fr.7$n.risk[1]+new.model.fr.13$n.risk[1]+new.model.fr.19$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.1<-p2.i*new.model.fr.1[quantiles.states,1]+p3.i*new.model.fr.7[quantiles.states,1]+p4.i*new.model.fr.13[quantiles.states,1]+p5.i*new.model.fr.19[quantiles.states,1]

t.together.1<-rep(t.1,n)

###state 2####
n.all.i<-new.model.fr.8$n.risk[1]+new.model.fr.14$n.risk[1]+new.model.fr.20$n.risk[1]+new.model.fr.25$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.2<-p2.i*new.model.fr.8[quantiles.states,1]+p3.i*new.model.fr.14[quantiles.states,1]+p4.i*new.model.fr.20[quantiles.states,1]+p5.i*new.model.fr.25[quantiles.states,1]

t.together.2<-rep(t.2,n)
### state 3####
n.all.i<-new.model.fr.2$n.risk[1]+new.model.fr.15$n.risk[1]+new.model.fr.21$n.risk[1]+new.model.fr.26$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.3<-p2.i*new.model.fr.2[quantiles.states,1]+p3.i*new.model.fr.15[quantiles.states,1]+p4.i*new.model.fr.21[quantiles.states,1]+p5.i*new.model.fr.26[quantiles.states,1]

t.together.3<-rep(t.3,n)

###state 4###

n.all.i<-new.model.fr.3$n.risk[1]+new.model.fr.9$n.risk[1]+new.model.fr.22$n.risk[1]+new.model.fr.27$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


t.4<-p2.i*new.model.fr.3[quantiles.states,1]+p3.i*new.model.fr.9[quantiles.states,1]+p4.i*new.model.fr.22[quantiles.states,1]+p5.i*new.model.fr.27[quantiles.states,1]

t.together.4<-rep(t.4,n)


########################
n.times<-11 ### ##

id.set<-unique(heart.sym$id)
###############################
pseudo.data.55<-matrix(0,n.times*n,4)
ps.time<-c(new.model.fr.4[quantiles.states,2])
pseudo.data.55[,1]<-sort(rep(id.set,n.times))
pseudo.data.55[,2]<-rep(ps.time,n)
####
pseudo.data.77<-matrix(0,n.times*n,4)
pseudo.data.77[,1]<-sort(rep(id.set,n.times))
pseudo.data.77[,2]<-rep(ps.time,n)
####
pseudo.data.66<-matrix(0,n.times*n,4)
pseudo.data.66[,1]<-sort(rep(id.set,n.times))
pseudo.data.66[,2]<-rep(ps.time,n)
###
pseudo.data.11<-matrix(0,n.times*n,4)
pseudo.data.11[,1]<-sort(rep(id.set,n.times))
pseudo.data.11[,2]<-rep(ps.time,n)
####
pseudo.data.22<-matrix(0,n.times*n,4)
pseudo.data.22[,1]<-sort(rep(id.set,n.times))
pseudo.data.22[,2]<-rep(ps.time,n)
####
pseudo.data.33<-matrix(0,n.times*n,4)
pseudo.data.33[,1]<-sort(rep(id.set,n.times))
pseudo.data.33[,2]<-rep(ps.time,n)
###

pseudo.data.44<-matrix(0,n.times*n,4)
pseudo.data.44[,1]<-sort(rep(id.set,n.times))
pseudo.data.44[,2]<-rep(ps.time,n)
###

for (i in 1: n){
id.to.exclude<-i
cat("i=",i)

new.data<-heart.sym[heart.sym$id!=id.to.exclude,]
new.model <- etm(new.data, c("1", "2", "3","4","5","6","7"), tra, NULL, 0)##we start from 0##

new.model.fr.4<-as.data.frame(summary(new.model)[[4]])
new.model.fr.10<-as.data.frame(summary(new.model)[[10]])
new.model.fr.16<-as.data.frame(summary(new.model)[[16]])
new.model.fr.28<-as.data.frame(summary(new.model)[[28]])

new.model.fr.6<-as.data.frame(summary(new.model)[[6]])
new.model.fr.12<-as.data.frame(summary(new.model)[[12]])
new.model.fr.18<-as.data.frame(summary(new.model)[[18]])
new.model.fr.24<-as.data.frame(summary(new.model)[[24]])
###state 6
new.model.fr.5<-as.data.frame(summary(new.model)[[5]])
new.model.fr.11<-as.data.frame(summary(new.model)[[11]])
new.model.fr.17<-as.data.frame(summary(new.model)[[17]])
new.model.fr.23<-as.data.frame(summary(new.model)[[23]])
####state 1
new.model.fr.1<-as.data.frame(summary(new.model)[[1]])
new.model.fr.7<-as.data.frame(summary(new.model)[[7]])
new.model.fr.13<-as.data.frame(summary(new.model)[[13]])
new.model.fr.19<-as.data.frame(summary(new.model)[[19]])
####state 4###
new.model.fr.3<-as.data.frame(summary(new.model)[[3]])
new.model.fr.9<-as.data.frame(summary(new.model)[[9]])
new.model.fr.22<-as.data.frame(summary(new.model)[[22]])
new.model.fr.27<-as.data.frame(summary(new.model)[[27]])
###state 3####
new.model.fr.2<-as.data.frame(summary(new.model)[[2]])
new.model.fr.15<-as.data.frame(summary(new.model)[[15]])
new.model.fr.21<-as.data.frame(summary(new.model)[[21]])
new.model.fr.26<-as.data.frame(summary(new.model)[[26]])
###state 2###
new.model.fr.8<-as.data.frame(summary(new.model)[[8]])
new.model.fr.14<-as.data.frame(summary(new.model)[[14]])
new.model.fr.20<-as.data.frame(summary(new.model)[[20]])
new.model.fr.25<-as.data.frame(summary(new.model)[[25]])

####

quantiles.state.rem<-rep(NA,11)

for( k in 1:11){
your.number<-time.quantiles[k]
quantiles.state.rem[k]<-which(abs(new.model.fr.4$time-your.number)==min(abs(new.model.fr.4$time-your.number)))[1]
}
###


n.all.i<-new.model.fr.4$n.risk[1]+new.model.fr.10$n.risk[1]+new.model.fr.16$n.risk[1]+new.model.fr.28$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.4$n.risk[1]/tall.i
p3.i<-new.model.fr.10$n.risk[1]/tall.i
p4.i<-new.model.fr.16$n.risk[1]/tall.i
p5.i<-new.model.fr.28$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.4$n.risk[1]
p3.i<-new.model.fr.10$n.risk[1]
p4.i<-new.model.fr.16$n.risk[1]
p5.i<-new.model.fr.28$n.risk[1] }



pseudo.data.55[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.4[quantiles.state.rem,1]+p3.i*new.model.fr.10[quantiles.state.rem,1]+
p4.i*new.model.fr.16[quantiles.state.rem,1]+p5.i*new.model.fr.28[quantiles.state.rem,1]


n.all.i<-new.model.fr.6$n.risk[1]+new.model.fr.12$n.risk[1]+new.model.fr.18$n.risk[1]+new.model.fr.24$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.6$n.risk[1]/tall.i
p3.i<-new.model.fr.12$n.risk[1]/tall.i
p4.i<-new.model.fr.18$n.risk[1]/tall.i
p5.i<-new.model.fr.24$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.6$n.risk[1]
p3.i<-new.model.fr.12$n.risk[1]
p4.i<-new.model.fr.18$n.risk[1]
p5.i<-new.model.fr.24$n.risk[1] }


pseudo.data.77[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.6[quantiles.state.rem,1]+p3.i*new.model.fr.12[quantiles.state.rem,1]+
p4.i*new.model.fr.18[quantiles.state.rem,1]+p5.i*new.model.fr.24[quantiles.state.rem,1]


n.all.i<-new.model.fr.5$n.risk[1]+new.model.fr.11$n.risk[1]+new.model.fr.17$n.risk[1]+new.model.fr.23$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.5$n.risk[1]/tall.i
p3.i<-new.model.fr.11$n.risk[1]/tall.i
p4.i<-new.model.fr.17$n.risk[1]/tall.i
p5.i<-new.model.fr.23$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.5$n.risk[1]
p3.i<-new.model.fr.11$n.risk[1]
p4.i<-new.model.fr.17$n.risk[1]
p5.i<-new.model.fr.23$n.risk[1] }


pseudo.data.66[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.5[quantiles.state.rem,1]+p3.i*new.model.fr.11[quantiles.state.rem,1]+
p4.i*new.model.fr.17[quantiles.state.rem,1]+p5.i*new.model.fr.23[quantiles.state.rem,1]

n.all.i<-new.model.fr.1$n.risk[1]+new.model.fr.7$n.risk[1]+new.model.fr.13$n.risk[1]+new.model.fr.19$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.1$n.risk[1]/tall.i
p3.i<-new.model.fr.7$n.risk[1]/tall.i
p4.i<-new.model.fr.13$n.risk[1]/tall.i
p5.i<-new.model.fr.19$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.1$n.risk[1]
p3.i<-new.model.fr.7$n.risk[1]
p4.i<-new.model.fr.13$n.risk[1]
p5.i<-new.model.fr.19$n.risk[1] }


pseudo.data.11[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.1[quantiles.state.rem,1]+p3.i*new.model.fr.7[quantiles.state.rem,1]+
p4.i*new.model.fr.13[quantiles.state.rem,1]+p5.i*new.model.fr.19[quantiles.state.rem,1]


n.all.i<-new.model.fr.8$n.risk[1]+new.model.fr.14$n.risk[1]+new.model.fr.20$n.risk[1]+new.model.fr.25$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.8$n.risk[1]/tall.i
p3.i<-new.model.fr.14$n.risk[1]/tall.i
p4.i<-new.model.fr.20$n.risk[1]/tall.i
p5.i<-new.model.fr.25$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.8$n.risk[1]
p3.i<-new.model.fr.14$n.risk[1]
p4.i<-new.model.fr.20$n.risk[1]
p5.i<-new.model.fr.25$n.risk[1] }


pseudo.data.22[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.8[quantiles.state.rem,1]+p3.i*new.model.fr.14[quantiles.state.rem,1]+
p4.i*new.model.fr.20[quantiles.state.rem,1]+p5.i*new.model.fr.25[quantiles.state.rem,1]

n.all.i<-new.model.fr.2$n.risk[1]+new.model.fr.15$n.risk[1]+new.model.fr.21$n.risk[1]+new.model.fr.26$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.2$n.risk[1]/tall.i
p3.i<-new.model.fr.15$n.risk[1]/tall.i
p4.i<-new.model.fr.21$n.risk[1]/tall.i
p5.i<-new.model.fr.26$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.2$n.risk[1]
p3.i<-new.model.fr.15$n.risk[1]
p4.i<-new.model.fr.21$n.risk[1]
p5.i<-new.model.fr.26$n.risk[1] }


pseudo.data.33[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.2[quantiles.state.rem,1]+p3.i*new.model.fr.15[quantiles.state.rem,1]+
p4.i*new.model.fr.21[quantiles.state.rem,1]+p5.i*new.model.fr.26[quantiles.state.rem,1]

n.all.i<-new.model.fr.3$n.risk[1]+new.model.fr.9$n.risk[1]+new.model.fr.22$n.risk[1]+new.model.fr.27$n.risk[1]
tall.i<-n.all.i
if(n.all.i>0){
p2.i<-new.model.fr.3$n.risk[1]/tall.i
p3.i<-new.model.fr.9$n.risk[1]/tall.i
p4.i<-new.model.fr.22$n.risk[1]/tall.i
p5.i<-new.model.fr.27$n.risk[1]/tall.i
} else{

p2.i<-new.model.fr.3$n.risk[1]
p3.i<-new.model.fr.9$n.risk[1]
p4.i<-new.model.fr.22$n.risk[1]
p5.i<-new.model.fr.27$n.risk[1] }


pseudo.data.44[(11*(i-1)+1):(11*(i-1)+11),3]<-p2.i*new.model.fr.3[quantiles.state.rem,1]+p3.i*new.model.fr.9[quantiles.state.rem,1]+
p4.i*new.model.fr.22[quantiles.state.rem,1]+p5.i*new.model.fr.27[quantiles.state.rem,1]

}



ps.time<-time.quantiles

pseudo.data.55[,2]<-rep(ps.time,n)
pseudo.data.55[,4]<-t.together.5
pseudo.fr.5<-as.data.frame(pseudo.data.55)
names(pseudo.fr.5)<-c("id","time","pseudoi","pseudo")

pseudo.data.77[,2]<-rep(ps.time,n)
pseudo.data.77[,4]<-t.together.7
pseudo.fr.7<-as.data.frame(pseudo.data.77)
names(pseudo.fr.7)<-c("id","time","pseudoi","pseudo")

pseudo.data.66[,2]<-rep(ps.time,n)
pseudo.data.66[,4]<-t.together.6
pseudo.fr.6<-as.data.frame(pseudo.data.66)
names(pseudo.fr.6)<-c("id","time","pseudoi","pseudo")

pseudo.data.11[,2]<-rep(ps.time,n)
pseudo.data.11[,4]<-t.together.1
pseudo.fr.1<-as.data.frame(pseudo.data.11)
names(pseudo.fr.1)<-c("id","time","pseudoi","pseudo")

pseudo.data.22[,2]<-rep(ps.time,n)
pseudo.data.22[,4]<-t.together.2
pseudo.fr.2<-as.data.frame(pseudo.data.22)
names(pseudo.fr.2)<-c("id","time","pseudoi","pseudo")

pseudo.data.33[,2]<-rep(ps.time,n)
pseudo.data.33[,4]<-t.together.3
pseudo.fr.3<-as.data.frame(pseudo.data.33)
names(pseudo.fr.3)<-c("id","time","pseudoi","pseudo")

pseudo.data.44[,2]<-rep(ps.time,n)
pseudo.data.44[,4]<-t.together.4
pseudo.fr.4<-as.data.frame(pseudo.data.44)
names(pseudo.fr.4)<-c("id","time","pseudoi","pseudo")

###calculate pseudo-values###

pseudoja<-pseudo.fr.5$pseudo*n-pseudo.fr.5$pseudoi*(n-1)
pseudo.fr.5$pseudoja<-pseudoja

pseudoja<-pseudo.fr.7$pseudo*n-pseudo.fr.7$pseudoi*(n-1)
pseudo.fr.7$pseudoja<-pseudoja

pseudoja<-pseudo.fr.6$pseudo*n-pseudo.fr.6$pseudoi*(n-1)
pseudo.fr.6$pseudoja<-pseudoja

pseudoja<-pseudo.fr.1$pseudo*n-pseudo.fr.1$pseudoi*(n-1)
pseudo.fr.1$pseudoja<-pseudoja

pseudoja<-pseudo.fr.2$pseudo*n-pseudo.fr.2$pseudoi*(n-1)
pseudo.fr.2$pseudoja<-pseudoja

pseudoja<-pseudo.fr.3$pseudo*n-pseudo.fr.3$pseudoi*(n-1)
pseudo.fr.3$pseudoja<-pseudoja

pseudoja<-pseudo.fr.4$pseudo*n-pseudo.fr.4$pseudoi*(n-1)
pseudo.fr.4$pseudoja<-pseudoja


##adding group###

group.ps<-matrix(0,length(pseudo.fr.5$id),1)
for ( i in 1: length(pseudo.fr.5$id)){
#cat("i=",i)
id.current<-pseudo.fr.5$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.5$group<-group.ps

#####
group.ps<-matrix(0,length(pseudo.fr.7$id),1)
for ( i in 1: length(pseudo.fr.7$id)){
#cat("i=",i)
id.current<-pseudo.fr.7$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.7$group<-group.ps

#####

group.ps<-matrix(0,length(pseudo.fr.6$id),1)
for ( i in 1: length(pseudo.fr.6$id)){
#cat("i=",i)
id.current<-pseudo.fr.6$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.6$group<-group.ps

###

group.ps<-matrix(0,length(pseudo.fr.1$id),1)
for ( i in 1: length(pseudo.fr.1$id)){
#cat("i=",i)
id.current<-pseudo.fr.1$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.1$group<-group.ps

###
group.ps<-matrix(0,length(pseudo.fr.2$id),1)
for ( i in 1: length(pseudo.fr.2$id)){
#cat("i=",i)
id.current<-pseudo.fr.2$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.2$group<-group.ps
###

group.ps<-matrix(0,length(pseudo.fr.3$id),1)
for ( i in 1: length(pseudo.fr.3$id)){
#cat("i=",i)
id.current<-pseudo.fr.3$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.3$group<-group.ps
###
group.ps<-matrix(0,length(pseudo.fr.4$id),1)
for ( i in 1: length(pseudo.fr.4$id)){
#cat("i=",i)
id.current<-pseudo.fr.4$id[i]
group.ps[i]<-heart.sym$group[heart.sym$id==id.current][1]
}

pseudo.fr.4$group<-group.ps




###fitting GEE

##3univariate GEE models###

library(geepack)


########

##true ones in chosen time points###

##state 5##

quantiles.true<-matrix(0,1,11)


for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.4.0$time-your.number)==min(abs(new.model.4.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.4.1$time-your.number)==min(abs(new.model.4.1$time-your.number)))[1]

}


p5.0.c<-p2.i.0*new.model.4.0[quantiles.true,1]+p3.i.0*new.model.10.0[quantiles.true,1]+p4.i.0*new.model.16.0[quantiles.true,1]+p5.i.0*new.model.28.0[quantiles.true,1]
p5.1.c<-p2.i.1*new.model.4.1[quantiles.true.2,1]+p3.i.1*new.model.10.1[quantiles.true.2,1]+p4.i.1*new.model.16.1[quantiles.true.2,1]+p5.i.1*new.model.28.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p5.0.c)),rep(1,length(p5.1.c)))
true.ones<- c(p5.0.c,p5.1.c)
true.data.5<-as.data.frame(cbind(group.true,true.ones))
true.data.5$time<-rep(time.quantiles,2)

fitgee.5<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.5,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")
summary(fitgee.5)

true.5<-lm(true.data.5$true.ones~true.data.5$time+true.data.5$group)

bias.5[iter]<-summary(fitgee.5)[[1]]$estimate[3]-summary(true.5)$coefficients[3]
stand.se.5[iter]<-summary(fitgee.5)[[1]]$san.se[3]
jack.se.5[iter]<-summary(fitgee.5)[[1]]$ajs.se[3]
stand.se.5[iter]<-summary(fitgee.5)[[1]]$san.se[3]
stand.cov.5[iter]<-as.numeric(summary(fitgee.5)[[1]]$p[3]<0.05)

######state 7  
quantiles.true<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.6.0$time-your.number)==min(abs(new.model.6.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.6.1$time-your.number)==min(abs(new.model.6.1$time-your.number)))[1]

}


p7.0.c<-p2.i.0*new.model.6.0[quantiles.true,1]+p3.i.0*new.model.12.0[quantiles.true,1]+p4.i.0*new.model.18.0[quantiles.true,1]+p5.i.0*new.model.24.0[quantiles.true,1]
p7.1.c<-p2.i.1*new.model.6.1[quantiles.true.2,1]+p3.i.1*new.model.12.1[quantiles.true.2,1]+p4.i.1*new.model.18.1[quantiles.true.2,1]+p5.i.1*new.model.24.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p7.0.c)),rep(1,length(p7.1.c)))
true.ones<- c(p7.0.c,p7.1.c)
true.data.7<-as.data.frame(cbind(group.true,true.ones))
true.data.7$time<-rep(time.quantiles,2)

fitgee.7<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.7,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")
summary(fitgee.7)
true.7<-lm(true.data.7$true.ones~true.data.7$time+true.data.7$group)

bias.7[iter]<-summary(fitgee.7)[[1]]$estimate[3]-summary(true.7)$coefficients[3]
stand.se.7[iter]<-summary(fitgee.7)[[1]]$san.se[3]
jack.se.7[iter]<-summary(fitgee.7)[[1]]$ajs.se[3]
stand.se.7[iter]<-summary(fitgee.7)[[1]]$san.se[3]
stand.cov.7[iter]<-as.numeric(summary(fitgee.7)[[1]]$p[3]<0.05)


####state 6

quantiles.true<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.5.0$time-your.number)==min(abs(new.model.5.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.5.1$time-your.number)==min(abs(new.model.5.1$time-your.number)))[1]

}


p6.0.c<-p2.i.0*new.model.5.0[quantiles.true,1]+p3.i.0*new.model.11.0[quantiles.true,1]+p4.i.0*new.model.17.0[quantiles.true,1]+p5.i.0*new.model.23.0[quantiles.true,1]
p6.1.c<-p2.i.1*new.model.5.1[quantiles.true.2,1]+p3.i.1*new.model.11.1[quantiles.true.2,1]+p4.i.1*new.model.17.1[quantiles.true.2,1]+p5.i.1*new.model.23.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p6.0.c)),rep(1,length(p6.1.c)))
true.ones<- c(p6.0.c,p6.1.c)
true.data.6<-as.data.frame(cbind(group.true,true.ones))
true.data.6$time<-rep(time.quantiles,2)

fitgee.6<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.6,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")
summary(fitgee.6)
true.6<-lm(true.data.6$true.ones~true.data.6$time+true.data.6$group)

bias.6[iter]<-summary(fitgee.6)[[1]]$estimate[3]-summary(true.6)$coefficients[3]
stand.se.6[iter]<-summary(fitgee.6)[[1]]$san.se[3]
jack.se.6[iter]<-summary(fitgee.6)[[1]]$ajs.se[3]
stand.se.6[iter]<-summary(fitgee.6)[[1]]$san.se[3]
stand.cov.6[iter]<-as.numeric(summary(fitgee.6)[[1]]$p[3]<0.05)

####### state  1

quantiles.true<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.1.0$time-your.number)==min(abs(new.model.1.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.1.1$time-your.number)==min(abs(new.model.1.1$time-your.number)))[1]

}


p1.0.c<-p2.i.0*new.model.1.0[quantiles.true,1]+p3.i.0*new.model.7.0[quantiles.true,1]+p4.i.0*new.model.13.0[quantiles.true,1]+p5.i.0*new.model.19.0[quantiles.true,1]
p1.1.c<-p2.i.1*new.model.1.1[quantiles.true.2,1]+p3.i.1*new.model.7.1[quantiles.true.2,1]+p4.i.1*new.model.13.1[quantiles.true.2,1]+p5.i.1*new.model.19.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p1.0.c)),rep(1,length(p1.1.c)))
true.ones<- c(p1.0.c,p1.1.c)
true.data.1<-as.data.frame(cbind(group.true,true.ones))
true.data.1$time<-rep(time.quantiles,2)

fitgee.1<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.1,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")
summary(fitgee.1)
true.1<-lm(true.data.1$true.ones~true.data.1$time+true.data.1$group)

bias.1[iter]<-summary(fitgee.1)[[1]]$estimate[3]-summary(true.1)$coefficients[3]
stand.se.1[iter]<-summary(fitgee.1)[[1]]$san.se[3]
jack.se.1[iter]<-summary(fitgee.1)[[1]]$ajs.se[3]
stand.se.1[iter]<-summary(fitgee.1)[[1]]$san.se[3]
stand.cov.1[iter]<-as.numeric(summary(fitgee.1)[[1]]$p[3]<0.05)


####state 4###

quantiles.true<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.3.0$time-your.number)==min(abs(new.model.3.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.3.1$time-your.number)==min(abs(new.model.3.1$time-your.number)))[1]

}


p4.0.c<-p2.i.0*new.model.3.0[quantiles.true,1]+p3.i.0*new.model.9.0[quantiles.true,1]+p4.i.0*new.model.22.0[quantiles.true,1]+p5.i.0*new.model.27.0[quantiles.true,1]
p4.1.c<-p2.i.1*new.model.3.1[quantiles.true.2,1]+p3.i.1*new.model.9.1[quantiles.true.2,1]+p4.i.1*new.model.22.1[quantiles.true.2,1]+p5.i.1*new.model.27.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p4.0.c)),rep(1,length(p4.1.c)))
true.ones<- c(p4.0.c,p4.1.c)
true.data.4<-as.data.frame(cbind(group.true,true.ones))
true.data.4$time<-rep(time.quantiles,2)

fitgee.4<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.4,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")
summary(fitgee.4)
true.4<-lm(true.data.4$true.ones~true.data.4$time+true.data.4$group)

bias.4[iter]<-summary(fitgee.4)[[1]]$estimate[3]-summary(true.4)$coefficients[3]
stand.se.1[iter]<-summary(fitgee.4)[[1]]$san.se[3]
jack.se.1[iter]<-summary(fitgee.4)[[1]]$ajs.se[3]
stand.se.1[iter]<-summary(fitgee.4)[[1]]$san.se[3]
stand.cov.1[iter]<-as.numeric(summary(fitgee.4)[[1]]$p[3]<0.05)

###state 3####

quantiles.true<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.2.0$time-your.number)==min(abs(new.model.2.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.2.1$time-your.number)==min(abs(new.model.2.1$time-your.number)))[1]

}


p3.0.c<-p2.i.0*new.model.2.0[quantiles.true,1]+p3.i.0*new.model.15.0[quantiles.true,1]+p4.i.0*new.model.21.0[quantiles.true,1]+p5.i.0*new.model.26.0[quantiles.true,1]
p3.1.c<-p2.i.1*new.model.2.1[quantiles.true.2,1]+p3.i.1*new.model.15.1[quantiles.true.2,1]+p4.i.1*new.model.21.1[quantiles.true.2,1]+p5.i.1*new.model.26.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p3.0.c)),rep(1,length(p3.1.c)))
true.ones<- c(p3.0.c,p3.1.c)
true.data.3<-as.data.frame(cbind(group.true,true.ones))
true.data.3$time<-rep(time.quantiles,2)

fitgee.3<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.3,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")
summary(fitgee.3)

true.3<-lm(true.data.4$true.ones~true.data.3$time+true.data.3$group)

bias.3[iter]<-summary(fitgee.3)[[1]]$estimate[3]-summary(true.3)$coefficients[3]
stand.se.3[iter]<-summary(fitgee.3)[[1]]$san.se[3]
jack.se.3[iter]<-summary(fitgee.3)[[1]]$ajs.se[3]
stand.se.3[iter]<-summary(fitgee.3)[[1]]$san.se[3]
stand.cov.3[iter]<-as.numeric(summary(fitgee.3)[[1]]$p[3]<0.05)



###state 2###


quantiles.true<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true[i]<-which(abs(new.model.8.0$time-your.number)==min(abs(new.model.8.0$time-your.number)))[1]

}

quantiles.true.2<-matrix(0,1,11)

for (i in 1:11){

your.number<-time.quantiles[i]
quantiles.true.2[i]<-which(abs(new.model.8.1$time-your.number)==min(abs(new.model.8.1$time-your.number)))[1]

}


p2.0.c<-p2.i.0*new.model.8.0[quantiles.true,1]+p3.i.0*new.model.14.0[quantiles.true,1]+p4.i.0*new.model.20.0[quantiles.true,1]+p5.i.0*new.model.25.0[quantiles.true,1]
p2.1.c<-p2.i.1*new.model.8.1[quantiles.true.2,1]+p3.i.1*new.model.14.1[quantiles.true.2,1]+p4.i.1*new.model.20.1[quantiles.true.2,1]+p5.i.1*new.model.25.1[quantiles.true.2,1]

group.true<-c(rep(0,length(p2.0.c)),rep(1,length(p2.1.c)))
true.ones<- c(p2.0.c,p2.1.c)
true.data.2<-as.data.frame(cbind(group.true,true.ones))
true.data.2$time<-rep(time.quantiles,2)

fitgee.2<-geese(pseudoja~time+as.factor(group),data=pseudo.fr.2,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")

summary(fitgee.2)
true.2<-lm(true.data.4$true.ones~true.data.2$time+true.data.2$group)

bias.2[iter]<-summary(fitgee.2)[[1]]$estimate[3]-summary(true.2)$coefficients[3]
stand.se.2[iter]<-summary(fitgee.2)[[1]]$san.se[3]
jack.se.2[iter]<-summary(fitgee.2)[[1]]$ajs.se[3]
stand.se.2[iter]<-summary(fitgee.2)[[1]]$san.se[3]
stand.cov.2[iter]<-as.numeric(summary(fitgee.2)[[1]]$p[3]<0.05)


 ####multivariate pseudo-values####

  pseudo.fr.5$respind<-5
  pseudo.fr.4$respind<-4
  pseudo.fr.3$respind<-3
  pseudo.fr.2$respind<-2
  pseudo.fr.1$respind<-1
  pseudo.fr.6$respind<-6
  pseudo.fr.7$respind<-7 

  pseudo.all<-rbind(pseudo.fr.1,pseudo.fr.2,pseudo.fr.3,pseudo.fr.4,pseudo.fr.5,pseudo.fr.6,pseudo.fr.7)

  fitgee.mult<-geese(pseudoja~as.factor(respind)*(time+as.factor(group)),data=pseudo.all,id=id,scale.fix=TRUE,family=gaussian,jack=TRUE,mean.link="identity",corstr="unstructured")

  summary(fitgee.mult)

    bias.1.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[9]-summary(true.1)$coefficients[3]
    bias.2.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[16]-summary(true.2)$coefficients[3]
    bias.3.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[17]-summary(true.3)$coefficients[3]
    bias.4.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[18]-summary(true.4)$coefficients[3]
    bias.5.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[19]-summary(true.5)$coefficients[3]
    bias.6.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[20]-summary(true.6)$coefficients[3]
    bias.7.mult[iter]<-summary(fitgee.mult)[[1]]$estimate[21]-summary(true.7)$coefficients[3]

        stand.se.1.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[9]
        jack.se.1.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[9]
        stand.cov.1.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[9]<0.05)

        stand.se.2.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[16]
        jack.se.2.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[16]
        stand.cov.2.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[16]<0.05)

        stand.se.3.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[17]
        jack.se.3.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[17]
        stand.cov.3.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[17]<0.05)

        stand.se.4.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[18]
        jack.se.4.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[18]
        stand.cov.4.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[18]<0.05)

        stand.se.5.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[19]
        jack.se.5.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[19]
        stand.cov.5.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[19]<0.05)

        stand.se.6.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[20]
        jack.se.6.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[20]
        stand.cov.6.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[20]<0.05)

        stand.se.7.mult[iter]<-summary(fitgee.mult)[[1]]$san.se[21]
        jack.se.7.mult[iter]<-summary(fitgee.mult)[[1]]$ajs.se[21]
        stand.cov.7.mult[iter]<-as.numeric(summary(fitgee.mult)[[1]]$p[21]<0.05)
 }
iter= 1

plot of chunk unnamed-chunk-1 plot of chunk unnamed-chunk-1

i= 1i= 2i= 3...
For example for state Transplanted we have the following univariate GEE model:
 
summary(fitgee.7)[1]

$mean
                     estimate      san.se      ajs.se         wald            p
(Intercept)       11.57883090 0.008515335 0.008032129 1.848953e+06 0.0000000000
time               6.02980522 0.004670620 0.004395196 1.666700e+06 0.0000000000
as.factor(group)1  0.05214013 0.014346656 0.013527635 1.320817e+01 0.0002787307
         bias.uni<-list(bias.1=bias.1,bias.2=bias.2,bias.3=bias.3,bias.4=bias.4,bias.5=bias.5,bias.6=bias.6,bias.7=bias.7)
         stand.se.uni<-list(stand.se.1=stand.se.1,stand.se.2=stand.se.2,
         stand.se.3=stand.se.3,stand.se.4=stand.se.4,stand.se.5=stand.se.5,stand.se.6=stand.se.6,stand.se.7=stand.se.7)
         jack.se.uni<-list(jack.se.1=jack.se.1,jack.se.2=jack.se.2,jack.se.3=jack.se.3,
         jack.se.4=jack.se.4,jack.se.5=jack.se.5,jack.se.6=jack.se.6,jack.se.7=jack.se.7)
         cov.uni<-list(stand.cov.1=stand.cov.1,stand.cov.2=stand.cov.2,stand.cov.3=stand.cov.3,stand.cov.4=stand.cov.4,
         stand.cov.5=stand.cov.5,stand.cov.6=stand.cov.6,stand.cov.7=stand.cov.7)

        bias.mult<-list(bias.1.mult=bias.1.mult,bias.2.mult=bias.2.mult,bias.3.mult=bias.3.mult,bias.4.mult=bias.4.mult,
        bias.5.mult=bias.5.mult,bias.6.mult=bias.6.mult,bias.7.mult=bias.7.mult)
        stand.se.mult<-list(stand.se.1.mult=stand.se.1.mult,stand.se.2.mult=stand.se.2.mult,stand.se.3.mult=stand.se.3.mult,
        stand.se.4.mult=stand.se.4.mult,stand.se.5.mult=stand.se.5.mult,stand.se.6.mult=stand.se.6.mult,stand.se.7.mult=stand.se.7.mult)
        jack.se.mult<-list(jack.se.1.mult=jack.se.1.mult,jack.se.2.mult=jack.se.2.mult,jack.se.3.mult=jack.se.3.mult,
        jack.se.4.mult=jack.se.4.mult,jack.se.5.mult=jack.se.5.mult,jack.se.6.mult=jack.se.6.mult,jack.se.7.mult=jack.se.7.mult)
        cov.mult<-list(stand.cov.1.mult=stand.cov.1.mult,stand.cov.2.mult=stand.cov.2.mult,stand.cov.3.mult=stand.cov.3.mult,stand.cov.4.mult=stand.cov.4.mult,
        stand.cov.5.mult=stand.cov.5.mult,stand.cov.6.mult=stand.cov.6.mult,stand.cov.7.mult=stand.cov.7.mult)

         save(bias.uni,  file = file.path("biasUni_nonMarkov_group_2.RData"))
         save(bias.mult,  file = file.path("biasMult_nonMarkov_group_2.RData"))
         save(stand.se.uni,  file = file.path("standSEUni_nonMarkov_group_2.RData"))
         save(stand.se.mult,  file = file.path("standSEMult_nonMarkov_group_2.RData"))
         save(jack.se.uni,  file = file.path("jackSEUni_nonMarkov_group_2.RData"))
         save(jack.se.mult,  file = file.path("jackSEMult_nonMarkov_group_2.RData"))
         save(cov.uni,  file = file.path("covUni_nonMarkov_group_2.RData"))
         save(cov.mult,  file = file.path("covMult_nonMarkov_group_2.RData"))