First things first
Let’s load the package and get our hands on theTasmaniadata set to look at the effects of disturbance treatment on invertebrateabundances. Note that theTasmania data set is a listobject. We will only look at thecopepods data frame forthis walk-through. Thecopepods data frame can be accessedusingTasmania$copepods or theattach()function which will make the contents ofTasmaniasearchable.
library(mvabund)data(Tasmania)attach(Tasmania)#> The following objects are masked from reveg:#>#> abund, treatmentskimr::skim(copepods)# Great function to get an overview of the data
Data summary| Name | copepods |
| Number of rows | 16 |
| Number of columns | 12 |
| _______________________ | |
| Column type frequency: | |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|
| Ameira | 0 | 1 | 55.44 | 46.60 | 4 | 6.75 | 58.5 | 92.25 | 142 | ▇▂▃▃▂ |
| Adopsyllus | 0 | 1 | 0.69 | 1.25 | 0 | 0.00 | 0.0 | 1.00 | 4 | ▇▂▁▁▁ |
| Ectinosoma | 0 | 1 | 1.31 | 2.30 | 0 | 0.00 | 0.0 | 1.25 | 7 | ▇▁▁▁▁ |
| Ectinosomat | 0 | 1 | 4.50 | 4.46 | 0 | 1.00 | 3.5 | 5.50 | 15 | ▇▃▂▁▂ |
| Haloschizo | 0 | 1 | 0.12 | 0.50 | 0 | 0.00 | 0.0 | 0.00 | 2 | ▇▁▁▁▁ |
| Lepta.A | 0 | 1 | 40.69 | 47.15 | 0 | 3.00 | 28.0 | 57.25 | 151 | ▇▂▁▂▁ |
| Lepta.B | 0 | 1 | 1.44 | 2.92 | 0 | 0.00 | 0.0 | 1.25 | 11 | ▇▁▁▁▁ |
| Lepta.C | 0 | 1 | 12.75 | 44.73 | 0 | 0.00 | 0.0 | 1.50 | 180 | ▇▁▁▁▁ |
| Mictyricola | 0 | 1 | 1.31 | 2.33 | 0 | 0.00 | 0.0 | 1.50 | 8 | ▇▁▁▁▁ |
| Parevansula | 0 | 1 | 0.25 | 0.58 | 0 | 0.00 | 0.0 | 0.00 | 2 | ▇▁▁▁▁ |
| Quin | 0 | 1 | 0.44 | 0.96 | 0 | 0.00 | 0.0 | 0.00 | 3 | ▇▁▁▁▁ |
| Rhizothrix | 0 | 1 | 0.81 | 2.04 | 0 | 0.00 | 0.0 | 0.00 | 6 | ▇▁▁▁▁ |
Visualise the multivariate data
We first need to turn our data into amvabund object sofunctions for this package and work with the data
copepod_abund<-mvabund(copepods)
Now lets take a look at abundance for each species across ourtreatment sites (Disturbed vs. Undisturbed). Observations were collectedusing a spatially blocked design where researchers took four samples ateach block (2 per treatment). We can set the colour (col)of the points to represent that four sampling blocks
plot(copepod_abund~treatment,col = block)#> Overlapping points were shifted along the y-axis to make them visible.#>#> PIPING TO 2nd MVFACTOR
Fitting Predictive Models
It was hypothesised that the abundance ofAmeira andEctinosoma was reduced in Disturbed sites, whereas theabundance ofMictyricola may have increased. Lets test thishypothesis using themanyglm() function. This function fitsa generalised linear model for each species. We specifiedfamily = "negative.binomial" as count data tends to followa negative binomial distribution. Other distributions are available too!See?manyglm()
cope.nb<-manyglm(copepods~ treatment*block,family ="negative.binomial")
Checking Model Assumptions
Before we look at the model output, we should check on the modelresiduals. What we want to see islittle pattern as thisimplies that our choice of negative binomial distribution isappropriate.
Now, lets proceed to check on themean-variancerelationship. We want to to see if themean-variancerelationship of our data adheres to that of a negative binomialdistribution which tends to be quadratic rather than linear. Themeanvar.plot() function plots the sample variance againstthe sample mean for each species within each factor level of(tr.block). A quadratic relationship seems appropriate forour sample mean and variance.
meanvar.plot(copepods~tr.block,col = treatment)#> START SECTION 2#> Plotting if overlay is TRUE#> using grouping variable tr.block 46 mean values were 0 and could#> not be included in the log-plot#> using grouping variable tr.block 49 variance values were 0 and could not#> be included in the log-plot
#> FINISHED SECTION 2
Hypothesis Testing
To test whethertreatment andblock had aneffect on the abundances of copepods we can use theanova()function. This function returns a Analysis of Deviance table which teststhe significance of each model term. Settingp.uni = "adjusted" allows for our p-values to be adjustedfor multiple testing of different species.
anova(cope.nb,p.uni ="adjusted")#> Time elapsed: 0 hr 0 min 34 sec#> Analysis of Deviance Table#>#> Model: copepods ~ treatment * block#>#> Multivariate test:#> Res.Df Df.diff Dev Pr(>Dev)#> (Intercept) 15#> treatment 14 1 32.52 0.047 *#> block 11 3 151.87 0.001 ***#> treatment:block 8 3 37.41 0.081 .#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#>#> Univariate Tests:#> Ameira Adopsyllus Ectinosoma#> Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)#> (Intercept)#> treatment 6.257 0.150 0.032 0.978 7.026 0.144#> block 12.615 0.153 17.691 0.054 15.361 0.072#> treatment:block 6.349 0.510 1.295 0.868 0.836 0.868#> Ectinosomat Haloschizo Lepta.A#> Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)#> (Intercept)#> treatment 0.366 0.943 1.497 0.805 0.281 0.943#> block 10.92 0.186 3.742 0.246 20.756 0.022#> treatment:block 1.122 0.868 0 0.868 13.072 0.162#> Lepta.B Lepta.C Mictyricola#> Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)#> (Intercept)#> treatment 0.575 0.926 2.141 0.748 7.093 0.144#> block 8.8 0.186 17.419 0.056 11.415 0.186#> treatment:block 7.402 0.453 0.34 0.868 5.268 0.510#> Parevansula Quin Rhizothrix#> Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)#> (Intercept)#> treatment 0 0.996 5.385 0.202 1.869 0.748#> block 6.1 0.186 9.44 0.186 17.613 0.056#> treatment:block 1.726 0.835 0 0.868 0 0.868#> Arguments:#> Test statistics calculated assuming uncorrelated response (for faster computation)#> P-value calculated using 999 iterations via PIT-trap resampling.
The test statistics in this table are by default calculated bysumming the change in deviance across all responses. You could use adifferent type of test statistic by changing thetest andcor.type arguments. We can see that there is asignificant effect of the treatment factor meaning thattreatment has a significant multiplicative effect on mean abundance. Theinteraction between blocks and treatments isnotsignificant, meaning that the multiplicative treatment effectis consistent across blocks.
If you do not have a specific hypothesis in mind that you want totest, and are instead interested in which model terms are statisticallysignificant, then thesummary() function will come inhandy. However results aren’t quite as trustworthy as foranova(). The reason is that re-samples are taken under thealternative hypothesis forsummary(), where there is agreater chance of fitted values being zero, especially for rarer taxa(e.g. if there is a treatment combination in which a taxon is neverpresent). Abundances don’t re-sample well if their predicted mean iszero.
summary(cope.nb)#>#> Test statistics:#> wald value Pr(>wald)#> (Intercept) 18.493 0.001 ***#> treatmentUndisturbed 3.330 0.238#> block2 4.710 0.052 .#> block3 6.650 0.001 ***#> block4 3.435 0.156#> treatmentUndisturbed:block2 2.747 0.266#> treatmentUndisturbed:block3 1.615 0.462#> treatmentUndisturbed:block4 4.262 0.043 *#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#>#> Test statistic: 14.81, p-value: 0.002#> Arguments:#> Test statistics calculated assuming response assumed to be uncorrelated#> P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).
If obtaining predicted values from the model is the goal, you may usethepredict() function. Note thattype = response will produce values on the scale of theresponse variable (i.e. counts)
predict(cope.nb,type ="response")#> Ameira Adopsyllus Ectinosoma Ectinosomat Haloschizo Lepta.A#> B1D1 53.0 3.678794e-07 3.678794e-07 8.0 3.678794e-07 63.5#> B1D2 53.0 3.678794e-07 3.678794e-07 8.0 3.678794e-07 63.5#> B1U1 114.5 3.678794e-07 1.000000e+00 7.0 1.000000e+00 134.0#> B1U2 114.5 3.678794e-07 1.000000e+00 7.0 1.000000e+00 134.0#> B2D1 4.5 3.678794e-07 3.678794e-07 9.0 3.678794e-07 31.0#> B2D2 4.5 3.678794e-07 3.678794e-07 9.0 3.678794e-07 31.0#> B2U1 74.0 3.678794e-07 3.678794e-07 4.5 3.678794e-07 51.5#> B2U2 74.0 3.678794e-07 3.678794e-07 4.5 3.678794e-07 51.5#> B3D1 6.5 3.678794e-07 3.678794e-07 2.5 3.678794e-07 2.0#> B3D2 6.5 3.678794e-07 3.678794e-07 2.5 3.678794e-07 2.0#> B3U1 35.0 5.000000e-01 2.500000e+00 2.5 3.678794e-07 1.5#> B3U2 35.0 5.000000e-01 2.500000e+00 2.5 3.678794e-07 1.5#> B4D1 37.0 2.500000e+00 5.000000e-01 1.0 3.678794e-07 38.0#> B4D2 37.0 2.500000e+00 5.000000e-01 1.0 3.678794e-07 38.0#> B4U1 119.0 2.500000e+00 6.500000e+00 1.5 3.678794e-07 4.0#> B4U2 119.0 2.500000e+00 6.500000e+00 1.5 3.678794e-07 4.0#> Lepta.B Lepta.C Mictyricola Parevansula Quin#> B1D1 6.000000e+00 3.678794e-07 3.678794e-07 3.678794e-07 3.678794e-07#> B1D2 6.000000e+00 3.678794e-07 3.678794e-07 3.678794e-07 3.678794e-07#> B1U1 3.678794e-07 3.678794e-07 3.678794e-07 5.000000e-01 2.500000e+00#> B1U2 3.678794e-07 3.678794e-07 3.678794e-07 5.000000e-01 2.500000e+00#> B2D1 1.500000e+00 3.678794e-07 5.500000e+00 1.000000e+00 3.678794e-07#> B2D2 1.500000e+00 3.678794e-07 5.500000e+00 1.000000e+00 3.678794e-07#> B2U1 3.500000e+00 3.678794e-07 3.678794e-07 5.000000e-01 3.678794e-07#> B2U2 3.500000e+00 3.678794e-07 3.678794e-07 5.000000e-01 3.678794e-07#> B3D1 3.678794e-07 9.500000e+01 5.000000e-01 3.678794e-07 3.678794e-07#> B3D2 3.678794e-07 9.500000e+01 5.000000e-01 3.678794e-07 3.678794e-07#> B3U1 3.678794e-07 5.000000e+00 5.000000e-01 3.678794e-07 3.678794e-07#> B3U2 3.678794e-07 5.000000e+00 5.000000e-01 3.678794e-07 3.678794e-07#> B4D1 5.000000e-01 2.000000e+00 4.000000e+00 3.678794e-07 3.678794e-07#> B4D2 5.000000e-01 2.000000e+00 4.000000e+00 3.678794e-07 3.678794e-07#> B4U1 3.678794e-07 3.678794e-07 3.678794e-07 3.678794e-07 1.000000e+00#> B4U2 3.678794e-07 3.678794e-07 3.678794e-07 3.678794e-07 1.000000e+00#> Rhizothrix#> B1D1 5.000000e-01#> B1D2 5.000000e-01#> B1U1 6.000000e+00#> B1U2 6.000000e+00#> B2D1 3.678794e-07#> B2D2 3.678794e-07#> B2U1 3.678794e-07#> B2U2 3.678794e-07#> B3D1 3.678794e-07#> B3D2 3.678794e-07#> B3U1 3.678794e-07#> B3U2 3.678794e-07#> B4D1 3.678794e-07#> B4D2 3.678794e-07#> B4U1 3.678794e-07#> B4U2 3.678794e-07