# example script for using nlmixr
# Created by Nathan Teuscher
# February 2025

# clear all objects
rm(list = ls())

# Load packages
library(nlmixr2)
library(tidyverse)

# List datsets in nlmixr2data package
data(package = "nlmixr2data")

# use this to view the sample analysis dataset in R Studio
View(theo_sd)

# plot data by subject

# filter to only observations and plot by ID
ggplot(data = theo_sd %>% filter(EVID == 0), aes(x = TIME, y = DV, color = as.factor(ID))) +
  geom_point() + geom_line()

# same plot with logarithmic y-axis
ggplot(data = theo_sd %>% filter(EVID == 0), aes(x = TIME, y = DV, color = as.factor(ID))) +
  geom_point() + geom_line() +
  scale_y_log10()

## The basic model consists of an ini block that has initial estimates
one.compartment <- function() {
  ini({
    tka <- log(1.57); label("Ka")
    tcl <- log(2.72); label("Cl")
    tv <- log(31.5); label("V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  # and a model block with the error specification and model specification
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
    cp ~ add(add.sd)
  })
}

# Fit the data to the model using SAEM
fit <- nlmixr2(one.compartment, theo_sd, est = "saem", saemControl(print=0))

# Print the model output to the Console
print(fit)

# Fit the data to the model using FOCEI
fit2 <- nlmixr2(one.compartment, theo_sd, est = "focei", foceiControl(print = 0))

# Print the model output to the Console
print(fit2)

# Generate diagnostic plots for SAEM model
pdf(file = "fit_plots.pdf")
plot(fit)
dev.off()

# Generate diagnostic pltos for FOCEI model
pdf(file = "fit2_plots.pdf")
plot(fit2)
dev.off()