using Pumas
Poisson Regression
This lesson is really comprehensive and has a wide audience in mind. We recommend the reader to briefly navigate the Table of Contents and choose which topics to dive in. If you are in doubt feel free to read it from top to bottom.
In this lesson we’ll dive into Poisson regression. We’ll start by defining what it is and, more importantly, when to use it.
1 What is Poisson Regression?
Poisson regression is a generalized linear model (GLM). GLMs uses a linear regression under the hood for inference but combines that with a non-linear transformation of the dependent variable (also called response).
This allows us to extend the interface of linear regression to more use cases that it would not be possible to apply it directly. This happens because of the nature of some dependent variables, being not linear and, subsequently, not continuous. Such instances occurs in binary data, discrete data, count data, or ordinal data.
Poisson Regression is a GLM that can handle count data. That is, when our response variable takes positive discrete values only, without an upper limit and with a natural lower limit of 0. It handles that by using the Poisson distrisbution to condition the values of our dependent variable.
1.1 Poisson Distribution
The Poisson distribution is named after mathematician Siméon Denis Poisson, whom had as advisor no less than Laplace. Poisson made important contributions not only to probability and statistics but also partial differential equations, fluid dynamics and more.
The Poisson distribution is used to model a random variable that takes positive discrete values. It accepts only a single parameter named \(\lambda\) which represents the constant mean rate of events in a single unit of time. In other words, \(\lambda\) represents the expected value of the distribution.
The parameter \(\lambda\) from Poisson distribution controls both the mean (sometimes referred as location parameter) and variance (sometimes referred as scale parameter). This makes the Poisson distribution not as flexible to model overdispersion in count data, since it has only one parameter to model both mean and variance. We can deal with overdispersion by introducing an extra parameter and, thus, using a different distribution. This distribution is named negative binomial and we have a separated tutorial for that. Don’t forget to also check it.
Take a look at some Poisson distribution plots for different values of \(\lambda\):
using Distributions
using CairoMakie
using AlgebraOfGraphics
= 0:0.1:10.0
x = [1, 3, 5]
λs = [Poisson(λ) for λ in λs]
dists = mapreduce(d -> pdf.(d, x), vcat, dists)
pdfs =
plt data((; λs = repeat(λs; inner = length(x)), x = repeat(x; outer = length(λs)), pdfs)) *
mapping(
:x => "outcome",
:pdfs => "PMF";
= :λs => nonnumeric => L"\mathrm{parameter}~\lambda",
color = :λs => nonnumeric,
col *
) visual(BarPlot)
draw(plt; axis = (; xticks = 0:10), legend = (; position = :bottom, titleposition = :left))
As you can see, the higher the value of \(\lambda\), higher the mean but also higher the variance of the outcome value.
1.2 Link Function \(e^x\)
The Poisson distribution has one caveat: \(\lambda\) has to be positive real number higher than 0; i.e. \(\lambda > 0\). This means that any parameter (or combinations of parameters) that we condition our Poisson distribution have to be higher than 0. We can accomplish this restriction with a non-linear transformation.
The exponential function is an ideal candidate to accomplish that. This is the exponential function:
\[\mathrm{exp}(x) = e^x\]
Let’s take a look at the graph of our exponential function:
data((; x = 0:0.1:10, y = exp.(0:0.1:10))) *
mapping(:x, :y => L"e^x") *
visual(Lines; color = :steelblue, linewidth = 3) |> draw
Let’s define our dependent variable as \(y\) which follows a Poisson Distribution. We can use the exponential function as a link function for the parameters of the model, such that \(y\) is parameterized by the exponential transformation of the parameter values:
\[y \sim \mathrm{Poisson}(e^\alpha)\]
That allows us to parametrize \(y\) with a real-valued parameters, while restraining \(y\) to being positive.
In this case, \(\alpha\) takes the role of an intercept, i.e. it is the basal rate of the values of \(y\).
1.3 Covariate \(x\)
Now with the dependent variable and intercept \(\alpha\) explained, let’s focus on covariate(s) \(x\).
Generally in regression problems, we attribute coefficients to act as weights to our covariate \(x\). These coefficients measure the strength of association between \(x\) and \(y\), where high positive coefficients means high positive association while high negative coefficients indicates high negative association.
Coefficients are usually denominated as the greek letter \(\beta\). Thus, we can parametrize \(y\) as:
\[y \sim \mathrm{Poisson}\left( e^\left(\alpha + \beta \cdot x \right) \right)\]
Almost everything that we want to infer or guess in statistics has greek letter such as \(\beta\) and \(\alpha\), while what is observed has roman letters such as \(x\) and \(y\).
We can express it in matrix form:
\[\mathbf{y} = \mathrm{Poisson}\left( e^\left(\boldsymbol{\alpha} + \mathbf{X} \cdot \boldsymbol{\beta} \right) \right)\]
where \(\mathbf{y}\), \(\boldsymbol{\alpha}\) and \(\boldsymbol{\beta}\) are vectors and \(\mathbf{X}\) is the data matrix where each row is an observation and each column a covariate.
Since \(\mathbf{y}\) and \(\mathbf{X}\) are known, all we have to do is to find the best value of our coefficient \(\beta\). That process is called model fit and we will deal with it in details further in this lesson.
1.3.1 How to Interpret Coefficient \(\beta\)?
Now, suppose we have found our ideal value for \(\beta\). How we would interpret our \(\beta\) estimated value?
First, to recap, \(\beta\) measures the strength of association between the covariate \(x\) and dependent variable \(y\). But, \(\beta\) is nested inside a non-linear transformation called exponential function:
\[\mathrm{linear~predictor} = \mathrm{exponential}(\alpha + \beta \cdot x) = e^{\left( \alpha + \beta \cdot x \right)}\]
where the linear predictor is the rate parameter (\(\lambda\)) of the Poisson distribution that we condition our dependent variable \(y\).
You can see that if \(\beta\) takes negative values, the association between covariate \(x\) with the linear predictor is negative, i.e. higher values of \(x\) means lower values of \(y\). Also, the converse is true. If \(\beta\) takes positive values, the association between covariate \(x\) with the linear predictor is positive, i.e. higher values of \(x\) means higher values of \(y\).
Thus, if you want just to analyze the raw directional effect of covariate \(x\) you can leave \(\beta\) in its original scale. However, if you want to have a deeper interpretation of the effect of covariate \(x\) you would need to transform \(\beta\).
First, we need to make \(\beta\) comparable to our dependent variable \(y\). Since \(y\) is modeled as the exponential of a linear combination of parameters and covariates, we need to also apply the exponential function to \(\beta\).
Second, notice that the exponential function has the following property:
\[e^{(a + b)} = e^a \cdot e^b\]
Hence, we can substitute \(a\) for our intercept \(\alpha\) and \(b\) for our covariate \(x\) multiplied by our coefficient \(\beta\):
\[e^{(\alpha + \beta \cdot x)} = e^\alpha \cdot e^{\beta \cdot x}\]
which means that the \(\beta\) has a compounding effect into the basal value (intercept \(\alpha\)), instead of a linear effect as we saw in other discrete models.
Now, in order to interpret \(\beta\) we need to first exponentiate all other parameters in the linear predictor term, such as the intercept \(\alpha\). \(\alpha\) is the basal value take \(y\) takes if all of the covariates \(x\) are set to zero. Then, for each increase of 1 unit of covariate \(x\) we expect to see a compounding effect of \(\beta\) into the basal value of \(y\), i.e. the intercept \(\alpha\).
Here’s an example. Suppose that \(\alpha = 1\) and \(\beta = 0.05\). This makes the basal value of \(y\):
\[e^\alpha = e^1 = e \approx 2.7182\]
The effect of the covariate \(x\) on \(y\) is:
\[e^{\beta \cdot x} = e^{0.05x} \approx 1.05x\]
Thus every increase of 1 unit in \(x\) increases \(y\) by 1.05, i.e. 5% increase.
2 When to use Poisson Regression?
Now that we know what Poisson regression is, let’s cover when to use it. The simple answer is whenever your dependent variable is discrete positive, i.e. have discrete positive outcome as values. Then you can use a Poisson likelihood inside your regression model.
3 seizurecount
Dataset
Let’s use the seizurecount
dataset from PharmaDatasets.jl
to run our Poisson regression example in this lesson:
using DataFrames
using DataFramesMeta
using PharmaDatasets
= dataset("seizurecount")
seizurecount first(seizurecount, 5)
Row | ID | Age | Gender | isF | TRT | NSeizure |
---|---|---|---|---|---|---|
Int64 | Int64 | String7 | Int64 | String7 | Int64 | |
1 | 1 | 43 | Female | 1 | Placebo | 6 |
2 | 2 | 55 | Male | 0 | Placebo | 2 |
3 | 3 | 76 | Female | 1 | Placebo | 5 |
4 | 4 | 93 | Male | 0 | Placebo | 0 |
5 | 5 | 95 | Female | 1 | Placebo | 2 |
We can see that we have 200 observations and also some key variables:
:NSeizure
(dependent variable): discrete positive outcome values are in the range (0, 16).:Age
: Age in years.:Gender
: patient gender with 2 levels – Female, Male.:isF
: discrete binary if patient is Female – 1, 0.:TRT
: treatment variable with 2 levels – Placebo, DrugA.
3.1 Exploratory Data Analysis
Let’s see some summary statistics with describe()
:
describe(seizurecount)
Row | variable | mean | min | median | max | nmissing | eltype |
---|---|---|---|---|---|---|---|
Symbol | Union… | Any | Union… | Any | Int64 | DataType | |
1 | ID | 75.5 | 1 | 75.5 | 150 | 0 | Int64 |
2 | Age | 48.4733 | 18 | 45.0 | 112 | 0 | Int64 |
3 | Gender | Female | Male | 0 | String7 | ||
4 | isF | 0.54 | 0 | 1.0 | 1 | 0 | Int64 |
5 | TRT | DrugA | Placebo | 0 | String7 | ||
6 | NSeizure | 3.39333 | 0 | 2.0 | 16 | 0 | Int64 |
As we can see :Age
has mean 48.47 and median 45.0 which signals for a bell-shaped (normally-distributed) distribution. We have weight values ranging from 18 and 112.
To conclude our summary statistics, :NSeizure
, our dependent variable, has mean 3.39 and median 2.0 which implies a long-tailed distribution (or assimetric shifted to the right).
It is also important to check the categorical columns to see what are the number of unique values (:nunique
) along with first and last values:
describe(seizurecount, :first, :last, :nunique; cols = [:Gender, :TRT])
Row | variable | first | last | nunique |
---|---|---|---|---|
Symbol | String7 | String7 | Int64 | |
1 | Gender | Female | Female | 2 |
2 | TRT | Placebo | DrugA | 2 |
Ok, no surprises in our categorical covariates.
@by seizurecount [:TRT, :isF] :count = length(:isF)
Row | TRT | isF | count |
---|---|---|---|
String7 | Int64 | Int64 | |
1 | Placebo | 0 | 24 |
2 | Placebo | 1 | 26 |
3 | DrugA | 0 | 45 |
4 | DrugA | 1 | 55 |
3.2 Visualizations
Now we turn to some visualizations to see what is going on beyond the summary statistics in our seizurecount
dataset.
First, we’ll explore our dependent variable: :NSeizure
.
data(seizurecount) *
mapping(:Gender, :NSeizure; color = :TRT, dodge = :TRT) *
expectation() |> draw AlgebraOfGraphics.
Female patients have more seizures than male patients. And male patients seems to not benefit from out treatment "DrugA"
.
data(seizurecount) *
mapping(:NSeizure; color = :TRT, layout = :Gender) *
density() |> draw AlgebraOfGraphics.
We can definitely see a difference between the densities of our dependent variable :NSeizures
if we stratify by :TRT
and :Gender
.
Now let’s turn our attention to our covariates :Age
and :TRT
.
We can see clearily a negative relationship between :Age
and our dependent variable :NSeizure
:
data(seizurecount) *
mapping(:Age, :NSeizure; color = :TRT) *
visual(Scatter; alpha = 0.7) + AlgebraOfGraphics.linear() * visual(; linewidth = 4)) |>
( draw
Also, this negative relationship does not appear to have any interaction effect with :TRT
since the slopes of the different :TRT
arms in the scatter plot above do not differ.
However, if we stratify by :Gender
, we can see that the interaction exists in males. The slopes have different values for different :TRT
arms in this specific stratum:
data(seizurecount) *
mapping(:Age, :NSeizure; color = :TRT, layout = :Gender) *
visual(Scatter; alpha = 0.7) + AlgebraOfGraphics.linear() * visual(; linewidth = 4)) |>
( draw
4 Pumas modeling
Finally, we’ll create a Pumas model for Poisson regression.
4.1 Pumas’ Workflow
Pumas’ workflow
- Define a model. It can be a
PumasModel
or aPumasEMModel
. - Define a
Subject
andPopulation
. - Fit your model.
- Perform inference on the model’s population-level parameters (also known as fixed effects).
- Predict from a fitted model using the empirical Bayes estimate or Simulate random observations from a fitted model using the sampling distribution (also known as likelihood).
- Model diagnostics and visual predictive checks.
Pumas modeling is a highly iterative process. Eventually you’ll probably go back and forth in the stages of workflow. This is natural and expected in Pumas’ modeling. Just remember to have a good code organization and version control for all of the workflow iterations.
Before we define the model, we need to convert the :TRT
column into a dummy variable. Currently :TRT
is a string:
unique(seizurecount.TRT)
2-element Vector{InlineStrings.String7}:
"Placebo"
"DrugA"
@rtransform! seizurecount :DrugA = ifelse(:TRT == "DrugA", 1, 0)
4.1.1 Defining a model in Pumas
To define a model in Pumas, you can use either the macros or the function interfaces. We will cover the macros interface in this lesson.
To define a model using the macro interface, you start with a begin .. end
block of code with @model
:
= @model begin end my_model
PumasModel
Parameters:
Random effects:
Covariates:
Dynamical system variables:
Dynamical system type: No dynamical model
Derived:
Observed:
This creates an empty model. Now we need to populate it with model components. These are additional macros that we can include inside the @model
definition. We won’t be covering all the possible macros that we can use in this lesson, but here is the full list:
@param
, fixed effects specifications.@random
, random effects specifications.@covariates
, covariate names.@pre
, pre-processing variables for the dynamic system and statistical specification.@dosecontrol
, specification of any dose control parameters present in the model.@vars
, shorthand notation.@init
, initial conditions for the dynamic system.@dynamics
, dynamics of the model.@derived
, statistical modeling of dependent variables.@observed
, model information to be stored in the model solution.
4.1.1.1 Model parameter with @param
Ok, let’s start with our poisson_model
using the macro @model
interface.
First, we’ll define our parameters with the @parameters
macro:
@param begin
∈ Domain(...)
ParameterName ...
end
To use the “in” (∈
) operator in Julia, you can either replace ∈
for in
or type \in <TAB>
for the Unicode character.
By using the begin ... end
block we can specify one parameter per line. We will just name the parameters the same name as it’s respective covariate, while also adding the β
prefix to make clear that this is a coefficient.
Regarding the Domain(...)
, Pumas has several types of domains for you to specify in your @param
block. Here is a list:
RealDomain
for scalar parametersVectorDomain
for vectorsPDiagDomain
for positive definite matrices with diagonal structurePSDDomain
for general positive semi-definite matrices
We will only use scalar parameters in our Poisson regression model, so we just need the RealDomain
.
If you don’t specify any arguments inside the domain constructor it will either error (for some domains that have required arguments) or will use the defaults. In the case of the RealDomain()
without arguments it just uses the following arguments:
RealDomain(; lower = -∞, upper = ∞, init = 0)
@param begin
∈ RealDomain()
α ∈ RealDomain()
βAge ∈ RealDomain()
βisF ∈ RealDomain()
βDrugA end
4.1.1.2 Model covariates with @covariates
We also have to specify our model covariates using the macro @covariates
. This can be done using a one-liner syntax, without begin ... end
block:
@covariates Age isF DrugA... # and so on
Or using the begin ... end
block and specifying one covariate per line. This is what we will do, since our model has several covariates, which makes easy to visualize them and also to comment out a covariate to perform iterative model exploration:
@covariates begin
Age
isF
DrugAend
4.1.1.3 Pre-processing variables with @pre
We can specify all the necessary variable and statistical pre-processing with the @pre
macro.
The @pre
block is traditionally used to specify the inputs of the Ordinary D ifferential Equations (ODE) system used for non-linear mixed-effects PK/PD models (NLME).
But we will use @pre
to specify our variable transformations instead.
Here, we are defining all the coefficients-covariates multiplication operation here. Note that we are adding the prefix _
to the resulting operation. Finally, we are also defining the linear_pred
as the linear predictors that should be added to the intercept and passed to the non-linear exponential function transformation:
@pre begin
= Age * βAge
_Age = isF * βisF
_isF = DrugA * βDrugA
_DrugA = _Age + _isF + _DrugA
linear_pred end
4.1.1.4 Statistical modeling of dependent variables with @derived
Our final block, @derived
, is used to specify all the assumed distributions of observed variables that are derived from the blocks above. This is where we include our dependent variable: NSeizures
.
Note that NSeizures
is being declared as following a Poisson
distribution. This is why we use the tilde notation ~
. It means (much like the mathematical model notation) that NSeizures
follows a Poisson
distribution. Since NSeizures
is a vector of values, we need to broadcast, i.e. vectorize, the operation with the dot .
operator:
~ Poisson.(exp.(λ)) NSeizures
where λ
is the single rate parameter that parametrizes the Poisson distribution. In some cases we would need to do a lot of broadcasting in the righthand side of the ~
. For example, since NSeizures
depends on the exponential function of all the linear predictors, i.e. the exponential function of α + linear_pred
, we would need to vectorize everything:
~ Poisson.(exp.(α .+ linear_pred)) NSeizures
This is cumbersome, so we can use the @.
macro which tells Julia to apply the .
in every operator and function call after it:
~ @. Poisson(exp(α + linear_pred)) NSeizures
We are using the @.
macro which tells Julia to vectorize (add the “dot syntax”) to all operations and function calls to the right of it.
Much better!
In this block we can use all variables defined in the previous blocks, in our example the @param
, @covariates
and @pre
blocks.
This is our @derived
block:
@derived begin
~ @. Poisson(exp(α + linear_pred))
NSeizure end
Here is our full model, once we combined all the macro blocks we just covered:
= @model begin
poisson_model @param begin
∈ RealDomain()
α ∈ RealDomain()
βAge ∈ RealDomain()
βisF ∈ RealDomain()
βDrugA end
@covariates begin
Age
isF
DrugAend
@pre begin
= Age * βAge
_Age = isF * βisF
_isF = DrugA * βDrugA
_DrugA = _Age + _isF + _DrugA
linear_pred end
@derived begin
~ @. Poisson(exp(α + linear_pred))
NSeizure end
end
PumasModel
Parameters: α, βAge, βisF, βDrugA
Random effects:
Covariates: Age, isF, DrugA
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: NSeizure
Observed: NSeizure
4.1.2 Define a Subject
and Population
Once we have our model defined we have to specify a Subject
and a Population
.
In Pumas, subjects (i.e. observations) are represented by the Subject
type and collections of subjects are represented as vectors of Subject
s are defined as Population
.
Subjects
can be constructed with the Subject()
constructor, for example:
Subject()
Subject
ID: 1
We just constructed a Subject
that has ID
equal to 1
and no extra information. In order to provide extra information, you would need to pass keyword arguments to the Subject()
constructor.
You can use the following keyword arguments inside Subject()
:
id
: identifier.observations
: a named tuple of the dependent variables.covariates
: a named tuple containing the covariates, ornothing
.events
: a vector ofEvent
s.time
: a vector of time stamps for the observations.
For this lesson, we will just use the covariates
keyword argument.
This is an example of a Subject
with ID
equal to 42
and with some covariates
specified as a named tuple:
Subject(id = 42, covariates = (; Age = 35, isF = 1, DrugA = 1))
Subject
ID: 42
Covariates: Age, isF, DrugA
Since a Population
is just a vector of Subjects
, we can use a simple map()
function over the ID
s present at our seizurecount
dataset:
= map(
pop -> Subject(
i = i,
id = (;
covariates = seizurecount.Age,
Age = seizurecount.isF,
isF = seizurecount.DrugA,
DrugA
),
),1:maximum(seizurecount.ID),
)
Population
Subjects: 150
Covariates: Age, isF, DrugA
Observations:
4.1.2.1 Reading Subject
s directly from a DataFrame
Another option is to use the read_pumas()
function to read Subject
s (and Population
) directly from a DataFrame
. This is more convenient than the map()
with the Subject()
constructor inside.
The read_pumas()
function accepts as first argument a DataFrame
followed by the following keyword options (which are very similar to the Subject()
constructor):
observations
: dependent variables specified by a vector of column names, i.e.[:NSeizure]
.covariates
: covariates specified by a vector of column names, i.e.[:Age, :isF, :DrugA]
.id
: specifies theID
column of theDataFrame
.time
: specifies thetime
column of theDataFrame
.event_data
: toggles assertions applicable to event data, eithertrue
or `false.
Our seizurecount
dataset do not have any event data, so we should set event_data = false
. Also, we don’t have a :time
column, but Pumas need us to specify a column for this keyword argument. So let’s create a dummy column :time
where all values are 0.0
:
@rtransform! seizurecount :time = 0.0
Now, we can create our population with read_pumas()
:
= read_pumas(
seizurecount_pop
seizurecount;= [:NSeizure],
observations = [:Age, :isF, :DrugA],
covariates = :ID,
id = :time,
time = false,
event_data )
Population
Subjects: 150
Covariates: Age, isF, DrugA
Observations: NSeizure
4.1.3 Fit your model
Now we are ready to fit our model! We already have a model specified, poisson_model
, along with a Population
: seizurecount_pop
. We can proceed with model fitting.
Model fiting in Pumas has the purpose of estimate parameters and is done by calling the fit()
function with the following positional arguments:
- Pumas model.
- a
Population
. - a named tuple of the initial parameter’s values.
- an inference algorithm.
4.1.3.1 Initial Parameter’s Values
We already covered model and Population
, now let’s talk about initial parameter’s values. It is the 3rd positional argument inside fit()
.
You can specify you initial parameters as a named tuple. For instance, if you want to have a certain parameter, β1
, as having an initial value as 3.14
, you can do so by passing it inside a named tuple in the 3rd positional argument of fit()
:
fit(model, population, (; β1 = 3.14))
You can also use the helper function init_params()
which will recover all the initial parameters we specified inside the model’s @param
block.
Since we used the RealDomain()
defaults, all of our example’s initial parameters are set to 0.0
.
4.1.3.2 Inference Algorithms
Finally, our last (4th) positional argument is the choice of inference algorithm.
Pumas has the following available inference algorithms:
Marginal Likelihood Estimation:
NaivePooled()
: first order approximation without random-effects.FO()
: first-order approximation.FOCE()
: first-order conditional estimation with automatic interaction detection.LaplaceI()
: second-order Laplace approximation.
Bayesian with Markov Chain Monte Carlo (MCMC):
BayesMCMC()
: MCMC using No-U-Turn Sampler (NUTS).
We can also use a Maximum A Posteriori (MAP) estimation procedure for any marginal likelihood estimation algorithm. You just need to call the MAP()
constructor with the desired marginal likelihood algorithm inside, for instance:
fit(model, population, init_params(model), MAP(FOCE()))
Ok, we are ready to fit our model. Let’s use the NaivePooled()
since we are not using random-effects in this lesson:
=
poisson_fit fit(poisson_model, seizurecount_pop, init_params(poisson_model), Pumas.NaivePooled())
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 7.084862e+02 1.520100e+04
* time: 0.03451895713806152
1 5.542122e+02 9.418562e+03
* time: 0.7631101608276367
2 5.389460e+02 4.901604e+03
* time: 0.7634921073913574
3 5.282896e+02 1.137961e+03
* time: 0.7638380527496338
4 5.263390e+02 4.599655e+02
* time: 0.7641961574554443
5 5.239327e+02 1.328781e+03
* time: 0.7645449638366699
6 5.112082e+02 4.376380e+03
* time: 0.7649240493774414
7 4.946020e+02 6.231816e+03
* time: 0.7652931213378906
8 4.738243e+02 5.619874e+03
* time: 0.7656619548797607
9 4.608481e+02 2.532336e+03
* time: 0.7660810947418213
10 4.576253e+02 3.761918e+02
* time: 0.7665300369262695
11 4.571941e+02 1.579012e+02
* time: 0.766977071762085
12 4.566451e+02 6.052517e+02
* time: 0.7674460411071777
13 4.551647e+02 1.328673e+03
* time: 0.7678649425506592
14 4.518011e+02 2.255583e+03
* time: 0.7682750225067139
15 4.447513e+02 3.202994e+03
* time: 0.7686269283294678
16 4.344023e+02 3.431150e+03
* time: 0.7689881324768066
17 4.262636e+02 2.411983e+03
* time: 0.7693359851837158
18 4.223312e+02 9.614849e+02
* time: 0.7696380615234375
19 4.213184e+02 2.425184e+01
* time: 0.7699191570281982
20 4.212106e+02 6.211133e+01
* time: 0.7701930999755859
21 4.210555e+02 1.521072e+02
* time: 0.770453929901123
22 4.206892e+02 2.831647e+02
* time: 0.7708189487457275
23 4.198621e+02 4.494193e+02
* time: 0.7712860107421875
24 4.184215e+02 5.611914e+02
* time: 0.7717220783233643
25 4.168899e+02 4.531058e+02
* time: 0.771996021270752
26 4.161840e+02 1.660627e+02
* time: 0.7722749710083008
27 4.160731e+02 9.336645e+00
* time: 0.7725400924682617
28 4.160686e+02 3.322691e+00
* time: 0.772813081741333
29 4.160685e+02 4.550751e-01
* time: 0.7730910778045654
30 4.160685e+02 1.900354e-02
* time: 0.7733631134033203
31 4.160685e+02 2.946794e-04
* time: 0.7736270427703857
FittedPumasModel
Successful minimization: true
Likelihood approximation: NaivePooled
Likelihood Optimizer: BFGS
Dynamical system type: No dynamical model
Log-likelihood value: -416.06855
Number of subjects: 150
Number of parameters: Fixed Optimized
0 4
Observation records: Active Missing
NSeizure: 150 0
Total: 150 0
---------------------
Estimate
---------------------
α 2.1503
βAge -0.018032
βisF 0.22675
βDrugA -0.35713
---------------------
4.1.4 Perform inference on the model’s population-level parameters
Once the model is fitted, we can analyze our inference and estimates. But first, let’s discuss what the is the meaning of the estimated coefficients in a Poisson regression model.
4.1.4.1 Interpreting the Coefficients
We can see that after the model is fitted, it prints a result with a summary and a table of the parameter estimates.
We can also recover the estimates as a named tuple with coef()
:
coef(poisson_fit)
(α = 2.150255882229772,
βAge = -0.018031843808920766,
βisF = 0.22675059376714593,
βDrugA = -0.35713447425647227,)
Or as a DataFrame
with coeftable()
:
coeftable(poisson_fit)
Row | parameter | estimate |
---|---|---|
String | Float64 | |
1 | α | 2.15026 |
2 | βAge | -0.0180318 |
3 | βisF | 0.226751 |
4 | βDrugA | -0.357134 |
Remember that the Poisson regression’s coefficients are the log scale, i.e. the natural logarithm of the effect size.
Log scales provide the facility that anything that is negative pulls down the base values of our dependent value \(y\). Anything positive pulls it up; and zero is a neutral effect. But since it is in a log scale it is not easily interpretable. This is why, in order for interpret it, often we undo the log transformation by exponentiating the effects, i.e. exponentiating the coefficients.
We can do that easily with the DataFrame
returned by coeftable()
:
@rtransform coeftable(poisson_fit) :estimate_exp = exp(:estimate)
Row | parameter | estimate | estimate_exp |
---|---|---|---|
String | Float64 | Float64 | |
1 | α | 2.15026 | 8.58706 |
2 | βAge | -0.0180318 | 0.98213 |
3 | βisF | 0.226751 | 1.25452 |
4 | βDrugA | -0.357134 | 0.699678 |
Now let us interpret the :estimate_exp
column.
Our basal value α
is the estimated number of seizures when all the other covariates are zero.
Let us now interpret the βAge
, the covariate Age
’s coefficient. It’s value is approximate 0.98, which means that for every increase in a single unit of Age
we can expect a decrease in 2% in the number of seizures. Numerically this is: \(\alpha \cdot 0.98 \cdot \mathrm{Age}\). For instance, a 30 increase in Age
would reduce the basal value \(8.58\) to \((8.58 - ( 8.58 \cdot 0.98)) \cdot 30 = 5.15\). Note that this association has the assumption that all other covariates are fixed/controlled.
The βisF
, which represents the effect for the dummy variable isF
. The exponentiated coefficient value is 1.25, which means that for every unit of increase in the covariate isF
we can expected on average an increase of 25% in the number of seizures. Since, isF
is a dummy variable taking either 0 or 1 as values, this means that female observations have, on average, 25% more seizures events than male observations.
Our DrugA
covariate, which is a binary representation of the observation being either on the placebo or treatment for drug A arm, has an exponentiated coefficient of 0.70. This means that observations treated with drug A, when all other covariates are controlled, have 30% less seizures than when treated with placebo.
We can also inspect our estimated coefficients’ standard errors (SE
) along with the 95% confidence intervals with the infer()
function:
infer(poisson_fit)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: NaivePooled
Likelihood Optimizer: BFGS
Dynamical system type: No dynamical model
Log-likelihood value: -416.06855
Number of subjects: 150
Number of parameters: Fixed Optimized
0 4
Observation records: Active Missing
NSeizure: 150 0
Total: 150 0
------------------------------------------------------------------------
Estimate SE 95.0% C.I.
------------------------------------------------------------------------
α 2.1503 0.26682 [ 1.6273 ; 2.6732 ]
βAge -0.018032 0.0049478 [-0.027729; -0.0083343]
βisF 0.22675 0.16322 [-0.093163; 0.54666 ]
βDrugA -0.35713 0.16018 [-0.67108 ; -0.043187 ]
------------------------------------------------------------------------
Also if you prefer other confidence interval band, you can choose with the keyword argument level
inside infer()
.
For instance, one common band for the confidence intervals are 89%:
infer(poisson_fit; level = 0.89)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: NaivePooled
Likelihood Optimizer: BFGS
Dynamical system type: No dynamical model
Log-likelihood value: -416.06855
Number of subjects: 150
Number of parameters: Fixed Optimized
0 4
Observation records: Active Missing
NSeizure: 150 0
Total: 150 0
----------------------------------------------------------------------
Estimate SE 89.0% C.I.
----------------------------------------------------------------------
α 2.1503 0.26682 [ 1.7238 ; 2.5767 ]
βAge -0.018032 0.0049478 [-0.025939; -0.010124]
βisF 0.22675 0.16322 [-0.034113; 0.48761 ]
βDrugA -0.35713 0.16018 [-0.61313 ; -0.10114 ]
----------------------------------------------------------------------
4.1.5 Predict from a fitted model or Simulate random observations from a fitted model
Now that we understand our model estimates, we can use it to make predictions or simulate random observations.
4.1.5.1 Predictions with predict()
In order to calculate prediction from a fitted model, we need to use the predict()
function and passing two positional arguments:
- a fitted Pumas model.
- either a single
Subject
or a collection of them with aPopulation
.
For example to make predictions onto a single Subject
, we can:
= Subject(id = 42, covariates = (; Age = 35, isF = 1, DrugA = 1)) my_subject
Subject
ID: 42
Covariates: Age, isF, DrugA
predict()
also needs a keyword argument obstimes
which represents the observation times that we want to predict the dependent variable. The value must be a vector. Since we used a dummy 0.0
value as :time
variable in our read_pumas
function, we can pass the same 0.0
but as a vector to the predict()
function.
predict(poisson_fit, my_subject; obstimes = [0.0])
SubjectPrediction
ID: 42
Predictions: NSeizure: (n=1)
Covariates: Age, isF, DrugA
Time: [0.0]
We can also generate prediction for a Population
:
predict(poisson_fit, seizurecount_pop; obstimes = [0.0])
Prediction
Subjects: 150
Predictions: NSeizure
Covariates: Age, isF, DrugA
4.1.5.2 Simulations with simobs()
We can generate random observations from our fitted Pumas model with the simobs()
function. You need to supply at least three positional arguments:
- a non-fitted Pumas model, in our case
poisson_model
. - either a single
Subject
or a collection of them with aPopulation
. - a named tuple of parameter values.
The 1st and 2nd positional arguments we already covered in our previous examples. Let’s focus on the 3rd positional argument: param
. It can be a named tuple with the parameter values. For instance, we can use our estimated parameters from the fitted model with coef(poisson_fit)
:
coef(poisson_fit)
(α = 2.150255882229772,
βAge = -0.018031843808920766,
βisF = 0.22675059376714593,
βDrugA = -0.35713447425647227,)
simobs(poisson_model, seizurecount_pop, coef(poisson_fit));
You can also pass a single parameter as a Pumas fitted model to simobs()
.
= simobs(poisson_fit); seizurecont_sim
4.1.6 Model diagnostics
Finally, our last step is to assess model diagnostics.
4.1.6.1 Assessing model diagnostics
To assess model diagnostics we can use the inspect()
function in our fitted Pumas model:
= DataFrame(inspect(poisson_fit))
my_inspect first(my_inspect, 5)
[ Info: Calculating predictions.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.
Row | id | time | evid | NSeizure | Age | isF | DrugA | NSeizure_pred | NSeizure_ipred | _Age | _isF | _DrugA | linear_pred |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
String | Float64 | Int64 | Float64? | Int64? | Int64? | Int64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | |
1 | 1 | 0.0 | 0 | 6.0 | 43 | 1 | 0 | 4.96115 | 4.96115 | -0.775369 | 0.226751 | -0.0 | -0.548619 |
2 | 10 | 0.0 | 0 | 5.0 | 50 | 1 | 0 | 4.37285 | 4.37285 | -0.901592 | 0.226751 | -0.0 | -0.674842 |
3 | 100 | 0.0 | 0 | 6.0 | 42 | 1 | 1 | 3.53437 | 3.53437 | -0.757337 | 0.226751 | -0.357134 | -0.887721 |
4 | 101 | 0.0 | 0 | 2.0 | 43 | 0 | 1 | 2.76697 | 2.76697 | -0.775369 | 0.0 | -0.357134 | -1.1325 |
5 | 102 | 0.0 | 0 | 0.0 | 55 | 0 | 1 | 2.22859 | 2.22859 | -0.991751 | 0.0 | -0.357134 | -1.34889 |
4.1.6.1.1 AIC
We can also check for Akaike Information Criterion (AIC):
aic(poisson_fit)
840.1370972852749
4.1.6.1.2 BIC
Besides AIC, there is also the Bayesian Information Criterion (BIC):
bic(poisson_fit)
852.1796384616599
4.1.6.2 Visual predictive checks
To conclude, we can inspect visual predictive checks with the vpc_plot()
function. But first, we need to generate a VPC
object with the vpc()
function.
By default, vpc
will use :time
as the independent variable for VPC. We need to change that, because or :time
variable is just a vector of 0
s in the Poisson regression.
Let’s use :Age
:
= vpc(poisson_fit; covariates = [:Age], bandwidth = 85); vpc_poisson
Now, we need to use the vpc_plot
function into our newly created VPC
object. For this, we need to import the package PumasUtilities
:
using PumasUtilities
Since our interval is not on \([0,1]\), we need to tell vpc_plot
to not use it with unit_yaxis=false
. Additionally, we change our ylabel
to represent a count instead of a proportion.
vpc_plot(vpc_poisson; unit_yaxis = false, axis = (; ylabel = "Count NSeizure"))
5 Conclusion
I hope you enjoyed this tutorial on Poisson regression. Whenever you have count data, this should be your first approach. We’ll explore a more robust and complex model named negative binomial regression next. If you are experiencing problems fitting your count data with Poisson regression you would also want to try it.