Merge branch 'master' of github.com:JuliaStats/MixedModels.jl into gl…

…mm_dispersion_deviance

.github/workflows/MKL.yml

@@ -0,0 +1,22 @@
+name: MKL
+on:
+  workflow_dispatch:
+jobs:
+  ci:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        julia-version: [1.5]
+        julia-arch: [x64]
+        os: [ubuntu-18.04]
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.julia-version }}
+      - uses: julia-actions/julia-buildpkg@v0.1
+      - run:  julia -e 'using Pkg; Pkg.add("MKL"); Pkg.build("MKL"; verbose=true)'
+      - uses: julia-actions/julia-runtest@v0.1
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

Artifacts.toml

@@ -5,21 +5,21 @@
 [TestData]
 # compute this using
 # using Tar, Inflate, SHA
-# filename = "download?version=4" # I just used wget for the URL below and this is how it saved it
+# filename = "download?version=5" # I just used wget for the URL below and this is how it saved it
 # println("sha256: ", bytes2hex(open(sha256, filename)))
 # println("git-tree-sha1: ", Tar.tree_hash(IOBuffer(inflate_gzip(filename))))
 # from https://julialang.github.io/Pkg.jl/dev/artifacts/
-git-tree-sha1 = "75260d131b693f26e5834adf855f4c35d627348d"
-lazy = true
+git-tree-sha1 = "91132469677f725c2e4097493ae8b1d566f90a3f"
+lazy = false
 
     [[TestData.download]]
-    # this is the SHA from https://osf.io/djaqb/download?version=4
-    sha256 = "d179cadfb27fc638fd8a06b0d58f5f2916f786c5719ab22de93fedb6129726f4"
+    # this is the SHA from https://osf.io/djaqb/download?version=5
+    sha256 = "8040933246179d4b46cf37e4dd1076752102dfcb4bd937d5f12c0f724775e119"
     # when updating this, make sure to change to change the version number,
     # because if the version number isn't included, it will always point to the
     # latest version, which means it will break existing users when we update
     # between releases.
-    url = "https://osf.io/djaqb/download?version=4"
+    url = "https://osf.io/djaqb/download?version=5"
 
     # for future work on using xz-compressed data:
     # Julia invokes wget without using HTTP metadata, so we need the link

NEWS.md

