deprecate lmm and glmm, update docs

benchmark/benchmarks.jl

@@ -72,7 +72,7 @@ const mods = Dict{Symbol,Vector{Expr}}(
     :star => []                   # not sure it is worthwhile working with these data
     );
 
-fitbobyqa(rhs::Expr, dsname::Symbol) = fit!(lmm(Formula(:Y, rhs), dat[dsname]))
+fitbobyqa(rhs::Expr, dsname::Symbol) = fit(LinearMixedModel, Formula(:Y, rhs), dat[dsname])
 compactstr(ds,rhs) = replace(string(ds, ':', rhs), ' ', "")
 
 SUITE["simplescalar"] = BenchmarkGroup(["single", "simple", "scalar"])

docs/jmd/MultipleTerms.jmd

@@ -7,7 +7,7 @@ In this chapter we consider models with multiple simple, scalar random-effects t
 
 ```{julia;term=true}
 using DataFrames, Distributions, FreqTables, Gadfly, MixedModels, RData
-using Gadfly.Geom: density, histogram, point
+using Gadfly.Geom: density, histogram, line, point
 using Gadfly.Guide: xlabel, ylabel
 const dat = convert(Dict{Symbol,DataFrame}, load(Pkg.dir("MixedModels", "test", "dat.rda")));
 const ppt250 = inv(500) : inv(250) : 1.;
@@ -81,7 +81,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;term=true}
-penm = fit!(lmm(@formula(Y ~ 1 + (1|G) + (1|H)), dat[:Penicillin]))
+penm = fit(LinearMixedModel, @formula(Y ~ 1 + (1|G) + (1|H)), dat[:Penicillin])
 ```
 
 This model display indicates that the sample-to-sample variability has the greatest contribution, then plate-to-plate variability and finally the “residual” variability that cannot be attributed to either the sample or the plate. These conclusions are consistent with what we see in the data plot (Fig. [fig:Penicillindot]).
@@ -406,7 +406,7 @@ it is sufficiently diffuse to warrant treating it as if it were a continuous res
 At this point we will fit models that have random effects for student, instructor, and department (or the combination) to these data. In the next chapter we will fit models incorporating fixed-effects for instructor and department to these data.
 
 ```{julia;term=true}
-@time instm = fit!(lmm(@formula(Y ~ 1 + A + (1|G) + (1|H) + (1|I)), dat[:InstEval]))
+@time instm = fit(LinearMixedModel, @formula(Y ~ 1 + A + (1|G) + (1|H) + (1|I)), dat[:InstEval])
 ```
 
 (Fitting this complex model to a moderately large data set takes a few seconds on a modest desktop computer. Although this is more time than required for earlier model fits, it is a remarkably short time for fitting a model of this size and complexity. In some ways it is remarkable that such a model can be fit at all on such a computer.)

docs/jmd/SimpleLMM.jmd

@@ -69,9 +69,9 @@ The data are described in Davies (), the fourth edition of the book mentioned ab
 
 First attach the packages to be used
 ```{julia;term=true}
-using DataFrames, Distributions, Gadfly, GLM, MixedModels, RData
+using DataFrames, Distributions, Gadfly, GLM, MixedModels, RData, RCall
 ```
-and allow for unqualified names for some graphics functions.
+and allow for unqualified names for some graphics functions
 ```{julia;term=true}
 using Gadfly.Geom: point, line, histogram, density, vline
 using Gadfly.Guide: xlabel, ylabel, yticks
@@ -84,11 +84,11 @@ Access the `Dyestuff` data
 ```{julia;term=true}
 const dat = convert(Dict{Symbol,DataFrame}, load(Pkg.dir("MixedModels", "test", "dat.rda")));
 dyestuff = dat[:Dyestuff];
-dump(dyestuff)
+describe(dyestuff)
 ```
 and plot it
 ```{julia;echo=false;fig_cap="Yield versus Batch for the Dyestuff data"; fig_width=8;}
