Adam,

Your FAQ contains rather a lot of bold and unsubstantiated claims. I'll
offer some examples with evidence.

It was well known for many years that at least some SAS procedures could fit
larger-than-memory problems whereas R's algorithms were limited to the size
of memory.  See
https://stat.ethz.ch/pipermail/r-help/2001-November/016869.html for a
typical discussion.  Progress has been made since then with some packages
like sqldf.

You correctly state that total development time should be considered, not
just processor time.  You appear to assume that development time is short,
which is typically not the case.  For example, suppose that you have code
written in nlme and then you decide to change all random factors to fixed
factors.  In PROC MIXED this is a trivial change.  In nlme it requires
changing from lme to lm or gls and then you cross your fingers that all the
extractor functions that you are using to process the results won't also
need re-writing.

There are a wealth of examples for PROC MIXED that can be used as a
reference since MIXED has had a fairly stable API, much more so than R.
Here is an example from Crawley's book "Statistical Computing (S-Plus
version)" page 673.

dat1 <- structure(list(y = c(12.82269837, 14.240145, 11.20959257,
                             11.63947395, 13.58207992, 11.85206633,
                             11.11016166, 10.36968933, 11.68580947,
10.34666665, 8.553623229, 11.25313156, 12.39250812, 9.709655149,
9.987958107, 14.01817061, 15.81948348, 14.62259798, 6.442767072,
8.731639184, 8.450912311, 8.947481868, 9.387362706, 8.810636629,
8.429724632, 12.72193612, 7.775626438, 9.190762911, 11.27483279,
10.14748798), fac = structure(as.integer(c(1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5,
5, 5, 5)), .Label = c("1", "2", "3", "4", "5"), class = "factor")), .Names =
c("y",
"fac"), row.names = c("1", "2", "3", "4", "5", "6", "7", "8",
"9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19",
"20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30"
), class = "data.frame")

library(nlme)
model1 <- lme(y~1, dat1,random=~1)

This code worked fine in S-Plus, but R says:

Error in getGroups.data.frame(dataMix, groups) :
  Invalid formula for groups

So the published example becomes almost useless as a reference.  There have
also been numerous API breakages with lme4 (just Google: lme4 broken), which
is normal for alpha versions, but people have come to use lme4 as though it
were stable, even publishing books based on it! (Google books search for
lme4).

Dealing with numerous examples like this can cause even devout open-source
fanatics a great deal of soul-searching.

There's definitely an opportunity for a FAQ, but it needs to be based on
evidence.

Kevin


