Hello,
Is there something similar to the SAME option in NONMEM to estimate one BOV, instead of a BOV per occasion? If so, could I get an example of how it’s done in Pumas? Thank you.
In Pumas, we can flexibly specify the random effect in the random block, as you probably know. We don’t use special ordering or syntax such as SAME just a name and the distribution:
@random begin
random_effect ~ Dist(Input...)
end
The above defines a random effect called random_effect that is assumed to have the distribution Dist parameterized by the inputs Input....
In your case, I assume you want a random effect per occasion that is normally distributed. Each occasion should have its own random effect, but they should all follow the same exact distribution: mean 0 and variance ω². So let us say that you have three occasions. This case can be programmed in different ways, but let me show you two. Once is using the scalar (univariate, a number) distribution Normal and the other uses the vector (multivariate) distribution MvNormal that you’re used to.
The scalar version would look like this
@param begin
ω ∈ RealDomain(lower = 0.00001)
end
@random begin
κ₁ ~ Normal(0, ω)
κ₂ ~ Normal(0, ω)
κ₃ ~ Normal(0, ω)
end
and then you can use the kappas where appropriate in your other blocks. Notice, that the scalar distribution Normal uses mean and standard deviations as parameters.
The slightly more concise version, especially so if you have many occasions, is to use the multivariate distribution MvNormal. As you probably realized, we often use the PSDDomain and PDiagDomain here. This is of course not appropriate because the elements of the matrix parameters in these domains are allowed to vary freely (as long as they are positive semi-definite or positive definite and diagonal). However, the PSDDomain and PDiagDomainparameter domains are only meant to be an easy way to maintain proper definiteness of the variance-covariance matrices, we do not in any way require that the variance-covariance matrix in the MvNormal random effect specifications are defined using parameters from these domains. For example, a PDiagDomain can be written like this
@param begin
Ω ∈ PDiagDomain(2)
end
@random begin
η ~ MvNormal(Ω)
end
but you could also have written
@param begin
ω₁² ∈ RealDomain(lower=1e-5)
ω₂² ∈ RealDomain(lower=1e-5)
end
@random begin
η ~ MvNormal(Diagonal([ω₁², ω₂²]))
end
Remember, a Diagonal is a special matrix type in Julia whose constructor accepts a vector of diagonal elements and returns an object that acts as a matrix with only diagonal elements, but saves on storage by not explicitly storing the off-diagonal elements.
julia> Diagonal([0.0, 3.0, 5.0])
3×3 Diagonal{Float64, Vector{Float64}}:
0.0 ⋅ ⋅
⋅ 3.0 ⋅
⋅ ⋅ 5.0
So, to complete a long answer to a short question, the multivariate distribution way of specifying BOV random effects is
@param begin
ω² ∈ RealDomain(lower=1e-5)
end
@random begin
κ ~ MvNormal(Diagonal([ω², ω², ω²]))
end
or, using the fill(value, number_of_elements) function
@param begin
ω² ∈ RealDomain(lower=1e-5)
end
@random begin
κ ~ MvNormal(Diagonal(fill(ω², 3))
end
Since this creates a vector random variable/effect κ where each element is univariately normal, but they share the variance. The vector version is especially useful if you code your occasion covariate OCC as values 1, 2, and 3. Then you can do the following in your pre-block to pick the appropriate occasion effect
@covariates OCC
@pre begin
CL = tvcl * exp(κ[OCC])
end
By the way, since occasions often change when an observation occurs, it’s often useful/correct to use read_pumas(df; covariates_direction=:left) instead of the default read_pumas(df; covariates_direction=:right) to make sure that the occasion effect is activated at the right time.
I hope this was helpful. Otherwise, please let me know.
2 Likes