-plot(dyestuff, x = "Y", y = "G", point, xlabel("Yield of dyestuff (g)"), ylabel("Batch"))
+plot(dyestuff, x = :Y, y = :G, point, xlabel("Yield of dyestuff (g)"), ylabel("Batch"))
 ```
 
 In the dotplot we can see that there is considerable variability in yield, even for preparations from the same batch, but there is also noticeable batch-to-batch variability.
@@ -99,11 +99,9 @@ In a plot, however, the order of the levels influences the perception of the pat
 Rather than providing an arbitrary pattern it is best to order the levels according to some criterion for the plot.
 In this case a good choice is to order the batches by increasing mean yield, which can be easily done in R.
 
-(Note: at present this plot fails because of the ongoing DataFrames conversion.)
-
 ```{julia;term=true}
-#dyestuff = rcopy("within(Dyestuff, Batch <- reorder(Batch, Yield, mean))");
-#plot(dyestuff, x="Y", y="G", point, xlabel("Yield of dyestuff (g)"))
+dyestuffR = rcopy(R"within(lme4::Dyestuff, Batch <- reorder(Batch, Yield, mean))");
+plot(dyestuffR, x = :Yield, y = :Batch, point, xlabel("Yield of dyestuff (g)"), ylabel("Batch"))
 ```
 
 In Sect. [sec:DyestuffLMM] we will use mixed models to quantify the variability in yield between batches.
@@ -122,11 +120,11 @@ The data are simulated data presented in Box and Tiao (1973), where the authors
 The structure and summary are intentionally similar to those of the `Dyestuff` data.
 ```{julia;term=true}
 dyestuff2 = dat[:Dyestuff2];
-dump(dyestuff2)
+describe(dyestuff2)
 ```
 As can be seen in a data plot
 ```{julia;echo=false;fig_width=8}
-plot(dyestuff2, x = "Y", y = "G", point, xlabel("Simulated response"), ylabel(""))
+plot(dyestuff2, x = :Y, y = :G, point, xlabel("Simulated response"), ylabel(""))
 ```
 the batch-to-batch variability in these data is small compared to the within-batch variability.
 In some approaches to mixed models it can be difficult to fit models to such data.
@@ -144,7 +142,7 @@ The structure of the formula will be explained after showing the example.
 
 A model allowing for an overall level of the `Yield` and for an additive random effect for each level of `Batch` can be fit as
 ```{julia;term=true}
-mm1 = fit!(lmm(@formula(Y ~ 1 + (1 | G)), dyestuff))
+mm1 = fit(LinearMixedModel, @formula(Y ~ 1 + (1 | G)), dyestuff)
 ```
 
 As shown in the summary of the model fit, the default estimation criterion is maximum likelihood.
@@ -187,7 +185,7 @@ The standard error of the intercept estimate is 17.69 g.
 Fitting a similar model to the `dyestuff2` data produces an estimate $\widehat{\sigma_1^2}=0$.
 
 ```{julia;term=true}
-mm2 = fit!(lmm(@formula(Y ~ 1 + (1 | G)), dyestuff2))
+mm2 = fit(LinearMixedModel, @formula(Y ~ 1 + (1 | G)), dyestuff2)
 ```
 
 An estimate of `0` for $\sigma_1$ does not mean that there is no variation between the groups.
@@ -335,7 +333,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!(lmm(@formula(Y ~ 1 + (1 | G)), dyestuff), true);
+mm1 = fit!(LinearMixedModel(@formula(Y ~ 1 + (1 | G)), dyestuff), 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$.
@@ -349,7 +347,7 @@ mm1.optsum
 The full list of fields in a `LinearMixedModel` object is
 
 ```{julia;term=true}
-fieldnames(LinearMixedModel)
+showcompact(fieldnames(LinearMixedModel))
 ```
 
 The `formula` field is a copy of the model formula
@@ -520,12 +518,12 @@ First set the random number seed for reproducibility.
 
 ```{julia;term=true}
 srand(1234321);
-mm1bstp = bootstrap(10000, mm1);
+mm1bstp = bootstrap(100000, mm1);
 size(mm1bstp)
 ```
 
 ```{julia;term=true}
