No data for 1 interaction combination: problem in R not in Genstat
On Tue, Apr 12, 2011 at 8:47 AM, Ben Bolker <bbolker at gmail.com> wrote:
Roger Humphry <roger.w.humphry at ...> writes:
I've found that when I do a linear mixed model either with lme or lmer with interactions amongst some fixed effects that if there are any interaction levels unoccupied with data that the model fails. "Error in MEEM(object, conLin, control$niterEM) : ? ?Singularity in backsolve at level 0, block 1" for lme and "Error in mer_finalize(ans) : Downdated X'X is not positive definite, 6." in lmer. Strangely in Genstat the model does produce output. My solutions will be either to use Genstat or to create a single new factor variable which has the levels of the interaction that *are* represented in the interaction. I haven't found anything about this when searching. I guess that this may be deliberate (e.g. Genstat uses a fudge considered inappropriate) but please could anybody advise me? I can provide a simple fictitious example that I invented in which the interaction would be of interest but can't be modelled. (I haven't posted it here because I'm unsure about the etiquette of posting large e-mails containing data onto the discussion list).
?I don't know that the fudge is inappropriate, it may just not have been implemented. ?For example, lm() in base R is sensible about setting empty factor levels to NA (if Genstat gives you an *estimate* for the interaction term with the missing levels, that would be weird).
The way that lm and glm handle the missing cells in a model with interactions is through pivoting on columns of the model matrix, but in a way that numerical analysts didn't anticipate and hence do not provide for in standard numerical linear algebra software. Most linear algebra operations in R are done using Lapack code except for the most important one - solving least squares problems - which uses modified Linpack code. (This, by the way, is why accelerated BLAS don't improve calculation speed for lm, glm, etc. The Linpack code only uses level-1 BLAS and there is not much you can do with clever programming to make level-1 BLAS run faster.) Most of the time rank deficiency in least squares problems is avoided by using contrasts for terms involving factors. However, missing cells in a model with an interaction term are one case where rank deficiency cannot (easily) be caught during the symbolic analysis. The methods used in lm perform a QR decomposition with the "pivot on apparent singularity" scheme that produces a set of columnst that are full rank. The methods used in lmer don't include a check for rank deficiency. Martin and I discussed the possibility of doing that immediately after the model matrix for the fixed effects was generated and performing the pivoting and rank reduction at that time but that modification remains on the To Do list. The fact that lm inserts NA for the estimates of the coefficients doesn't really relate to the computational method. That is done, by reversing the column permutation, after the other parameter estimates have been calculated.
?One way of posting compact examples is to make them from synthesized data, for example: d <- expand.grid(f1=LETTERS[1:3],f2=letters[1:3],rep=1:5) d2 <- subset(d, !(f1=="A" & f2=="a")) with(d2,table(f1,f2)) ? f2 f1 ?a b c ?A 0 5 5 ?B 5 5 5 ?C 5 5 5 d3 <- data.frame(d2,y=runif(nrow(d2))) lm(y~f1*f2,data=d3) Call: lm(formula = y ~ f1 * f2, data = d3) Coefficients: (Intercept) ? ? ? ? ?f1B ? ? ? ? ?f1C ? ? ? ? ?f2b ? ? ? ? ?f2c ? ? ?f1B:f2b ? ?0.60964 ? ? -0.32352 ? ? -0.15239 ? ? ?0.13997 ? ? -0.11326 ? ? ?0.10615 ? ?f1C:f2b ? ? ?f1B:f2c ? ? ?f1C:f2c ? -0.03831 ? ? ?0.47481 ? ? ? ? ? NA ?showing in this case that lm() does a sensible thing, setting one confounded parameter to NA. ?If you have long examples you can also try posting the data somewhere and supplying a URL. ?Your solutions sound sensible.
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models