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