Title

This is an R Markdown document. Markdown is a simple formatting syntax for authoring web pages (click the MD toolbar button for help on Markdown).

When you click the Knit HTML button a web page will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

n.iter <- 2

p7.0.list <- list()
p7.1.list <- list()

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



    #
    # sym.data.heart.nozero<-sym.data.heart.nozero[sym.data.heart.nozero$from!=sym.data.heart.nozero$to,]


    # sym.data.heart.nozero$time<-sym.data.heart.nozero$time

    # heart.sym2<-heart.sym

    # heart.sym2$id<-seq(1:length(heart.sym2$id))

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

    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.i - 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=2,xlim=c(0,10),ylim=c(0,1),col=1)

    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=2,xlim=c(0,10),ylim=c(0,1),col=2)

    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.i - 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=2,xlim=c(0,10),ylim=c(0,1),col=2)

    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,7),ylim=c(0,1),lty=1)
    # lines(new.model.6.1[,2],p7.1,xlab=list('Time',cex=1.5),type='l',lwd=2,xlim=c(0,7),ylim=c(0,1),col=1,lty=1)

    #
    # lines(new.model.6.0[,2],p7.0,xlab=list('Time',cex=1.5),type='l',lwd=2,xlim=c(0,7),ylim=c(0,1),col=1,lty=2)

    p7.1.list[[iter]] <- p7.1
    p7.0.list[[iter]] <- p7.0
}
## iter= 1
## Warning: package 'etm' was built under R version 2.15.3
## Loading required package: lattice
## Loading required package: survival
## Loading required package: splines
## iter= 2




p7.1.list.m <- list()
p7.0.list.m <- list()

n.iter <- 50



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]
                  time35[i, j] <- time35[i, j]
                  time45[i, j] <- time45[i, j]

                }

                if (prev2.state == 4) {

                  time24[i, j] <- time24[i, j]
                  time34[i, j] <- time34[i, j]
                  time54[i, j] <- time54[i, j]

                }

                if (prev2.state == 3) {

                  time23[i, j] <- time23[i, j]
                  time43[i, j] <- time43[i, j]
                  time53[i, j] <- time53[i, j]

                }
            }




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



    #
    # sym.data.heart.nozero<-sym.data.heart.nozero[sym.data.heart.nozero$from!=sym.data.heart.nozero$to,]


    # sym.data.heart.nozero$time<-sym.data.heart.nozero$time

    # heart.sym2<-heart.sym

    # heart.sym2$id<-seq(1:length(heart.sym2$id))

    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=2,xlim=c(0,10),ylim=c(0,1),col=2)

    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.i - 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=2,xlim=c(0,10),ylim=c(0,1),col=2)

    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,7),ylim=c(0,1),lty=1)
    # lines(new.model.6.1[,2],p7.1,xlab=list('Time',cex=1.5),type='l',lwd=2,xlim=c(0,7),ylim=c(0,1),col=2,lty=1)

    #
    # lines(new.model.6.0[,2],p7.0,xlab=list('Time',cex=1.5),type='l',lwd=2,xlim=c(0,7),ylim=c(0,1),col=2,lty=2)

    p7.1.list.m[[iter]] <- p7.1
    p7.0.list.m[[iter]] <- p7.0
}
## iter= 1iter= 2iter= 3iter= 4iter= 5iter= 6iter= 7iter= 8iter= 9iter= 10iter= 11iter= 12iter= 13iter= 14iter= 15iter= 16iter= 17iter= 18iter= 19iter= 20iter= 21iter= 22iter= 23iter= 24iter= 25iter= 26iter= 27iter= 28iter= 29iter= 30iter= 31iter= 32iter= 33iter= 34iter= 35iter= 36iter= 37iter= 38iter= 39iter= 40iter= 41iter= 42iter= 43iter= 44iter= 45iter= 46iter= 47iter= 48iter= 49iter= 50

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

new.table <- do.call(rbind, p7.1.list)
## Warning: number of columns of result is not a multiple of vector length
## (arg 1)
new.table.m <- do.call(rbind, p7.1.list.m)
## Warning: number of columns of result is not a multiple of vector length
## (arg 1)


new.table.0 <- do.call(rbind, p7.0.list)
## Warning: number of columns of result is not a multiple of vector length
## (arg 1)
new.table.0.m <- do.call(rbind, p7.0.list.m)
## Warning: number of columns of result is not a multiple of vector length
## (arg 1)

quantiles.1 <- apply(new.table, 2, function(x) {
    qx <- quantile(x, c(0.025, 0.5, 0.975), na.rm = TRUE)
    qx
})

quantiles.1.m <- apply(new.table.m, 2, function(x) {
    qx <- quantile(x, c(0.025, 0.5, 0.975), na.rm = TRUE)
    qx
})

quantiles.0 <- apply(new.table.0, 2, function(x) {
    qx <- quantile(x, c(0.025, 0.5, 0.975), na.rm = TRUE)
    qx
})