-show(names(mm1bstp))
+showcompact(names(mm1bstp))
 ```
 
 #### Histograms, kernel density plots and quantile-quantile plots
@@ -539,15 +537,15 @@ Finally, the extent to which the distribution of a sample can be approximated by
 The [`Gadfly`](https://github.com/GiovineItalia/Gadfly.jl) package for Julia uses a "grammar of graphics" specification, similar to the [`ggplot2`](http://ggplot2.org/) package for R.  A histogram or a kernel density plot are describes as *geometries* and specified by `Geom.histogram` and `Geom.density`, respectively.
 
 ```{julia;term=true}
-plot(mm1bstp, x = :β₁, Geom.histogram)
+plot(mm1bstp, x = :β₁, histogram)
 ```
 
 ```{julia;term=true}
-plot(mm1bstp, x = :σ, Geom.histogram)
+plot(mm1bstp, x = :σ, histogram)
 ```
 
 ```{julia;term=true}
-plot(mm1bstp, x = :σ₁, Geom.histogram)
+plot(mm1bstp, x = :σ₁, histogram)
 ```
 
 The last two histograms show that, even if the models are defined in terms of variances, the variance is usually not a good scale on which to assess the variability of the parameter estimates.  The standard deviation or, in some cases, the logarithm of the standard deviation is a more suitable scale.
@@ -561,19 +559,19 @@ length(mm1bstp[:θ₁]) - countnz(mm1bstp[:θ₁])
 That is, nearly 1/10 of the `theta1` values are zeros.  Because such a spike or pulse will be spread out or diffused in a kernel density plot,
 
 ```{julia;term=true}
-plot(mm1bstp, x = :θ₁, Geom.density)
+plot(mm1bstp, density, x = :θ₁)
 ```
 
 such a plot is not suitable for a sample of a bounded parameter that includes values on the boundary.
 
 The density of the estimates of the other two parameters, $\beta_1$ and $\sigma$, are depicted well in kernel density plots.
 
 ```{julia;term=true}
-plot(mm1bstp, x = :β₁, Geom.density)
+plot(mm1bstp, density, x = :β₁)
 ```
 
 ```{julia;term=true}
-plot(mm1bstp, x = :σ, Geom.density)
+plot(mm1bstp, density, x = :σ)
 ```
 
 The standard approach of summarizing a sample by its mean and standard deviation, or of constructing a confidence interval using the sample mean, the standard error of the mean and quantiles of a *t* or normal distribution, are based on the assumption that the sample is approximately normal (also called Gaussian) in shape.  A *normal probability plot*, which plots sample quantiles versus quantiles of the standard normal distribution, $\mathcal{N}(0,1)$, can be used to assess the validity of this assumption.  If the points fall approximately along a straight line, the assumption of normality should be valid.  Systematic departures from a straight line are cause for concern.
@@ -597,34 +595,34 @@ 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), Geom.line,
-    Guide.xlabel("Standard Normal Quantiles"), Guide.ylabel("β₁"))
+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.
 
 ```{julia;echo=false;fig_width=8}
-plot(x = zquantiles, y = quantile(mm1bstp[:σ], ppt250), Geom.line,
-    Guide.xlabel("Standard Normal quantiles"), Guide.ylabel("σ"))
+plot(x = zquantiles, y = quantile(mm1bstp[:σ], ppt250), line,
+     xlabel("Standard Normal quantiles"), ylabel("σ"))
 ```
 
 The normal probability plot of $\sigma_1$ has a flat section at $\sigma_1 = 0$.
 
 ```{julia;echo=false;fig_width=8}
-plot(x = zquantiles, y = quantile(mm1bstp[:σ₁], ppt250), Geom.line,
-    Guide.xlabel("Standard Normal Quantiles"), Guide.ylabel("σ₁"))
+plot(x = zquantiles, y = quantile(mm1bstp[:σ₁], ppt250), line,
+     xlabel("Standard Normal Quantiles"), ylabel("σ₁"))
 ```
 
 In terms of the variances, $\sigma^2$ and $\sigma_1^2$, the normal probability plots are
 
 ```{julia;echo=false}
