Custom structures

To make your own covariance structure first you should make struct that <: AbstractCovarianceType:

Example:

struct YourCovarianceStruct <: AbstractCovarianceType end

You also can specify additional field if you need to use them inside gmat!/rmat! functions.

Then you can make function for construction random effect matrix (gmat!) and repeated effect (rmat!). Only upper triangular can be updated.

Function gmat! have 3 arguments: mx - zero matrix, θ - theta vector for this effect, and your custom structure object.

Next this function used to make "random" part of variance-covariance matrix: V' = Z * G * Z'

function Metida.gmat!(mx, θ, ::YourCovarianceStruct)
    @inbounds @simd for i = 1:size(mx, 1)
        mx[i, i] = θ[i] ^ 2
    end
    nothing
end

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_, ::Int)
    l     = size(rz, 2)
    vec   = rz * (θ[1:l])
    s   = size(mx, 1)
    if s > 1 && ct.p > 1
        for m = 1:s - 1
            for n = m + 1:(m + ct.p - 1 > s ? s : m + ct.p - 1)
                @inbounds  mx[m, n] += vec[m] * vec[n] * θ[n - m + l]
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] += vec[m] * vec[m]
    end
    nothing
end

One more function you shoud make is covstrparam, this function need to know how many parameters included in theta vector for optimization. Function returns number of variance parameters and rho parameters for this structure. Where t - number of columns in individual Z matrix for random effect or number of columns in repeated effect matrix (rZ).

Example for Heterogeneous Autoregressive and Heterogeneous Compound Symmetry structures:

function Metida.covstrparam(ct::Union{ARH_, CSH_}, t::Int)::Tuple{Int, Int}
    return (t, 1)
end

For better printing you can add:

function Metida.rcoefnames(s, t, ct::YourCovarianceStruct)
    return ["σ² ", "γ ", "ρ "]
end

Where, s - effect schema, t - number of parameters, this function returns names for your covariance structure for printing in LMM output.

Add this method for better printing:

function Base.show(io::IO, ct::YourCovarianceStruct)
    print(io, "YourCovarianceStruct")
end

Then just make model and fit it:

lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;
  random = Metida.VarEffect(Metida.@covstr(factor|subject&factor), YourCovarianceStruct()),
  repeated = Metida.VarEffect(Metida.@covstr(1|subject&factor), YourCovarianceStruct()),
  )
  Metida.fit!(lmm)

Example:

using Metida, DataFrames, CSV, CategoricalArrays

spatdf       = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "spatialdata.csv"); types = [Int, Int, String, Float64, Float64, Float64, Float64, Float64]) |> DataFrame
ftdf = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame
df0 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv",  "df0.csv"); types = [String, String, String, String,Float64, Float64, Float64]) |> DataFrame

struct CustomCovarianceStructure <: Metida.AbstractCovarianceType end
function Metida.covstrparam(ct::CustomCovarianceStructure, t::Int)::Tuple{Int, Int}
    return (t, 1)
end
function Metida.gmat!(mx, θ, ct::CustomCovarianceStructure)
    s = size(mx, 1)
    @inbounds @simd for m = 1:s
        mx[m, m] = θ[m]
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = mx[m, m] * mx[n, n] * θ[end]
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] = mx[m, m] * mx[m, m]
    end
    nothing
end

lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), CustomCovarianceStructure()),
)
Metida.fit!(lmm)

# for R matrix

function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure, ::Int)
    vec = Metida.tmul_unsafe(rz, θ)
    rn    = size(mx, 1)
    if rn > 1
        for m = 1:rn - 1
            @inbounds @simd for n = m + 1:rn
                mx[m, n] += vec[m] * vec[n] * θ[end]
            end
        end
    end
        @inbounds  for m ∈ axes(mx, 1)
        mx[m, m] += vec[m] * vec[m]
    end
    nothing
end

lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
repeated = Metida.VarEffect(Metida.@covstr(period|subject), CustomCovarianceStructure()),
)
Metida.fit!(lmm)

Linear Mixed Model: var ~ sequence + period + formulation
Random 1: 
   No
Repeated: 
    Model: period|subject
    Type: CustomCovarianceStructure (5)
Blocks: 5, Maximum block size: 4
Status: converged (No Errors)
    -2 logREML: 8.7401    BIC: 22.603

    Fixed-effects parameters:
───────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)      1.31285     0.28779     4.56    <1e-05
sequence: 2      0.283193    0.328706    0.86    0.3889
period: 2        0.152412    0.0881943   1.73    0.0840
period: 3        0.08        0.136257    0.59    0.5571
period: 4        0.152412    0.0986418   1.55    0.1223
formulation: 2  -0.0379392   0.0717112  -0.53    0.5968
───────────────────────────────────────────────────────
    Variance components:
    θ vector: [0.43779, 0.426869, 0.630382, 0.49783, 0.899017]
  Residual   Val    var   0.19166
  Residual   Val    var   0.182217
  Residual   Val    var   0.397381
  Residual   Val    var   0.247835
  Residual   Val    rho   0.899017

Custom distance estimation for spatial structures

If you want to use coordinates or some other structures for distance estimation you can define method Metida.edistance to calculate distance:

function Metida.edistance(mx::AbstractMatrix{<:CartesianIndex}, i::Int, j::Int)
    return sqrt((mx[i, 1][1] - mx[j, 1][1])^2 + (mx[i, 1][2] - mx[j, 1][2])^2)
end

For example this method returns distance between two vectors represented as CartesianIndex.

Make vector of CartesianIndex:

spatdf.ci = map(x -> CartesianIndex(x[:x], x[:y]), eachrow(spatdf))

200-element Vector{CartesianIndex{2}}:
 CartesianIndex(8, 94)
 CartesianIndex(9, 40)
 CartesianIndex(40, 62)
 CartesianIndex(49, 24)
 CartesianIndex(20, 89)
 CartesianIndex(57, 14)
 CartesianIndex(44, 68)
 CartesianIndex(64, 90)
 CartesianIndex(30, 81)
 CartesianIndex(79, 32)
 ⋮
 CartesianIndex(100, 45)
 CartesianIndex(55, 59)
 CartesianIndex(11, 22)
 CartesianIndex(37, 25)
 CartesianIndex(57, 8)
 CartesianIndex(81, 99)
 CartesianIndex(68, 83)
 CartesianIndex(89, 17)
 CartesianIndex(62, 79)

Then use new column as "raw" variable with Metida.RawCoding contrast and fit the model:

lmm = Metida.LMM(@formula(r2 ~ f), spatdf;
    repeated = Metida.VarEffect(Metida.@covstr(ci|1), Metida.SPEXP; coding = Dict(:ci => Metida.RawCoding())),
    )
Metida.fit!(lmm)

Linear Mixed Model: r2 ~ f
Random 1: 
   No
Repeated: 
    Model: ci|1
    Type: SPEXP (2)
Blocks: 1, Maximum block size: 200
Status: converged (No Errors)
    -2 logREML: 1985.34    BIC: 1995.92

    Fixed-effects parameters:
────────────────────────────────────────────────
               Coef.  Std. Error     z  Pr(>|z|)
────────────────────────────────────────────────
(Intercept)  6.26619     8.841    0.71    0.4785
f: 1         3.85137     4.12172  0.93    0.3501
────────────────────────────────────────────────
    Variance components:
    θ vector: [46.8999, 8.28173]
  Residual   σ²    var     2199.6
  Residual   θ     theta   8.28173