ADM Assignment 6

Code:

Erweiterung des selbstgeschriebenen Proxel-Programms mit den Tests und der Symptom-Sequenz. (siehe Code-Button rechts)

max.ts = 12
delta = 1
startproxel=c(1,0,0,1,0)
RG=matrix(c(1,1,0,1,
            0,1,1,1,
            0,0,1,1,
            0,0,0,1),ncol = 4,byrow = T)
FEVER=data.frame(nf=c(0.9,0.5,0.2,1),f=c(0.1,0.5,0.8,0))
BLOOD=data.frame(Ba=c(.7,.5,.2),Bb=c(.2,.3,.3),Bc=c(.1,.2,.5))
URINE=data.frame(Ua=c(.8,.6,.5),Ub=c(.1,.3,.3),Uc=c(.1,.1,.2))

Proxel_algo=function(RG,startproxel,delta,max.ts,Print=T,RETURN,FEVER=NA,BLOOD=NA,URINE=NA,SEQ=NA){
  
  SF=subset(SEQ,SEQ[,2] %in% c("F","nF"))
  SB=subset(SEQ,SEQ[,2] %in% c("Ba","Bb","Bc"))
  SU=subset(SEQ,SEQ[,2] %in% c("Ua","Ub","Uc"))
  
  unihrf <- function(age, a, b) {
    if (a <= age && age < b) {
      return (1.0 / (b - age))
    } else {
      return (0)
    }
  }
  
  exphrf <- function(age, l) {
    return (l)
  }
    
  proxel_df=data.frame(tau=numeric(),S1=numeric(),S2=numeric(),S3=numeric(),H=numeric())
  
  proxel_list=list()
  proxel_list[[max.ts+1]]=list()
  proxel_list[[1]]=list(timestep=c(0),proxel=data.frame(state=startproxel[1],tau=startproxel[2],tau_k=startproxel[3],prob=startproxel[4],prev=startproxel[5]),
                        HnMM=tidyr::tibble(pos.path=list(c(1)),WKT=FEVER[1,1]))
  
  if(Print) print("initial proxel:")
  if(Print) print(proxel_list[[1]])
  
  proxel_df[1,]=c(0,proxel_list[[1]]$proxel$prob[proxel_list[[1]]$proxel$state == 1],0,0,0)
  
  for (timestep in 1:max.ts) {
    proxel_list[[timestep+1]]=list(timestep=c(timestep),proxel=data.frame(state=numeric(),tau=numeric(),tau_k=numeric(),prob=numeric(),prev=numeric()),
                                   HnMM=tidyr::tibble(pos.path=list(),WKT=numeric()))     # erstelle schonmal eine neue leere liste für den nächsten timestep
    
    for (l in 1:nrow(proxel_list[[timestep]]$proxel)) {
      currprox = proxel_list[[timestep]]$proxel[l,]                                     # aktueller proxel zu dem folgeproxels erzeugt werden sollen
      folger = RG[as.numeric(currprox[1]),]                                             # folgestates von aktuellem proxel
      
      for (i in 1:4) {
        if (folger[i]==1) {                                                             # erstelle folgeproxel falls pfeil dahin gibt
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel)+1,1] = i                       # state vom nächsten proxel
          
          if(currprox$state==i) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),2] = currprox$tau + delta
            else proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),2] = 0                    # tau vom nächsten proxel
          
          if(currprox$state %in% c(1,2,3) & i!=4) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = currprox$tau_k + delta
            else if(currprox$state==4 & i==4) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = currprox$tau_k - delta
              else proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = 0                   # tau_k vom nächsten proxel
            
          if(currprox$state==1 & i==2) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (unihrf(currprox$tau_k,2,4) * delta) * currprox$prob}                 # transition von S1 zu S2
          if(currprox$state==2 & i==3) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (exphrf(currprox$tau_k,0.1) * delta) * currprox$prob}                 # transition von S2 zu S3
          if(currprox$state!=4 & i==4) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (unihrf(currprox$tau_k,7,14)* delta) * currprox$prob}                 # transtion zu heal von S1,S2 oder S3
          if(currprox$state==1 & i==1) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (unihrf(currprox$tau_k,2,4) * delta) - (unihrf(currprox$tau_k,7,14)* delta)) * currprox$prob}
          if(currprox$state==2 & i==2) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (exphrf(currprox$tau_k,0.1) * delta) - (unihrf(currprox$tau_k,7,14)* delta)) * currprox$prob}
          if(currprox$state==3 & i==3) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (unihrf(currprox$tau_k,7,14)* delta)) * currprox$prob}
          if(currprox$state==4 & i==4) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = 1 }
              
           proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),5] = l   
           
        }
      } 
    }
    
    if(length(which(proxel_list[[timestep+1]]$proxel[,4]<=0))!=0) proxel_list[[timestep+1]]$proxel=proxel_list[[timestep+1]]$proxel[-which(proxel_list[[timestep+1]]$proxel[,4]<=0),]
    if(length(which(proxel_list[[timestep+1]]$proxel[,3]< 0))!=0) proxel_list[[timestep+1]]$proxel=proxel_list[[timestep+1]]$proxel[-which(proxel_list[[timestep+1]]$proxel[,3]< 0),]
    
    for (prox in 1:nrow(proxel_list[[timestep+1]]$proxel)) {
      
      proxel_list[[timestep+1]]$HnMM[prox,1][[1]] = list(c(unlist(proxel_list[[timestep]]$HnMM[proxel_list[[timestep+1]]$proxel$prev[prox],1]),proxel_list[[timestep+1]]$proxel$state[prox]))
      
      if(timestep+1 <= nrow(SF)){
        if(SF[timestep+1,2]=="F") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$proxel$prob[prox] * FEVER$f[proxel_list[[timestep+1]]$proxel$state[prox]] *
            proxel_list[[timestep]]$HnMM$WKT[proxel_list[[timestep+1]]$proxel$prev[prox]]
        if(SF[timestep+1,2]=="nF") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$proxel$prob[prox] * FEVER$nf[proxel_list[[timestep+1]]$proxel$state[prox]] *
            proxel_list[[timestep]]$HnMM$WKT[proxel_list[[timestep+1]]$proxel$prev[prox]]
      }
      
      if((timestep+1)%%2 == 0){
        if(SU[SU[,1]==timestep+1,2]=="Ua") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            URINE$Ua[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SU[SU[,1]==timestep+1,2]=="Ub") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            URINE$Ub[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SU[SU[,1]==timestep+1,2]=="Uc") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            URINE$Uc[proxel_list[[timestep+1]]$proxel$state[prox]]
      }
      
      if((timestep+1)%%3 == 0){
        if(SB[SB[,1]==timestep+1,2]=="Ba") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            BLOOD$Ba[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SB[SB[,1]==timestep+1,2]=="Bb") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            BLOOD$Bb[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SB[SB[,1]==timestep+1,2]=="Bc") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            BLOOD$Bc[proxel_list[[timestep+1]]$proxel$state[prox]]
      }
    }
    
    if(Print) print(paste("Time step:",timestep))
    if(Print) print(proxel_list[[timestep+1]])
    
    proxel_df[timestep+1,1]=timestep*delta
    for (v in 2:4) {
      if(length(proxel_list[[timestep+1]]$proxel$prob[proxel_list[[timestep+1]]$proxel$state == v-1]) != 0) {
        proxel_df[timestep+1,v] = round(sum(proxel_list[[timestep+1]]$proxel$prob[proxel_list[[timestep+1]]$proxel$state == v-1]),5)}
          else proxel_df[timestep+1,v] = 0
    }
    proxel_df[timestep+1,5]=round(1-sum(proxel_df[timestep+1,2:4]),5)
      
  }
  
  if(RETURN=="TREE") return(proxel_list)
  if(RETURN=="DF") return(proxel_df)
}
proxel_list=Proxel_algo(RG,startproxel,delta,max.ts,F,"TREE",FEVER,BLOOD,URINE,patient1)

