[R-sig-eco] mvabund-interprating the uni.test
David Warton
david.warton at unsw.edu.au
Fri Oct 24 00:50:07 CEST 2014
Hi Kendra,
Which values are in uni.test depends which test statistic you choose to use (in the test argument) - the default is log-likelihood ratio statistics (not sum-of-LR, but the original LR's that later get summed to make a multivariate statistic).
For something in the output that tells you the direction of the effect, try model coefficients, easy enough to access using something like:
coef(glm.spid)
I'd advise plotting too! Apart from using a residual plot to check your assumptions
plot(glm.spid)
you could for example plot a subset corresponding to the spp with largest test statistics. Something like
uniSorted = sort(an$uni.test, index.return=TRUE, decreasing=TRUE)
Then plot
spiddat[ , uniSorted$ix[1:10] ]
all the best
David
-----Original Message-----
From: Maas, Kendra [mailto:kendrami at mail.ubc.ca]
Sent: Friday, 24 October 2014 9:27 AM
To: David Warton; 'r-sig-ecology at r-project.org'
Subject: RE: [R-sig-eco] mvabund-interprating the uni.test
thanks for responding. To clarify, the values in uni.test are the sum-of-LR (I didn't specify which test to use, default is LR)? So that value gives me the relative influence of treatment on each species but not the direction, correct. Is there something in the output or either manyglm or the anova that gives the direction? My treatments are a gradient from control to harshest, so negative intercept (positive slope) are species that increase with treatment-that is what I'm trying to find for the ones that have a significant p. For my datasets, I have ~80 species that are responding so would like to find a way to assess direction of impact other than plotting each one.
thanks
Kendra
--
Kendra Maas, Ph.D.
Post Doctoral Research Fellow
University of British Columbia
604-822-5646
________________________________________
From: David Warton [david.warton at unsw.edu.au]
Sent: Thursday, October 23, 2014 2:57 PM
To: 'r-sig-ecology at r-project.org'; Maas, Kendra
Subject: RE: [R-sig-eco] mvabund-interprating the uni.test
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