Skip to content

Commit

Permalink
v 0.2.0
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Oct 23, 2019
1 parent 0ba9c97 commit b10b5a7
Show file tree
Hide file tree
Showing 6 changed files with 152 additions and 70 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ Generate DataFrame with random multivariate data. Where:

## Structures

Struct information see [here](https://github.com/PharmCat/ReplicateBE.jl/blob/master/docs/src/struct.md).
Struct information see [here](https://pharmcat.github.io/ReplicateBE.jl/build/struct/).

## Acknowledgments

Expand Down
107 changes: 60 additions & 47 deletions docs/build/api/index.html

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion docs/build/search_index.js

Large diffs are not rendered by default.

35 changes: 33 additions & 2 deletions src/randrbeds.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,42 @@ mutable struct RandRBEDS
end

"""
```julia
randrbeds(;n=24, sequence=[1,1],
design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2],
intercept = 0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], formcoef = [0.0, 0.0],
dropsubj = 0.0, dropobs::Int = 0, seed::Int = 0)
```
Random dataset generation for bioequivalence.
# Parameters
* n: number of subjects
* sequence: distribution in sequences [1,1] means 1:1, [1,3] - 1:4 etc.
* design: desin matrix, each line is a sequence, each column - periods, cell - formulation id
* inter: inter-subject variance vector for G matrix (length 3)
* intra: intra-subject variance vector for R matrix (length 2)
* intercept: intercept effect value
* seqcoef: coefficients of sequences (length(sequence) == length(seqcoef) == size(design, 1))
* periodcoef: coefficients of periods (length(periodcoef) == size(design, 2))
* formcoef: coeficients of formulations
* dropobs: number of randomly dropped subjects
* seed: seed for random
Multivariate normal disribution:
```math
f(\\mathbf{x}; \\boldsymbol{\\mu}, \\boldsymbol{V}) = \\frac{1}{(2 \\pi)^{d/2} |\\boldsymbol{V}|^{1/2}}
\\exp \\left( - \\frac{1}{2} (\\mathbf{x} - \\boldsymbol{\\mu})^T V^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) \\right)
```
Where V:
```math
V_{i} = Z_{i}GZ_i'+R_{i}
```
Random dataset for bioequivalence.
"""
function randrbeds(;n=24, sequence=[1,1],
design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
Expand All @@ -32,13 +61,15 @@ function randrbeds(;n=24, sequence=[1,1],
end

"""
```julia
randrbeds(n::Int, sequence::Vector,
design::Matrix,
θinter::Vector, θintra::Vector,
intercept::Real, seqcoef::Vector, periodcoef::Vector, formcoef::Vector,
dropsubj::Float64, dropobs::Int, seed::Int)
```
Random dataset for bioequivalence.
Simple interface.
"""
function randrbeds(n::Int, sequence::Vector,
Expand Down
42 changes: 31 additions & 11 deletions src/rbe.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
# GENARAL REPLICATE BIOEQUIVALENCE STRUCTURE
#
"""
```julia
struct RBE
model::ModelFrame #Model frame
rmodel::ModelFrame #Random effect model
Expand All @@ -27,6 +28,7 @@
preoptim::Union{Optim.MultivariateOptimizationResults, Nothing} #Pre-optimization result object
optim::Optim.MultivariateOptimizationResults #Optimization result object
end
```
Replicate bioequivalence structure.
Expand Down Expand Up @@ -72,13 +74,16 @@ Mixed model fitting function for replicate bioequivalence without data preparati
Mixed model in matrix form:
<p align="center">
``y = X\\beta+Zu+\\epsilon``
</p>
```math
y = X\\beta+Zu+\\epsilon
```math
with covariance matrix for each subject:
``V_{i} = Z_{i}GZ_i'+R_{i}``
```math
V_{i} = Z_{i}GZ_i'+R_{i}
```
"""
function rbe(df; dvar::Symbol,
Expand Down Expand Up @@ -262,11 +267,13 @@ end #END OF rbe()
"""
This function apply following code for each factor before executing:
```julia
categorical!(df, subject);
categorical!(df, formulation);
categorical!(df, period);
categorical!(df, sequence);
sort!(df, [subject, formulation, period])
```
It can takes more time, but can help to avoid some errors like: "ERROR: type ContinuousTerm has no field contrasts".
"""
Expand Down Expand Up @@ -314,9 +321,12 @@ end
Returm -2logREML for rbe model
``logREML(\\theta,\\beta) = -\\frac{N-p}{2} - \\frac{1}{2}\\sum_{i=1}^nlog|V_{i}|-``
```math
logREML(\\theta,\\beta) = -\\frac{N-p}{2} - \\frac{1}{2}\\sum_{i=1}^nlog|V_{i}|-
``-\\frac{1}{2}log|\\sum_{i=1}^nX_i'V_i^{-1}X_i|-\\frac{1}{2}\\sum_{i=1}^n(y_i - X_{i}\\beta)'V_i^{-1}(y_i - X_{i}\\beta)``
-\\frac{1}{2}log|\\sum_{i=1}^nX_i'V_i^{-1}X_i|-\\frac{1}{2}\\sum_{i=1}^n(y_i - X_{i}\\beta)'V_i^{-1}(y_i - X_{i}\\beta)
```
"""
function reml2(rbe::RBE)
Expand All @@ -338,7 +348,9 @@ end
Return the standard errors for the coefficients of the model.
``se = \\sqrt{LCL'}``
```math
se = \\sqrt{LCL'}
```
"""
function StatsBase.stderror(rbe::RBE)
return collect(rbe.fixed.se)
Expand All @@ -348,7 +360,9 @@ end
Return model coefficients.
``\\beta = {(\\sum_{i=1}^n X_{i}'V_i^{-1}X_{i})}^{-1}(\\sum_{i=1}^n X_{i}'V_i^{-1}y_{i})``
```math
\\beta = {(\\sum_{i=1}^n X_{i}'V_i^{-1}X_{i})}^{-1}(\\sum_{i=1}^n X_{i}'V_i^{-1}y_{i})
```
"""
function StatsBase.coef(rbe::RBE)
return collect(rbe.fixed.est)
Expand All @@ -368,6 +382,8 @@ end
Compute confidence intervals for coefficients, with confidence level ```level``` (by default 95%).
# Arguments
```expci = true```: return exponented CI.
```inv = true```: return ```-estimate ± t*se```
Expand All @@ -376,6 +392,10 @@ Compute confidence intervals for coefficients, with confidence level ```level```
```df = :df3``` or ```df = :cont```: DF (contain) = N - rank(ZX).
```math
CI = stimate ± t(alpha, df)*se
```
"""
function StatsBase.confint(obj::RBE; level::Real=0.95, expci::Bool = false, inv::Bool = false, df = :sat)
confint(obj, 1 - level; expci = expci, inv = inv, df = df)
Expand Down Expand Up @@ -425,7 +445,7 @@ end
"""
theta(rbe::RBE)
Return theta vector (vector of variation parameters from optimization procedure).
Return theta (θ) vector (vector of variation parameters from optimization procedure).
"""
function theta(rbe::RBE)
return collect(rbe.θ)
Expand All @@ -443,7 +463,7 @@ end
Return design information object, where:
```julia
struct Design
obs::Int # Number of observations
subj::Int # Number of statistica independent subjects
Expand All @@ -457,7 +477,7 @@ Return design information object, where:
df3::Int # obs - rankxz (Contain DF for sequence and period)
df4::Int # obs - rankxz + p
end
```
"""
function design(rbe::RBE)::Design
return rbe.design
Expand Down
34 changes: 26 additions & 8 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,15 @@
Return contrast table for L matrix.
``F = \\frac{\\beta'L'(LCL')^{-1}L\\beta}{rank(LCL')}``
```math
F = \\frac{\\beta'L'(LCL')^{-1}L\\beta}{rank(LCL')}
```
DF for one-dimetion case:
``df = \\frac{2(LCL')^{2}}{g'Ag}``
```math
df = \\frac{2(LCL')^{2}}{g'Ag}
```
where ``A = 2H``
Expand Down Expand Up @@ -47,25 +51,39 @@ end
Return estimate table for L 1xp matrix.
``estimate = L\\beta``
```math
estimate = L\\beta
```
``se = \\sqrt{LCL'}``
```math
se = \\sqrt{LCL'}
```
``t = estimate/se``
```math
t = estimate/se
```
For ```df = :sat```:
``df = \\frac{2(LCL')^{2}}{g'Ag}``
```math
df = \\frac{2(LCL')^{2}}{g'Ag}
```
where ``A = 2H``
where ``g = \\triangledown _{\\theta}(LC^{-1}L')``
For ```df = :cont``` (contain):
``df = N - rank(ZX)``
```math
df = N - rank(ZX)
```
CI estimate is ``CI = stimate ± t(alpha, df)*se ``
CI estimate is:
```math
CI = stimate ± t(alpha, df)*se
```
"""
function estimate(rbe::RBE, L::Matrix; df = :sat, name = "Estimate", memopt = true, alpha = 0.05)
lcl = L*rbe.C*L'
Expand Down

2 comments on commit b10b5a7

@PharmCat
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

Release notes:

  • v0.2.0

    • bugfix
    • optimization
    • documentation
    • validation tests
    • memory caching
    • post optimization (can be done if Newton() failed)

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/4660

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if Julia TagBot is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.0 -m "<description of version>" b10b5a741f3df680a4de0e8e60fe162bc0ab9128
git push origin v0.2.0

Please sign in to comment.