You are right that both the individual predictions \mathrm{E}(y|\eta=\eta^*) (for any choice of \eta^*) and the population predictions \mathrm{E}(y) are supposed to be positive when y models drug concentrations. However, E(y) is a complicated integral and typically we approximate it by linearizing the conditional expectation \mathrm{E}(y|\eta) since \mathrm{E}(y)=\mathrm{E}(\mathrm{E}(y|\eta)). So if we denote f(\eta) = \mathrm{E}(y|\eta) then we have

\mathrm{E}(y) = \mathrm{E}(\mathrm{E}(y|\eta)) = \mathrm{E}(f(\eta)) \approx \mathrm{E}(f(\eta^*) + f'(\eta^*)(\eta - \eta^*)) = f(\eta^*) - f'(\eta^*)\eta^*

i.e. the first order approximation of the population prediction is the individual prediction minus a first order term and this quantity might sometimes be negative even though \mathrm{E}(y) is positive.

Negative approximate population predictions suggest that a first-order approximation might be too rough and you might want to consider a Monte Carlo approximation of the integral instead. You can currently do that in Pumas per subject with

```
Pumas.epredict(fpm, subject, nsim)
```

where `fpm`

is your model fit and `nsim`

is the number of simulations used in the Monte Carlo integration. For a large enough `nsim`

you shouldn’t see any negative predictions. (Be warned, though, that this API will most likely change soon.)