@@ -69,7 +69,7 @@ Availability of test data sets
 * Several data sets from the literature were previously saved in `.rda` format
   in the `test` directory and read using the `RData` package. These are now available
   in an `Artifact` in the [`Arrow`](https://github.com/JuliaData/Arrow.jl.git) format [#382].
-* Call `MixedModel.datasets()` to get a listing of the names of available datasets
+* Call `MixedModels.datasets()` to get a listing of the names of available datasets
 * To load, e.g. the `dyestuff` data, use `MixedModels.dataset(:dyestuff)`
 * Data sets are loaded using `Arrow.Table` which returns a column table.  Wrap the
   call in `DataFrame` if you prefer a `DataFrame`.

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>", "Jose Bayoan Santiago Calderon <jbs3hp@virginia.edu>"]
-version = "3.1.0"
+version = "3.1.2"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
@@ -24,15 +24,15 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
-Arrow = "0.3, 0.4"
-BlockArrays = "0.11, 0.12"
-DataAPI = "1.1, 1.2, 1.3"
+Arrow = "1"
+BlockArrays = "0.11, 0.12, 0.13"
+DataAPI = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24"
 GLM = "1"
 NLopt = "0.5, 0.6"
 PooledArrays = "0.5"
 ProgressMeter = "1"
-StaticArrays = "0.11, 0.12"
+StaticArrays = "0.11, 0.12, 1"
 StatsBase = "0.31, 0.32, 0.33"
 StatsFuns = "0.8, 0.9"
 StatsModels = "0.6"

benchmark/benchmarks.jl

@@ -1,171 +1,110 @@
-using BenchmarkTools, DataFrames, MixedModels, RData, Tables
+using BenchmarkTools, MixedModels, StatsModels
+using MixedModels: dataset
 
 const SUITE = BenchmarkGroup()
 
-const dat = Dict(Symbol(k) => v for (k, v) in load(joinpath(
-    dirname(pathof(MixedModels)),
-    "..",
-    "test",
-    "dat.rda",
-)));
-
-categorical!(dat[:ml1m], [:G,:H]);  # forgot to convert these grouping factors
-
-const mods = Dict{Symbol,Vector{Expr}}(
-    :Alfalfa => [:(1 + A * B + (1 | G)), :(1 + A + B + (1 | G))],
-    :Animal => [:(1 + (1 | G) + (1 | H))],
-    :Arabidopsis => [],              # glmm and rename variables
-    :Assay => [:(1+A+B*C+(1|G)+(1|H))],
-    :AvgDailyGain => [:(1 + A * U + (1 | G)), :(1 + A + U + (1 | G))],
-    :BIB => [:(1 + A * U + (1 | G)), :(1 + A + U + (1 | G))],
-    :Bond => [:(1 + A + (1 | G))],
-    :Chem97 => [:(1 + (1 | G) + (1 | H)), :(1 + U + (1 | G) + (1 | H))],
-    :Contraception => [],            # glmm and rename variables
-    :Cultivation => [:(1 + A * B + (1 | G)), :(1 + A + B + (1 | G)), :(1 + A + (1 | G))],
-    :Demand => [:(1 + U + V + W + X + (1 | G) + (1 | H))],
-    :Dyestuff => [:(1 + (1 | G))],
-    :Dyestuff2 => [:(1 + (1 | G))],
-    :Early => [:(1 + U + U & A + (1 + U | G))],
-    :Exam => [:(1 + A * U + B + (1 | G)), :(1 + A + B + U + (1 | G))],
-    :Gasoline => [:(1 + U + (1 | G))],
-    :Gcsemv => [:(1 + A + (1 | G))], # variables must be renamed
-    :Genetics => [:(1 + A + (1 | G) + (1 | H))],
-    :HR => [:(1 + A * U + V + (1 + U | G))],
-    :Hsb82 => [:(1 + A + B + C + U + (1 | G))],
-    :IncBlk => [:(1 + A + U +  + W + Z + (1 | G))],
-    :InstEval => [:(1 + A + (1 | G) + (1 | H) + (1 | I)), :(1 + A * I + (1 | G) + (1 | H))],
-    :KKL => [],                      # variables must be renamed
-    :KWDYZ => [],                    # variables must be renamed
-    :Mississippi => [:(1 + A + (1 | G))],
-    :Mmmec => [],                    # glmm (and offset) and variables renamed
-    :Multilocation => [:(1 + A + (0 + A | G) + (1 | H))],
-    :Oxboys => [:(1 + U + (1 + U | G))],
-    :PBIB => [:(1 + A + (1 | G))],
-    :Pastes => [:(1 + (1 | G) + (1 | H))],
-    :Penicillin => [:(1 + (1 | G) + (1 | H))],
-    :Pixel => [:(1 + U + V + (1 + U | G) + (1 | H))],  # variables must be renamed
-    :Poems => [:(1 + U + V + W + (1 | G) + (1 | H) + (1 | I))],
-    :Rail => [:(1 + (1 | G))],
-    :SIMS => [:(1 + U + (1 + U | G))],
-    :ScotsSec => [:(1 + A + U + V + (1 | G) + (1 | H))],
-    :Semi2 => [:(1 + A + (1 | G) + (1 | H))],
-    :Semiconductor => [:(1 + A * B + (1 | G))],
-    :Socatt => [],                   # variables must be renamed - binomial glmm?
-    :TeachingII => [:(1 + A + T + U + V + W + X + Z + (1 | G))],
-    :VerbAgg => [:(1 + A + B + C + U + (1 | G) + (1 | H))], # Bernoulli glmm and rename variables
-    :Weights => [:(1 + A * U + (1 + U | G))],
-    :WWheat => [:(1 + U + (1 + U | G))],
-    :bdf => [],                      # rename variables and look up model
-    :bs10 => [:(1 + U + V + W + ((1 + U + V + W) | G) + ((1 + U + V + W) | H))],
-    :cake => [:(1 + A * B + (1 | G))],
-    :cbpp => [:(1 + A + (1 | G))],   # Binomial glmm, create and rename variables
+const global contrasts = Dict(
+    :mrk17_exp1 => merge(Dict(n => HelmertCoding() for n in (:F, :P, :Q, :lQ, :lT)),
+     Dict(n => Grouping() for n in (:item, :subj))),
+)
+const global fms = Dict(
+    :dyestuff => [
+        @formula(yield ~ 1 + (1|batch)),
+        ],
+    :dyestuff2 => [
+        @formula(yield ~ 1 + (1|batch)),
+        ],
     :d3 => [
-        :(1 + U + (1 | G) + (1 | H) + (1 | I)),
-        :(1 + U + (1 + U | G) + (1 + U | H) + (1 + U | I)),
-    ],
-    :dialectNL => [:(1 + A + T + U + V + W + X + (1 | G) + (1 | H) + (1 | I))],
-    :egsingle => [:(1 + A + U + V + (1 | G) + (1 | H))],
-    :epilepsy => [],                 # unknown origin
-    :ergoStool => [:(1 + A + (1 | G))],
-    :gb12 => [:(1 + S + T + U + V + W + X + Z + ((1 + S + U + W) | G) +
-                ((1 + S + T + V) | H))],
-    :grouseticks => [],              # rename variables, glmm needs formula
-    :guImmun => [],                  # rename variables, glmm needs formula
-    :guPrenat => [],                 # rename variables, glmm needs formula
+        @formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i)),
+        ],
+    :insteval => [
+        @formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)),
+        @formula(y ~ 1 + service*dept + (1|s) + (1|d)),
+        ],
     :kb07 => [
-        :(1 + S + T + U + V + W + X + Z + ((1 + S + T + U + V + W + X + Z) | G) +
-          ((1 + S + T + U + V + W + X + Z) | H)),
-        :(1 + S + T + U + V + W + X + Z +
-          zerocorr((1 + S + T + U + V + W + X + Z) | G) +
-          zerocorr((1 + S + T + U + V + W + X + Z) | H)),
-    ],
-    :ml1m => [:(1 + (1 | G) + (1 | H))],
-    :paulsim => [:(1 + S + T + U + (1 | H) + (1 | G))],  # names of H and G should be reversed
-    :sleepstudy => [:(1 + U + (1 + U | G)), :(1 + U + zerocorr(1 + U | G))],
-    :s3bbx => [],                    # probably drop this one
-    :star => [],                     # not sure it is worthwhile working with these data
-);
+        @formula(rt_trunc ~ 1+spkr+prec+load+(1|subj)+(1|item)),
+        @formula(rt_trunc ~ 1+spkr*prec*load+(1|subj)+(1+prec|item)),
+        @formula(rt_trunc ~ 1+spkr*prec*load+(1+spkr+prec+load|subj)+(1+spkr+prec+load|item)),
+        ],
+    :machines => [
+        @formula(score ~ 1 + (1|Worker) + (1|Machine)),
+        ],
+    :ml1m => [
+        @formula(y ~ 1 + (1|g) + (1|h)),
+        ],
+    :mrk17_exp1 => [
+        @formula(1000/rt ~ 1+F*P*Q*lQ*lT + (1|item) + (1|subj)),
+        @formula(1000/rt ~ 1+F*P*Q*lQ*lT + (1+P+Q+lQ+lT|item) + (1+F+P+Q+lQ+lT|subj)),
+        ],
+    :pastes => [
+        @formula(strength ~ 1 + (1|batch&cask)),
+        @formula(strength ~ 1 + (1|batch/cask)),
+        ],
+    :penicillin => [
+        @formula(diameter ~ 1 + (1|plate) + (1|sample)),
+        ],
+    :sleepstudy => [
+        @formula(reaction ~ 1 + days + (1|subj)),
+        @formula(reaction ~ 1 + days + zerocorr(1+days|subj)),
+        @formula(reaction ~ 1 + days + (1|subj) + (0+days|subj)),
+        @formula(reaction ~ 1 + days + (1+days|subj)),
+        ],
+)
 
