pharmcat / metida.jl Goto Github PK

View Code? Open in Web Editor NEW

24.0 2.0 0.0 4.6 MB

Julia package for fitting mixed-effects models with flexible random/repeated covariance structure.

Home Page: http://metidajl.org

License: MIT License

Julia 100.00%

variance-component mixed-models statistics julia mixed-effects statistical-models regression mixed-effects-models

metida.jl's Introduction

Metida

This program comes with absolutely no warranty. No liability is accepted for any loss and risk to public health resulting from use of this software.

Version	pkgeval	Cover	Build	Docs

Metida.jl is Julia package for fitting mixed-effects models with flexible covariance structure.

Install:

import Pkg; Pkg.add("Metida")

Using:

using Metida, CSV, DataFrames, CategoricalArrays
df = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","df0.csv")) |> DataFrame
transform!(df, :subject => categorical, renamecols=false)
transform!(df, :period => categorical, renamecols=false)
transform!(df, :sequence => categorical, renamecols=false)
transform!(df, :formulation => categorical, renamecols=false)

lmm = LMM(@formula(var~sequence+period+formulation), df;
random = VarEffect(@covstr(formulation|subject), CSH),
repeated = VarEffect(@covstr(formulation|subject), DIAG),
)

fit!(lmm)

# Or you can use macro @lmmformula

lmm = LMM(@lmmformula(var~sequence+period+formulation,
    random = formulation|subject:CSH,
    repeated = formulation|subject:DIAG),
    df0)
fit!(lmm)

Also you can use this package with MatidaNLopt.jl and MetidaCu.jl.

See also MixedModels.jl: powerful package for mixed models.

metida.jl's People

Contributors

Stargazers

Watchers

metida.jl's Issues

Rank deficient fixed effect matrix

Now it impossible to fit model with rank deficient fixed effect matrix. There are two ways to fix it:

drop compatibility for StatsModels
make re-implementation of StatsModels

First lead to lost compatibility with packages like Effects.jl
Second take tonnes of "hardcoding".

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Add to model result check for positive definite G matrix

Check rank before variance-covariance matrix of β calc

Multiple repeated effects

There is no need to use multiple repeated effects in most cases, but in some rare cases, it can be useful.

Remove ModelFrame from model object

as deprecated

1.6 test error on experimental functions

fitting LMM with given kinship matrix

Hi. Is there any detailed documetation about the parameters used in Metida? Like, if I have a kinship matrix, which commonly used in genetic studies reflects the pairwised genetic relatedness, and i want to include it into LMM as a random effect, how could I possibly do it?

Option for AR(1)

Is there an option for using the Auto Regressive structure (1)?
It seems to default to AR(2).

Error specifying formula and model fit

Hi, new to Julia coming from R, I have a question about an error message I get trying to apply a mixed effects model using Metida on my data.

using Metida, CSV, DataFrames, CategoricalArrays
df = CSV.File("../data/mydata.csv") |> DataFrame;

transform!(df, :Peptide => categorical, renamecols=false);
transform!(df, :Protein => categorical, renamecols=false);
transform!(df, :Treatment => categorical, renamecols=false);
transform!(df, :Subject => categorical, renamecols=false);

df
24,104 rows × 5 columns
Peptide | Protein | Measurement | Treatment | Subject
Cat...  | Cat...  | Float64?..  | Cat...    | Cat...

AACAQLNDFLQEYGTQGC | P0C0L4 | missing | G | 6
AACAQLNDFLQEYGTQGC | P0C0L5 | 22.3928 | G | 6
...

# This is the full model I want to specify
lmm = LMM(@formula(Measurement ~ Treatment * Peptide + Protein), df;
random = VarEffect(@covstr(1 | Subject / Protein), CSH));

MethodError: no method matching fulldummycodingdict(::FunctionTerm{typeof(/), var"#68#70", (:Subject, :Protein)})
...

# Removing the nested term lets me specify the formula
lmm = LMM(@formula(Measurement ~ Treatment * Peptide + Protein), df;
random = VarEffect(@covstr(1 | Subject), CSH));

# But the fit throws an error
fit!(lmm)

BoundsError: attempt to access 20683×1109 Matrix{Float64} at index [UInt32[0x00004dc9, 0x00004dca, 0x00004dcb, 0x00004dcc, 0x00004dcd, 0x00004dce, 0x00004dcf, 0x00004dd0, 0x00004dd1, 0x00004dd2  …  0x000051d7, 0x000051d8, 0x000051d9, 0x000051da, 0x000051db, 0x000051dc, 0x000051dd, 0x000051de, 0x000051df, 0x000051e0], 1:1109]

Stacktrace:
  [1] throw_boundserror(A::Matrix{Float64}, I::Tuple{Vector{UInt32}, Base.Slice{Base.OneTo{Int64}}})
...


# This works in MixedModels.jl, perhaps the nested terms have different syntax i Metida?
fm = @formula(Measurement ~ 1 + Treatment*Peptide + Protein + (1 | Subject/Protein))
mm = fit(MixedModel, fm, df, REML=true)

However in MixedModels.jl I am unable to get degrees of freedom and confidence intervals for my model.

Ultimately I want to obtain pairwise contrasts from my model. In R using lme4 and emmeans I would do this:

lmm = lmer(Measurement ~ Treatment  *  Peptide + Protein + (1  |  Subject / Protein), mydata)

emm <- emmeans(lmm, pairwise ~ Treatment | Peptide)
emm <- confint(emm)
emm <- as.data.frame(emm$contrasts)

Any help appreciated