ANOVA error in comparing two models

[John Benning]

Hello all!

I'm trying to test my block effect with the following code:

####################################################################

## Full Aster model

####################################################################

# full fixed aster model

t1a.aster.full <- aster(resp ~ varb + 

                         total_seeds:(site*Seed*Caged) + 

                         total_seeds:site/block,

                         pred1a, fam1a, varb, id, root, data = transplant_yr1a_aster)

summary(t1a.aster.full, show.graph = T)

####################################################################

## test block

####################################################################

t1a.aster.full_noBlock <- aster(resp ~ varb + 

                          total_seeds:(site*Seed*Caged),

                          pred1a, fam1a, varb, id, root, data = transplant_yr1a_aster)

anova(t1a.aster.full_noBlock, t1a.aster.full)

But receive the error message: 

Error in anovaAsterOrReasterList(allargs, tolerance) :
model matrices for fixed effects not nested
model 1 and model 2

Any idea what could be causing this? They seem properly nested to me...

Thanks,

John Benning

PhD Student, Moeller Lab, UMN PMB

[John Benning]

I got it to run by specifying the tolerance exponent = 0.5:

anova(t1a.aster.full_noBlock, t1a.aster.full, tolerance = .Machine$double.eps ^ 0.5)

I think the default exponent is 0.75. What exactly is this tolerance controlling?

John

[John Benning]

Also works when tolerance exponent = 0.7. But even testing my simpler models (eg, main effects only), I have to increase the tolerance. 

John

[geyer]

What the tolerance is for is that the two model matrices (for the big model and the little model) are expressed in the computer's so-called real numbers, which have about 16 significant figure.  So we do not expect the column space of the little model to be exactly contained in the column space of the big model.  Every comparison involving the computer's so-called real numbers cannot be expected to be exact and must involve a tolerance.  To tell whether the model matrices have nested column spaces, we regress the columns of the little model on those of the big model.  This involves a lot of calculation (QR decomposition of the model matrix of the big model) and inaccuracies due to inexactness of computer arithmetic build up.  Depending on how well- or ill-conditioned these model matrices are we can need a large tolerance.

I perhaps see the problem. total_seeds:site/block isn't supposed to work with aster models.  It is an lme-ism AFAIK.  The / operator isn't documented in ?formula.


So having no idea what that model is supposed to be or what the corresponding model matrix is, I would guess that the models are not nested.  So anova.asterOrReaster has done you a big favor.





On Friday, December 21, 2018 at 1:48:12 PM UTC-6, John Benning wrote:

[John Benning]

Thanks for the quick response! And apologies for not posting the full code and data. Data are here and code is attached.

Good to know exactly what the tolerance argument is in reference to (I assumed it was computer arithmetic but didn't know exactly how that related to comparing nesting of models).

A brief summary of the experiment:

These are data from a transplant experiment with Clarkia xantiana in southern California, similar in design to the experiment Ryan Briscoe Runquist was analyzing with Aster previously. I had six sites (with six blocks per site), three seed sources (populations), and a caging treatment to test for effects of biotic interactions across the species' geographic range. The aster LHS graph is: germination, early survival, late survival, survival to flowering, any fruits produced, total seeds produced.

Re: nesting terms:

I believe R / aster may be handling the "/" operator correctly -- as I understand it, its functionality for nesting works in base R (eg, lm()). In this case, I think it just translates to site:block, but keeps each set of blocks within each site, since the blocks are not uniquely coded (i.e., blocks are simply labeled 1-6 in each site). If my blocks were uniquely coded (eg, s1.1, s1.2,...s4.1, s4.2,...), I could simply include site + site:block_unique to nest block properly. I include an example in the attached file. Indeed, the models using site/block and site:block_unique appear to be the same model. 

Interestingly, however, while the "/" operator seems to "alert" R / aster that it shouldn't try to include block interactions across sites (ie, it keeps blocks within their sites), when I use unique block IDs, the "nonsense" predictor variables (eg, site1 : site2.block1, site2 : site3.block4) are identified as aliased and dropped. It seems aster tries to include them but then drops them (but maybe that's what R does behind the scenes and doesn't tell us).

Re: aliased terms

I think I understand why total_seeds:siteWB would be aliased, as there were 0 seeds collected from that site, so it would be aliased with the total_seeds:siteWB:Seed variables in that site, which, obviously, had zero seeds collected. I guess I'm a little unclear when we should care much about aliased terms (beyond examples like those in the aster tutorial) -- is it sometimes simply unavoidable?

Re: tolerance, I still have to adjust the argument exponent to 0.7 even when testing main effects in simple models like,

t1a.mains <- aster(resp ~ varb + 

                    total_seeds:site +      

                    total_seeds:Seed + 

                    total_seeds:Caged,

                    pred1a, fam1a, varb, id, root, data = transplant_yr1a_aster)

through backward step-wise deletion.

John

[John Benning]

Shoot data link didn't go through. Here they are.