-plot(x = zquantiles, y = quantile(abs2.(mm1bstp[:σ]), ppt250), Geom.line,
-    Guide.xlabel("Standard Normal quantiles"), Guide.ylabel("σ²"))
+plot(x = zquantiles, y = quantile(abs2.(mm1bstp[:σ]), ppt250), line,
+     xlabel("Standard Normal quantiles"), ylabel("σ²"))
 ```
 
 ```{julia;echo=false}
-plot(x = zquantiles, y = quantile(abs2.(mm1bstp[:σ₁]), ppt250), Geom.line,
-    Guide.xlabel("Standard Normal Quantiles"), Guide.ylabel("σ₁²"))
+plot(x = zquantiles, y = quantile(abs2.(mm1bstp[:σ₁]), ppt250), line,
+     xlabel("Standard Normal Quantiles"), ylabel("σ₁²"))
 ```
 
 ### Confidence intervals based on bootstrap samples
@@ -741,7 +739,7 @@ The empirical cumulative distribution function (ecdf) of a sample maps the range
 plot(layer(x = quantile(mm1bstp[:σ₁], ppt250), y = ppt250, line),
      layer(xintercept = quantile(mm1bstp[:σ₁], [0.1, 0.9]), vline(color = colorant"orange")),
      layer(xintercept = hpdinterval(mm1bstp[:σ₁], 0.8), vline(color=colorant"red")),
-     ylabel(""), Guide.xlabel("σ₁"), yticks(ticks=[0.0, 0.1, 0.9, 1.0])
+     ylabel(""), xlabel("σ₁"), yticks(ticks=[0.0, 0.1, 0.9, 1.0])
 )
 ```
 

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!(lmm(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), true)
+sleepm = fit!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), 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.
@@ -274,7 +274,7 @@ show(names(oxboys))
 ```
 
 ```{julia;term=true}
-oxboysm = fit!(lmm(@formula(height ~ 1 + age + (1+age | Subject)), oxboys))
+oxboysm = fit(LinearMixedModel, @formula(height ~ 1 + age + (1+age | Subject)), oxboys)
 ```
 
 ```{julia;term=true}
@@ -383,7 +383,7 @@ When the time origin is the beginning of the treatment there is not generally a
 
 ```{julia;term=true}
 early = rcopy(R"subset(Early, select = c(cog, tos, id, trt, trttos))");
-earlym = fit!(lmm(@formula(cog ~ 1 + tos + trttos + (1 + tos | id)), early))
+earlym = fit(LinearMixedModel, @formula(cog ~ 1 + tos + trttos + (1 + tos | id)), early)
 ```
 
 The model converges to a singular covariance matrix for the random effects.

docs/jmd/SubjectItem.jmd

