[R-sig-ME] Error en mer_finalize(ans) : Downdated X'X is not positive definite, 1. What is wrong with my model?
PALACIO BLASCO, SARA
s.palacio at ipe.csic.es
Tue Mar 5 12:05:32 CET 2013
Below you can see the contingency table for Treatment*Bud_type with
the observations in each combination. As you can see the design is not
balanced but there are no combinations with particularly few
observations. I have checked this table per species and (taking into
account that each species only has one Bud_type), all treatment levels
have a representative number of observations...
Treatment
Bud_type-80 -34.7 -30.2 -26.3 -22.1 -17.2 -13.1 -6.6 4 Total
hy 42 45 47 47 54 71 79 64 39 488
na 24 31 30 34 30 30 30 30 24 263
sc 62 71 72 68 71 82 69 78 62 635
Total 128 147 149 149 155 183 178 172 125 1386
The only issue I can see is that of each species only having one Bud_type...
Thanks for your help!
sara palacio
Quoting Emmanuel Curis <emmanuel.curis at parisdescartes.fr>:
> Then my hint is that you are missing some combinations with for
> instance some species not having some values for the treatment or
> something like that. If this is the cause, the model without treatment
> or without the treatment:fBud_type interaction may not have the
> message. And to adjust the treatment*fBud_type part of the model, some
> coefficients are « stolen » from the fBud_type:Species term, hence the
> NA appearing on last species (and when removing these species,
> appearing on other species).
>
> To check this, you can also build contingency tables
> with treatment and fBud_type, and see if some combinations have only a
> few values, less than the species in the corresponding fBud_type for
> instance... Or 3-D contingency tables, but more difficult to read...
>
> Best regards,
>
> On Tue, Mar 05, 2013 at 09:47:50AM +0100, PALACIO BLASCO, SARA wrote:
> « You are right! This is weird since when I check the data table
> « "species" I can see values of the variable Bud_type for both
> « species... Vu has "sc" and Vv has "hy"...
> «
> « Then, if I follow your suggestion and try to run the model in glmer
> « without these two species it still gives the same error:
> «
> « Error en mer_finalize(ans) : Downdated X'X is not positive definite, 1.
> «
> « If I then run the model in glm to see where the NAs are, I get this
> « output, where, surprisingly, the species "En" that in the previous
> « run had But_type=hy, now has NAs for all the levels of Bud_type!:
> «
> « Call:
> « glm(formula = Dead ~ Treatment * fBud_type + fBud_type:Species,
> « family = binomial, data = species)
> «
> « Deviance Residuals:
> « Min 1Q Median 3Q Max
> « -5.4628 -0.1953 0.0599 0.3561 2.6043
> «
> « Coefficients: (14 not defined because of singularities)
> « Estimate Std. Error z value Pr(>|z|)
> « (Intercept) -5.89770 0.81116 -7.271 3.58e-13 ***
> « Treatment -0.38498 0.04952 -7.774 7.63e-15 ***
> « fBud_typena -1.68437 1.36330 -1.236 0.21664
> « fBud_typesc 2.61123 0.94512 2.763 0.00573 **
> « Treatment:fBud_typena 0.01934 0.07001 0.276 0.78234
> « Treatment:fBud_typesc 0.15635 0.05528 2.828 0.00468 **
> « fBud_typehy:SpeciesEc 1.08250 0.38154 2.837 0.00455 **
> « fBud_typena:SpeciesEc NA NA NA NA
> « fBud_typesc:SpeciesEc NA NA NA NA
> « fBud_typehy:SpeciesEn NA NA NA NA
> « fBud_typena:SpeciesEn NA NA NA NA
> « fBud_typesc:SpeciesEn NA NA NA NA
> « fBud_typehy:SpeciesLp NA NA NA NA
> « fBud_typena:SpeciesLp 1.21835 0.49763 2.448 0.01435 *
> « fBud_typesc:SpeciesLp NA NA NA NA
> « fBud_typehy:SpeciesRf NA NA NA NA
> « fBud_typena:SpeciesRf NA NA NA NA
> « fBud_typesc:SpeciesRf 1.18017 0.44178 2.671 0.00755 **
> « fBud_typehy:SpeciesRh NA NA NA NA
> « fBud_typena:SpeciesRh NA NA NA NA
> « fBud_typesc:SpeciesRh -0.08258 0.40105 -0.206 0.83686
> « fBud_typehy:SpeciesVm NA NA NA NA
> « fBud_typena:SpeciesVm NA NA NA NA
> « fBud_typesc:SpeciesVm NA NA NA NA
> « ---
> « Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> «
> « (Dispersion parameter for binomial family taken to be 1)
> «
> « Null deviance: 1313.62 on 1015 degrees of freedom
> « Residual deviance: 522.38 on 1006 degrees of freedom
> « AIC: 542.38
> «
> « Number of Fisher Scoring iterations: 8
> «
> « I really don't know what is going on, but thanks heaps for your help!!
> «
> « Sara Palacio
> «
> « Quoting Emmanuel Curis <emmanuel.curis at parisdescartes.fr>:
> «
> « >Just looking quickly, is seems strange that for the last two species
> « >(Vu and Vv), _all_ coefficients are to NA? If you try without these
> « >two species, does it work better?
> « >
> « >After that, I'm not specialist enough to traceback why for these two
> « >species there are only NAs --- may be only one observation only for
> « >each of them? or associated to another fBud_type not used in the
> « >analysis for some reason? --- if unused levels have been removed...
> « >
> « >Hope this helps,
> « >
> « >Best regards,
> « >
> « >On Tue, Mar 05, 2013 at 08:29:45AM +0100, PALACIO BLASCO, SARA wrote:
> « >« Dear Ben
> « >«
> « >« This is what the summary(M_bud_type0) says. As expected, there are
> « >« plenty of NAs in the interactions between (uncrossed) levels of the
> « >« interaction between the nested factors (fBud_type:Species):
> « >«
> « >«
> « >« Call:
> « >« glm(formula = Dead ~ Treatment * fBud_type + fBud_type:Species,
> « >« family = binomial, data = species)
> « >«
> « >« Deviance Residuals:
> « >« Min 1Q Median 3Q Max
> « >« -5.8281 -0.2220 0.0703 0.3323 2.3882
> « >«
> « >« Coefficients: (18 not defined because of singularities)
> « >« Estimate Std. Error z value Pr(>|z|)
> « >« (Intercept) -7.77657 0.85126 -9.135 < 2e-16 ***
> « >« Treatment -0.31190 0.03200 -9.747 < 2e-16 ***
> « >« fBud_typena 0.19449 1.38718 0.140 0.88850
> « >« fBud_typesc 5.36751 0.91869 5.843 5.14e-09 ***
> « >« Treatment:fBud_typena -0.05374 0.05892 -0.912 0.36172
> « >« Treatment:fBud_typesc 0.06949 0.03837 1.811 0.07012 .
> « >« fBud_typehy:SpeciesEc 3.96261 0.52793 7.506 6.10e-14 ***
> « >« fBud_typena:SpeciesEc NA NA NA NA
> « >« fBud_typesc:SpeciesEc NA NA NA NA
> « >« fBud_typehy:SpeciesEn 3.01308 0.48926 6.158 7.35e-10 ***
> « >« fBud_typena:SpeciesEn NA NA NA NA
> « >« fBud_typesc:SpeciesEn NA NA NA NA
> « >« fBud_typehy:SpeciesLp NA NA NA NA
> « >« fBud_typena:SpeciesLp 1.21835 0.49753 2.449 0.01433 *
> « >« fBud_typesc:SpeciesLp NA NA NA NA
> « >« fBud_typehy:SpeciesRf NA NA NA NA
> « >« fBud_typena:SpeciesRf NA NA NA NA
> « >« fBud_typesc:SpeciesRf 0.14214 0.39921 0.356 0.72180
> « >« fBud_typehy:SpeciesRh NA NA NA NA
> « >« fBud_typena:SpeciesRh NA NA NA NA
> « >« fBud_typesc:SpeciesRh -1.18370 0.37535 -3.154 0.00161 **
> « >« fBud_typehy:SpeciesVm NA NA NA NA
> « >« fBud_typena:SpeciesVm NA NA NA NA
> « >« fBud_typesc:SpeciesVm -1.09756 0.37513 -2.926 0.00344 **
> « >« fBud_typehy:SpeciesVu NA NA NA NA
> « >« fBud_typena:SpeciesVu NA NA NA NA
> « >« fBud_typesc:SpeciesVu NA NA NA NA
> « >« fBud_typehy:SpeciesVv NA NA NA NA
> « >« fBud_typena:SpeciesVv NA NA NA NA
> « >« fBud_typesc:SpeciesVv NA NA NA NA
> « >« ---
> « >« Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> « >«
> « >« (Dispersion parameter for binomial family taken to be 1)
> « >«
> « >« Null deviance: 1797.06 on 1385 degrees of freedom
> « >« Residual deviance: 736.43 on 1374 degrees of freedom
> « >« AIC: 760.43
> « >«
> « >« Number of Fisher Scoring iterations: 7
> « >«
> « >«
> « >« ### If I try your second suggestion and run the model in glm, the
> « >« number of NAs goes down, but there are still a few:
> « >«
> « >« Call:
> « >« glm(formula = Dead ~ Treatment * fBud_type + budspecies, family
> = binomial,
> « >« data = species)
> « >«
> « >« Deviance Residuals:
> « >« Min 1Q Median 3Q Max
> « >« -5.8281 -0.2220 0.0703 0.3323 2.3882
> « >«
> « >« Coefficients: (2 not defined because of singularities)
> « >« Estimate Std. Error z value Pr(>|z|)
> « >« (Intercept) -7.77657 0.85126 -9.135 < 2e-16 ***
> « >« Treatment -0.31190 0.03200 -9.747 < 2e-16 ***
> « >« fBud_typena 0.19449 1.38718 0.140 0.88850
> « >« fBud_typesc 5.36751 0.91869 5.843 5.14e-09 ***
> « >« budspecieshy.Ec 3.96261 0.52793 7.506 6.10e-14 ***
> « >« budspecieshy.En 3.01308 0.48926 6.158 7.35e-10 ***
> « >« budspeciesna.Lp 1.21835 0.49753 2.449 0.01433 *
> « >« budspeciessc.Rf 0.14214 0.39921 0.356 0.72180
> « >« budspeciessc.Rh -1.18370 0.37535 -3.154 0.00161 **
> « >« budspeciessc.Vm -1.09756 0.37513 -2.926 0.00344 **
> « >« budspeciessc.Vu NA NA NA NA
> « >« budspecieshy.Vv NA NA NA NA
> « >« Treatment:fBud_typena -0.05374 0.05892 -0.912 0.36172
> « >« Treatment:fBud_typesc 0.06949 0.03837 1.811 0.07012 .
> « >« ---
> « >« Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> « >«
> « >« (Dispersion parameter for binomial family taken to be 1)
> « >«
> « >« Null deviance: 1797.06 on 1385 degrees of freedom
> « >« Residual deviance: 736.43 on 1374 degrees of freedom
> « >« AIC: 760.43
> « >«
> « >« Number of Fisher Scoring iterations: 7
> « >«
> « >« I also don't know how to include the new factor with droplevels in
> « >« the glmer model... should this new factor replace the nested one?
> « >«
> « >« Cheers,
> « >«
> « >« Sara
> « >«
> « >«
> « >«
> « >«
> « >« Quoting Ben Bolker <bbolker at gmail.com>:
> « >«
> « >«
> « >« >
> « >« > Did you try to fit
> « >« >
> « >« >M_bud_type0 = glm(Dead~Treatment* fBud_type +
> « >« > fBud_type:Species, family=binomial, data=species)
> « >« >
> « >« >as suggested in the FAQ to see where the rank-deficiencies are
> « >« >(i.e. are there NA-valued coefficients?)
> « >« >
> « >« > It's not immediately obvious to me that the fBud_type:Species
> « >« >interaction should be causing trouble, because lme4 internally
> « >« >drops unused levels of factors. You could *try*
> « >« >
> « >« >species$budspecies <- with(species,
> « >« > droplevels(interaction(fBud_type,Species)))
> « >« >
> « >« >just to check that, but I don't think it will help.
> « >« >
> « >« > Using Species as a random effect does *not* mean you "will
> not be able
> « >« >to know its effect" -- you just won't be able to test hypotheses about
> « >« >differences between particular species/combinations of species.
> « >« >You can still use ranef() to get a value (technically not an
> "estimate")
> « >« >for the conditional mode of each species.
> « >« >
> « >« >
> « >« >>
> « >« >>Quoting Ben Bolker <bbolker at gmail.com>:
> « >« >>
> « >« >>>PALACIO BLASCO, SARA <s.palacio at ...> writes:
> « >« >>>
> « >« >>>[snip]
> « >« >>>
> « >« >>>>I am trying to run the following model in glmer:
> « >« >>>>
> « >« >>>>> M_bud_type1=glmer(Dead~Treatment* fBud_type + fBud_type:Species +
> « >« >>>>> (1|fRep), family=binomial, data=species)
> « >« >>>>
> « >« >>>>where:
> « >« >>>>- Dead is a binomial response variable
> « >« >>>>- fBud_type is a fixed factor with 3 levels
> « >« >>>>- Species is a fixed factor with 9 levels nested within
> fBud_type and
> « >« >>>>- fRep is a random factor with 27 levels nested within Species
> « >« >>>>
> « >« >>>>I have 1386 observations.
> « >« >>>>The error message I receive reads:
> « >« >>>>
> « >« >>>>Error en mer_finalize(ans) : Downdated X'X is not positive
> definite, 1.
> « >« >>>>
> « >« >>>
> « >« >>> Did you already read the
> http://glmm.wikidot.com/faq#errors section?
> « >« >>>
> « >« >>> It sounds like all your predictors are categorical
> (although we don't
> « >« >>>know about Treatment), so centering isn't really as important/as
> « >« >>>practical
> « >« >>>an option (you can use sum-to-zero contrasts, but it probably won't
> « >« >>>make a big difference).
> « >« >>>
> « >« >>> Ben Bolker
> « >« >>>
> « >« >>>_______________________________________________
> « >« >>>R-sig-mixed-models at r-project.org mailing list
> « >« >>>https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> « >« >>
> « >« >>
> « >« >>
> « >«
> « >« _______________________________________________
> « >« R-sig-mixed-models at r-project.org mailing list
> « >« https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> « >
> « >--
> « > Emmanuel CURIS
> « > emmanuel.curis at parisdescartes.fr
> « >
> « >Page WWW: http://emmanuel.curis.online.fr/index.html
> «
> «
>
> --
> Emmanuel CURIS
> emmanuel.curis at parisdescartes.fr
>
> Page WWW: http://emmanuel.curis.online.fr/index.html
More information about the R-sig-mixed-models
mailing list