-fitbobyqa(rhs::Expr, dsname::Symbol) =
-    fit(MixedModel, @eval(@formula(Y ~ $rhs)), dat[dsname])
-compactstr(ds, rhs) = replace(string(ds, ':', rhs), ' ' => "")
+function fitbobyqa(dsname::Symbol, index::Integer)
+    fit(
+        MixedModel, 
+        fms[dsname][index], 
+        dataset(dsname), 
+        contrasts=get!(contrasts, dsname, Dict{Symbol,StatsModels.AbstractContrasts}()),
+        )
+end
 
 SUITE["simplescalar"] = BenchmarkGroup(["single", "simple", "scalar"])
-for ds in [
-    :Alfalfa,
-    :AvgDailyGain,
-    :BIB,
-    :Bond,
-    :cake,
-    :Cultivation,
-    :Dyestuff,
-    :Dyestuff2,
-    :ergoStool,
-    :Exam,
-    :Gasoline,
-    :Hsb82,
-    :IncBlk,
-    :Mississippi,
-    :PBIB,
-    :Rail,
-    :Semiconductor,
-    :TeachingII,
+for (ds, i) in [  
+    (:dyestuff, 1,),
+    (:dyestuff2, 1,),
+    (:pastes, 1),
+    (:sleepstudy, 1,),
 ]
-    for rhs in mods[ds]
-        SUITE["simplescalar"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
-            $(QuoteNode(rhs)),
-            $(QuoteNode(ds)),
-        )
-    end
+    SUITE["simplescalar"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
 end
 
 SUITE["singlevector"] = BenchmarkGroup(["single", "vector"])
-for ds in [:Early, :HR, :Oxboys, :SIMS, :sleepstudy, :Weights, :WWheat]
-    for rhs in mods[ds]
-        SUITE["singlevector"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
-            $(QuoteNode(rhs)),
-            $(QuoteNode(ds)),
-        )
-    end
+for (ds, i) in [  
+    (:sleepstudy, 2,),
+    (:sleepstudy, 3,),
+    (:sleepstudy, 4,),
+]
+    SUITE["singlevector"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
 end
 
 SUITE["nested"] = BenchmarkGroup(["multiple", "nested", "scalar"])
-for ds in [:Animal, :Chem97, :Genetics, :Pastes, :Semi2]
-    for rhs in mods[ds]
-        SUITE["nested"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
-            $(QuoteNode(rhs)),
-            $(QuoteNode(ds)),
-        )
-    end
+for (ds, i) in [  
+    (:pastes, 2,),
+]
+    SUITE["nested"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
 end
 
 SUITE["crossed"] = BenchmarkGroup(["multiple", "crossed", "scalar"])
-
-for ds in [
-    :Assay,
-    :Demand,
-    :InstEval,
-    :Penicillin,
-    :ScotsSec,
-    :dialectNL,
-    :egsingle,
-    :ml1m,
-    :paulsim,
+for (ds, i) in [
+    (:insteval, 1),
+    (:insteval, 2),
+    (:kb07, 1),
+    (:machines, 1),
+    (:ml1m, 1),
+    (:mrk17_exp1, 1),
+    (:penicillin, 1),
 ]
-    for rhs in mods[ds]
-        SUITE["crossed"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
-            $(QuoteNode(rhs)),
-            $(QuoteNode(ds)),
-        )
-    end
+    SUITE["crossed"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
 end
 
 SUITE["crossedvector"] = BenchmarkGroup(["multiple", "crossed", "vector"])
-for ds in [:bs10, :d3, :gb12, :kb07]
-    for rhs in mods[ds]
-        SUITE["crossedvector"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa(
-            $(QuoteNode(rhs)),
-            $(QuoteNode(ds)),
-        )
-    end
+for (ds, i) in [
+    (:d3, 1),
+    (:kb07, 2),
+    (:kb07, 3),
+    (:mrk17_exp1, 2),
+]
+    SUITE["crossedvector"][string(ds, ':', i)] = @benchmarkable fitbobyqa($(QuoteNode(ds)), $(QuoteNode(i)))
 end

docs/src/constructors.md

@@ -167,7 +167,7 @@ end
 
 ### Simplifying the random effect correlation structure
 
-MixedEffects.jl estimates not only the *variance* of the effects for each random effect level, but also the *correlation* between the random effects for different predictors.
+MixedModels.jl estimates not only the *variance* of the effects for each random effect level, but also the *correlation* between the random effects for different predictors.
 So, for the model of the *sleepstudy* data above, one of the parameters that is estimated is the correlation between each subject's random intercept (i.e., their baseline reaction time) and slope (i.e., their particular change in reaction time per day of sleep deprivation).
 In some cases, you may wish to simplify the random effects structure by removing these correlation parameters.
 This often arises when there are many random effects you want to estimate (as is common in psychological experiments with many conditions and covariates), since the number of random effects parameters increases as the square of the number of predictors, making these models difficult to estimate from limited data.

src/MixedModels.jl

@@ -18,6 +18,14 @@ using StatsBase
 using StatsModels
 using Tables
 
+# When we move to 1.6 as the support lower minimum, we should change Artifact.toml to be lazy
+# and add LazyArtifacts to our dependencies
+# @static if VERSION > v"1.6.0-DEV.1588" # the actual bound may be lower
+#     @warn """Loading LazyArtifacts
+#              This will generate a dependency warning until compatibility with Julia 1.4+1.5 is removed"""
+#     using LazyArtifacts
+# end
+
 using LinearAlgebra: BlasFloat, BlasReal, HermOrSym, PosDefException, copytri!
 using Base: Ryu
 using GLM: Link, canonicallink

src/blockdescription.jl

@@ -25,11 +25,11 @@ function BlockDescription(m::LinearMixedModel)
     ALtypes = fill("", k, k)
     Ltypes = fill(Nothing, k, k)
     for i in 1:k, j in 1:i
-        ALtypes[i, j] = shorttype(A[Block(i,j)], L[Block(i,j)])
+        ALtypes[i, j] = shorttype(getblock(A, i,j), getblock(L, i,j))
     end
     BlockDescription(
         blknms,
-        [size(A[Block(i, 1)], 1) for i in 1:k],
+        [size(getblock(A, i, 1), 1) for i in 1:k],
         ALtypes,
      )
 end