@@ -5,7 +5,7 @@ const dat = convert(Dict{Symbol,DataFrame}, load(Pkg.dir("MixedModels", "test",
 ```
 
 ```{julia;term=true}
-mm1 = fit!(lmm(@formula(Y ~ 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)), dat[:kb07]))
+mm1 = fit(LinearMixedModel, @formula(Y ~ 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)), dat[:kb07])
 ```
 
 ```{julia;term=true}

docs/jmd/bootstrap.jmd

@@ -22,63 +22,36 @@ parameter, `θ`, that defines the variance-covariance matrices of the random eff
 For example, a simple linear mixed-effects model for the `Dyestuff` data in the [`lme4`](http://github.com/lme4/lme4)
 package for [`R`](https://www.r-project.org) is fit by
 ```{julia;term=true}
-using DataFrames, Gadfly, MixedModels, RData
+using DataFrames, MixedModels, RData, Gadfly
 ```
 ```{julia;echo=false;results="hidden"}
 const dat = convert(Dict{Symbol,DataFrame}, load(Pkg.dir("MixedModels", "test", "dat.rda")));
 ```
 ```{julia;term=true}
 ds = names!(dat[:Dyestuff], [:Batch, :Yield])
-m1 = fit!(lmm(@formula(Yield ~ 1 + (1 | Batch)), ds))
+m1 = fit(LinearMixedModel, @formula(Yield ~ 1 + (1 | Batch)), ds)
 ```
 
-
-## Using the `bootstrap!` function
-
-This quick explanation is provided for those who only wish to use the `bootstrap!` method and do not need
-detailed explanations of how it works.
-The three arguments to `bootstrap!` are the matrix that will be overwritten with the results, the model to bootstrap,
-and a function that overwrites a vector with the results of interest from the model.
-
-Suppose the objective is to obtain 100,000 parametric bootstrap samples of the estimates of the "variance
-components", `σ²` and `σ₁²`, in this model.  In many implementations of mixed-effects models the
-estimate of `σ₁²`, the variance of the scalar random effects, is reported along with a
-standard error, as if the estimator could be assumed to have a Gaussian distribution.
-Is this a reasonable assumption?
-
-A suitable function to save the results is
-```{julia;term=true}
-function saveresults!(v, m)
-    v[1] = varest(m)
-    v[2] = abs2(getθ(m)[1]) * v[1]
-end
-```
-The `varest` extractor function returns the estimate of `σ²`.  As seen above, the estimate of the
-`σ₁` is the product of `Θ` and the estimate of `σ`.  The expression `abs2(getΘ(m)[1])` evaluates to
-`Θ²`. The `[1]` is necessary because the value returned by `getθ` is a vector and a scalar is needed
-here.
-
-As with any simulation-based method, it is advisable to set the random number seed before calling
-`bootstrap!` for reproducibility.
+Now bootstrap the model parameters
 ```{julia;term=true;}
-srand(1234321);
-```
-```{julia;term=true;}
-results = bootstrap!(zeros(2, 100000), m1, saveresults!);
+results = bootstrap(100_000, m1);
+showcompact(names(results))
 ```
 The results for each bootstrap replication are stored in the columns of the matrix passed in as the first
-argument.  A density plot of the first row using the [`Gadfly`](https://github.com/dcjones/Gadfly.jl) package
-is created as
+argument.  A density plot of the bootstrapped values of `σ` is created as
 ```{julia;eval=false;term=true}
-plot(x = view(results, 1, :), Geom.density(), Guide.xlabel("Parametric bootstrap estimates of σ²"))
+plot(results, x = :σ, Geom.density, Guide.xlabel("Parametric bootstrap estimates of σ"))
+```
+```{julia;echo=false;fig_cap="Density of parametric bootstrap estimates of σ from model m1"; fig_width=8;}
+plot(results, x = :σ, Geom.density, Guide.xlabel("Parametric bootstrap estimates of σ"))
 ```
-```{julia;echo=false;fig_cap="Density of parametric bootstrap estimates of σ² from model m1"; fig_width=8;}
-plot(x = view(results, 1, :), Geom.density(), Guide.xlabel("Parametric bootstrap estimates of σ²"))
+```{julia;echo=false;fig_cap="Density of parametric bootstrap estimates of σ₁ from model m1"; fig_width=8;}
+plot(results, x = :σ₁, Geom.density, Guide.xlabel("Parametric bootstrap estimates of σ₁"))
 ```
-```{julia;echo=false;fig_cap="Density of parametric bootstrap estimates of σ₁² from model m1"; fig_width=8;}
-plot(x = view(results, 2, :), Geom.density(), Guide.xlabel("Parametric bootstrap estimates of σ₁²"))
+```{julia;echo=false;fig_cap="Histogram of parametric bootstrap estimates of σ₁ from model m1"; fig_width=8;}
+plot(results, x = :σ₁, Geom.histogram, Guide.xlabel("Parametric bootstrap estimates of σ₁"))
 ```
 
-The distribution of the bootstrap samples of `σ²` is a bit skewed but not terribly so.  However, the
-distribution of the bootstrap samples of the estimate of `σ₁²` is highly skewed and has a spike at
+The distribution of the bootstrap samples of `σ` is a bit skewed but not terribly so.  However, the
+distribution of the bootstrap samples of the estimate of `σ₁` is highly skewed and has a spike at
 zero.