The first step to explain the RxODE model plots is to setup:

- An RxODE model
- An event table
- A solved object

```
## RxODE 1.1.1 using 1 threads (see ?getRxThreads)
## no cache: create with `rxCreateCache()`
```

```
## Model from RxODE tutorial
m1 <- RxODE({
KA <- 2.94E-01
CL <- 1.86E+01
V2 <- 4.02E+01
Q <- 1.05E+01
V3 <- 2.97E+02
Kin <- 1
Kout <- 1
EC50 <- 200
## Added modeled bioavaiblity, duration and rate
fdepot <- 1
durDepot <- 8
rateDepot <- 1250
C2 <- centr / V2
C3 <- peri / V3
d/dt(depot) <- -KA * depot
f(depot) <- fdepot
dur(depot) <- durDepot
rate(depot) <- rateDepot
d/dt(centr) <- KA * depot - CL * C2 - Q * C2 + Q * C3
d/dt(peri) <- Q * C2 - Q * C3
d/dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff
eff(0) <- 1
})
ev <- et(timeUnits = "hr") %>%
et(amt = 10000, ii = 12, until = 24) %>%
et(seq(0, 24, length.out = 100))
s <- rxSolve(m1, ev)
```

The default plot in `RxODE`

is to plot time versus `state`

or calculated `lhs`

values:

`plot(s)`

Things to note:

First, these plots are actually `ggplot2`

plots and can be modified by the standard `ggplot2`

grammar of graphics

Also note that the `plot`

method for `rxSolve`

objects can currently use some of the arguments of the default plot method:

`log`

a character string which contains`"x"`

if the x axis is to be logarithmic,`"y"`

if the y axis is to be logarithmic and`"xy"`

or`"yx"`

if both axes are to be logarithmic. If available, this will use the`xgxr`

`xgx_scale_x_log10()`

and`xgx_scale_y_log10()`

functions.`xlab`

This is a bit different than simply supplying the x-label through`ggplot2`

, it will retain the unit information originally supplied in the event table. Also note, if available, this unit information is used will be used to have more meaningful ticks with`xgx_scale_x_time_units()`

`ylab`

This supplies the y label

The following are currently unsupported: a. `type`

b. `xlim`

c. `ylim`

d. `main`

e. `sub`

f. `ann`

g. `axes`

While this is very interesting, it is often more useful to subset plots based on values you are more interested in; For example, it is more interesting in this plot to `eff`

and `C2`

. To do this, you simply specify the interesting items after the plot function. For instance:

`plot(s, C2, eff)`

Semi-log plots of PK concentrations are very common; To do this you simply need to use `log="y"`

`plot(s, C2, log="y", ylab="Concentration")`

`## Warning: Transformation introduced infinite values in continuous y-axis`

This uses the `xgx_scale_y_log10()`

when available to more clearly show the semi-log nature of the plot.

If you have multi-subject plots you can easily plot applying the same principles as above. To illustrate this plot, lets expand an event table to include a plot of 4 subjects with lognormal random variability between subjects;

```
# Setup the new problem
m2 <- RxODE({
KA <- 2.94E-01
CL <- 1.86E+01 * exp(eta.Cl)
V2 <- 4.02E+01
Q <- 1.05E+01
V3 <- 2.97E+02
Kin <- 1
Kout <- 1
EC50 <- 200
## Added modeled bioavaiblity, duration and rate
fdepot <- 1
durDepot <- 8
rateDepot <- 1250
C2 <- centr / V2
C3 <- peri / V3
d / dt(depot) <- -KA * depot
f(depot) <- fdepot
dur(depot) <- durDepot
rate(depot) <- rateDepot
d / dt(centr) <- KA * depot - CL * C2 - Q * C2 + Q * C3
d / dt(peri) <- Q * C2 - Q * C3
d / dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff
eff(0) <- 1
})
# Add the between subject variability:
omega <- lotri(eta.Cl ~ 0.4^2)
# Create the event table
ev <- et(timeUnits = "hr") %>%
et(amt = 10000, until = units::set_units(3, days), ii = 12) %>% # loading doses
et(seq(0, 48, length.out = 200)) %>%
et(id = 1:4)
s <- rxSolve(m2, ev, omega=omega)
```

Once that is complete, you may plot it by the same method:

`plot(s, C2, eff)`

Notice that this is colored by each individual and labeled with a legend.

If you are only interested in the concentration, it produces a similar plot:

`plot(s, C2, log="y", ylab="Concentration")`

`## Warning: Transformation introduced infinite values in continuous y-axis`

Notice that this plot each individual is labeled by an attached id to let you know where the individual comes from. This is done by `ggrepel`

if it is available, otherwise a legend is retained.

This can get a bit much when there are many subjects that are solved (in this example lets simulate 100)

```
ev <- et(timeUnits = "hr") %>%
et(amt = 10000, until = units::set_units(3, days), ii = 12) %>% # loading doses
et(seq(0, 48, length.out = 200)) %>%
et(id = 1:100) # 100 subjects
s <- rxSolve(m2, ev, omega=omega)
plot(s, C2, log="y", ylab="Concentration")
```

`## Warning: Transformation introduced infinite values in continuous y-axis`

In this case, all the individuals are put on the plot in transparent greyscale and plot on the same pane. This allows the places where more subjects are present to be darker.

The number of individuals where the plots switch from legend to grayscale is controlled by changing the options for `RxODE.spaghetti`

. ie. `options(RxODE.spaghetti=7)`

. 7 individuals is the default value when the plotting changes from individual to grayscale spaghetti plots.

You can also create a confidence interval of these simulations with `confint`

:

`s2 <- confint(s, parm="C2")`

`## ! in order to put confidence bands around the intervals, you need at least 2500 simulations`

`## summarizing data...done`

And plot this with `plot`

`plot(s2)`