What is the probability of the given sequence of symptoms (patient1.seq)?

Prob_dis1 = sum(proxel_list[[13]]$HnMM$WKT)

• Prob = 1.228379^{-21}

What are the 5 most likely generating paths of that sequence? & What are the actual portions of time spent in each disease stage for each of these paths?

ord=order(proxel_list[[13]]$HnMM$WKT,decreasing = T)

as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[1]]])
proxel_list[[13]][["HnMM"]][["WKT"]][[ord[1]]]
table(as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[1]]]))

as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[2]]])
proxel_list[[13]][["HnMM"]][["WKT"]][[ord[2]]]
as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[2]]]) |> table()

as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[3]]])
proxel_list[[13]][["HnMM"]][["WKT"]][[ord[3]]]
as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[3]]]) |> table()

as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[4]]])
proxel_list[[13]][["HnMM"]][["WKT"]][[ord[4]]]
as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[4]]]) |> table()

as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[5]]])
proxel_list[[13]][["HnMM"]][["WKT"]][[ord[5]]]
as.vector(proxel_list[[13]][["HnMM"]][["pos.path"]][[ord[5]]]) |> table()

1. Path: 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2
   • Prob = 5.6322607^{-22}
   • days in stage…
     • 1: 4
     • 2: 9
2. Path: 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2
   • Prob = 3.9137598^{-22}
   • days in stage…
     • 1: 3
     • 2: 10
