Example 1 - Continuous and categorical predictors
using Metida, CSV, DataFrames, CategoricalArrays, MixedModels;
rds = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame
Metida result:
lmm = LMM(@formula(response ~1 + factor*time), rds;
random = VarEffect(@covstr(1 + time|subject&factor), CSH),
)
fit!(lmm)
Linear Mixed Model: response ~ 1 + factor + time + factor & time
Random 1:
Model: :(1 + time)|:(subject & factor)
Type: CSH (3), Subjects: 40
Repeated:
Residual only
Blocks: 40, Maximum block size: 10
Status: converged (No Errors)
-2 logREML: 1300.18 BIC: 1324.11
Fixed-effects parameters:
──────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────────────
(Intercept) 22.1331 0.304339 72.73 <1e-99
factor: 1.0 2.00049 0.430401 4.65 <1e-05
time 1.11853 0.0264038 42.36 <1e-99
factor: 1.0 & time 0.404971 0.0373406 10.85 <1e-26
──────────────────────────────────────────────────────────
Variance components:
θ vector: [1.2401, 0.101122, -0.0427208, 0.985061]
Random 1 σ² (Intercept) var 1.53785
Random 1 σ² time var 0.0102257
Random 1 ρ rho -0.0427208
Residual σ² var 0.970345
MixedModels result:
fm = @formula(response ~ 1 + factor*time + (1 + time|subject&factor))
mm = fit(MixedModel, fm, rds, REML=true)
Minimizing 2 Time: 0:00:00 ( 0.22 s/it)
objective: 1426.9231858044589
Minimizing 68 Time: 0:00:00 ( 6.97 ms/it)
Linear mixed model fit by REML
response ~ 1 + factor + time + factor & time + (1 + time | subject & factor)
REML criterion at convergence: 1300.180759818896
Variance components:
Column Variance Std.Dev. Corr.
subject & factor (Intercept) 1.537850 1.240101
time 0.010226 0.101122 -0.04
Residual 0.970345 0.985061
Number of obs: 400; levels of grouping factors: 40
Fixed-effects parameters:
──────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────────────
(Intercept) 22.1331 0.304339 72.73 <1e-99
factor: 1.0 2.00049 0.4304 4.65 <1e-05
time 1.11853 0.0264037 42.36 <1e-99
factor: 1.0 & time 0.404971 0.0373405 10.85 <1e-26
──────────────────────────────────────────────────────────
Example 2 - Two random factors (Penicillin data)
Metida:
df = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "Penicillin.csv"); types = [String, Float64, String, String]) |> DataFrame
df.diameter = float.(df.diameter)
lmm = LMM(@formula(diameter ~ 1), df;
random = [VarEffect(@covstr(1|plate), SI), VarEffect(@covstr(1|sample), SI)]
)
fit!(lmm)
Linear Mixed Model: diameter ~ 1
Random 1:
Model: 1|plate
Type: SI (1), Subjects: 24
Random 2:
Model: 1|sample
Type: SI (1), Subjects: 6
Repeated:
Residual only
Blocks: 1, Maximum block size: 144
Status: converged (No Errors)
-2 logREML: 330.861 BIC: 345.749
Fixed-effects parameters:
─────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
─────────────────────────────────────────────────
(Intercept) 22.9722 0.808573 28.41 <1e-99
─────────────────────────────────────────────────
Variance components:
θ vector: [0.846704, 1.93156, 0.549923]
Random 1 σ² var 0.716908
Random 2 σ² var 3.73092
Residual σ² var 0.302415
MixedModels:
fm2 = @formula(diameter ~ 1 + (1|plate) + (1|sample))
mm = fit(MixedModel, fm2, df, REML=true)
Minimizing 2 Time: 0:00:00 ( 0.29 s/it)
objective: 366.5302937135882
Minimizing 46 Time: 0:00:00 (12.69 ms/it)
Linear mixed model fit by REML
diameter ~ 1 + (1 | plate) + (1 | sample)
REML criterion at convergence: 330.86058899096145
Variance components:
Column Variance Std.Dev.
plate (Intercept) 0.716908 0.846704
sample (Intercept) 3.730907 1.931556
Residual 0.302415 0.549923
Number of obs: 144; levels of grouping factors: 24, 6
Fixed-effects parameters:
─────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
─────────────────────────────────────────────────
(Intercept) 22.9722 0.808572 28.41 <1e-99
─────────────────────────────────────────────────
Example 3 - Repeated ARMA/AR/ARH
rds = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "1freparma.csv"); types = [String, String, Float64, Float64]) |> DataFrame
ARMA:
lmm = LMM(@formula(response ~ 1 + factor*time), rds;
random = VarEffect(@covstr(factor|subject&factor), DIAG),
repeated = VarEffect(@covstr(1|subject&factor), ARMA),
)
fit!(lmm)
Linear Mixed Model: response ~ 1 + factor + time + factor & time
Random 1:
Model: factor|:(subject & factor)
Type: DIAG (2), Subjects: 24
Repeated:
Model: 1|:(subject & factor)
Type: ARMA (3)
Blocks: 24, Maximum block size: 10
Status: converged See warnings in log.
-2 logREML: 348.11 BIC: 375.429
Fixed-effects parameters:
──────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────────────
(Intercept) 44.3198 0.718497 61.68 <1e-99
factor: 1.0 5.60925 0.979922 5.72 <1e-07
time 0.514843 0.0316349 16.27 <1e-58
factor: 1.0 & time 0.267151 0.0447385 5.97 <1e-08
──────────────────────────────────────────────────────────
Variance components:
θ vector: [2.39579, 2.20751, 0.773677, 0.621233, -0.814062]
Random 1 σ² factor: 0.0 var 5.73981
Random 1 σ² factor: 1.0 var 4.8731
Residual σ² var 0.598576
Residual γ rho 0.621233
Residual ρ rho -0.814062
AR:
lmm = Metida.LMM(@formula(response ~ 1 + factor*time), rds;
random = VarEffect(@covstr(factor|subject&factor), DIAG),
repeated = VarEffect(@covstr(1|subject&factor), AR),
)
fit!(lmm)
Linear Mixed Model: response ~ 1 + factor + time + factor & time
Random 1:
Model: factor|:(subject & factor)
Type: DIAG (2), Subjects: 24
Repeated:
Model: 1|:(subject & factor)
Type: AR (2)
Blocks: 24, Maximum block size: 10
Status: converged (No Errors)
-2 logREML: 710.096 BIC: 731.952
Fixed-effects parameters:
──────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────────────
(Intercept) 44.3915 1.34686 32.96 <1e-99
factor: 1.0 5.39757 1.60271 3.37 0.0008
time 0.508075 0.0456524 11.13 <1e-28
factor: 1.0 & time 0.290317 0.0645622 4.50 <1e-05
──────────────────────────────────────────────────────────
Variance components:
θ vector: [4.54797, 2.82342, 1.05771, 0.576979]
Random 1 σ² factor: 0.0 var 20.684
Random 1 σ² factor: 1.0 var 7.9717
Residual σ² var 1.11875
Residual ρ rho 0.576979
ARH:
lmm = Metida.LMM(@formula(response ~ 1 + factor*time), rds;
random = VarEffect(@covstr(factor|subject&factor), DIAG),
repeated = VarEffect(@covstr(1|subject&factor), ARH),
)
fit!(lmm)
Linear Mixed Model: response ~ 1 + factor + time + factor & time
Random 1:
Model: factor|:(subject & factor)
Type: DIAG (2), Subjects: 24
Repeated:
Model: 1|:(subject & factor)
Type: ARH (2)
Blocks: 24, Maximum block size: 10
Status: converged (No Errors)
-2 logREML: 710.096 BIC: 731.952
Fixed-effects parameters:
──────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────────────
(Intercept) 44.3915 1.34686 32.96 <1e-99
factor: 1.0 5.39757 1.60271 3.37 0.0008
time 0.508075 0.0456524 11.13 <1e-28
factor: 1.0 & time 0.290317 0.0645622 4.50 <1e-05
──────────────────────────────────────────────────────────
Variance components:
θ vector: [4.54797, 2.82342, 1.05771, 0.576979]
Random 1 σ² factor: 0.0 var 20.684
Random 1 σ² factor: 1.0 var 7.9717
Residual σ² (Intercept) var 1.11875
Residual ρ rho 0.576979
Example 4 - SAS relation
Model 1
df0 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "df0.csv")) |> DataFrame
lmm = LMM(@formula(var ~ sequence + period + formulation), df0;
random = VarEffect(@covstr(formulation|subject), CSH),
repeated = VarEffect(@covstr(formulation|subject), DIAG),
)
fit!(lmm)
SAS code:
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM formulation/TYPE=CSH SUB=subject G V;
REPEATED/GRP=formulation SUB=subject R;
RUN;
Model 2
lmm = LMM(
@formula(var ~ sequence + period + formulation), df0;
random = VarEffect(@covstr(formulation|subject), SI),
repeated = VarEffect(@covstr(formulation|subject), DIAG),
)
fit!(lmm)
SAS code:
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM formulation/TYPE=VC SUB=subject G V;
REPEATED/GRP=formulation SUB=subject R;
RUN;
Model 3
lmm = LMM(@formula(var ~ sequence + period + formulation), df0;
random = VarEffect(@covstr(subject|1), SI)
)
fit!(lmm)
SAS code:
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM subject/TYPE=VC G V;
RUN;
Example 5 - Working with Effects.jl
using Effects, StatsModels
lmm = LMM(@formula(var ~ sequence + period + formulation), df0;
random = VarEffect(@covstr(subject|1), SI)
)
fit!(lmm)
table_model = StatsModels.TableRegressionModel(lmm, lmm.mf, lmm.mm)
emmeans(tm)
effects(Dict(:period => ["1", "2", "3", "4"]), tm)