quantiles.0.m <- apply(new.table.0.m, 2, function(x) {
    qx <- quantile(x, c(0.025, 0.5, 0.975), na.rm = TRUE)
    qx
})

You can also embed plots, for example:


plot(new.model.6.1[1:400, 2], quantiles.1[2, 1:400], type = "l", col = 1, lwd = 2, 
    lty = 1, ylim = c(0, 0.7), xlim = c(0, 1.2), xlab = "Days", ylab = "Occupation probability", 
    cex.lab = 1.5)
lines(new.model.6.1[1:400, 2], quantiles.1[1, 1:400], type = "l", col = 1, lty = 2, 
    lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.1[3, 1:400], type = "l", col = 1, lty = 3, 
    lwd = 2)


lines(new.model.6.1[1:400, 2], quantiles.1.m[2, 1:400], type = "l", col = 2, 
    lty = 1, lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.1.m[1, 1:400], type = "l", col = 2, 
    lty = 3, lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.1.m[3, 1:400], type = "l", col = 2, 
    lty = 3, lwd = 2)

lines(new.model.6.1[1:400, 2], quantiles.0[2, 1:400], type = "l", col = 3, lty = 1, 
    lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.0[1, 1:400], type = "l", col = 3, lty = 3, 
    lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.0[3, 1:400], type = "l", col = 3, lty = 3, 
    lwd = 2)

lines(new.model.6.1[1:400, 2], quantiles.0.m[2, 1:400], type = "l", col = 4, 
    lty = 1, lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.0.m[1, 1:400], type = "l", col = 4, 
    lty = 3, lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.0.m[3, 1:400], type = "l", col = 4, 
    lty = 3, lwd = 2)

legend(0.05, 1.2, c("Group 1 non-Markov", "Group 1 non-Markov 25%", "Group 1 non-Markov 97.5%", 
    "Group 1 Markov", "Group 1 Markov 25%", "Group 1 Markov 97.5%", "Group 0 non-Markov", 
    "Group 0 non-Markov 25%", "Group 0 non-Markov 97.5%", "Group 0 Markov", 
    "Group 0 Markov 25%", "Group 0 Markov 97.5%"), col = c(1, 1, 1, 2, 2, 2, 
    3, 3, 3, 4, 4, 4), lty = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3), lwd = c(2, 
    1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1), merge = TRUE, bg = "gray90")
title("Aalen-Johansen estimator for OP for state Transplanted")

plot of chunk unnamed-chunk-2


plot(new.model.6.1[1:400, 2], quantiles.1[2, 1:400], type = "l", col = 1, lwd = 2, 
    lty = 1, ylim = c(0, 1.2), xlim = c(0, 1.2), xlab = "Days", ylab = "Occupation probability")
lines(new.model.6.1[1:400, 2], quantiles.1[1, 1:400], type = "l", col = 1, lty = 2)
lines(new.model.6.1[1:400, 2], quantiles.1[3, 1:400], type = "l", col = 1, lty = 3)


lines(new.model.6.1[1:400, 2], quantiles.1.m[2, 1:400], type = "l", col = 2, 
    lty = 1, lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.1.m[1, 1:400], type = "l", col = 2, 
    lty = 2)
lines(new.model.6.1[1:400, 2], quantiles.1.m[3, 1:400], type = "l", col = 2, 
    lty = 3)

lines(new.model.6.1[1:400, 2], quantiles.0[2, 1:400], type = "l", col = 3, lty = 1, 
    lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.0[1, 1:400], type = "l", col = 3, lty = 2)
lines(new.model.6.1[1:400, 2], quantiles.0[3, 1:400], type = "l", col = 3, lty = 3)

lines(new.model.6.1[1:400, 2], quantiles.0.m[2, 1:400], type = "l", col = 4, 
    lty = 1, lwd = 2)
lines(new.model.6.1[1:400, 2], quantiles.0.m[1, 1:400], type = "l", col = 4, 
    lty = 2)
lines(new.model.6.1[1:400, 2], quantiles.0.m[3, 1:400], type = "l", col = 4, 
    lty = 3)

legend(0.05, 1.2, c("Group 1 non-Markov", "Group 1 non-Markov 25%", "Group 1 non-Markov 97.5%", 
    "Group 1 Markov", "Group 1 Markov 25%", "Group 1 Markov 97.5%", "Group 0 non-Markov", 
    "Group 0 non-Markov 25%", "Group 0 non-Markov 97.5%", "Group 0 Markov", 
    "Group 0 Markov 25%", "Group 0 Markov 97.5%"), col = c(1, 1, 1, 2, 2, 2, 
    3, 3, 3, 4, 4, 4), lty = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3), lwd = c(2, 
    1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1), merge = TRUE, bg = "gray90")
title("Aalen-Johansen estimator for OP for state Transplanted")

plot of chunk unnamed-chunk-2