Merge pull request #48 from PharmCat/dev

Dev

.github/workflows/Documenter.yml

@@ -23,12 +23,11 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 45
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: "1.8"
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-docdeploy@v1
+      - uses: julia-actions/cache@v1
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

.github/workflows/Tier1.yml

@@ -39,15 +39,26 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
+      - uses: actions/cache@v4
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
+        if: ${{ matrix.os == 'ubuntu-latest' && matrix.version == '1' && matrix.arch == 'x64' }}
         with:
-          files: lcov.info
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}

Project.toml

@@ -3,7 +3,7 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2"
 keywords = ["lenearmodel", "mixedmodel"]
 desc = "Mixed-effects models with flexible covariance structure."
 authors = ["Vladimir Arnautov <[email protected]>"]
-version = "0.15.1"
+version = "0.16.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

docs/src/api.md

@@ -124,6 +124,12 @@ Metida.ToeplitzParameterized
 Metida.Unstructured
 ```
 
+### Metida.ScaledWeightedCov
+
+```@docs
+Metida.ScaledWeightedCov
+```
+
 ### Methods
 
 ### Metida.caic

docs/src/custom.md

@@ -26,10 +26,10 @@ function Metida.gmat!(mx, θ, ::YourCovarianceStruct)
 end
 ```
 
