Since I am reading the Joel and Jill's textbook, I tried to convert the code example in the book to the R-based code.

Here it is:

# Make a same model with the p.31 of Joel and Jill's book (Intro to popPK/PD Anal. with NONMEM)

####     The NONMEM Code in the book
####     
####     $PROBLEM Example NM-TRAN Control Stream for a 1-cmt Model
####     $DATA /home/pathtomydataset/datafilename.dat  IGNORE=#
####     $INPUT ID TIME AMT DV CMT EVID MDV
####     $SUBROUTINES ADVAN2 TRANS2
####     $PK
####     
####     TVCL = THETA(1)
####     CL = TVCL * EXP(ETA(1))
####     TVV = THETA(2)
####     V = TVV * EXP(ETA(2))
####     TVKA = THETA(3)
####     KA = TVKA
####     
####     ; scale predictions based on dose (mg) and Cp (ng/mL)
####     S2 = V/1000
####     
####     $ERROR
####     
####     Y = F*(1 + EPS(1))
####     
####     $THETA    (0, 3.9)    ; Clearance (L/h)
####               (0, 52.1)   ; Volume (L)
####               (0, 0.987)  ; Absorption Rate (1/h)
####     
####     $OMEGA    (0.09)      ; IIV on CL (use 30%CV as init est)
####               (0.16)      ; IIV on V  (use 40%CV as init est)
####     
####     $SIGMA    (0.04)      ; RV (use 20%CV as init est)
####     
####     ; Use conditional estimation with interaction
####     $EST METHOD=1 INTER MAXEVAL=9999 PRINT=5 MSFO=example.msf
####     $COVAR PRINT=E
####     
####     $TABLE ID IME AMT CMT ONEHEADER NOPRINT FILE=example.tbl

# Loading NLMIXR2 library
library(nlmixr2)
library(xpose)
library(xpose.nlmixr2)
library(ggplot2)

# Load data: Here, I will use the example data. Instead of the sample, someone might be import data using read_csv method of Tidyverse library.
# Header of the imported data would be specified by header element of read_csv method. In this sample case, the header is already marked.
data <- theo_sd # This is not the same data with the one of textbook author. So, this code would not results in a well fitted parameters. This is just a demo code for the converting.

# I didn't found $SUBROUTINE alternatives in the NLMIXR2 manual. linCmt() pseudo-function is doing similar things.

# Defining the PK model
one.cmt <- function() {
  ini({
    TVCL <- c(0, 3.9); label("CL")       # Clearance (L/h); c(lower bound, initial estimate)
    TVV  <- c(0, 52.1); label("V")       # Volume (L); c(lower bound, initial estimate)
    TVKA <- c(0, 0.987); label("Ka")     # Absorption Rate (1/h); c(lower bound, initial estimate)
      
    Eta.CL ~ 0.09             # IIV on CL (use 30%CV as init est)
    Eta.V  ~ 0.16             # IIV on V (use 40%CV as init est)
    
    Eps.Cp <- 0.04            # RV (use 20%CV as init est)
  })
  model({
    CL <- TVCL * exp(Eta.CL)
    V  <- TVV  * exp(Eta.V)
    KA <- TVKA
    
    # scale predictions based on dose (mg) and Cp (ng/mL)
    S2 <- V/1000
    
    linCmt() ~ prop(Eps.Cp)
  })
}

fitF <- nlmixr2(one.cmt, theo_sd, "focei", control=list(print=0))
#fitS <- nlmixr(one.cmt, theo_sd, "saem")

fitF  <- fitF %>% addCwres() %>% addNpde()
plot(fitF)

################################################################################
## Xpose plots; Need to print otherwise running a script won't
## show xpose plots
################################################################################
xpdb <- xpose_data_nlmixr(fitF) ## first convert to nlmixr object

print(dv_vs_pred(xpdb) +
        ylab("Observed Concentrations (Unknown unit)") +
        xlab("Population Predicted Concentrations (Unknown unit)"))

print(dv_vs_ipred(xpdb) +
        ylab("Observed Concentrations (Unknown unit)") +
        xlab("Individual Predicted Concentrations (Unknown unit)"))

print(res_vs_pred(xpdb) +
        ylab("Conditional Weighted Residuals") +
        xlab("Population Predicted Concentrations (Unknown unit)"))

print(res_vs_idv(xpdb) +
        ylab("Conditional Weighted Residuals") +
        xlab("Time (Unknown unit)"))

print(absval_res_vs_idv(xpdb, res = 'IWRES') +
        ylab("Individual Weighted Residuals") +
        xlab("Time (Unknown unit)"))

print(absval_res_vs_pred(xpdb, res = 'IWRES') +
        ylab("Individual Weighted Residuals") +
        xlab("Population Predicted Concentrations (Unknown unit)"))

print(ind_plots(xpdb, nrow=3, ncol=4) +
        ylab("Predicted and Observed concentrations (Unknown unit)") +
        xlab("Time (Unknown unit)"))

print(res_distrib(xpdb) +
        ylab("Density") +
        xlab("Conditional Weighted Residuals"))

################################################################################
##Visual Predictive Checks
################################################################################
vpcPlot(fitF,n=1000, show=list(obs_dv=T),
        bins = c(-0.5, 0, 0.1, 0.4, 0.7, 1.5, 3, 4.5, 6, 8, 10, 15, 33),
        pred_corr = T, # T for pcVPC,
        log_y = F, # y-scale
        smooth = F, # Smoooth the VPC (connect bin mid points)
        ylab = "Concentrations (Unknown unit)", xlab = "Time (Unknown unit)")