Discrete Response Models

Vijay Ivaturi, Chris Rackauckas
July 19th, 2019

Introduction

In this tutorial we will go over the simulation of discrete responses. Many pharmacometrics scenarios have observables, such as pain scores or counts, which necessarily have to be discrete. Handling this discreteness can be paramount to getting an appropriate data fit and to properly understand the variation.

Luckily, in Pumas, discrete outputs are handled no differently from the rest of the Pumas toolchain. In Pumas, to have a discrete distribution as output, simply have that your derived or observed variables come from a discrete distribution like a Poisson process.

Binary Response Example

First, let's take a look at a binary response. A binary response is a model which gives an output of 0 or 1 with a probability p. In Pumas, this is represented by a Bernoulli(p) distribution.

To get started, first let's load Pumas and read in some example data:

using Pumas, StatsFuns
data = read_pumas(example_nmtran_data("pain_remed"),
    cvs = [:arm, :dose, :conc, :painord,:remed],event_data=false)
Population
  Subjects: 160
  Covariates: arm, dose, conc, painord, remed
  Observables: dv

Next let's implement a model with a binary response. Here, we do not have an @dynamics porition. Pumas will automatically handle this. In our derived, we define a logistic from StatsFuns

import StatsFuns.logistic
binary_model = @model begin
    @param begin
        intercept  RealDomain(init=0.001)
        tvslope  RealDomain(init=0.0001)
        Ω  VectorDomain(1)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates arm dose

    @pre begin
        rx = dose > 0 ? 1 : 0
        slope = tvslope*rx
        logit = intercept + slope + η[1]
    end

    @derived begin
        pain = logistic(logit)
        dv ~ Bernoulli(logistic(logit))
    end
end
PumasModel
  Parameters: intercept, tvslope, Ω
  Random effects: η
  Covariates: arm, dose
  Dynamical variables: 
  Derived: pain, dv
  Observed: pain, dv

Note that we could have alternatively defined our @pre like:

@pre begin
    logit = intercept + tvslope*(dose > 0 ? 1 : 0) + η[1]
    logit = intercept + tvslope*Int(dose > 0) + η[1]
end

more directly instead of using logistic.

Now let's fit our model to the data:

param = (
    intercept = 0.001,
    tvslope = 0.0001,
    Ω = [1.0]
    )
res = fit(binary_model,data,param,Pumas.LaplaceI())
FittedPumasModel

Successful minimization:                true

Likelihood approximation:     Pumas.LaplaceI
Objective function value:            1083.01
Total number of observation records:    1920
Number of active observation records:   1920
Number of subjects:                      160

-------------------
          Estimate
-------------------
intercept -1.3085
tvslope    1.7386
Ω₁         1.5374
-------------------

and simulate some outputs:

sim = simobs(binary_model,data,res.param)
simdf = DataFrame(sim, include_events=false)
first(simdf,6) # Print only the first 6 rows

6 rows × 9 columns

idtimepaindvarmdoseconcpainordremed
StringFloat64Float64BoolInt64Int64Float64Int64Int64
110.00.792461true2200.030
210.50.792461true2201.1557810
311.00.792461true2201.3721100
411.50.792461true2201.3005800
512.00.792461true2201.1919510
612.50.792461true2201.1360210

Note that now our simulation output for dv is true/false values pulled with probability given by logit dependent on the individual's random effects.

Poisson Response Example

Next let's use a Poisson counting process in our model. Here we generate a population where everyone is receiving the same doses as a covariate.

pop = Population(map(i -> Subject(id=i,cvs=(dose = i.*10,),time=[0.0]),1:10))
Population
  Subjects: 10
  Covariates: dose

Now we define our model without dynamics, and directly use the dose information to predict the count for some observable dv:

poisson_model = @model begin
  @param begin
    tvbase  RealDomain(init=3.0, lower=0.1)
    d50  RealDomain(init=50, lower=0.1)
    Ω   PSDDomain(fill(0.1, 1, 1))
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    baseline = tvbase*exp(η[1])
  end

  @covariates dose

  @derived begin
    dv ~ @. Poisson(baseline*(1-dose/(dose + d50)))
  end
end
PumasModel
  Parameters: tvbase, d50, Ω
  Random effects: η
  Covariates: dose
  Dynamical variables: 
  Derived: dv
  Observed: dv

and simulate runs from the model:

sim = simobs(poisson_model,pop)
simdf = DataFrame(sim, include_events=false)

10 rows × 4 columns

idtimedvdose
StringFloat64Int64Int64
110.0110
220.0220
330.0130
440.0140
550.0050
660.0160
770.0070
880.0380
990.0190
10100.00100

Here dv is an integer output probabilistically dependent on the dose.

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

Appendix