-Function `rmat!` have 4 arguments and add repeated effect to V': V = V' + R (so V = Z * G * Z' + R), `mx` - V' matrix, `θ` - theta vector for this effect, `rz` - subject effect matrix, `ct` - your covariance type object. For example, `rmat!` for Heterogeneous Toeplitz Parameterized structure is specified bellow (`TOEPHP_  <: AbstractCovarianceType`).
+Function `rmat!` have 5 arguments and add repeated effect to V': V = V' + R (so V = Z * G * Z' + R), `mx` - V' matrix, `θ` - theta vector for this effect, `rz` - subject effect matrix, `ct` - your covariance type object, `sb` = block number. For example, `rmat!` for Heterogeneous Toeplitz Parameterized structure is specified bellow (`TOEPHP_  <: AbstractCovarianceType`).
 
 ```
-function Metida.rmat!(mx, θ, rz, ct::TOEPHP_)
+function Metida.rmat!(mx, θ, rz, ct::TOEPHP_, ::Int)
     l     = size(rz, 2)
     vec   = rz * (θ[1:l])
     s   = size(mx, 1)
@@ -123,7 +123,7 @@ Metida.fit!(lmm)
 
 # for R matrix
 
-function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure)
+function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure, ::Int)
     vec = Metida.tmul_unsafe(rz, θ)
     rn    = size(mx, 1)
     if rn > 1

docs/src/details.md

@@ -13,6 +13,33 @@ In matrix notation a mixed effect model can be represented as:
 y = X\beta + Zu + \epsilon
 ```
 
+where:
+
+
+```math
+\epsilon \sim N(0, R)
+
+\\
+
+u \sim N(0, G)
+
+\\
+
+y \sim N(X\beta, V)
+
+```
+
+where V depends on covariance sructure and parameters ``\theta``:
+
+```math
+V = CovStruct(\theta)
+
+```
+
+The unknown parameters include the regression parameters in ``\beta`` and covariance parameters in ``\theta``.
+
+Estimation of these model parameters relies on the use of a Newton-Ralphson (by default) algorithm. When we use either algorithm for finding REML solutions, we need to compute ``V^{-1}`` and its derivatives with respect to ``\theta``, which are computationally difficult for large ``n``, therefor SWEEP (see https://github.com/joshday/SweepOperator.jl) algorithm used to meke oprtimization less computationaly expensive.
+
 #### V
 
 ```math
@@ -33,7 +60,7 @@ logREML(\theta,\beta) = -\frac{N-p}{2} - \frac{1}{2}\sum_{i=1}^nlog|V_{\theta, i
 -\frac{1}{2}log|\sum_{i=1}^nX_i'V_{\theta, i}^{-1}X_i|-\frac{1}{2}\sum_{i=1}^n(y_i - X_{i}\beta)'V_{\theta, i}^{-1}(y_i - X_{i}\beta)
 ```
 
-Actually ```L(\theta) = -2logREML = L_1(\theta) + L_2(\theta) + L_3(\theta) + c``` used for optimization, where:
+Actually `` L(\theta) = -2logREML = L_1(\theta) + L_2(\theta) + L_3(\theta) + c`` used for optimization, where:
 
 ```math
 L_1(\theta) = \frac{1}{2}\sum_{i=1}^nlog|V_{i}| \\
@@ -51,16 +78,36 @@ L_3(\theta) = \frac{1}{2}\sum_{i=1}^n(y_i - X_{i}\beta)'V_i^{-1}(y_i - X_{i}\bet
 \mathcal{H}\mathcal{L}(\theta) =  \mathcal{H}L_1(\theta)  + \mathcal{H}L_2(\theta) +  \mathcal{H} L_3(\theta)
 ```
 
-#### Weights
+#### [Weights](@id weights_header)
 
 If weights defined:
 
 ```math
+
+W_{i} = diag(wts_{i})
+
+\\
+
 V_{i} = Z_{i} G Z_i'+ W^{- \frac{1}{2}}_i R_{i} W^{- \frac{1}{2}}_i
 ```
 
+where ``W`` - diagonal matrix of weights.
+
+
+If wts is matrix then:
+
+
+```math
+
+W_{i} = wts_{i}
+
+\\
+
+V_{i} = Z_{i} G Z_i'+  R_{i} \circ W_{i}
+```
+
+where ``\circ`` - element-wise multiplication.
 
-where ```W``` - diagonal matrix of weights.
 
 
 #### Multiple random and repeated effects

docs/src/index.md

@@ -38,7 +38,7 @@ Implemented covariance structures:
 
 Actually Metida can fit datasets with wore than  160k observation and 40k subjects levels on PC with 64 GB RAM. This is not "hard-coded" limitation, but depends on your model and data structure. Fitting of big datasets can take a lot of time. Optimal dataset size is less than 100k observations with maximum length of block less than 400.
 
-!!! note
+!!! warning
 
     Julia v1.8 or higher required.
 

src/Metida.jl

@@ -17,8 +17,10 @@ import StatsBase: fit, fit!, coef, coefnames, confint, nobs, dof_residual, dof,
 import Base: show, rand, ht_keyindex, getproperty
 import Random: default_rng, AbstractRNG, rand!
 
+
 export @formula, @covstr, @lmmformula,
 SI, ScaledIdentity,
+SWC, ScaledWeightedCov,
 DIAG, Diag,
 AR, Autoregressive,
 ARH, HeterogeneousAutoregressive,

src/fvalue.jl

@@ -1,8 +1,5 @@
 # fvalue.jl
 
-#=
-Metida.fvalue(lmm, [0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0])
-=#
 """
     fvalue(lmm::LMM, l::Matrix)
 

src/gmat.jl

@@ -9,7 +9,7 @@
             gmat!(gt[r].data, view(θ, covstr.tr[r]), covstr.random[r].covtype.s)
         end
     end
-    gt
+    return gt
 end
 # Main
 @noinline function zgz_base_inc!(mx::AbstractArray, G, covstr, bi)
@@ -23,7 +23,7 @@ end
             end
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 ################################################################################
@@ -32,22 +32,22 @@ function gmat!(::Any, ::Any, ::AbstractCovarianceType)
     error("No gmat! method defined for thit structure!")
 end
 function gmat!(mx, ::Any, ::ZERO)
-    mx
+    return mx
 end
 #SI
 Base.@propagate_inbounds function gmat!(mx, θ, ::SI_)
     val = θ[1] ^ 2
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] = val
     end
-    mx
+    return mx
 end
 #DIAG
 function gmat!(mx, θ, ::DIAG_)
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] = θ[i] ^ 2
     end
-    mx
+    return mx
 end
 #AR
 function gmat!(mx, θ, ::AR_)
@@ -64,7 +64,7 @@ function gmat!(mx, θ, ::AR_)
             end
         end
     end
-    mx
+    return mx
 end
 #ARH
 function gmat!(mx, θ, ::ARH_)
@@ -84,7 +84,7 @@ function gmat!(mx, θ, ::ARH_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m] 
     end
-    mx
+    return mx
 end
 #CS
 function gmat!(mx, θ, ::CS_)
@@ -99,7 +99,7 @@ function gmat!(mx, θ, ::CS_)
             end
         end
     end
-    mx
+    return mx
 end
 #CSH
 function gmat!(mx, θ, ::CSH_)
@@ -118,7 +118,7 @@ function gmat!(mx, θ, ::CSH_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 ################################################################################
 #ARMA
@@ -136,7 +136,7 @@ function gmat!(mx, θ, ::ARMA_)
             end
         end
     end
-    mx
+    return mx
 end
 #TOEP
 function gmat!(mx, θ, ::TOEP_)
@@ -152,7 +152,7 @@ function gmat!(mx, θ, ::TOEP_)
             end
         end
     end
-    mx
+    return mx
 end
 function gmat!(mx, θ, ct::TOEPP_)
     de  = θ[1] ^ 2    #diagonal element
@@ -167,7 +167,7 @@ function gmat!(mx, θ, ct::TOEPP_)
             end
         end
     end
-    mx
+    return mx
 end
 #TOEPH
 function gmat!(mx, θ, ::TOEPH_)
@@ -186,7 +186,7 @@ function gmat!(mx, θ, ::TOEPH_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 #TOEPHP
 function gmat!(mx, θ, ct::TOEPHP_)
@@ -205,7 +205,7 @@ function gmat!(mx, θ, ct::TOEPHP_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 #UN
 function gmat!(mx, θ, ::UN_)
@@ -224,7 +224,7 @@ function gmat!(mx, θ, ::UN_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 
 function tpnum(m, n, s)
@@ -233,4 +233,5 @@ function tpnum(m, n, s)
         b += s - i
     end
     b -= s - n
+    return b
 end