[R-sig-eco] mvabund-interprating the uni.test
David Warton
david.warton at unsw.edu.au
Thu Oct 23 23:57:13 CEST 2014
Hi Kendra,
Thanks for the question. It is great to see you looking at the relative size/significance of univariate tests as a way to gauge which species most strongly respond to some environmental gradient (/treatment) - this is a much more reliable way to do things than using something like SIMPER, which muddles up effect size and residual variability.
The univariate P-values are found in uni.p, the second last argument listed in your str call. You can find out what else is in your anova.manyglm object by typing
?anova.manyglm
And looking at the section titled "Value", which lists the output arguments with brief descriptions of what they show.
All the best
David
David Warton
Associate Professor and Australian Research Council Future Fellow
School of Mathematics and Statistics and the Evolution & Ecology Research Centre
The University of New South Wales NSW 2052 AUSTRALIA
phone (61)(2) 9385-7031
fax (61)(2) 9385-7123
http://www.eco-stats.unsw.edu.au/
-------- Original Message --------
Subject: [R-sig-eco] mvabund-interprating the uni.test
Date: Wed, 22 Oct 2014 22:07:38 +0000
From: Maas, Kendra <kendrami at mail.ubc.ca>
To: r-sig-ecology at r-project.org <r-sig-ecology at r-project.org>
I'm using mvabund uni.test to select species that are changing according to a treatment gradient. I've pulled out species that have a significantly different intercept from the null but I'd also like to indicate if they are positive or negative with treatment and extract the "sum-of-LR" as Warton does in the paper showing the usefulness of this approach (http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2011.00127.x/abstract). I've run manyglm followed by anova on the spider data to try to match the mvabund documentation but I can't figure out where the results beyond the p-values are stored, here are the str for the manyglm and anova output. (the dput are huge)
thanks
Kendra
> str(glm.spid)
List of 41
$ coefficients : num [1:3, 1:12] 2.2353 0.0582 -1.3425 2.1552 -0.4229 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ var.coefficients : num [1:3, 1:12] 0.1068 0.0185 0.1351 0.252 0.0466 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ fitted.values : num [1:28, 1:12] 9.349 0.843 9.349 9.349 9.349 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ linear.predictor : num [1:28, 1:12] 2.24 -0.17 2.24 2.24 2.24 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ residuals : num [1:28, 1:12] 1.767 -0.778 0.638 -0.83 -0.943 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ PIT.residuals : num [1:28, 1:12] 0.934 0.534 0.79 0.238 0.101 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ sqrt.1_Hii : num [1:28, 1:12] 0.945 0.89 0.945 0.945 0.945 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ var.estimator : num [1:28, 1:12] 87.79 1.48 87.79 87.79 87.79 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ sqrt.weight : num [1:28, 1:12] 0.998 0.693 0.998 0.998 0.998 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ theta : Named num [1:12] 1.114 0.428 1.301 0.156 0.715 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ two.loglike : Named num [1:12] -121.5 -142.1 -78.9 -60.6 -44.5 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ deviance : Named num [1:12] 22.2 29.4 18.7 16 10.5 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ aic : Named num [1:12] 129.5 150.1 86.9 68.6 52.5 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ iter : Named num [1:12] 3 3 6 4 5 4 4 5 4 4 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ data :'data.frame': 28 obs. of 3 variables:
..$ spiddat : int [1:28, 1:12] 25 0 15 2 1 0 2 0 1 3 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
.. ..- attr(*, "class")= chr [1:2] "mvabund" "matrix"
..$ bare.sand : num [1:28] 0 0 0 0 0 ...
..$ fallen.leaves: num [1:28] 0 1.79 0 0 0 ...
..- attr(*, "terms")=Classes 'terms', 'formula' length 3 spiddat ~ bare.sand + fallen.leaves
.. .. ..- attr(*, "variables")= language list(spiddat, bare.sand, fallen.leaves)
.. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:3] "spiddat" "bare.sand" "fallen.leaves"
.. .. .. .. ..$ : chr [1:2] "bare.sand" "fallen.leaves"
.. .. ..- attr(*, "term.labels")= chr [1:2] "bare.sand" "fallen.leaves"
.. .. ..- attr(*, "order")= int [1:2] 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. .. ..- attr(*, "predvars")= language list(spiddat, bare.sand, fallen.leaves)
.. .. ..- attr(*, "dataClasses")= Named chr [1:3] "nmatrix.12" "numeric" "numeric"
.. .. .. ..- attr(*, "names")= chr [1:3] "spiddat" "bare.sand" "fallen.leaves"
$ stderr.coefficients: num [1:3, 1:12] 0.327 0.136 0.368 0.502 0.216 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ phi : Named num [1:12] 0.897 2.338 0.768 6.43 1.399 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ tol : num 1.49e-08
$ maxiter : num 25
$ maxiter2 : num 10
$ AIC : Named num [1:12] 129.5 150.1 86.9 68.6 52.5 ...
..- attr(*, "names")= chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ AICsum : num 1658
$ family : chr "negative.binomial"
$ K : num 1
$ theta.method : chr "PHI"
$ cor.type : chr "I"
$ shrink.param : num 0
$ call : language manyglm(formula = spiddat ~ bare.sand + fallen.leaves, family = "negative.binomial", data = X)
$ terms :Classes 'terms', 'formula' length 3 spiddat ~ bare.sand + fallen.leaves
.. ..- attr(*, "variables")= language list(spiddat, bare.sand, fallen.leaves)
.. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:3] "spiddat" "bare.sand" "fallen.leaves"
.. .. .. ..$ : chr [1:2] "bare.sand" "fallen.leaves"
.. ..- attr(*, "term.labels")= chr [1:2] "bare.sand" "fallen.leaves"
.. ..- attr(*, "order")= int [1:2] 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(spiddat, bare.sand, fallen.leaves)
.. ..- attr(*, "dataClasses")= Named chr [1:3] "nmatrix.12" "numeric" "numeric"
.. .. ..- attr(*, "names")= chr [1:3] "spiddat" "bare.sand" "fallen.leaves"
$ assign : int [1:3] 0 1 2
$ formula :Class 'formula' length 3 spiddat ~ bare.sand + fallen.leaves
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ rank : int 3
$ xlevels : Named list()
$ df.residual : int 25
$ x : num [1:28, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
$ y : num [1:28, 1:12] 25 0 15 2 1 0 2 0 1 3 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ qr :List of 4
..$ qr : num [1:28, 1:3] -5.292 0.189 0.189 0.189 0.189 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:28] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
..$ rank : int 3
..$ qraux: num [1:3] 1.19 1.11 1.18
..$ pivot: int [1:3] 1 2 3
..- attr(*, "class")= chr "qr"
$ show.coef : logi FALSE
$ show.fitted : logi FALSE
$ show.residuals : logi FALSE
$ offset : num [1:28, 1:12] 0 0 0 0 0 0 0 0 0 0 ...
- attr(*, "class")= chr [1:2] "manyglm" "mglm"
> str(an.spid)
List of 11
$ family : chr "negative.binomial"
$ p.uni : chr "adjusted"
$ test : chr "Dev"
$ cor.type : chr "I"
$ resamp : chr "pit.trap"
$ nBoot : num 1000
$ shrink.param: num [1:3] 0 0 0
$ n.bootsdone : num 999
$ table :'data.frame': 3 obs. of 4 variables:
..$ Res.Df : int [1:3] 27 26 25
..$ Df.diff : num [1:3] NA 1 1
..$ Dev : num [1:3] NA 70.4 70.3
..$ Pr(>Dev): num [1:3] NA 0.002 0.007
..- attr(*, "heading")= chr [1:3] "Analysis of Deviance Table\n" "Model: manyglm(formula = spiddat ~ bare.sand + fallen.leaves, family = \"negative.binomial\", " "Model: data = X)"
..- attr(*, "title")= chr "\nMultivariate test:\n"
$ uni.p : num [1:3, 1:12] NA 0.568 0.001 NA 0.485 0.988 NA 0.001 0.074 NA ...
..- attr(*, "title")= chr "\nUnivariate Tests:\n"
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
$ uni.test : num [1:3, 1:12] NA 1.22 27.64 NA 2.42 ...
..- attr(*, "title")= chr "\nUnivariate Tests:\n"
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "(Intercept)" "bare.sand" "fallen.leaves"
.. ..$ : chr [1:12] "Alopacce" "Alopcune" "Alopfabr" "Arctlute" ...
- attr(*, "class")= chr "anova.manyglm"
--
Kendra Maas, Ph.D.
Post Doctoral Research Fellow
University of British Columbia
604-822-5646
_______________________________________________
R-sig-ecology mailing list
R-sig-ecology at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
More information about the R-sig-ecology
mailing list