Admin — Apr 14, 2013, 9:53 PM
library(etm)
Warning: package 'etm' was built under R version 2.15.3
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 Transplantable")
#####
##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
i= 1i= 2i= 3i= 4i= 5i= 6i= 7i= 8i= 9i= 10i= 11i= 12i= 13i= 14i= 15i= 16i= 17i= 18i= 19i= 20i= 21i= 22i= 23i= 24i= 25i= 26i= 27i= 28i= 29i= 30i= 31i= 32i= 33i= 34i= 35i= 36i= 37i= 38i= 39i= 40i= 41i= 42i= 43i= 44i= 45i= 46i= 47i= 48i= 49i= 50i= 51i= 52i= 53i= 54i= 55i= 56i= 57i= 58i= 59i= 60i= 61i= 62i= 63i= 64i= 65i= 66i= 67i= 68i= 69i= 70i= 71i= 72i= 73i= 74i= 75i= 76i= 77i= 78i= 79i= 80i= 81i= 82i= 83i= 84i= 85i= 86i= 87i= 88i= 89i= 90i= 91i= 92i= 93i= 94i= 95i= 96i= 97i= 98i= 99i= 100i= 101i= 102i= 103i= 104i= 105i= 106i= 107i= 108i= 109i= 110i= 111i= 112i= 113i= 114i= 115i= 116i= 117i= 118i= 119i= 120i= 121i= 122i= 123i= 124i= 125i= 126i= 127i= 128i= 129i= 130i= 131i= 132i= 133i= 134i= 135i= 136i= 137i= 138i= 139i= 140i= 141i= 142i= 143i= 144i= 145i= 146i= 147i= 148i= 149i= 150i= 151i= 152i= 153i= 154i= 155i= 156i= 157i= 158i= 159i= 160i= 161i= 162i= 163i= 164i= 165i= 166i= 167i= 168i= 169i= 170i= 171i= 172i= 173i= 174i= 175i= 176i= 177i= 178i= 179i= 180i= 181i= 182i= 183i= 184i= 185i= 186i= 187i= 188i= 189i= 190i= 191i= 192i= 193i= 194i= 195i= 196i= 197i= 198i= 199i= 200i= 201i= 202i= 203i= 204i= 205i= 206i= 207i= 208i= 209i= 210i= 211i= 212i= 213i= 214i= 215i= 216i= 217i= 218i= 219i= 220i= 221i= 222i= 223i= 224i= 225i= 226i= 227i= 228i= 229i= 230i= 231i= 232i= 233i= 234i= 235i= 236i= 237i= 238i= 239i= 240i= 241i= 242i= 243i= 244i= 245i= 246i= 247i= 248i= 249i= 250i= 251i= 252i= 253i= 254i= 255i= 256i= 257i= 258i= 259i= 260i= 261i= 262i= 263i= 264i= 265i= 266i= 267i= 268i= 269i= 270i= 271i= 272i= 273i= 274i= 275i= 276i= 277i= 278i= 279i= 280i= 281i= 282i= 283i= 284i= 285i= 286i= 287i= 288i= 289i= 290i= 291i= 292i= 293i= 294i= 295i= 296i= 297i= 298i= 299i= 300i= 301i= 302i= 303i= 304i= 305i= 306i= 307i= 308i= 309i= 310i= 311i= 312i= 313i= 314i= 315i= 316i= 317i= 318i= 319i= 320i= 321i= 322i= 323i= 324i= 325i= 326i= 327i= 328i= 329i= 330i= 331i= 332i= 333i= 334i= 335i= 336i= 337i= 338i= 339i= 340i= 341i= 342i= 343i= 344i= 345i= 346i= 347i= 348i= 349i= 350i= 351i= 352i= 353i= 354i= 355i= 356i= 357i= 358i= 359i= 360i= 361i= 362i= 363i= 364i= 365i= 366i= 367i= 368i= 369i= 370i= 371i= 372i= 373i= 374i= 375i= 376i= 377i= 378i= 379i= 380i= 381i= 382i= 383i= 384i= 385i= 386i= 387i= 388i= 389i= 390i= 391i= 392i= 393i= 394i= 395i= 396i= 397i= 398i= 399i= 400i= 401i= 402i= 403i= 404i= 405i= 406i= 407i= 408i= 409i= 410i= 411i= 412i= 413i= 414i= 415i= 416i= 417i= 418i= 419i= 420i= 421i= 422i= 423i= 424i= 425i= 426i= 427i= 428i= 429i= 430i= 431i= 432i= 433i= 434i= 435i= 436i= 437i= 438i= 439i= 440i= 441i= 442i= 443i= 444i= 445i= 446i= 447i= 448i= 449i= 450i= 451i= 452i= 453i= 454i= 455i= 456i= 457i= 458i= 459i= 460i= 461i= 462i= 463i= 464i= 465i= 466i= 467i= 468i= 469i= 470i= 471i= 472i= 473i= 474i= 475i= 476i= 477i= 478i= 479i= 480i= 481i= 482i= 483i= 484i= 485i= 486i= 487i= 488i= 489i= 490i= 491i= 492i= 493i= 494i= 495i= 496i= 497i= 498i= 499i= 500
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"))