Merge pull request #163 from dmbates/REML

add REML estimation criterion

Author: dmbates

.travis.yml

@@ -7,6 +7,7 @@ os:
 julia:
   - 0.7
   - 1.0
+  - 1.1
   - nightly
 matrix:
   allow_failures:

appveyor.yml

@@ -2,6 +2,7 @@ environment:
   matrix:
   - julia_version: 0.7
   - julia_version: 1.0
+  - julia_version: 1.1
   - julia_version: nightly
 
 platform:

docs/jmd/SimpleLMM.jmd

@@ -336,7 +336,7 @@ For a linear mixed model, where all the conditional and unconditional distributi
 The optional second argument, `verbose`, in a call to `fit!` of a `LinearMixedModel` object produces output showing the progress of the iterative optimization of $\tilde{d}(\bf\theta|\bf y)$.
 
 ```{julia;term=true}
-mm1 = fit!(LinearMixedModel(@formula(Y ~ 1 + (1 | G)), dyestuff), true);
+mm1 = fit!(LinearMixedModel(@formula(Y ~ 1 + (1 | G)), dyestuff), verbose = true);
 ```
 
 The algorithm converges after 18 function evaluations to a profiled deviance of 327.32706 at $\theta=0.752581$. In this model the parameter $\theta$ is of length 1, the single element being the ratio $\sigma_1/\sigma$.
@@ -597,9 +597,9 @@ const ppt250 = ppoints(250)
 The kernel density estimate of $\sigma$ is more symmetric
 
 ```{julia;echo=false;fig_width=8}
-zquantiles = quantile(Normal(), ppt250);
-plot(x = zquantiles, y = quantile(mm1bstp[:β₁], ppt250), line)
-#    Guide.xlabel("Standard Normal Quantiles"), Guide.ylabel("β₁"))
+zquantiles = quantile.(Normal(), ppt250);
+plot(x = zquantiles, y = quantile(mm1bstp[:β₁], ppt250), line,
+    Guide.xlabel("Standard Normal Quantiles"), Guide.ylabel("β₁"))
 ```
 
 and the normal probability plot of $\sigma$ is also reasonably straight.

docs/jmd/SingularCovariance.jmd

@@ -25,7 +25,7 @@ As is customary (though not required) in Julia, a function whose name ends in `!
 An optional second argument of `true` in the call to `fit!` produces verbose output from the optimization.
 
 ```{julia;term=true}
-sleepm = fit!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), true)
+sleepm = fit!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), verbose=true)
 ```
 
 The variables in the optimization are the elements of a lower triangular matrix, $\Lambda$, which is the relative covariance factor of the random effects.
@@ -332,7 +332,7 @@ orthfm = fit!(LinearMixedModel(@formula(distance ~ 1 + age + (1 + age | Subject)
 ```
 
 ```{julia;term=true}
-srand(1234123)
+Random.seed!(1234123)
 orthfmbtstrp = bootstrap(10000, orthfm);
 ```
 

docs/jmd/constructors.jmd

@@ -28,14 +28,19 @@ Categorical covariates not suitable as grouping factors are named starting with
 The formula language in *Julia* is similar to that in *R* except that the formula must be enclosed in a call to the `@formula` macro.
 A basic model with simple, scalar random effects for the levels of `G` (the batch of an intermediate product, in this case) is declared and fit as
 ```{julia;term=true}
-fm1 = fit(LinearMixedModel, @formula(Y ~ 1 + (1|G)), dat[:Dyestuff])
+fm1 = fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff]))
 ```
 
 (If you are new to Julia you may find that this first fit takes an unexpectedly long time, due to Just-In-Time (JIT) compilation of the code.
 The second and subsequent calls to such functions are much faster.)
 
 ```{julia;term=true}
-@time fit(LinearMixedModel, @formula(Y ~ 1 + (1|G)), dat[:Dyestuff2]);
+@time fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff2]));
+```
+
+By default, the model fit is by maximum likelihood.  To use the `REML` criterion instead, add the optional named argument `REML = true` to the call to `fit!`
+```{julia;term=true}
+fm1R = fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff]), REML=true)
 ```
 
 ### Simple, scalar random effects

docs/jmd/optimization.jmd

@@ -190,10 +190,10 @@ Note that the first `ScalarFactorReTerm` in `fm4.trms` corresponds to grouping f
 
 ### Progress of the optimization
 