3. Path: 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3
   • Prob = 1.5902854^{-22}
   • days in stage…
     • 1: 4
     • 2: 8
     • 3: 1
4. Path: 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3
   • Prob = 1.1050616^{-22}
   • days in stage…
     • 1: 3
     • 2: 9
     • 3: 1
5. Path: 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3
   • Prob = 1.9005546^{-24}
   • days in stage…
     • 1: 4
     • 2: 7
     • 3: 2

Code 2:

Abänderung des HnMM-Programms für Disease 2

max.ts = 12
delta = 1
startproxel=c(1,0,0,1,0)
RG=matrix(c(1,1,0,1,
            0,1,1,1,
            0,0,1,1,
            0,0,0,1),ncol = 4,byrow = T)
FEVER2=data.frame(nf=c(0.7,0.3,0.1,1),f=c(0.3,0.7,0.9,0))
BLOOD2=data.frame(Ba=c(.7,.5,.2),Bb=c(.2,.3,.3),Bc=c(.1,.2,.5))
URINE2=data.frame(Ua=c(.8,.7,.7),Ub=c(.2,.2,.2),Uc=c(.0,.1,.1))

Proxel_algo2=function(RG,startproxel,delta,max.ts,Print=T,RETURN,FEVER=NA,BLOOD=NA,URINE=NA,SEQ=NA){
  
  SF=subset(SEQ,SEQ[,2] %in% c("F","nF"))
  SB=subset(SEQ,SEQ[,2] %in% c("Ba","Bb","Bc"))
  SU=subset(SEQ,SEQ[,2] %in% c("Ua","Ub","Uc"))
  ###########################################################################################################################
  normalpdf <- function(x, m, s) {
  z <- (x - m) / s
  return (exp(-z^2 / 2) / (sqrt(2 * pi) * s))
}

logGamma <- function(x) {
  coef <- c(76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.00120858003, -0.00000536382)
  step <- 2.50662827465
  fpf <- 5.5
  t <- x - 1
  tmp <- t + fpf
  tmp <- (t + 0.5) * log(tmp) - tmp
  ser <- 1
  for (i in 1:6) {
    t <- t + 1
    ser <- ser + coef[i] / t
  }
  return (tmp + log(step * ser))
}

gammaSeries <- function(x, a) {
  maxit <- 100
  eps <- 0.0000003
  sum <- 1.0 / a
  ap <- a
  gln <- logGamma(a)
  del <- sum

  for (n in 1:maxit) {
    ap <- ap + 1
    del <- del * x / ap
    sum <- sum + del
    if (abs(del) < abs(sum) * eps) break
  }
  return (sum * exp(-x + a * log(x) - gln))
}

gammaCF <- function(x, a) {
  maxit <- 100
  eps <- 0.0000003
  gln <- logGamma(a)
  g <- 0
  gold <- 0
  a0 <- 1
  a1 <- x
  b0 <- 0
  b1 <- 1
  fac <- 1
  an <- 0
  ana <- 0
  anf <- 0

  for (n in 1:maxit) {
    an <- 1.0 * n
    ana <- an - a
    a0 <- (a1 + a0 * ana) * fac
    b0 <- (b1 + b0 * ana) * fac
    anf <- an * fac
    a1 <- x * a0 + anf * a1
    b1 <- x * b0 + anf * b1
    if (a1 != 0) {
      fac <- 1.0 / a1
      g <- b1 * fac
      if (abs((g - gold) / g) < eps) break
      gold <- g
    }
  }
  return (exp(-x + a * log(x) - gln) * g)
}

gammacdf <- function(x, a) {
  if (x <= 0)
    return (0)
  else if (x < a + 1)
    return (gammaSeries(x, a))
  else
    return (1 - gammaCF(x, a))
}

normalcdf <- function(x, m, s) {
  z <- (x - m) / s
  if (z >= 0)
    return (0.5 + 0.5 * gammacdf(z^2 / 2, 0.5))
  else
    return (0.5 - 0.5 * gammacdf(z^2 / 2, 0.5))
}

normalhrf <- function(x, m, s) {
  return (normalpdf(x, m, s) / (1 - normalcdf(x, m, s)))
}
  
  unihrf <- function(age, a, b) {
    if (a <= age && age < b) {
      return (1.0 / (b - age))
    } else {
      return (0)
    }
  }
  
  exphrf <- function(age, l) {
    return (l)
  }
   #####################################################################################################################################################
  
  proxel_df=data.frame(tau=numeric(),S1=numeric(),S2=numeric(),S3=numeric(),H=numeric())
  
  proxel_list=list()
  proxel_list[[max.ts+1]]=list()
  proxel_list[[1]]=list(timestep=c(0),proxel=data.frame(state=startproxel[1],tau=startproxel[2],tau_k=startproxel[3],prob=startproxel[4],prev=startproxel[5]),
                        HnMM=tidyr::tibble(pos.path=list(c(1)),WKT=FEVER[1,1]))
  
  if(Print) print("initial proxel:")
  if(Print) print(proxel_list[[1]])
  
  proxel_df[1,]=c(0,proxel_list[[1]]$proxel$prob[proxel_list[[1]]$proxel$state == 1],0,0,0)
  
  for (timestep in 1:max.ts) {
    proxel_list[[timestep+1]]=list(timestep=c(timestep),proxel=data.frame(state=numeric(),tau=numeric(),tau_k=numeric(),prob=numeric(),prev=numeric()),
                                   HnMM=tidyr::tibble(pos.path=list(),WKT=numeric()))     # erstelle schonmal eine neue leere liste für den nächsten timestep
    
    for (l in 1:nrow(proxel_list[[timestep]]$proxel)) {
      currprox = proxel_list[[timestep]]$proxel[l,]                                     # aktueller proxel zu dem folgeproxels erzeugt werden sollen
      folger = RG[as.numeric(currprox[1]),]                                             # folgestates von aktuellem proxel
      
      for (i in 1:4) {
        if (folger[i]==1) {                                                             # erstelle folgeproxel falls pfeil dahin gibt
          proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel)+1,1] = i                       # state vom nächsten proxel
          
          if(currprox$state==i) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),2] = currprox$tau + delta
            else proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),2] = 0                    # tau vom nächsten proxel
          
          if(currprox$state %in% c(1,2,3) & i!=4) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = currprox$tau_k + delta
            else if(currprox$state==4 & i==4) proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = currprox$tau_k - delta
              else proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),3] = 0                   # tau_k vom nächsten proxel
            
          if(currprox$state==1 & i==2) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (normalhrf(currprox$tau_k,3,0.5) * delta) * currprox$prob}           # transition von S1 zu S2
          if(currprox$state==2 & i==3) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (exphrf(currprox$tau_k,0.3) * delta) * currprox$prob}           # transition von S2 zu S3
          if(currprox$state!=4 & i==4) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (normalhrf(currprox$tau_k,10,2)* delta) * currprox$prob}           # transtion zu heal von S1,S2 oder S3
          if(currprox$state==1 & i==1) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (normalhrf(currprox$tau_k,3,0.5) * delta) - (normalhrf(currprox$tau_k,10,2) * delta)) * currprox$prob}
          if(currprox$state==2 & i==2) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (exphrf(currprox$tau_k,0.3) * delta) - (normalhrf(currprox$tau_k,10,2) * delta)) * currprox$prob}
          if(currprox$state==3 & i==3) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = (1 - (normalhrf(currprox$tau_k,10,2) * delta)) * currprox$prob}
          if(currprox$state==4 & i==4) {
            proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),4] = 1 }
              
           proxel_list[[timestep+1]]$proxel[nrow(proxel_list[[timestep+1]]$proxel),5] = l   
           
        }
      } 
    }
    
    if(length(which(proxel_list[[timestep+1]]$proxel[,4]<=0))!=0) proxel_list[[timestep+1]]$proxel=proxel_list[[timestep+1]]$proxel[-which(proxel_list[[timestep+1]]$proxel[,4]<=0),]
    if(length(which(proxel_list[[timestep+1]]$proxel[,3]< 0))!=0) proxel_list[[timestep+1]]$proxel=proxel_list[[timestep+1]]$proxel[-which(proxel_list[[timestep+1]]$proxel[,3]< 0),]
    
    for (prox in 1:nrow(proxel_list[[timestep+1]]$proxel)) {
      
      proxel_list[[timestep+1]]$HnMM[prox,1][[1]] = list(c(unlist(proxel_list[[timestep]]$HnMM[proxel_list[[timestep+1]]$proxel$prev[prox],1]),proxel_list[[timestep+1]]$proxel$state[prox]))
      
      if(timestep+1 <= nrow(SF)){
        if(SF[timestep+1,2]=="F") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$proxel$prob[prox] * FEVER$f[proxel_list[[timestep+1]]$proxel$state[prox]] *
            proxel_list[[timestep]]$HnMM$WKT[proxel_list[[timestep+1]]$proxel$prev[prox]]
        if(SF[timestep+1,2]=="nF") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$proxel$prob[prox] * FEVER$nf[proxel_list[[timestep+1]]$proxel$state[prox]] *
            proxel_list[[timestep]]$HnMM$WKT[proxel_list[[timestep+1]]$proxel$prev[prox]]
      }
      
      if((timestep+1)%%2 == 0){
        if(SU[SU[,1]==timestep+1,2]=="Ua") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            URINE$Ua[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SU[SU[,1]==timestep+1,2]=="Ub") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            URINE$Ub[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SU[SU[,1]==timestep+1,2]=="Uc") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            URINE$Uc[proxel_list[[timestep+1]]$proxel$state[prox]]
      }
      
      if((timestep+1)%%3 == 0){
        if(SB[SB[,1]==timestep+1,2]=="Ba") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            BLOOD$Ba[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SB[SB[,1]==timestep+1,2]=="Bb") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            BLOOD$Bb[proxel_list[[timestep+1]]$proxel$state[prox]]
        if(SB[SB[,1]==timestep+1,2]=="Bc") proxel_list[[timestep+1]]$HnMM[prox,2] = proxel_list[[timestep+1]]$HnMM[prox,2] * proxel_list[[timestep+1]]$proxel$prob[prox] *
            BLOOD$Bc[proxel_list[[timestep+1]]$proxel$state[prox]]
      }
    }
    
    if(Print) print(paste("Time step:",timestep))
    if(Print) print(proxel_list[[timestep+1]])
    
    proxel_df[timestep+1,1]=timestep*delta
    for (v in 2:4) {
      if(length(proxel_list[[timestep+1]]$proxel$prob[proxel_list[[timestep+1]]$proxel$state == v-1]) != 0) {
        proxel_df[timestep+1,v] = round(sum(proxel_list[[timestep+1]]$proxel$prob[proxel_list[[timestep+1]]$proxel$state == v-1]),5)}
          else proxel_df[timestep+1,v] = 0
    }
    proxel_df[timestep+1,5]=round(1-sum(proxel_df[timestep+1,2:4]),5)
      
  }
  
  if(RETURN=="TREE") return(proxel_list)
  if(RETURN=="DF") return(proxel_df)
}
proxel_list2=Proxel_algo2(RG,startproxel,delta,max.ts,F,"TREE",FEVER2,BLOOD2,URINE2,patient1)

