testing/scratch2.R

In vsritter/sim.etm: What the Package Does (One Line, Title Case)

# Q[1,2,] <- c(1, 1)
# Q[1,3,] <- c(1, .2)
# Q[2,3,] <- c(1, 1)


K=5000;N=K;s=0
tau=600;a=160
a12=1;a13=0.2;a23=1
b12=0;b13=-.8;b23=0
set.seed(237985)
# seed.vector=sample(1:99999,1000)
set.seed(76564)
markov12=matrix(0,6,(2*K)+7)
nonmarkov12=matrix(0,6,K+6)
semimarkov12=matrix(0,6,K+6)
for(k in 1:K)
{
  if (k %% 500 == 0) cat(k, "\n")
  # set.seed(seed.vector[k])
  #Startbygeneratingdata:
  U12=runif(N,0,1)
  U13=runif(N,0,1)
  V23=runif(N,0,1)
  T12=(-((b12+1)/a12)*log(U12))^(1/(b12+1))
  T13=(-((b13+1)/a13)*log(U13))^(1/(b13+1))
  U=runif(N,0,a)+Inf
  C=rep(0,N)
  for(e in 1:N)
  {
    C[e]=min(tau,U[e])
  }
  T1=rep(0,N)
  delta=rep(0,N)
  for(e in 1:N)
  {
    T1[e]=min(T12[e],T13[e],C[e])
    if(T1[e]==T12[e])
    {delta[e]=1}
  }
  MT23=rep(0,N)
  for(e in 1:N)
  {
    MT23[e]=(-(b23+1)/a23*log(V23[e])+
                T1[e]^(b23+1))^(1/(b23+1))
  }
  MT23status=(MT23-T1)*delta
  T2=rep(0,N)
  for(e in 1:N)
  {
    if(C[e]>=T1[e]&T1[e]!=T13[e])
    {T2[e]=min(MT23status[e],C[e]-T1[e])}
    else
    {T2[e]=0}
  }
  time=rep(0,N)
  for(e in 1:N)
  {
    time[e]=min(C[e],T1[e]+T2[e])
  }
  status=rep(0,N)
  for(e in 1:N)
  {
    if(min(T12[e],T13[e],C[e])==T12[e])
    {
      if(T1[e]+MT23status[e]<=C[e])
      {status[e]=1}
    }
    if(min(T12[e],T13[e],C[e])==T13[e])
    {status[e]=1}
  }
}

#Generateddata:
dataMarkov=matrix(0,N,5)
dataMarkov[,1]=T1
dataMarkov[,2]=delta
dataMarkov[,3]=T2
dataMarkov[,4]=time
dataMarkov[,5]=status


dt_etm <- data.frame(
  id = 1:N,
  from = 'T1',
  to = ifelse(delta == 1, 'T2', 'T3'),
  entry = 0,
  exit = T1, stringsAsFactors = F
) %>% bind_rows(data.frame(
  id = c(1:N)[which(delta==1)],
  from = 'T2',
  to = 'T3',
  entry = T1[which(delta==1)],
  exit = time[which(delta==1)], stringsAsFactors = F
)) %>%
  arrange(id, exit)


#Calculatethetruetransitionprobabilities:
np=1000
g <- function(u) {
  exp(-((a12/(b12+1)*u^(b12+1))+(a13/(b13+1)*u^(b13+1)))
       +((a12/(b12+1)*s^(b12+1))+(a13/(b13+1)*s^(b13+1))))*
    a12*u^b12*exp(-(a23/(b23+1)*t^(b23+1))+
                     (a23/(b23+1)*u^(b23+1)))
}
P12=rep(0,np)
e=1
for(t in seq(0,5,length.out = np))
{
  if(t<=s)
  {P12[e]=0}
  else
  {P12[e]=integrate(g,s,t)$value}
  e=e+1
}


etm <- etm::etm(dt_etm, state_names, tra, "cens", s=0, covariance = F)
etm <- data.frame(p = as.vector(etm$est),
                    t = rep(etm$time, each = n_states^2),
                    trans = rep(trans, times = length(t)),
                    stringsAsFactors = F) %>%
  filter(trans %in% poss_trans)


plot(filter(p2, trans == 'T12')$t, filter(p2, trans == 'T12')$p, col=1, type='l')
lines(filter(p, trans == 'T12')$t, filter(p, trans == 'T12')$p, col=5, type='l')
lines(filter(etm_p, trans == 'T12')$t, filter(etm_p, trans == 'T12')$p, col=2, type='l')

lines(seq(0,5,length.out = np), P12, col = 3, type = 'l')
lines(filter(etm, trans == 'T12')$t, filter(etm, trans == 'T12')$p, col=4, type='l')