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:kendr...@mail.ubc.ca]
Sent: Friday, 24 October 2014 9:27 AM
To: David Warton; 'r-sig-ecology@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.war...@unsw.edu.au]
Sent: Thursday, October 23, 2014 2:57 PM
To: 'r-sig-ecology@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 <kendr...@mail.ubc.ca>
To: r-sig-ecology@r-project.org <r-sig-ecology@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@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
_______________________________________________
R-sig-ecology mailing list
R-sig-ecology@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