These tutorials are part of the PumasTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using PumasTutorials
PumasTutorials.weave_file("discrete","discrete_response_models.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = "C:\Users\accou\AppData\Local\atom\app-1.38.2\atom.exe"  -a
  JULIA_NUM_THREADS = 4

Package Information:

Status `C:\Users\accou\.julia\environments\v1.1\Project.toml`
[621f4979-c628-5d54-868e-fcf4e3e8185c] AbstractFFTs 0.4.1
[c52e3926-4ff0-5f6e-af25-54175e0327b1] Atom 0.8.8
[f0abef60-9ec0-11e9-27de-db6506a91768] AutoOffload 0.1.0
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[4ece37e6-a012-11e8-38cd-91247efc2c34] Bioequivalence 0.1.0
[336ed68f-0bac-5ca0-87d4-7b16caf5d00b] CSV 0.5.9
[c5f51814-7f29-56b8-a69c-e4d8f6be1fde] CUDAdrv 3.0.1
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.1
[49dc2e85-a5d0-5ad3-a950-438e2897f1b9] Calculus 0.5.0
[7057c7e9-c182-5462-911a-8362d720325c] Cassette 0.2.5
[34da2185-b29b-5c13-b0c7-acf172513d20] Compat 2.1.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.1.0
[667455a9-e2ce-5579-9412-b964f529a492] Cubature 1.4.0
[a93c6f00-e57d-5684-b7b6-d8193f3e46c0] DataFrames 0.18.4
[82cc6244-b520-54b8-b5a6-8a565e85f1d0] DataInterpolations 0.2.0
[31a5f54b-26ea-5ae9-a837-f05ce5417438] Debugger 0.5.0
[bcd4f6db-9728-5f36-b5f7-82caef46ccdb] DelayDiffEq 5.9.1
[2b5f629d-d688-5b77-993f-72d75c75574e] DiffEqBase 5.16.3
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.2.0
[31c91b34-3c75-11e9-0341-95557aab0344] DiffEqBenchmarks 0.1.0
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2+
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.13.0
[01453d9d-ee7c-5054-8395-0335cb756afa] DiffEqDiffTools 0.14.0
[aae7a2af-3d4f-5e19-a356-7da93b79d9d0] DiffEqFlux 0.5.2
[071ae1c0-96b5-11e9-1965-c90190d839ea] DiffEqGPU 0.1.0
[c894b116-72e5-5b58-be3c-e6d8d4ac2b12] DiffEqJump 6.1.1+
[8f2b45d5-b17b-5532-9e92-98ae0077e2e3] DiffEqMachineLearning 0.1.0
[78ddff82-25fc-5f2b-89aa-309469cbf16f] DiffEqMonteCarlo 0.15.1
[77a26b50-5914-5dd7-bc55-306e6241c503] DiffEqNoiseProcess 3.3.1
[9fdde737-9c7f-55bf-ade8-46b3f136cc48] DiffEqOperators 3.5.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.2.0
[a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9] DiffEqProblemLibrary 4.3.0
[41bf760c-e81c-5289-8e54-58b1f1f8abe2] DiffEqSensitivity 3.3.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.6.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[e30172f5-a6a5-5a46-863b-614d45cd2de4] Documenter 0.23.0
[587475ba-b771-5e3f-ad9e-33799f191a9c] Flux 0.8.3
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3+
[ba82f77b-6841-5d2e-bd9f-4daf811aec27] GPUifyLoops 0.2.5
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.1.0
[bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26] GitHub 5.1.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[42fd0dbc-a981-5370-80f2-aaf504508153] IterativeSolvers 0.8.1
[033835bb-8acc-5ee8-8aae-3f567f8a3819] JLD2 0.1.2
[e5e0dc1b-0480-54bc-9374-aad01c23163d] Juno 0.7.0
[2d691ee1-e668-5016-a719-b2531b85e0f5] LIBLINEAR 0.5.1
[7f56f5a3-f504-529b-bc02-0b1fe5e64312] LSODA 0.4.0
[6f1fad26-d15e-5dc8-ae53-837a1d7b8c9f] Libtask 0.3.0
[c7f686f2-ff18-58e9-bc7b-31028e88f75d] MCMCChains 0.3.10
[33e6dc65-8f57-5167-99aa-e5a354878fb2] MKL 0.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.5.0
[4886b29c-78c9-11e9-0a6e-41e1f4161f7b] MonteCarloIntegration 0.0.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[8faf48c0-8b73-11e9-0e63-2155955bfa4d] NeuralNetDiffEq 0.1.0
[09606e27-ecf5-54fc-bb29-004bd9f985bf] ODEInterfaceDiffEq 3.3.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.12.0
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.2.0
[14b8a8f1-9102-5b29-a752-f990bacb7fe1] PkgTemplates 0.6.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.3
[92933f4c-e287-5a05-a399-4b506db050ca] ProgressMeter 1.0.0
[d7b8c89e-ad89-52e0-b9fd-d0ed321fa021] Pumas 0.1.0
[b7b41870-aa11-11e9-048a-09266ec4a62f] PumasTutorials 0.0.1
[438e738f-606a-5dbb-bf0a-cddfbfd45ab0] PyCall 1.91.2
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[1fd47b50-473d-5c70-9696-f719f8f3bcdc] QuadGK 2.1.0
[612083be-0b0f-5412-89c1-4e7c75506a58] Queryverse 0.3.1
[6f49c342-dc21-5d91-9882-a32aef131414] RCall 0.13.3
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[37e2e3b7-166d-5795-8a7a-e32c996b4267] ReverseDiff 0.3.1
[295af30f-e4ad-537b-8983-00126c2a3abe] Revise 2.1.6
[2b6d1eac-7baa-5078-8adc-e6a3e659f14f] SingleFloats 0.1.3
[47a9eef4-7e08-11e9-0b38-333d64bd3804] SparseDiffTools 0.4.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[4c63d2b9-4356-54db-8cca-17b64c39e42c] StatsFuns 0.8.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[9672c7b4-1e72-59bd-8a11-6ac3964bc41f] SteadyStateDiffEq 1.5.0
[789caeaf-c7a9-5a7d-9973-96adeb23e2a0] StochasticDiffEq 6.6.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[fd094767-a336-5f1f-9728-57cf17d0bbfb] Suppressor 0.1.1
[6fc51010-71bc-11e9-0e15-a3fcc6593c49] Surrogates 0.1.0
[9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c] Tracker 0.2.2
[fce5fe82-541a-59a6-adf8-730c64b5f9a0] Turing 0.6.18
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.16.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.1
[e88e6eb3-aa80-5325-afca-941959d7151f] Zygote 0.3.2