Understanding errors when comparing procD.pgls outputs

[Ben Sulser]

Dear All - 

I'm encountering an error message when trying to compare outputs of the "procD.pgls()" function, and I'm hoping someone may be able to put me on the right track. Utilizing both the anova() and model.comparison() functions throw an error on my end. To recreate this: 

******************

> data(plethspecies)
> Y.gpa<-gpagen(plethspecies$land)

Performing GPA
  |==============================================================================================================| 100%

Making projections... Finished!
> gdf <- geomorph.data.frame(Y.gpa, phy = plethspecies$phy)
> pleth.pgls <- procD.pgls(coords ~ Csize, phy = phy, data = gdf, iter = 999)

Please be aware that printing progress slows down the analysis (perhaps slightly).

Preliminary Model Fit...

Sums of Squares calculations: 1000 permutations.
  |==============================================================================================================| 100%
> pleth.null <- procD.pgls(coords ~ 1, phy = phy, data = gdf, iter = 999)

Please be aware that printing progress slows down the analysis (perhaps slightly).

Preliminary Model Fit...

Sums of Squares calculations: 1000 permutations.
  |==============================================================================================================| 100%

No terms for ANOVA; only RSS calculated in each permutation
> model.comparison(pleth.null, pleth.pgls, type = "logLik")
Error in sqrt(w) : non-numeric argument to mathematical function
> anova(pleth.null, pleth.pgls)
Error in sqrt(refModel$LM$weights) :
  non-numeric argument to mathematical function

*************************************

This does not seem to occur with OLS models:

*************************************

> pleth.Lm <- procD.lm(coords ~ Csize, data = gdf, iter = 999)
> pleth.nullLm <- procD.lm(coords ~ 1, data = gdf, iter = 999)
> model.comparison(pleth.nullLm, pleth.Lm, type = "logLik")
        logLik residual.pc.no penalty       AIC
1     308.9513              8      88 -529.9027
Csize 119.6700              8     104 -135.3401
> anova(pleth.nullLm, pleth.Lm)

Sums of Squares calculations for 1 model: 1000 permutations.
  |==============================================================================================================| 100%

Analysis of Variance, using Residual Randomization
Permutation procedure: Randomization of null model residuals
Number of permutations: 1000
Estimation method: Ordinary Least Squares
Effect sizes (Z) based on F distributions

                  ResDf Df       RSS        SS        MS     Rsq      F       Z     P Pr(>F)
coords ~ 1 (Null)     8  1 0.0047685                     0.00000                            
coords ~ Csize        7  1 0.0041314 0.0006371 0.0006371 0.13361 1.0795 0.27376 0.395      
Total                 8    0.0047685          

******************************* 

Based on the error messages, in both cases it looks like a missing argument that should be stored in "Model_Name$LM$weights" is to blame (I do see the printed message about ANOVA terms during the calculation of the null model, but I'm having this issue with more complex analyses as well). 

Do I need to calculate/provide these weights manually to make such a comparison? Or perhaps I am missing a key concept in how GLS results should be compared and interpreted?

Thank you for your time and help, regardless!

Best regards, 
~Ben

[Mike Collyer]

Dear Ben,

There is an inherent problem with the anova.procD.lm and model.comparison functions to find some of the output it needs to do their jobs, which I will have to investigate.  This problem can be temporarily avoided, if you add verbose = TRUE as arguments to the model fits, e.g., 

pleth.null <- procD.pgls(coords ~ 1, phy = phy, data = gdf, iter = 999, verbose = TRUE)

pleth.pgls <- procD.pgls(coords ~ Csize, phy = phy, data = gdf, iter = 999, verbose = TRUE)

This assures all output is available.  I’ll have to see why the functions are not capable of generating the output when verbose = FALSE.  (Using verbose = FALSE means not appending a lot of results that are rarely used, a benefit when one has large data sets.)

Thanks for the alert!

Mike

--
You received this message because you are subscribed to the Google Groups "geomorph R package" group.
To unsubscribe from this group and stop receiving emails from it, send an email to geomorph-r-pack...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/cfdb50d2-fd48-4e54-b8c7-f1e46c9b630fn%40googlegroups.com.

[Ben Sulser]

Hi Mike- 

Just tested on my end with the verbose = TRUE, everything seems to be working. Thanks for the help & speedy response! 

Best, 
~Ben

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/BA85E48A-40AB-4CD2-AD34-A4CF619A0D3B%40gmail.com.

[Mike Collyer]

Hi Ben,

I’m glad it worked.  I have also found and eliminated the bugs that caused the issues, which were in the RRPP package rather than the geomorphology package.  An update can be readily installed via Github, e.g.,

devtools::install_github("mlcollyer/RRPP”)

I will likely update the package on CRAN in a week or two, after I check thoroughly for other issues.

Regards,

Mike

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/CAP7QoF036eCwatB1OCDGwWARtX13ZS37cHnnK1VYtvMe3zBzRQ%40mail.gmail.com.