-An optional `Bool` argument of `true` in the call to `fit!` of a `LinearMixedModel` causes printing of the objective and the $\theta$ parameter at each evaluation during the optimization.
+An optional named argument, `verbose=true`, in the call to `fit!` of a `LinearMixedModel` causes printing of the objective and the $\theta$ parameter at each evaluation during the optimization.
 ```{julia;term=true}
-fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff]), true);
-fit!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), true);
+fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff]), verbose=true);
+fit!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), verbose=true);
 ``` 
 
 A shorter summary of the optimization process is always available as an

docs/make.jl

@@ -1,7 +1,6 @@
 using Documenter, MixedModels, StatsBase
 
 makedocs(
-    format = :html,
     sitename = "MixedModels",
     pages = ["index.md",
              "constructors.md",

docs/src/GaussHermite.md

@@ -33,7 +33,7 @@ For a `k`th order normalized Gauss-Hermite rule the tridiagonal matrix has zeros
 julia> using LinearAlgebra, Gadfly
 
 julia> sym3 = SymTridiagonal(zeros(3), sqrt.(1:2))
-3×3 SymTridiagonal{Float64,Array{Float64,1}}:
+3×3 LinearAlgebra.SymTridiagonal{Float64,Array{Float64,1}}:
  0.0  1.0       ⋅     
  1.0  0.0      1.41421
   ⋅   1.41421  0.0    
@@ -73,11 +73,11 @@ julia> gausshermitenorm(3)
 
 The weights and positions are often shown as a *lollipop plot*.
 For the 9th order rule these are
-![Lollipop plot of 9th order normalized Gauss-Hermite rule](./assets//GaussHermite_4_1.svg)
+![Lollipop plot of 9th order normalized Gauss-Hermite rule](./assets/GaussHermite_4_1.svg)
 
 
 Notice that the magnitudes of the weights drop quite dramatically away from zero, even on a logarithmic scale
-![Lollipop plot of 9th order normalized Gauss-Hermite rule (logarithmic scale](./assets//GaussHermite_5_1.svg)
+![Lollipop plot of 9th order normalized Gauss-Hermite rule (logarithmic scale](./assets/GaussHermite_5_1.svg)
 
 
 
@@ -136,7 +136,7 @@ julia> const contra = @transform(dat[:Contraception],
      a2 = abs2.(:a), urbdist = string.(:urb, :d));
 
 julia> describe(contra)
-8×8 DataFrame. Omitted printing of 2 columns
+8×8 DataFrames.DataFrame. Omitted printing of 2 columns
 │ Row │ variable │ mean       │ min    │ median  │ max     │ nunique │
 │     │ Symbol   │ Union…     │ Any    │ Union…  │ Any     │ Union…  │
 ├─────┼──────────┼────────────┼────────┼─────────┼─────────┼─────────┤
@@ -156,7 +156,7 @@ julia> describe(contra)
 
 
 Because a smoothed scatterplot of contraception use versus age
-![Scatterplot smooth of contraception use versus age](./assets//GaussHermite_9_1.svg)
+![Scatterplot smooth of contraception use versus age](./assets/GaussHermite_9_1.svg)
 
 
 shows that the proportion of women using artificial contraception is approximately quadratic in age,
@@ -170,8 +170,8 @@ julia> m1 = fit!(GeneralizedLinearMixedModel(form1, contra,
     Bernoulli()), fast=true)
 Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
   Formula: use ~ 1 + a + a2 + l + urb + (1 | d)
-  Distribution: Bernoulli{Float64}
-  Link: LogitLink()
+  Distribution: Distributions.Bernoulli{Float64}
+  Link: GLM.LogitLink()
 
   Deviance: 2372.7844
 
@@ -232,7 +232,7 @@ This is primarily due to different sample sizes in the different districts.
 julia> using FreqTables
 
 julia> freqtable(contra, :d)'
-1×60 Named Adjoint{Int64,Array{Int64,1}}
+1×60 Named LinearAlgebra.Adjoint{Int64,Array{Int64,1}}
 ' ╲ d │   1    2    3    4    5    6    7  …   55   56   57   58   59   60   61
 ──────┼────────────────────────────────────────────────────────────────────────
 1     │ 117   20    2   30   39   65   18  …    6   45   27   33   10   32   42
@@ -266,12 +266,12 @@ end
 
 
 A plot of the deviance contribution versus $u_1$
-![Deviance contribution of u₁](./assets//GaussHermite_14_1.svg)
+![Deviance contribution of u₁](./assets/GaussHermite_14_1.svg)
 
 
 shows that the deviance contribution is very close to a quadratic.
 This is also true for $u_3$
-![Deviance contribution of u₃](./assets//GaussHermite_15_1.svg)
+![Deviance contribution of u₃](./assets/GaussHermite_15_1.svg)
 
 
 
@@ -292,7 +292,7 @@ As shown below, this is an estimate of the conditional standard deviations of th
 julia> const s = inv.(m1.LMM.L.data[Block(1,1)].diag);
 
 julia> s'
-1×60 Adjoint{Float64,Array{Float64,1}}:
+1×60 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
  0.406889  0.713511  0.952164  0.627135  …  0.839679  0.654965  0.60326
 
 ````
@@ -316,9 +316,9 @@ end
 
 
 
-![Scaled and shifted deviance contributions](./assets//GaussHermite_19_1.svg)
+![Scaled and shifted deviance contributions](./assets/GaussHermite_19_1.svg)
 
-![Scaled and shifted deviance contributions](./assets//GaussHermite_20_1.svg)
+![Scaled and shifted deviance contributions](./assets/GaussHermite_20_1.svg)
 
 
 
@@ -337,9 +337,9 @@ end
 
 
 
-![Scaled and shifted conditional density](./assets//GaussHermite_22_1.svg)
+![Scaled and shifted conditional density](./assets/GaussHermite_22_1.svg)
 
-![Scaled and shifted conditional density](./assets//GaussHermite_23_1.svg)
+![Scaled and shifted conditional density](./assets/GaussHermite_23_1.svg)
 
 
 and the function to be integrated with the normalized Gauss-Hermite rule is
@@ -357,6 +357,6 @@ end
 
 
 
-![Function to be integrated with normalized Gauss-Hermite rule](./assets//GaussHermite_25_1.svg)
+![Function to be integrated with normalized Gauss-Hermite rule](./assets/GaussHermite_25_1.svg)
 
-![Function to be integrated with normalized Gauss-Hermite rule](./assets//GaussHermite_26_1.svg)
+![Function to be integrated with normalized Gauss-Hermite rule](./assets/GaussHermite_26_1.svg)

docs/src/MultipleTerms.md

@@ -19,7 +19,7 @@ julia> const dat = Dict(Symbol(k)=>v for (k,v) in
 
 julia> const ppt250 = inv(500) : inv(250) : 1.;
 
-julia> const zquantiles = quantile(Normal(), ppt250);
+julia> const zquantiles = quantile.(Normal(), ppt250);
 
 julia> function hpdinterval(v, level=0.95)
     n = length(v)
@@ -62,7 +62,7 @@ The data are derived from Table 6.6, p. 144 of Davies (), where they are descr
 
 ````julia
 julia> describe(dat[:Penicillin])
-3×8 DataFrame. Omitted printing of 1 columns
+3×8 DataFrames.DataFrame. Omitted printing of 1 columns
 │ Row │ variable │ mean    │ min  │ median │ max  │ nunique │ nmissing │
 │     │ Symbol   │ Union…  │ Any  │ Union… │ Any  │ Union…  │ Nothing  │
 ├─────┼──────────┼─────────┼──────┼────────┼──────┼─────────┼──────────┤
@@ -122,7 +122,7 @@ Even when we apply each of the six samples to each of the 24 plates, something c
 A model incorporating random effects for both the plate and the sample is straightforward to specify — we include simple, scalar random effects terms for both these factors.
 
 ````julia
-julia> penm = fit(LinearMixedModel, @formula(Y ~ 1 + (1|G) + (1|H)), dat[:Penicillin])
+julia> penm = fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G) + (1|H)), dat[:Penicillin]))
 Linear mixed model fit by maximum likelihood
  Formula: Y ~ 1 + (1 | G) + (1 | H)
    logLik   -2 logLik     AIC        BIC    
@@ -157,7 +157,7 @@ In chapter [chap:ExamLMM] we saw that a model with a single, simple, scalar ran
 The relative covariance factor, $\Lambda_\theta$, (Fig. [fig:fm03LambdaLimage], left panel) is no longer a multiple of the identity. It is now block diagonal, with two blocks, one of size 24 and one of size 6, each of which is a multiple of the identity. The diagonal elements of the two blocks are $\theta_1$ and $\theta_2$, respectively. The numeric values of these parameters can be obtained as
 
 ````julia
-julia> show(getθ(penm))
+julia> show(penm.θ)
 [1.53758, 3.21975]
 ````
 
@@ -201,7 +201,7 @@ A bootstrap simulation of the model
 
 ````julia
 julia> @time penmbstp = bootstrap(10000, penm);
- 24.313179 seconds (83.13 M allocations: 2.532 GiB, 3.95% gc time)
+ 26.932617 seconds (77.19 M allocations: 2.076 GiB, 3.17% gc time)
 
 ````
 
@@ -211,11 +211,11 @@ julia> @time penmbstp = bootstrap(10000, penm);
 
 provides the density plots
 
-![](./assets//MultipleTerms_8_1.svg)
+![](./assets/MultipleTerms_8_1.svg)
 
-![](./assets//MultipleTerms_9_1.svg)
+![](./assets/MultipleTerms_9_1.svg)
 
-![](./assets//MultipleTerms_10_1.svg)
+![](./assets/MultipleTerms_10_1.svg)
 
 ````julia
 julia> plot(penmbstp, x = :σ₂, density, xlabel("σ₂"))
@@ -281,7 +281,7 @@ The structure and summary of the data object are
 
 ````julia
 julia> describe(dat[:Pastes])
-4×8 DataFrame. Omitted printing of 1 columns
+4×8 DataFrames.DataFrame. Omitted printing of 1 columns
 │ Row │ variable │ mean    │ min  │ median │ max  │ nunique │ nmissing │
 │     │ Symbol   │ Union…  │ Any  │ Union… │ Any  │ Union…  │ Nothing  │
 ├─────┼──────────┼─────────┼──────┼────────┼──────┼─────────┼──────────┤
@@ -403,13 +403,13 @@ confirm this impression in that all the prediction intervals for the random effe
 julia> Random.seed!(4321234);
 
 julia> @time pstsbstp = bootstrap(10000, pstsm);
- 17.509696 seconds (67.42 M allocations: 2.024 GiB, 4.36% gc time)
+ 21.205995 seconds (66.25 M allocations: 1.788 GiB, 3.44% gc time)
 
 ````
 
 
 
-![](./assets//MultipleTerms_16_1.svg)
+![](./assets/MultipleTerms_16_1.svg)
 
 ````julia
 julia> plot(pstsbstp, x = :σ, density, xlabel("σ"))
@@ -450,7 +450,7 @@ Plot(...)
 
 ````julia
 julia> count(x -> x < 1.0e-5, pstsbstp[:σ₂])
-3671
+3669
 
 ````
 
@@ -522,7 +522,7 @@ is compared to model `pstsm` with
 
 ````julia
 julia> MixedModels.lrt(pstsm1, pstsm)
-2×4 DataFrame
+2×4 DataFrames.DataFrame
 │ Row │ Df    │ Deviance │ Chisq    │ pval     │
 │     │ Int64 │ Float64  │ Float64  │ Float64  │
 ├─────┼───────┼──────────┼──────────┼──────────┤
@@ -545,7 +545,7 @@ A bootstrap sample
 
 ````julia
 julia> @time psts1bstp = bootstrap(10000, pstsm1);
-  6.902661 seconds (23.19 M allocations: 700.737 MiB, 4.32% gc time)
+  9.347764 seconds (24.71 M allocations: 704.028 MiB, 3.37% gc time)
 
 ````
 
@@ -555,13 +555,13 @@ julia> @time psts1bstp = bootstrap(10000, pstsm1);
 
 provides empirical density plots
 
-![](./assets//MultipleTerms_26_1.svg)
+![](./assets/MultipleTerms_26_1.svg)
 
 
 
 and
 
-![](./assets//MultipleTerms_27_1.svg)
+![](./assets/MultipleTerms_27_1.svg)
 
 
 
@@ -631,7 +631,7 @@ Although the response, `Y`, is on a scale of 1 to 5,
 
 ````julia
 julia> freqtable(dat[:InstEval][:Y])'
-1×5 Named Adjoint{Int64,Array{Int64,1}}
+1×5 Named LinearAlgebra.Adjoint{Int64,Array{Int64,1}}
 ' ╲ Dim1 │     1      2      3      4      5
 ─────────┼──────────────────────────────────
 1        │ 10186  12951  17609  16921  15754
@@ -648,7 +648,7 @@ At this point we will fit models that have random effects for student, instructo
 
 ````julia
 julia> @time instm = fit(LinearMixedModel, @formula(Y ~ 1 + A + (1|G) + (1|H) + (1|I)), dat[:InstEval])
-  3.911847 seconds (2.74 M allocations: 323.059 MiB, 8.79% gc time)
+  6.545960 seconds (680.61 k allocations: 218.683 MiB, 1.01% gc time)
 Linear mixed model fit by maximum likelihood
  Formula: Y ~ 1 + A + (1 | G) + (1 | H) + (1 | I)
      logLik        -2 logLik          AIC             BIC