Scenario IV (sim_non_Markov_group_multinomial.r)
library(nnet)
n.iter<-1
coef.l5 <- matrix(NA,n.iter,30)
z.l5<-matrix(NA,n.iter,30)
p.l5<-matrix(NA,n.iter,30)
sig.l5<-matrix(NA,n.iter,30)
coef.l5.5<-matrix(NA,n.iter,1)
coef.l5.7<-matrix(NA,n.iter,1)
coverage.l5.5<- matrix(NA,n.iter,1)
coverage.l5.7<- matrix(NA,n.iter,1)
coef.l4 <- matrix(NA,n.iter,30)
z.l4<-matrix(NA,n.iter,30)
p.l4<-matrix(NA,n.iter,30)
sig.l4<-matrix(NA,n.iter,30)
coef.l4.5<-matrix(NA,n.iter,1)
coef.l4.7<-matrix(NA,n.iter,1)
coverage.l4.5<- matrix(NA,n.iter,1)
coverage.l4.7<- matrix(NA,n.iter,1)
coef.l3 <- matrix(NA,n.iter,30)
z.l3<-matrix(NA,n.iter,30)
p.l3<-matrix(NA,n.iter,30)
sig.l3<-matrix(NA,n.iter,30)
coef.l3.5<-matrix(NA,n.iter,1)
coef.l3.7<-matrix(NA,n.iter,1)
coverage.l3.5<- matrix(NA,n.iter,1)
coverage.l3.7<- matrix(NA,n.iter,1)
coef.l2 <- matrix(NA,n.iter,30)
z.l2<-matrix(NA,n.iter,30)
p.l2<-matrix(NA,n.iter,30)
sig.l2<-matrix(NA,n.iter,30)
coef.l2.5<-matrix(NA,n.iter,1)
coef.l2.7<-matrix(NA,n.iter,1)
coverage.l2.5<- matrix(NA,n.iter,1)
coverage.l2.7<- matrix(NA,n.iter,1)
coef.tot <- matrix(NA,n.iter,36)
z.tot<-matrix(NA,n.iter,36)
p.tot<-matrix(NA,n.iter,36)
sig.tot<-matrix(NA,n.iter,36)
coef.tot.7<-matrix(NA,n.iter,1)
coverage.tot.7<- matrix(NA,n.iter,1)
for (iter in 1:n.iter){
cat("iter=", iter)
##MARKOV with a group effect######
n<-500
state<-matrix(0,n,10)
ile<-matrix(1,n,1)
time<-matrix(0,n,10)
group<-rbinom(n,1,0.5)
###
time21<-time23<-time24<-time25<-time26<-time27<-matrix(0,n,10)
time31<-time32<-time34<-time35<-time36<-time37<-matrix(0,n,10)
time41<-time42<-time43<-time45<-time46<-time47<-matrix(0,n,10)
time51<-time52<-time53<-time54<-time56<-time57<-matrix(0,n,10)
for ( i in 1: n){
state[i,1]<-5
time[i,1]<-0
for (j in 1:10){
time21[i,j]<-rexp(1,rate=exp(-1))
#time22[i,j]<-rexp(1,rate=exp(-1))
time23[i,j]<-rexp(1,rate=exp(-1))
time24[i,j]<-rexp(1,rate=exp(-1))
time25[i,j]<-rexp(1,rate=exp(-1))
time26[i,j]<-rexp(1,rate=exp(-1))
time27[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
time31[i,j]<-rexp(1,rate=exp(-1))
time32[i,j]<-rexp(1,rate=exp(-1))
#time33[i,j]<-rexp(1,rate=exp(-1))
time34[i,j]<-rexp(1,rate=exp(-1))
time35[i,j]<-rexp(1,rate=exp(-1))
time36[i,j]<-rexp(1,rate=exp(-1))
time37[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
time41[i,j]<-rexp(1,rate=exp(-1))
time42[i,j]<-rexp(1,rate=exp(-1))
time43[i,j]<-rexp(1,rate=exp(-1))
#time44[i,j]<-rexp(1,rate=exp(-1))
time45[i,j]<-rexp(1,rate=exp(-1))
time46[i,j]<-rexp(1,rate=exp(-1))
time47[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
time51[i,j]<-rexp(1,rate=exp(-1))
time52[i,j]<-rexp(1,rate=exp(-1))
time53[i,j]<-rexp(1,rate=exp(-1))
time54[i,j]<-rexp(1,rate=exp(-1))
#time55[i,j]<-rexp(1,rate=exp(-1))
time56[i,j]<-rexp(1,rate=exp(-1))
time57[i,j]<-rexp(1,rate=exp(-1))/(1+group[i])
}
}
for ( i in 1:n){
for (j in 2:10){
if (j>2){
prev2.state<-state[i,j-2]
if (prev2.state==5){
time25[i,j]<-time25[i,j]/(group[i]+1)
time35[i,j]<-time35[i,j]/(group[i]+1)
time45[i,j]<-time45[i,j]/(group[i]+1)
}
if (prev2.state==4){
time24[i,j]<-time24[i,j]/(group[i]+1)
time34[i,j]<-time34[i,j]/(group[i]+1)
time54[i,j]<-time54[i,j]/(group[i]+1)
}
if (prev2.state==3){
time23[i,j]<-time23[i,j]/(group[i]+1)
time43[i,j]<-time43[i,j]/(group[i]+1)
time53[i,j]<-time53[i,j]/(group[i]+1)
}
}
prev.state<-state[i,j-1]
if (prev.state==1){
state[i,j]<-1}
if (prev.state==6){
state[i,j]<-6}
if (prev.state==7){
state[i,j]<-7}
if (prev.state==5){
time[i,j]<-min(time51[i,j],time52[i,j],time53[i,j],time54[i,j],9999,time56[i,j],time57[i,j])
times<-c(time51[i,j],time52[i,j],time53[i,j],time54[i,j],9999,time56[i,j],time57[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time51[i,j],time52[i,j],time53[i,j],time54[i,j],9999,time56[i,j],time57[i,j]))
}
if (prev.state==4){
time[i,j]<-min(time41[i,j],time42[i,j],time43[i,j],9999,time45[i,j],time46[i,j],time47[i,j])
times<-c(time41[i,j],time42[i,j],time43[i,j],999,time45[i,j],time46[i,j],time47[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time41[i,j],time42[i,j],time43[i,j],9999,time45[i,j],time46[i,j],time47[i,j]))
}
if (prev.state==3){
time[i,j]<-min(time31[i,j],time32[i,j],9999,time34[i,j],time35[i,j],time36[i,j],time37[i,j])
times<-c(time31[i,j],time32[i,j],9999,time34[i,j],time35[i,j],time36[i,j],time37[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time31[i,j],time32[i,j],9999,time34[i,j],time35[i,j],time36[i,j],time37[i,j]))
}
if (prev.state==2){
time[i,j]<-min(time21[i,j],9999,time23[i,j],time24[i,j],time25[i,j],time26[i,j],time27[i,j])
times<-c(time21[i,j],9999,time23[i,j],time24[i,j],time25[i,j],time26[i,j],time27[i,j])
times==time[i,j]
state[i,j]<-which(times==min(time21[i,j],9999,time23[i,j],time24[i,j],time25[i,j],time26[i,j],time27[i,j]))
}
}#end for i
}# end for j
ile<-matrix(0,n,1)
for ( i in 1:n){
ile[i,1]<-sum(as.numeric(time[i,]>0))
}
##total nr trans##
total.trans<-sum(ile)
###
###DATA WITH GROUP####
id.data<-matrix(0,total.trans,1)
group.data<-matrix(0,total.trans,1)
from<-matrix(0,total.trans,1)
to<-matrix(0,total.trans,1)
tstart<-tstop<-matrix(0,total.trans,1)
ile.juz<-0
for (i in 1: n){
for ( j in 1: ile[i,]){
id.data[ile.juz+j]<-i
group.data[ile.juz+j]<-group[i]
if (j==1){
from[ile.juz+j]<-state[i,j]
to[ile.juz+j]<-state[i,j+1]
tstart[ile.juz+j]<-0
tstop[ile.juz+j]<-time[i,j+1]
}
if(j>1){
from[ile.juz+j]<-to[ile.juz+j-1]
to[ile.juz+j]<-state[i,j+1]
tstart[ile.juz+j]<-tstop[ile.juz+j-1]
tstop[ile.juz+j]<-tstop[ile.juz+j-1]+time[i,j]
}
}
ile.juz<-ile.juz+ile[i,]
}
heart.sym<-as.data.frame(cbind(id.data,group.data,from,to,tstart,tstop))
names(heart.sym)<-c("id","group","from","to","tstart","tstop")
heart.sym$time<-heart.sym$tstop
###
####history dependence ###
###creating history- a list####
heart.sym$hist<-list("N")
num<-1
for (i in 2:(length(heart.sym$id))){
current.id<-heart.sym$id[i]
id.prev<-heart.sym$id[i-1]
if(current.id!=id.prev){num<-1}
if (current.id==id.prev){
if(heart.sym$hist[i-1]=="N"){
num<-2
heart.sym$hist[[i]][num-1]<-as.character(heart.sym$from[i-1])}else{
num<-length(heart.sym$hist[[i-1]])
for( k in 1:num){
heart.sym$hist[[i]][k]<-as.character(heart.sym$hist[[i-1]][k])}
heart.sym$hist[[i]][num+1]<-as.character(heart.sym$from[i-1])
}
}}
#####creating history
rownames(heart.sym)<-seq(1:length(heart.sym$id))
###counting the length of the history###
count<-matrix(0,1,length(heart.sym$id))
for ( i in 1: length(heart.sym$id)){
count[,i]<-as.integer(length(heart.sym$hist[[i]]))
}
heart.sym$count<-t(count)
###
lasteq5<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="5")
lasteq5[i,]<-1}
heart.sym$lasteq5<-lasteq5
###last 2
lasteq2<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="2")
lasteq2[i,]<-1}
heart.sym$lasteq2<-lasteq2
##last 3###
lasteq3<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="3")
lasteq3[i,]<-1}
heart.sym$lasteq3<-lasteq3
##last 4###
lasteq4<-matrix(0,length(heart.sym$id),1)
for (i in 1: length(heart.sym$id)){
if (heart.sym$hist[[i]][heart.sym$count[i]]=="4")
lasteq4[i,]<-1}
heart.sym$lasteq4<-lasteq4
####history 5 selected####
heart.sym.l5<-heart.sym[heart.sym$lasteq5==1,]
#heart.sym.l5$id<-seq(1:length(heart.sym.l5$id))
id.set.l5<-unique(heart.sym.l5$id)
n.id.l5 <-length(id.set.l5)
jack.se.l5<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l5)
for (i in 1:n.id.l5) {
id.c<-heart.sym.l5$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l5[heart.sym.l5$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l5[r] <- sqrt(((n.id.l5 - 1)/n.id.l5) * sum((u - mean(u))^2))
}
coef.l5[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5))
z.l5[iter,]<-(coef.l5[iter,1:30])/jack.se.l5
p.l5[iter,]<-2-2*pnorm(abs(z.l5[iter,]))
sig.l5[iter,]<-as.numeric(p.l5[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l5.5[iter]<-coef.l5[iter,22]
coef.l5.7[iter]<-coef.l5[iter,24]
coverage.l5.5[iter]<- sig.l5[iter,22]
coverage.l5.7[iter]<- sig.l5[iter,24]
######we estimate transition times for the subgroup with pooling individuals###
###history 4 selected####
heart.sym.l4<-heart.sym[heart.sym$lasteq4==1,]
id.set.l4<-unique(heart.sym.l4$id)
n.id.l4 <-length(id.set.l4)
jack.se.l4<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l4)
for (i in 1:n.id.l4) {
id.c<-heart.sym.l5$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l4[heart.sym.l4$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l4[r] <- sqrt(((n.id.l4 - 1)/n.id.l4) * sum((u - mean(u))^2))
}
coef.l4[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l4))
z.l4[iter,]<-(coef.l4[iter,1:30])/jack.se.l4
p.l4[iter,]<-2-2*pnorm(abs(z.l4[iter,]))
sig.l4[iter,]<-as.numeric(p.l4[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT (7)#
coef.l4.5[iter]<-coef.l4[iter,22]
coef.l4.7[iter]<-coef.l4[iter,24]
coverage.l4.5[iter]<- sig.l4[iter,22]
coverage.l4.7[iter]<- sig.l4[iter,24]
#####last 3
heart.sym.l3<-heart.sym[heart.sym$lasteq3==1,]
id.set.l3<-unique(heart.sym.l3$id)
n.id.l3 <-length(id.set.l3)
jack.se.l3<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l3)
for (i in 1:n.id.l3) {
id.c<-heart.sym.l5$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l3[heart.sym.l3$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l3[r] <- sqrt(((n.id.l3 - 1)/n.id.l3) * sum((u - mean(u))^2))
}
coef.l3[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l3))
z.l3[iter,]<-(coef.l3[iter,1:30])/jack.se.l3
p.l3[iter,]<-2-2*pnorm(abs(z.l3[iter,]))
sig.l3[iter,]<-as.numeric(p.l3[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l3.5[iter]<-coef.l3[iter,22]
coef.l3.7[iter]<-coef.l3[iter,24]
coverage.l3.5[iter]<- sig.l3[iter,22]
coverage.l3.7[iter]<- sig.l3[iter,24]
###last 2
heart.sym.l2<-heart.sym[heart.sym$lasteq2==1,]
id.set.l2<-unique(heart.sym.l2$id)
n.id.l2 <-length(id.set.l2)
jack.se.l2<-rep(0,30)
for ( r in 1:30){
u <- rep(0,n.id.l2)
for (i in 1:n.id.l2) {
id.c<-heart.sym.l2$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym.l2[heart.sym.l2$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.l2[r] <- sqrt(((n.id.l2 - 1)/n.id.l2) * sum((u - mean(u))^2))
}
coef.l2[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l2))
z.l2[iter,]<-(coef.l2[iter,1:30])/jack.se.l2
p.l2[iter,]<-2-2*pnorm(abs(z.l2[iter,]))
sig.l2[iter,]<-as.numeric(p.l2[iter,]<0.05)
###we asses only bias and coverages for coming back to 5 and transition to TT#
coef.l2.5[iter]<-coef.l2[iter,22]
coef.l2.7[iter]<-coef.l2[iter,24]
coverage.l2.5[iter]<- sig.l2[iter,22]
coverage.l2.7[iter]<- sig.l2[iter,24]
###total to asses better 7####
id.set.tot<-unique(heart.sym$id)
n.id.tot <-length(id.set.tot)
jack.se.tot<-rep(0,36)
for ( r in 1:36){
u <- rep(0,n.id.tot)
for (i in 1:n.id.tot) {
id.c<-heart.sym$id[i]
u[i] <- coef(multinom(to~as.factor(from)+as.factor(group)+time, data=heart.sym[heart.sym$id!=id.c,],contrasts = NULL,maxit=3000))[r]
}
jack.se.tot[r] <- sqrt(((n.id.tot - 1)/n.id.tot) * sum((u - mean(u))^2))
}
coef.tot[iter,] <- coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym))
z.tot[iter,]<-(coef.tot[iter,1:36])/jack.se.tot
p.tot[iter,]<-2-2*pnorm(abs(z.tot[iter,]))
sig.tot[iter,]<-as.numeric(p.tot[iter,]<0.05)
coef.tot.7[iter]<-coef.tot[iter,30]
coverage.tot.7[iter]<- sig.tot[iter,30]
}
iter= 1# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 431.933646
iter 20 value 424.189870
iter 30 value 423.130775
iter 40 value 423.060801
final value 423.059918
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 433.356309
iter 20 value 424.703050
iter 30 value 423.666997
iter 40 value 423.570330
final value 423.569537
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 432.273103
iter 20 value 424.536869
iter 30 value 423.513515
iter 40 value 423.418895
final value 423.417887
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 431.257315
iter 20 value 424.137881
iter 30 value 423.154080
iter 40 value 423.066380
final value 423.065316
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 430.855419
iter 20 value 424.640749
iter 30 value 423.600830
iter 40 value 423.524824
final value 423.524152
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 432.375229
iter 20 value 424.217145
iter 30 value 422.993218
iter 40 value 422.895711
final value 422.894611
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 432.428738
iter 20 value 424.115228
iter 30 value 423.134155
iter 40 value 423.004464
final value 423.002782
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 434.204872
iter 20 value 424.584257
iter 30 value 423.343592
iter 40 value 423.240819
final value 423.240024
converged
# weights: 42 (30 variable)
initial value 476.747987
iter 10 value 431.199322
iter 20 value 424.696606
iter 30 value 423.700171
iter 40 value 423.612941
final value 423.612031
converged
...
Showing results for multinomial model with last state equal T:
names<-colnames(coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5)))
model.show<-coef(multinom(to~as.factor(from)+as.factor(group)+time,data=heart.sym.l5))
colnames(model.show)<-names
jack.se.show<-matrix(jack.se.l5,nrow=6,byrow=TRUE)
colnames(jack.se.show)<-names
model.show
(Intercept) as.factor(from)3 as.factor(from)4 as.factor(group)1 time
2 -15.2893928 14.7929774 15.7293827 -0.24395013 0.3325791
3 0.2919006 -16.6952682 0.4077411 -0.40012116 0.1487177
4 -0.5179829 0.1592948 -16.5820002 0.01348601 0.4403449
5 0.3346142 -0.3547416 0.7267854 -0.09079121 0.1717782
6 0.3856212 -0.6485112 0.1653082 -0.98617957 0.1711206
7 -0.2223938 -0.2254897 -0.1327517 0.70863071 0.2872732
Jackknife SE:
jack.se.show
(Intercept) as.factor(from)3 as.factor(from)4 as.factor(group)1 time
[1,] 2.9278207 0.5421481 0.7927232 0.5576469 0.5590972
[2,] 0.5506199 2.9393823 5.4270839 0.6940920 0.5796293
[3,] 0.6154395 0.5647047 2.9785828 0.6222915 6.9418369
[4,] 0.6510179 0.6447650 0.6376941 0.5658553 0.5417642
[5,] 0.6438745 0.4832308 0.5345415 0.4929255 0.3881445
[6,] 0.3926737 0.4770983 0.3703444 0.3865607 0.3625086
list.l5<-list( coef.l5=coef.l5,z.l5=z.l5,p.l5=p.l5,sig.l5=sig.l5,coef.l5.5=coef.l5.5,coef.l5.7=coef.l5.7,coverage.l5.5=coverage.l5.5,coverage.l5.7=coverage.l5.7)
list.l4<-list( coef.l4=coef.l4,z.l4=z.l4,p.l4=p.l4,sig.l4=sig.l4,coef.l4.5=coef.l4.5,coef.l4.7=coef.l4.7,coverage.l4.5=coverage.l4.5,coverage.l4.7=coverage.l4.7)
list.l3<-list( coef.l3=coef.l3,z.l3=z.l3,p.l3=p.l3,sig.l3=sig.l3,coef.l3.5=coef.l3.5,coef.l3.7=coef.l3.7,coverage.l3.5=coverage.l3.5,coverage.l3.7=coverage.l3.7)
list.l2<-list( coef.l2=coef.l2,z.l2=z.l2,p.l2=p.l2,sig.l2=sig.l2,coef.l2.5=coef.l2.5,coef.l2.7=coef.l2.7,coverage.l2.5=coverage.l2.5,coverage.l2.7=coverage.l2.7)
list.l.tot<-list( coef.tot=coef.tot,z.tot=z.tot,p.tot=p.tot,sig.tot=sig.tot,coef.tot.7=coef.tot.7,coverage.tot.7=coverage.tot.7)
save(list.l5, file = file.path("list.l5_nonMarkov_group.RData"))
save(list.l4, file = file.path("list.l4_nonMarkov_group.RData"))
save(list.l3, file = file.path("list.l3_nonMarkov_group.RData"))
save(list.l2, file = file.path("list.l2_nonMarkov_group.RData"))
save(list.l.tot, file = file.path("list.l.tot_nonMarkov_group.RData"))