Compute the probability of each disease.

Prob_dis2=sum(proxel_list2[[13]]$HnMM$WKT)
#Prob_dis2>Prob_dis1

• Prob(Disease 1) = 1.228379^{-21}
• Prob(Disease 2) = 2.9354376^{-16}

Which disease is the patient most likely suffering from?

• Prob(Disease 1) < Prob(Disease 2) –> suffering from disease 2

We have further measurements from patients (2-5) that were infected by our original subject, so they must be suffering from the same disease. What does your program say to that claim?

pat2_A=Proxel_algo(RG,startproxel,delta,9,F,"TREE",FEVER,BLOOD,URINE,patient2)
pat2_B=Proxel_algo2(RG,startproxel,delta,9,F,"TREE",FEVER2,BLOOD2,URINE2,patient2)
pat3_A=Proxel_algo(RG,startproxel,delta,7,F,"TREE",FEVER,BLOOD,URINE,patient3)
pat3_B=Proxel_algo2(RG,startproxel,delta,7,F,"TREE",FEVER2,BLOOD2,URINE2,patient3)
pat4_A=Proxel_algo(RG,startproxel,delta,11,F,"TREE",FEVER,BLOOD,URINE,patient4)
pat4_B=Proxel_algo2(RG,startproxel,delta,11,F,"TREE",FEVER2,BLOOD2,URINE2,patient4)
pat5_A=Proxel_algo(RG,startproxel,delta,6,F,"TREE",FEVER,BLOOD,URINE,patient5)
pat5_B=Proxel_algo2(RG,startproxel,delta,6,F,"TREE",FEVER2,BLOOD2,URINE2,patient5)

sum(pat2_A[[10]]$HnMM$WKT,na.rm = T) < sum(pat2_B[[10]]$HnMM$WKT,na.rm = T)
sum(pat3_A[[8]]$HnMM$WKT,na.rm = T) < sum(pat3_B[[8]]$HnMM$WKT,na.rm = T)
sum(pat4_A[[12]]$HnMM$WKT,na.rm = T) < sum(pat4_B[[12]]$HnMM$WKT,na.rm = T)
sum(pat5_A[[7]]$HnMM$WKT,na.rm = T) < sum(pat5_B[[7]]$HnMM$WKT,na.rm = T)

• Patient 2:
  • Prob(D1) = 3.6309174^{-13}
  • Prob(D2) = 9.6647641^{-24}
• Patient 3:
  • Prob(D1) = 6.9146422^{-9}
  • Prob(D2) = 2.132833^{-5}
• Patient 4:
  • Prob(D1) = 1.3978813^{-16}
  • Prob(D2) = 1.7138522^{-12}
• Patient 5:
  • Prob(D1) = 1.1893694^{-7}
  • Prob(D2) = 3.1100656^{-4}
• except patient 2 all the others also suffering from disease 2