Hello @Vaibhavdixit02

This will certainly allow you to get a scatter plot of the observations overlayed on the simulation results. There is another way to do this using only `sim_plot`

, so I will share an example of how you could do that.

We could use the `nlme_sample`

dataset from `PharmaDatasets`

for this example. You can load it and convert it to a population with the following code.

```
using Pumas, PumasUtilities
using PharmaDatasets # Get dataset
using CairoMakie # Customize plots and show legend
# Get dataset
pkdata = dataset("nlme_sample")
# Create population with read_pumas
population = read_pumas(
pkdata;
id = :ID,
time = :TIME,
amt = :AMT,
covariates = [:WT, :AGE, :SEX, :CRCL, :GROUP],
observations = [:DV],
cmt = :CMT,
evid = :EVID,
rate = :RATE,
)
```

Now we can define a 1 compartment model to be used in our simulation:

```
# Create the model
model = @model begin
@metadata begin
desc = "base model: 1comp"
timeu = u"hr"
end
@param begin
"Clearance (L/hr)"
tvcl ∈ RealDomain(; lower = 0.0001)
"Volume (L)"
tvvc ∈ RealDomain(; lower = 0.0001)
"""
- ΩCL
- ΩVc
"""
Ω ∈ PDiagDomain(2) #covariance matrix
"Additive RUV"
σ²_add ∈ RealDomain(; lower = 0.0001) # variance
"Proportional RUV"
σ²_prop ∈ RealDomain(; lower = 0.0001) # variance
end
@random begin
η ~ MvNormal(Ω)
end
@covariates WT AGE SEX CRCL GROUP
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
end
@dynamics Central1
@derived begin
CONC := @. Central / Vc
"""
DrugY Concentration (ng/mL)
"""
DV ~ @. Normal(CONC, sqrt(CONC^2 * σ²_prop + σ²_add)) # using variance
end
end
```

Also, we should define the initial values for the model parameters. In this case, we can use the following.

```
params = (;
tvvc = 5,
tvcl = 0.2,
Ω = Diagonal([0.01, 0.01]),
σ²_add = 0.01,
σ²_prop = 0.01
)
```

Let’s choose one of the subjects from our population. I will be choosing one using the `id`

(say, subject 8), but you could use different criteria depending on your needs.

```
subject = filter(i -> i.id == "8", population)
```

Now that we have our model, initial parameters, and subject, we can run the simulation using `simobs`

.

```
simulation = simobs(model, subject, params)
```

Having done that, we can generate the plot using `sim_plot`

.

```
sim_plot(
model,
simulation;
observations=[:DV],
axis=(;
xlabel="Time (hr)",
ylabel="Observed/Predicted \n Concentration (ng/mL)"
),
color=:red
)
```

This will return a scatter plot of the observations and the simulated values as a red line on top. However, this plot will not have a legend, which I think would be very useful for visualizing the results.

To get a legend, we will take advantage of the fact that `Makie.jl`

returns a tuple containing the figure, axis, and plot objects. We will store these values using the following syntax.

```
fig, ax, plt = sim_plot(
model,
simulation;
observations=[:DV],
axis=(;
xlabel="Time (hr)",
ylabel="Observed/Predicted \n Concentration (ng/mL)"
),
color=:red
)
```

Now we can display the legend by using `axislegend`

and then inspect `fig`

to see the final result. We also use the keyword argument `position`

to tell `Makie`

to place the legend on the upper left portion of the figure.

```
axislegend(ax, position=:lt)
fig
```

This will produce the following result: