# Item response theory (IRT) pharmacodynamic model for morphine
# analgesia in mechanically ventilated preterm neonates undergoing
# endotracheal suctioning. Extracted from DDMORE Foundation Model
# Repository entry DDMODEL00000247 (Valitalo et al. 2017,
# CPT Pharmacometrics Syst Pharmacol 6(4):239-248,
# doi:10.1002/psp4.12156).
#
# Source bundle note: the DDMORE bundle for DDMODEL00000247 ships an
# Output_real_OriginalModelCode.lst that is for a different model
# (Midazolam 2-cmt popPK in critically ill children, fit to
# Data20150701.csv). The bundle's .mod is the correct Valitalo 2017
# IRT morphine PD model; its inline $THETA / $OMEGA / $SIGMA values
# match Table 2 of the publication (paper PDF on disk at
# /home/bill/gitlab/nlmixr2lib_ingestion/from_people/literature_2018_search/
# Vlitalo_2017_Morphine_Pharmacodynamics_in_Mechanicall_966c39.pdf).
# Parameter VALUES below come from Table 2 of the paper and are
# cross-verified against the .mod's inline initials, which equal the
# final estimates. The mismatched .lst is documented in the vignette
# Errata.

Valitalo_2017_morphine <- function() {
  description <- paste(
    "Item response theory (IRT) pharmacodynamic model for morphine analgesia in mechanically ventilated preterm neonates undergoing endotracheal suctioning.",
    "A latent pain variable is driven by suctioning state (pre/during/after), a linear morphine plasma-concentration effect, and a linear study-time effect.",
    "The 8 ordinal pain-item scores (COMFORT-B alertness/calmness/respiratory/body movement/facial tension; PIPP brow bulge/eye squeeze/nasolabial furrow; NIPS total) and 1 continuous VAS score map onto the latent variable via a graded-response IRT model.",
    "Morphine plasma concentration is supplied as a time-varying covariate CP_MORPH_NGML from an upstream PK model (Knibbe 2009; see modellib('Knibbe_2009_morphine')).",
    "The IRT discrimination and difficulty parameters were estimated separately in a graded-response sub-model and are FIX'd here per the source NONMEM control stream.",
    "Population: 140 preterm neonates (mean weight 1.4 kg, mean postmenstrual age 211 days)."
  )
  reference <- paste(
    "Valitalo PA, Krekels EHJ, van Dijk M, Simons SHP, Tibboel D, Knibbe CAJ (2017).",
    "Morphine pharmacodynamics in mechanically ventilated preterm neonates undergoing endotracheal suctioning.",
    "CPT Pharmacometrics Syst Pharmacol 6(4):239-248.",
    "doi:10.1002/psp4.12156.",
    "Upstream PK structure adapted from Knibbe et al. (2009)",
    "Clin Pharmacokinet 48(6):371-385,",
    "doi:10.2165/00003088-200948060-00003;",
    "see modellib('Knibbe_2009_morphine').",
    "DDMORE Foundation Model Repository: DDMODEL00000247.",
    sep = " "
  )
  vignette <- "Valitalo_2017_morphine"
  ddmore_id <- "DDMODEL00000247"
  replicate_of <- NULL
  units <- list(
    time = "day",
    dosing = "(none; PD-only IRT model fed by an external morphine PK trajectory)",
    concentration = "(unitless latent pain score; derived VAS score on 0-10 cm scale; morphine plasma concentration CP_MORPH_NGML supplied as covariate in ng/mL)"
  )

  covariateData <- list(
    CP_MORPH_NGML = list(
      description        = "Instantaneous morphine plasma concentration (ng/mL) at the time of each pain assessment, supplied as a time-varying covariate from an upstream PK model.",
      units              = "ng/mL",
      type               = "continuous",
      reference_category = NULL,
      notes              = paste(
        "Time-varying per event row. In Valitalo 2017 this column was generated by simulating individual morphine concentrations from a Knibbe 2009 (doi:10.2165/00003088-200948060-00003) popPK model for the same population; see modellib('Knibbe_2009_morphine').",
        "Source NONMEM column name is CP.",
        "Enters the latent-pain equation as a linear reducer:",
        "pain = pain_state - e_morph_pain * CP_MORPH_NGML + e_time_pain * time.",
        "Observed individual predicted concentrations in Valitalo 2017 were mostly within 0-60 ng/mL (Figure 2a)."
      ),
      source_name        = "CP"
    ),
    MOMENT = list(
      description        = "Endotracheal suctioning state at the time of the pain assessment: 1 = before suctioning, 2 = during suctioning, 3 = after suctioning.",
      units              = "(categorical, 3 levels)",
      type               = "categorical",
      reference_category = "1 (before suctioning)",
      notes              = paste(
        "Time-varying per assessment. Selects between the three baseline pain typical values presuct / suct / aftsuct in the model body.",
        "Source NONMEM column name is MOMENT, preserved here."
      ),
      source_name        = "MOMENT"
    ),
    ITEM = list(
      description        = "Identifier of the pain-assessment item being scored at each observation row. Source NONMEM coding: 1 = COMFORT-B alertness; 2 = COMFORT-B calmness/agitation; 3 = COMFORT-B respiratory response; 5 = COMFORT-B body movement; 7 = COMFORT-B facial tension; 12 = VAS (cm, range 0-10); 25 = PIPP brow bulge; 26 = PIPP eye squeeze; 27 = PIPP nasolabial furrow; 28 = NIPS total.",
      units              = "(categorical, 9 active levels)",
      type               = "categorical",
      reference_category = NULL,
      notes              = paste(
        "Selects which IRT discrimination/difficulty set governs the item-grade probabilities at each observation row.",
        "The COMFORT-B muscle-tension item was excluded by Valitalo 2017 because it could not be assessed from video recordings."
      ),
      source_name        = "ITEM"
    ),
    OBSTYPE = list(
      description        = "Observer type for VAS observations: 1 = investigator (video-based assessment), 2 = bedside nurse.",
      units              = "(categorical, 2 levels)",
      type               = "categorical",
      reference_category = "1 (video)",
      notes              = paste(
        "Applies only to ITEM = 12 (VAS).",
        "Selects between the two VAS difficulty / discrimination THETAs (49 vs 50, 51 vs 52) and between the two additive residual error SDs (addSd_vas_video vs addSd_vas_bedside)."
      ),
      source_name        = "OBSTYPE"
    )
  )

  population <- list(
    species        = "human",
    n_subjects     = 140L,
    n_studies      = 1L,
    age_range      = "postmenstrual age mean 211 days (SD 24.4 days); 133/140 enrolled at postnatal age 1 day, remaining 7 at postnatal age 2-3 days",
    age_median     = "postmenstrual age ~211 days (~30 weeks PMA equivalent)",
    weight_range   = "mean birthweight 1400 g (SD 725 g)",
    weight_median  = "~1.4 kg",
    sex_female_pct = NA_real_,
    race_ethnicity = NULL,
    disease_state  = "Preterm neonates undergoing mechanical ventilation and endotracheal suctioning",
    dose_range     = paste(
      "Morphine arm (n=71): loading dose 100 mcg/kg IV, followed by continuous infusion 10 mcg/kg/h.",
      "Placebo arm (n=69): saline.",
      "Both arms allowed open-label rescue morphine (50 mcg/kg bolus + 5-10 mcg/kg/h infusion) as needed.",
      "The current PD model treats morphine concentration as a covariate (CP_MORPH_NGML) supplied by an upstream PK model, so dosing does not enter the PD model directly."
    ),
    regions        = "Netherlands (Erasmus MC Sophia Children's Hospital, Rotterdam; Leiden University)",
    notes          = paste(
      "Re-analysis of the Simons 2003 (JAMA 290:2419-2427, doi:10.1001/jama.290.18.2419) randomized double-blind double-dummy trial.",
      "29/140 neonates were 'small for gestational age' (birthweight < 2.5th percentile by gestational age).",
      "Ten patients of the original 150 in Simons 2003 were excluded because they were extubated early or had other reasons for early study termination.",
      "16,257 item-level COMFORT-B / PIPP / NIPS / VAS recordings were analysed across the 140 neonates.",
      "Demographics from Valitalo 2017 Methods 'Study design' paragraph and Table 1.",
      "Sex split not reported in the available trimmed source."
    )
  )

  ini({
    # --------------------------------------------------------------
    # IRT graded-response sub-model parameters (FIX'd in source).
    #
    # All discrimination (discr_*) and difficulty (diff_*_k) parameters
    # below come from the FIX'd $THETA(1..52) block of the bundle .mod
    # (Executable_OriginalModelCode.mod). These were estimated in a
    # prior graded-response sub-model (Valitalo 2017 Methods "Graded
    # response model for pain"; Supplementary Electronic Material S2)
    # and held fixed for the PD modelling stage. The values are not
    # individually tabulated in the paper main text -- they are in the
    # paper's online supplement S2 which is not on disk; the .mod
    # carries them as the authoritative copy.
    #
    # Discrimination THETAs are positive (`(0,value)` in $THETA);
    # difficulty THETAs can be negative (lowest-grade boundary is
    # often negative). All wrapped in fixed().
    #
    # --- Item 1: COMFORT-B alertness (5 grades, indices 1..5) ---
    discr_cmt_alert  <- fixed(3.6)
    label("IRT discrimination for COMFORT-B alertness item (unitless)")
    # .mod $THETA(1) = (0, 3.6) FIX

    diff_cmt_alert_1 <- fixed(-0.761)
    label("IRT difficulty 1 (boundary 1->2) for COMFORT-B alertness")
    # .mod $THETA(2) = -0.761 FIX
    diff_cmt_alert_2 <- fixed(1.39)
    label("IRT difficulty 2 increment for COMFORT-B alertness")
    # .mod $THETA(3) = (0, 1.39) FIX
    diff_cmt_alert_3 <- fixed(0.887)
    label("IRT difficulty 3 increment for COMFORT-B alertness")
    # .mod $THETA(4) = (0, 0.887) FIX
    diff_cmt_alert_4 <- fixed(1.01)
    label("IRT difficulty 4 increment for COMFORT-B alertness")
    # .mod $THETA(5) = (0, 1.01) FIX

    # --- Item 2: COMFORT-B calmness/agitation (5 grades) ---
    discr_cmt_calm   <- fixed(4.43)
    label("IRT discrimination for COMFORT-B calmness/agitation")
    # .mod $THETA(6) = (0, 4.43) FIX
    diff_cmt_calm_1  <- fixed(-0.162)
    label("IRT difficulty 1 for COMFORT-B calmness/agitation")
    # .mod $THETA(7) = -0.162 FIX
    diff_cmt_calm_2  <- fixed(0.992)
    label("IRT difficulty 2 increment for COMFORT-B calmness/agitation")
    # .mod $THETA(8) = (0, 0.992) FIX
    diff_cmt_calm_3  <- fixed(0.808)
    label("IRT difficulty 3 increment for COMFORT-B calmness/agitation")
    # .mod $THETA(9) = (0, 0.808) FIX
    diff_cmt_calm_4  <- fixed(0.533)
    label("IRT difficulty 4 increment for COMFORT-B calmness/agitation")
    # .mod $THETA(10) = (0, 0.533) FIX

    # --- Item 3: COMFORT-B respiratory response (5 grades) ---
    discr_cmt_resp   <- fixed(2.91)
    label("IRT discrimination for COMFORT-B respiratory response")
    # .mod $THETA(11) = (0, 2.91) FIX
    diff_cmt_resp_1  <- fixed(-2.67)
    label("IRT difficulty 1 for COMFORT-B respiratory response")
    # .mod $THETA(12) = -2.67 FIX
    diff_cmt_resp_2  <- fixed(3.2)
    label("IRT difficulty 2 increment for COMFORT-B respiratory response")
    # .mod $THETA(13) = (0, 3.2) FIX
    diff_cmt_resp_3  <- fixed(1.72)
    label("IRT difficulty 3 increment for COMFORT-B respiratory response")
    # .mod $THETA(14) = (0, 1.72) FIX
    diff_cmt_resp_4  <- fixed(1.09)
    label("IRT difficulty 4 increment for COMFORT-B respiratory response")
    # .mod $THETA(15) = 1.09 FIX

    # --- Item 5: COMFORT-B body movement (5 grades) ---
    discr_cmt_body   <- fixed(2.89)
    label("IRT discrimination for COMFORT-B body movement")
    # .mod $THETA(16) = (0, 2.89) FIX
    diff_cmt_body_1  <- fixed(-1.74)
    label("IRT difficulty 1 for COMFORT-B body movement")
    # .mod $THETA(17) = -1.74 FIX
    diff_cmt_body_2  <- fixed(0.777)
    label("IRT difficulty 2 increment for COMFORT-B body movement")
    # .mod $THETA(18) = (0, 0.777) FIX
    diff_cmt_body_3  <- fixed(1.94)
    label("IRT difficulty 3 increment for COMFORT-B body movement")
    # .mod $THETA(19) = (0, 1.94) FIX
    diff_cmt_body_4  <- fixed(1.94)
    label("IRT difficulty 4 increment for COMFORT-B body movement")
    # .mod $THETA(20) = (0, 1.94) FIX

    # --- Item 7: COMFORT-B facial tension (5 grades) ---
    discr_cmt_face   <- fixed(4.04)
    label("IRT discrimination for COMFORT-B facial tension")
    # .mod $THETA(21) = (0, 4.04) FIX
    diff_cmt_face_1  <- fixed(-1.49)
    label("IRT difficulty 1 for COMFORT-B facial tension")
    # .mod $THETA(22) = -1.49 FIX
    diff_cmt_face_2  <- fixed(1.76)
    label("IRT difficulty 2 increment for COMFORT-B facial tension")
    # .mod $THETA(23) = (0, 1.76) FIX
    diff_cmt_face_3  <- fixed(1.43)
    label("IRT difficulty 3 increment for COMFORT-B facial tension")
    # .mod $THETA(24) = (0, 1.43) FIX
    diff_cmt_face_4  <- fixed(0.334)
    label("IRT difficulty 4 increment for COMFORT-B facial tension")
    # .mod $THETA(25) = (0, 0.334) FIX

    # --- Item 25: PIPP brow bulge (4 grades, indices 0..3) ---
    discr_pipp_brow  <- fixed(2.21)
    label("IRT discrimination for PIPP brow bulge")
    # .mod $THETA(26) = (0, 2.21) FIX
    diff_pipp_brow_1 <- fixed(0.467)
    label("IRT difficulty 1 for PIPP brow bulge")
    # .mod $THETA(27) = 0.467 FIX
    diff_pipp_brow_2 <- fixed(1.08)
    label("IRT difficulty 2 increment for PIPP brow bulge")
    # .mod $THETA(28) = (0, 1.08) FIX
    diff_pipp_brow_3 <- fixed(0.937)
    label("IRT difficulty 3 increment for PIPP brow bulge")
    # .mod $THETA(29) = (0, 0.937) FIX
    # NB: .mod $THETA(30) = 100000 FIX is a sentinel placeholder for
    # the unused 5th-grade boundary on a 4-grade item; not exposed
    # here because the PIPP brow item has only 4 grades (0..3).

    # --- Item 26: PIPP eye squeeze (4 grades, indices 0..3) ---
    discr_pipp_eye   <- fixed(2.14)
    label("IRT discrimination for PIPP eye squeeze")
    # .mod $THETA(31) = (0, 2.14) FIX
    diff_pipp_eye_1  <- fixed(0.343)
    label("IRT difficulty 1 for PIPP eye squeeze")
    # .mod $THETA(32) = 0.343 FIX
    diff_pipp_eye_2  <- fixed(1.17)
    label("IRT difficulty 2 increment for PIPP eye squeeze")
    # .mod $THETA(33) = (0, 1.17) FIX
    diff_pipp_eye_3  <- fixed(0.98)
    label("IRT difficulty 3 increment for PIPP eye squeeze")
    # .mod $THETA(34) = (0, 0.98) FIX
    # .mod $THETA(35) = 100000 FIX sentinel placeholder, not exposed.

    # --- Item 27: PIPP nasolabial furrow (4 grades, indices 0..3) ---
    discr_pipp_naso  <- fixed(2.66)
    label("IRT discrimination for PIPP nasolabial furrow")
    # .mod $THETA(36) = (0, 2.66) FIX
    diff_pipp_naso_1 <- fixed(0.442)
    label("IRT difficulty 1 for PIPP nasolabial furrow")
    # .mod $THETA(37) = 0.442 FIX
    diff_pipp_naso_2 <- fixed(0.829)
    label("IRT difficulty 2 increment for PIPP nasolabial furrow")
    # .mod $THETA(38) = (0, 0.829) FIX
    diff_pipp_naso_3 <- fixed(1.27)
    label("IRT difficulty 3 increment for PIPP nasolabial furrow")
    # .mod $THETA(39) = (0, 1.27) FIX
    # .mod $THETA(40) = 100000 FIX sentinel placeholder, not exposed.

    # --- Item 28: NIPS total (8 grades, indices 0..7) ---
    discr_nips       <- fixed(4.37)
    label("IRT discrimination for NIPS total")
    # .mod $THETA(41) = (0, 4.37) FIX
    diff_nips_1      <- fixed(-0.082)
    label("IRT difficulty 1 for NIPS total")
    # .mod $THETA(42) = -0.082 FIX
    diff_nips_2      <- fixed(0.332)
    label("IRT difficulty 2 increment for NIPS total")
    # .mod $THETA(43) = (0, 0.332) FIX
    diff_nips_3      <- fixed(0.239)
    label("IRT difficulty 3 increment for NIPS total")
    # .mod $THETA(44) = (0, 0.239) FIX
    diff_nips_4      <- fixed(0.216)
    label("IRT difficulty 4 increment for NIPS total")
    # .mod $THETA(45) = (0, 0.216) FIX
    diff_nips_5      <- fixed(0.211)
    label("IRT difficulty 5 increment for NIPS total")
    # .mod $THETA(46) = (0, 0.211) FIX
    diff_nips_6      <- fixed(0.249)
    label("IRT difficulty 6 increment for NIPS total")
    # .mod $THETA(47) = (0, 0.249) FIX
    diff_nips_7      <- fixed(0.492)
    label("IRT difficulty 7 increment for NIPS total")
    # .mod $THETA(48) = (0, 0.492) FIX

    # --- Item 12: VAS continuous score (cm, 0-10) ---
    # Two thetas each for difficulty and discrimination depending on
    # OBSTYPE (1=video investigator, 2=bedside nurse).
    diff_vas_video    <- fixed(1.907720)
    label("VAS difficulty (latent value giving 50% VAS) for video investigator")
    # .mod $THETA(49) = 1.907720 FIX
    diff_vas_bedside  <- fixed(2.576450)
    label("VAS difficulty (latent value giving 50% VAS) for bedside nurse")
    # .mod $THETA(50) = 2.576450 FIX
    discr_vas_video   <- fixed(1.160130)
    label("VAS discrimination (steepness) for video investigator")
    # .mod $THETA(51) = (0, 1.160130) FIX
    discr_vas_bedside <- fixed(0.663043)
    label("VAS discrimination (steepness) for bedside nurse")
    # .mod $THETA(52) = (0, 0.663043) FIX

    # --------------------------------------------------------------
    # Estimated parameters of the PD layer (Valitalo 2017 Table 2).
    #
    # Baseline pain typical values by suctioning state (latent-variable
    # scale, no log transform -- the latent variable is N(0,1) by
    # convention so baselines can be negative).
    presuct  <- -0.27
    label("Typical baseline pain, before suctioning (latent units)")
    # Valitalo 2017 Table 2, h_bl,1 = -0.27 (SE 0.0615);
    # .mod $THETA(53) = -0.27.
    suct     <- 1.15
    label("Typical baseline pain, during suctioning (latent units)")
    # Valitalo 2017 Table 2, h_bl,2 = 1.15 (SE 0.0581);
    # .mod $THETA(54) = 1.15.
    aftsuct  <- -0.393
    label("Typical baseline pain, after suctioning (latent units)")
    # Valitalo 2017 Table 2, h_bl,3 = -0.393 (SE 0.0602);
    # .mod $THETA(55) = -0.393.

    # Morphine effect slope (linear), units (ng/mL)^-1.
    # NONMEM lower bound is 0 (`(0,0.0091)` would be the typical
    # encoding but the .mod actually has just `0.0091`, no bound).
    e_morph_pain <- 0.0091
    label("Morphine effect slope on latent pain ((ng/mL)^-1) -- subtracted from baseline")
    # Valitalo 2017 Table 2, h_morph = 0.0091 (SE 0.0016, 5-95% CI 0.00567-0.0116);
    # .mod $THETA(56) = 0.0091.

    # Study-time effect slope (linear), units 1/day.
    e_time_pain  <- 0.0476
    label("Study-time effect slope on latent pain (1/day) -- added to baseline")
    # Valitalo 2017 Table 2, h_time = 0.0476 (SE 0.0185, 5-95% CI 0.018-0.0782);
    # .mod $THETA(57) = (0, 0.0476).

    # --------------------------------------------------------------
    # IIV: 3x3 OMEGA BLOCK on (pre, intra, post) suctioning baselines.
    # Source $OMEGA BLOCK(3) values (lower-triangle, row-major):
    #   0.311
    #   0.172  0.34
    #   0.216  0.153  0.317
    # Diagonals reproduce Valitalo 2017 Table 2 BSV SDs:
    #   sqrt(0.311) =  0.558 -> x_bl,1 = 0.558 (Table 2 SE 0.0271);
    #   sqrt(0.34)  =  0.583 -> x_bl,2 = 0.583 (Table 2 SE 0.0237);
    #   sqrt(0.317) =  0.563 -> x_bl,3 = 0.563 (Table 2 SE 0.0231).
    # Correlations reported in Valitalo 2017 Results "Pharmacodynamic
    # model" paragraph: 54% (pre vs intra), 44% (intra vs post),
    # 69% (pre vs post).
    etapresuct + etasuct + etaaftsuct ~ c(0.311,
                                    0.172, 0.34,
                                    0.216, 0.153, 0.317)
    # .mod $OMEGA BLOCK(3); paper Table 2 SDs + correlations paragraph.

    # --------------------------------------------------------------
    # Residual error for the VAS observation (additive, two SDs by
    # observer type). NONMEM .mod ships SIGMAs as variances; nlmixr2
    # `add()` expects standard deviations.
    addSd_vas_video   <- sqrt(0.755)
    label("VAS additive residual SD (video investigator, 0-10 cm scale)")
    # .mod $SIGMA EPS(1) = 0.755 (variance); SD = sqrt(0.755) ~= 0.869.
    addSd_vas_bedside <- sqrt(2.55)
    label("VAS additive residual SD (bedside nurse, 0-10 cm scale)")
    # .mod $SIGMA EPS(2) = 2.55 (variance); SD = sqrt(2.55) ~= 1.597.
  })

  model({
    # Time in units$time = "day". Source NONMEM has TIME in minutes
    # and computes COVEFF = THETA(57) * TIME / (24 * 60); here `time`
    # is already in days so the conversion is implicit.
    days <- time

    # --- Baseline pain by individual and suctioning state ----------
    # Three correlated etas (3x3 OMEGA block above). Source NONMEM:
    #   MU_1 = THETA(53); presuctpain = MU_1 + ETA(1)
    #   MU_2 = THETA(54); suctpain   = MU_2 + ETA(2)
    #   MU_3 = THETA(55); aftsuctpain = MU_3 + ETA(3)
    presuctpain  <- presuct + etapresuct
    suctpain     <- suct    + etasuct
    aftsuctpain  <- aftsuct + etaaftsuct

    # Select the baseline by MOMENT covariate (1 = before, 2 = during,
    # 3 = after). Source NONMEM uses:
    #   PAIN = presuctpain; IF (moment.EQ.2) PAIN=suctpain;
    #   IF (moment.EQ.3) PAIN=aftsuctpain.
    pain_state_typ <-
      presuct  * (MOMENT == 1) +
      suct     * (MOMENT == 2) +
      aftsuct  * (MOMENT == 3)
    pain_state <-
      presuctpain * (MOMENT == 1) +
      suctpain    * (MOMENT == 2) +
      aftsuctpain * (MOMENT == 3)

    # --- Morphine effect and study-time effect on the latent --------
    # Source NONMEM: EFF = THETA(56)*CP; COVEFF = THETA(57)*TIME/(24*60).
    # Here time is days already, so no conversion is needed.
    eff <- e_morph_pain * CP_MORPH_NGML
    cov_eff <- e_time_pain * days

    # PAIN -= EFF (morphine reduces pain) ; PAIN += COVEFF (study time
    # increases pain). Source NONMEM: PAIN = PAIN - EFF + COVEFF.
    pain      <- pain_state     - eff + cov_eff
    pain_typ  <- pain_state_typ - eff + cov_eff

    # --------------------------------------------------------------
    # IRT graded-response probabilities for each item.
    # Source NONMEM dispatches via $ABBREVIATED INTEGER IND with IND
    # = 1, 6, 11, 16, 21, 26, 31, 36, 41 keyed off ITEM, then
    # `discr = THETA(IND); diff1..diff7 = THETA(IND+1..IND+7)`. Here
    # we expand the dispatch with an explicit `(ITEM == k) * ...`
    # multiplexer so all 8 categorical items contribute to derived
    # state outputs.
    #
    # The graded-response probability that the observed grade exceeds
    # level k is `pGTk = expit(discr * (pain - sum(diff_1..diff_k)))`.

    discr <-
      discr_cmt_alert * (ITEM == 1) +
      discr_cmt_calm  * (ITEM == 2) +
      discr_cmt_resp  * (ITEM == 3) +
      discr_cmt_body  * (ITEM == 5) +
      discr_cmt_face  * (ITEM == 7) +
      discr_pipp_brow * (ITEM == 25) +
      discr_pipp_eye  * (ITEM == 26) +
      discr_pipp_naso * (ITEM == 27) +
      discr_nips      * (ITEM == 28)

    diff1 <-
      diff_cmt_alert_1 * (ITEM == 1) +
      diff_cmt_calm_1  * (ITEM == 2) +
      diff_cmt_resp_1  * (ITEM == 3) +
      diff_cmt_body_1  * (ITEM == 5) +
      diff_cmt_face_1  * (ITEM == 7) +
      diff_pipp_brow_1 * (ITEM == 25) +
      diff_pipp_eye_1  * (ITEM == 26) +
      diff_pipp_naso_1 * (ITEM == 27) +
      diff_nips_1      * (ITEM == 28)

    diff2 <-
      diff_cmt_alert_2 * (ITEM == 1) +
      diff_cmt_calm_2  * (ITEM == 2) +
      diff_cmt_resp_2  * (ITEM == 3) +
      diff_cmt_body_2  * (ITEM == 5) +
      diff_cmt_face_2  * (ITEM == 7) +
      diff_pipp_brow_2 * (ITEM == 25) +
      diff_pipp_eye_2  * (ITEM == 26) +
      diff_pipp_naso_2 * (ITEM == 27) +
      diff_nips_2      * (ITEM == 28)

    diff3 <-
      diff_cmt_alert_3 * (ITEM == 1) +
      diff_cmt_calm_3  * (ITEM == 2) +
      diff_cmt_resp_3  * (ITEM == 3) +
      diff_cmt_body_3  * (ITEM == 5) +
      diff_cmt_face_3  * (ITEM == 7) +
      diff_pipp_brow_3 * (ITEM == 25) +
      diff_pipp_eye_3  * (ITEM == 26) +
      diff_pipp_naso_3 * (ITEM == 27) +
      diff_nips_3      * (ITEM == 28)

    diff4 <-
      diff_cmt_alert_4 * (ITEM == 1) +
      diff_cmt_calm_4  * (ITEM == 2) +
      diff_cmt_resp_4  * (ITEM == 3) +
      diff_cmt_body_4  * (ITEM == 5) +
      diff_cmt_face_4  * (ITEM == 7) +
      diff_nips_4      * (ITEM == 28)
    # Items 25/26/27 (PIPP) have only 3 active difficulty levels (4
    # grades 0..3); items 1/2/3/5/7 (COMFORT-B) have 4 active levels
    # (5 grades 1..5); item 28 (NIPS) has 7 active levels (8 grades
    # 0..7). diff5..diff7 contribute only for NIPS.

    diff5 <- diff_nips_5 * (ITEM == 28)
    diff6 <- diff_nips_6 * (ITEM == 28)
    diff7 <- diff_nips_7 * (ITEM == 28)

    diffGT1 <- diff1
    diffGT2 <- diff1 + diff2
    diffGT3 <- diff1 + diff2 + diff3
    diffGT4 <- diff1 + diff2 + diff3 + diff4
    diffGT5 <- diff1 + diff2 + diff3 + diff4 + diff5
    diffGT6 <- diff1 + diff2 + diff3 + diff4 + diff5 + diff6
    diffGT7 <- diff1 + diff2 + diff3 + diff4 + diff5 + diff6 + diff7

    pGT1 <- expit(discr * (pain - diffGT1))
    pGT2 <- expit(discr * (pain - diffGT2))
    pGT3 <- expit(discr * (pain - diffGT3))
    pGT4 <- expit(discr * (pain - diffGT4))
    pGT5 <- expit(discr * (pain - diffGT5))
    pGT6 <- expit(discr * (pain - diffGT6))
    pGT7 <- expit(discr * (pain - diffGT7))

    # Per-grade probabilities (source NONMEM `peq1..peq8` block).
    # For COMFORT-B / PIPP items, only peq1..peq5 are populated;
    # for NIPS, peq1..peq8 are populated. Negative probabilities
    # (numerical artifact) are floored to 0 to match the source.
    peq1 <- 1 - pGT1
    peq2 <- pGT1 - pGT2
    peq3 <- pGT2 - pGT3
    peq4 <- pGT3 - pGT4
    # peq5 / peq6 / peq7 / peq8 are item-dependent: NIPS uses 5..8 as
    # internal differences; non-NIPS items use peq5 = pGT4.
    peq5_nips <- pGT4 - pGT5
    peq5_other <- pGT4
    peq5 <- peq5_nips * (ITEM == 28) + peq5_other * (ITEM != 28)
    peq6 <- (pGT5 - pGT6) * (ITEM == 28)
    peq7 <- (pGT6 - pGT7) * (ITEM == 28)
    peq8 <- pGT7         * (ITEM == 28)

    # --------------------------------------------------------------
    # VAS predicted value (item 12, OBSTYPE-dependent).
    # Source NONMEM:
    #   VASDIFF   = THETA(49); if OBSTYPE==2: VASDIFF   = THETA(50).
    #   VASDISCR  = THETA(51); if OBSTYPE==2: VASDISCR  = THETA(52).
    #   VASPRED   = 10 * expit(VASDISCR * (PAIN - VASDIFF)).
    vas_diff_typ  <- diff_vas_video   * (OBSTYPE == 1) +
      diff_vas_bedside * (OBSTYPE == 2)
    vas_discr_typ <- discr_vas_video  * (OBSTYPE == 1) +
      discr_vas_bedside * (OBSTYPE == 2)

    vas_pred     <- 10 * expit(vas_discr_typ * (pain     - vas_diff_typ))
    vas_pred_typ <- 10 * expit(vas_discr_typ * (pain_typ - vas_diff_typ))

    # Observation: continuous VAS (cm, 0-10) with additive residual
    # error switching on OBSTYPE. Matches source NONMEM
    # `Y = VASPRED + EPS(1)*(2-OBSTYPE) + EPS(2)*(OBSTYPE-1)`.
    vas_add <- addSd_vas_video   * (OBSTYPE == 1) +
      addSd_vas_bedside * (OBSTYPE == 2)
    vas_pred ~ add(vas_add)
  })
}
