## ---- echo=FALSE, message=FALSE, warning=FALSE---------------------------
library(mlxR)
theme_set(theme_bw())
htmltools::includeHTML("../macros.txt")
knitr::opts_chunk$set(warning=FALSE, message=FALSE)

## ------------------------------------------------------------------------
project.file <- 'monolixRuns/theophylline_project.mlxtran'

## ------------------------------------------------------------------------
sim.res1  <- simulx(project = project.file)

## ------------------------------------------------------------------------
print(ggplot(data=sim.res1$y1) + 
        geom_point(aes(x=time, y=y1, colour=id)) +
        geom_line(aes(x=time, y=y1, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))


## ------------------------------------------------------------------------
sim.param <- c(b=0)
out  <- list(name = 'Cc', time = seq(0, 25, by=0.1))
sim.res2  <- simulx(project   = project.file,
                    output    = out,
                    parameter = sim.param)

print(ggplot() + 
        geom_point(data=sim.res2$y1, aes(x=time, y=y1, colour=id)) +
        geom_line(data=sim.res2$Cc, aes(x=time, y=Cc, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))

## ------------------------------------------------------------------------
N <- 50
sim.res3  <- simulx(project = project.file,
                    group   = list(size = N))

print(ggplot(data=sim.res3$y1) + 
        geom_point(aes(x=time, y=y1, colour=id)) +
        geom_line(aes(x=time, y=y1, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration") +
        theme(legend.position="none"))


## ------------------------------------------------------------------------
out  <- list(name = 'y1', time = (0:12))

sim.res4  <- simulx(project = project.file,
                    output = out)

print(ggplot(data=sim.res4$y1) + 
        geom_point(aes(x=time, y=y1, colour=id)) +
        geom_line(aes(x=time, y=y1, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))

## ------------------------------------------------------------------------
out1  <- list(name = 'y1', time = seq(0, 24, by=2))
out2  <- list(name = 'Cc', time = seq(0, 24, by=0.1))
out3  <- list(name = c('V', 'Cl', 'WEIGHT'))

sim.res5  <- simulx(project = project.file,
                    output  = list(out1, out2, out3))

print(sim.res5$parameter)

print(ggplot() + 
        geom_point(data=sim.res5$y1,aes(x=time, y=y1, colour=id)) +
        geom_line(data=sim.res5$Cc,aes(x=time, y=Cc, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))

## ------------------------------------------------------------------------
adm   <- list(time = c(0,6), amount = c(320, 320))

sim.res6  <- simulx(project   = project.file,
                    treatment = adm,
                    output    = list(out1, out2))

print(ggplot() + 
        geom_point(data=sim.res6$y1,aes(x=time, y=y1, colour=id)) +
        geom_line(data=sim.res6$Cc,aes(x=time, y=Cc, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))

## ------------------------------------------------------------------------
N   	<- 100
weight <- list( name     = 'WEIGHT', 
                colNames = c('id', 'WEIGHT'),
                value    = cbind(c(1:N),c(rep(50, N/2), rep(90, N/2))))

adm   <- list(time = 0, amount = 500)

sim.res7 <- simulx(project = project.file, 
                   output = list(out1, out2),
                   treatment = adm,
                   parameter = weight)

sim.res7$Cc$weight <- 50
sim.res7$Cc$weight[as.numeric(sim.res7$Cc$id)>N/2] <-90
sim.res7$Cc$weight <- as.factor(sim.res7$Cc$weight)

print(ggplot() + 
        geom_line(data=sim.res7$Cc,aes(x=time, y=Cc, group=id, colour=weight)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))

## ------------------------------------------------------------------------
library(plyr, warn.conflicts = FALSE)

#-- Define the number of patients to simulate
N <- 100

#-- Generate individual weights by sampling a lognormal distribution
sim.weight <- data.frame(id=1:N, WEIGHT=rlnorm(N, log(70), 0.2))

outy  <- list(name = 'y1', time = 120)
Dose.amount <- c(50, 100, 250, 500, 1000)
Dose.times  <- seq(from = 0, to = 120, by = 12)

#-- Use the same patients for each of the dose levels
s  <- list(seed = 123456)

#-- run simulx with each dose level
sim.data <- NULL
for(n.c in 1:length(Dose.amount)) {
  adm   <- list(time = Dose.times, amount = Dose.amount[n.c])

	tmp <- simulx(project   = project.file,
	              output    = outy,
	              treatment = adm,
	              parameter = sim.weight,
	              settings  = s)
	
	tmp2 			<- tmp$y1
	tmp2['Dose'] 	<- Dose.amount[n.c]
	sim.data <- rbind(sim.data, tmp2)
}

#-- Compute statistics
sim.data.stat <- ddply(sim.data, .(time, Dose), summarize,
	median = median(y1),
	p05  = quantile(y1, 0.05),
	p95  = quantile(y1, 0.95),
	Dose = Dose[1])

print(ggplot(data=sim.data.stat) +
  geom_line(aes(x=Dose, y=median)) +
  geom_point(aes(x=Dose, y=median)) +
  geom_ribbon(aes(x=Dose, ymin=p05, ymax=p95), alpha = 0.3) +
  scale_x_continuous("Dose") + scale_y_continuous("Concentration"))

## ------------------------------------------------------------------------
#-- Define the output as summary statistics of the simulated concentration at t=120
out.y <- list(name='y1', time=120, FUN = "quantile", probs = c(0.05, 0.50, 0.95))

Dose.amount <- c(100, 250)
Dose.times  <- seq(from = 0, to = 120, by = 12)

#-- Define the number of trial simulations
N.trial <- 50

#-- Define the number of patients to simulate
vN <- c(10, 30, 50, 100)

sim.data <- NULL
for(n.c in 1:length(Dose.amount)){	
  adm   <- list(time = Dose.times, amount = Dose.amount[n.c])
  for(N in vN){

    tmp1 <- simulx(project   = project.file,
                   nrep = N.trial,
                   output    = out.y,
                   treatment = adm,
                   group = list(size=N))
    
    tmp2 <- tmp1$y1
    tmp2['Dose'] <- factor(Dose.amount[n.c])
    tmp2['N.patients']<- N
    sim.data <- rbind(sim.data, tmp2)   
  }
}
head(sim.data)

# Distribution of the median
print(ggplot(data=sim.data) + geom_point(aes(x=N.patients, y=y1.p50, color=Dose)))

## ------------------------------------------------------------------------
WEIGHT = 70
sim.param <- c(omega_V=0, omega_ka=0, omega_Cl=0, WEIGHT=WEIGHT)  # define a typical patient

out.Cc  <- list(name = 'Cc', time = seq(0,100,by=0.2))

Dose.times  <- seq(from = 0, to = 72, by = 6)
nb.times <- length(Dose.times)

Dose.perkg <- c(1.5, 2, 2.5)
nb.group <- length(Dose.perkg)

g <- list()
for (k in (1:nb.group)) {
  trt.k <- list(time = Dose.times, amount = round(Dose.perkg[k]*WEIGHT))
  g[[k]] <- list(treatment = trt.k, size=1)
}

res <- simulx(project   = project.file, 
              group     = g, 
              parameter = sim.param,
              output    = out.Cc)
pl <- ggplot(data=res$Cc, ) + geom_line(aes(x=time, y=Cc, colour=group)) +  scale_y_continuous("predicted concentration") +
  scale_colour_discrete(name  ="Dose", labels=c("1.5 mg/kg", "2 mg/kg", "2.5 mg/kg")) 
print(pl)

## ------------------------------------------------------------------------
N <- 200
sim.weight <- data.frame(id=1:N, WEIGHT=rlnorm(N, log(70), 0.2))

out.Cc  <- list(name = 'Cc', time = 0:100)
Dose.times  <- seq(from = 0, to = 72, by = 6)
nb.times <- length(Dose.times)
Dose.perkg <- c(1, 2)
nb.group <- length(Dose.perkg)
g <- list()
for (k in (1:nb.group)) {
  trt.k <- data.frame(id     = rep(1:N,each=nb.times),
                      time   = rep(Dose.times,N),
                      amount = rep(round(Dose.perkg[k]*sim.weight$WEIGHT),each=nb.times))
  g[[k]] <- list(treatment = trt.k, size=N)
}

res <- simulx(project   = project.file, 
              group     = g, 
              parameter = sim.weight,
              output    = out.Cc)
pl <- prctilemlx(res$Cc, number = 2, level = 90, labels = c("1 mg/kg", "2 mg/kg")) 
print(pl)