>
> --Adam
>
> On Wed, 23 Sep 2009, Douglas Bates wrote:
>
>> I forgot to cc: the list on this reply.
>>
>> ---------- Forwarded message ----------
>> From: Douglas Bates <bates@stat.wisc.edu>
>> Date: Wed, Sep 23, 2009 at 11:29 AM
>> Subject: Re: [R-sig-ME] More naive questions: Speed comparisons? what
>> is a "stack imbalance" in lmer? does lmer center variables?
>> To: Paul Johnson <pauljohn32@gmail.com>
>>
>>
>>
>> On Wed, Sep 23, 2009 at 1:36 AM, Paul Johnson <pauljohn32@gmail.com>
>> wrote:
>>>
>>> Sent this to r-sig-debian by mistake the first time.  Depressing.
>>
>>> 1.  One general question for general discussion:
>>
>>> Is HLM6 faster than lmer? If so, why?
>>
>>> I'm always advocating R to students, but some faculty members are
>>> skeptical.  A colleague compared the commercial HLM6 software to lmer.
>>>  HLM6 seems to fit the model in 1 second, but lmer takes 60 seconds.
>>
>>> If you have HLM6 (I don't), can you tell me if you see similar
>>> differences?
>>
>> Honestly, Paul, I don't see a point in discussing vague generalities
>> and rumors.  If your colleague can make the data (or even simuiated
>> data with a similar structure) available and describe the model being
>> fit then we can see why lmer appears to be so much slower on it.  But
>> do bear in mind that the folks marketing HLM6 have the ability -
>> indeed they have the right - to run lmer on their test cases and to
>> examine every single line of code in lmer so they can determine
>> exactly what it is doing.  They have the right to do that because
>> everyone has the right to do that.  If we want to find out how HLM6
>> performs on a test case we are reduced to your approach of asking if
>> someone has a copy of the software (on one of or perhaps the only
>> operating system under which it runs) and all that can be done is to
>> quote the results.  The internal details of HLM6 are proprietary.
>>
>> I can tell you my goals in developing the lme4 package are:
>>  1. Provide the ability to fit general forms of mixed-effects models.
>>  2. To the best of my ability, provide reliable estimates.
>>  3. Allow for models to be fit to very large data sets.
>>  4. Take advantage of the R environment to allow for straightforward
>> model specification and for very general data manipulation and
>> graphical capabilities.
>>  5. Subject to objectives 1-4 provide the fastest software that I can
>>
>>> My first thought was that LM6 uses PQL by default, and it would be
>>> faster.  However, in the output, HLM6 says:
>>
>>> Method of estimation: restricted maximum likelihood
>>
>>> But that doesn't tell me what quadrature approach they use, does it?
>>
>> To me this is very confusing.  In my world REML is only meaningful for
>> linear mixed models.  (I know, at one time Mary Lindstrom and I wrote
>> about REML estimates for nonlinear mixed models but, as I have said
>> before, I consider that a youthful indiscretion.)  For linear mixed
>> models there is no need to use PQL or to discuss quadrature because
>> the evaluation of the profiled deviance or the profiled REML criterion
>> is a direct calculation.
>>
>> Again, without a more complete description of the data and the model
>> to be fit it is impossible to discuss these issues meaningfully.
>>
>>> Another explanation for the difference in time might be the way HLM6
>>> saves the results of some matrix calculations and re-uses them behind
>>> the scenes.  If every call to lmer is re-calculating some big matrix
>>> results, I suppose that could explain it.
>>
>> But you don't need to speculate about what lmer does.  It is Open
>> Source so you can check for yourself.
>>
>> However, this does bring up another point which is the need to compare
>> apples with apples when you are benchmarking software.  If the data
>> import and model specification stages in HLM6 create the necessary
>> matrix structures for performing the iterative fit then does the time
>> to fit the model consist solely of the optimization and summary
>> stages?  Using an expression like
>>
>> system.time(fm1 <- lmer(...))
>>
>> is assessing the time to take the original data, which could be in a
>> very general form, create all those internal structures and perform
>> the optimization.
>>
>> You should bear in mind that a lot of that construction of the model
>> structures is written in R exactly so that it is capable of fitting
>> very general model specifications.  The code in HLM6 is, I imagine,
>> compiled code, which is possible because it targets a very specific
>> task, and compiled code is always going to be much faster than
>> interpreted code.
>>
>>> There are comparisons from 2006 here
>>
>>>
>>>
http://www.cmm.bristol.ac.uk/learning-training/multilevel-m-software/tables.shtml
>>
>>> that indicate that lme was much slower than HLM, but that doesn't help
>>> me understand *why* there is a difference.
>>
>> I believe that those examples are available for testing with lmer in
>> the mlmRev package for R.  However, the results can't be compared
>> meaningfully to the tabled results because today's hardware is so
>> different from what was used for those comparisons.
>>
>> Besides, none of those examples is really challenging.  I was just
>> doing some timing tests on a model fit to 1.7 million observations
>> with 750,000 random effects associated with four different, non-nested
>> grouping factors and 40 fixed-effects parameters.  To me that is a
>> meaningful test because it took from 1 to 2 hours to fit on a fast
>> server computer.  Unfortunately, I can't make that data available
>> because of confidentiality restrictions but, even if I could, I don't
>> think the model could be fit in HLM6 or MLWin or SAS PROC MIXED or
>> Stata because there are 750,000 random effects in a non-nested
>> configuration.
>>
>> I can make a similar but smaller test case available using the star
>> (Student-Teacher Achievement Ratio) data from the mlmRev package.  One
>> model could be
>>
>>> system.time(fm1 <- lmer(math ~ gr + cltype + sx + eth + (1|id) + (1|tch)
>>> + (1|sch), star))
>>
>>   user  system elapsed
>>  9.885   0.044  10.508
>>>
>>> print(fm1, corr = FALSE)
>>
>> Linear mixed model fit by REML
>> Formula: math ~ gr + cltype + sx + eth + (1 | id) + (1 | tch) + (1 | sch)
>>   Data: star
>>    AIC    BIC  logLik deviance REMLdev
>>  239236 239365 -119602   239245  239204
>> Random effects:
>>  Groups   Name        Variance Std.Dev.
>>  id       (Intercept) 1001.60  31.648
>>  tch      (Intercept)  295.03  17.176
>>  sch      (Intercept)  104.78  10.236
>>  Residual              397.32  19.933
>> Number of obs: 24578, groups: id, 10732; tch, 1374; sch, 80
>>
>> Fixed effects:
>>            Estimate Std. Error t value
>> (Intercept) 560.8868     1.5745   356.2
>> gr.L         95.6126     1.0080    94.9
>> gr.Q         -4.6783     0.9857    -4.7
>> gr.C         -3.2320     0.9729    -3.3
>> cltypereg    -7.7943     1.3322    -5.9
>> cltypereg+A  -7.0055     1.3265    -5.3
>> sxF           2.8766     0.6836     4.2
>> ethB        -21.9197     1.2219   -17.9
>> ethA          3.1940     6.7037     0.5
>> ethH          4.5411     9.6653     0.5
>> ethI        -28.1618    13.5359    -2.1
>> ethO          3.4356     8.3133     0.4
>>
>> Unfortunately, because that model is based on only 25000 observations
>> and 11000 random effects it is difficult to get a meaningful timing.
>> Timings under 10 seconds, like this, are subject to too much
>> variability.
>>
>> Nevertheless you can take that data and try to fit that model in any
>> of the other software systems.  Note that student (the "id" factor) is
>> not nested within teacher ("tch") or school, as will almost always
>> happen when fitting longitudinal data so the random effects are
>> partially crossed.  In other words, they are not nested so the
>> structure is not hierarchical or multi-level in the sense of HLM or
>> MLWin.  I would invite those who have access to such software to fit
>> the model using lme4 and to fit an equivalent model in other systems
>> and tell us how they compare.  I can't do that because I choose not to
>> use proprietary software.
>>
>>> 2. What does "stack imbalance in .Call" mean in lmer?
>>
>> It is a sign of a bug in the C code underlying the lme4 package.  It
>> should be reported as a bug if it occurred in the currently released
>> version of lme4.
>>
>>> Here's why I ask.  Searching for comparisons of lmer and HLM,  I went
>>> to CRAN &  I checked this document:
>>>
>>> http://cran.r-project.org/web/packages/mlmRev/vignettes/MlmSoftRev.pdf
>>>
>>> I *think* these things are automatically generated.  The version
>>> that's up there at this moment  (mlmRev edition 0.99875-1)  has pages
>>> full of the error message:
>>>
>>> stack imbalance in .Call,
>>>
>>> Were those always there?  I don't think so.   What do they mean?
>>
>> No they shouldn't be there.  I'll update the mlmRev package.
>>
>>> 3. In the HLM6 output, there is a message at the end of the variable
>>> list:
>>>
>>> '%' - This level-1 predictor has been centered around its grand mean.
>>> '$' - This level-2 predictor has been centered around its grand mean.
>>>
>>> What effect does that have on the estimates?  I believe it should have
>>> no effect on the fixed effect slope estimates, but it seems to me the
>>> estimates of the variances of random parameters would be
>>> changed.  In order to make the estimates from lmer as directly
>>> comparable as possible, should I manually center all of the variables
>>> before fitting the model?   I'm a little stumped on how to center a
>>> multi-category factor before feeding it to lmer.  Know what I mean?
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>

	[[alternative HTML version deleted]]


