QK Box 16.3

Parkinson et al. (2020) set tested experimentally the effect of artificial light at night (ALAN) and fish presence on abundances of terrestrial and aquatic invertebrates in the littoral zone of a small lake. They established 20 mesocosms (1 m2 plastic mesh walls and bottom submerged with part of the walls above the water surface) with ten provided with ALAN (solar lights in each corner) and ten controls. Half the mesocosms in each group had fish added, but there were no effects of fish on any variables, so Parkinson et al. ignored fish in the subsequent analyses. After about six weeks, a single pan trap was placed in each mesocosm overnight, and invertebrates were identified to family and counted. We will use these data to do an ordination of mesocosms based on invertebrate family abundances and then fit models to examine differences between the two treatment groups (ALAN vs. no-ALAN).

The paper is here

Parkinson, E., Lawson, J. & Tiegs, S. D. (2020). Artificial light at night at the terrestrial-aquatic interface: Effects on predators and fluxes of insect prey. PLoS One, 15, e0240138.

Preliminaries

# Load standard set of graphics packages and other libraries
source("../R/appearance.R")
source("../R/libraries.R")

Load packages specific to this example: vegan,mvabund

library(vegan)
library(mvabund)

Read in data file (park2020.csv) and remove labels

park <- read.csv("../data/park2020.csv")
park1 <- park[,-(1:3)]

MDS on invert abundances

park1.bc <- vegdist(park1,'bray')
park1.mds <- metaMDS(park1.bc,k=2,autotransform=FALSE)

Run 0 stress 0.0314779 
Run 1 stress 0.03147804 
... Procrustes: rmse 0.000362395  max resid 0.001017012 
... Similar to previous best
Run 2 stress 0.03171596 
... Procrustes: rmse 0.004781443  max resid 0.01452972 
Run 3 stress 0.04438 
Run 4 stress 0.04025036 
Run 5 stress 0.03147831 
... Procrustes: rmse 0.0005104671  max resid 0.001430591 
... Similar to previous best
Run 6 stress 0.03426509 
Run 7 stress 0.03148572 
... Procrustes: rmse 0.001484178  max resid 0.004134368 
... Similar to previous best
Run 8 stress 0.0370026 
Run 9 stress 0.04025026 
Run 10 stress 0.03147912 
... Procrustes: rmse 0.0007710645  max resid 0.002162996 
... Similar to previous best
Run 11 stress 0.03147796 
... Procrustes: rmse 0.0003048454  max resid 0.0008546811 
... Similar to previous best
Run 12 stress 0.03150071 
... Procrustes: rmse 0.002858939  max resid 0.008187161 
... Similar to previous best
Run 13 stress 0.03147832 
... Procrustes: rmse 0.0005116255  max resid 0.001435861 
... Similar to previous best
Run 14 stress 0.03148999 
... Procrustes: rmse 0.002114141  max resid 0.006009197 
... Similar to previous best
Run 15 stress 0.03148076 
... Procrustes: rmse 0.001111898  max resid 0.003117315 
... Similar to previous best
Run 16 stress 0.03148167 
... Procrustes: rmse 0.001250987  max resid 0.003506942 
... Similar to previous best
Run 17 stress 0.03170596 
... Procrustes: rmse 0.003903063  max resid 0.0123791 
Run 18 stress 0.03426473 
Run 19 stress 0.03147825 
... Procrustes: rmse 0.0004821467  max resid 0.001352596 
... Similar to previous best
Run 20 stress 0.04150753 
*** Best solution repeated 11 times

stressplot(park1.mds, main="Shepard plot")

park1.mds

Call:
metaMDS(comm = park1.bc, k = 2, autotransform = FALSE) 

global Multidimensional Scaling using monoMDS

Data:     park1.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.0314779 
Stress type 1, weak ties
Best solution was repeated 11 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

ordiplot(park1.mds$points,display="sites",type="n")
points(park1.mds$points,pch=as.numeric(park$symb))

With standardisation (wisconsin)

park1s <- wisconsin(park1)
park1s.bc <- vegdist(park1s,'bray')   
park1s.mds <- metaMDS(park1s.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1448505 
Run 1 stress 0.1448503 
... New best solution
... Procrustes: rmse 0.0005742005  max resid 0.001327975 
... Similar to previous best
Run 2 stress 0.1553992 
Run 3 stress 0.1547493 
Run 4 stress 0.1448503 
... Procrustes: rmse 0.0004012248  max resid 0.0009023985 
... Similar to previous best
Run 5 stress 0.1537683 
Run 6 stress 0.1448505 
... Procrustes: rmse 0.0005829848  max resid 0.001384117 
... Similar to previous best
Run 7 stress 0.1631977 
Run 8 stress 0.1552571 
Run 9 stress 0.1459782 
Run 10 stress 0.1553015 
Run 11 stress 0.1537685 
Run 12 stress 0.1459783 
Run 13 stress 0.1448505 
... Procrustes: rmse 0.0006084014  max resid 0.001457572 
... Similar to previous best
Run 14 stress 0.1448505 
... Procrustes: rmse 0.000559664  max resid 0.001322976 
... Similar to previous best
Run 15 stress 0.157652 
Run 16 stress 0.1448506 
... Procrustes: rmse 0.0006319359  max resid 0.001482624 
... Similar to previous best
Run 17 stress 0.1550382 
Run 18 stress 0.1459783 
Run 19 stress 0.1609314 
Run 20 stress 0.1459786 
Run 21 stress 0.1448503 
... New best solution
... Procrustes: rmse 0.0003238484  max resid 0.000775218 
... Similar to previous best
Run 22 stress 0.1448502 
... New best solution
... Procrustes: rmse 0.0001464993  max resid 0.0003554879 
... Similar to previous best
Run 23 stress 0.1459782 
Run 24 stress 0.1459782 
Run 25 stress 0.1547494 
Run 26 stress 0.1448503 
... Procrustes: rmse 0.0001443632  max resid 0.0002833306 
... Similar to previous best
Run 27 stress 0.1554602 
Run 28 stress 0.1459782 
Run 29 stress 0.153768 
Run 30 stress 0.1552822 
Run 31 stress 0.1459784 
Run 32 stress 0.1548378 
Run 33 stress 0.1448504 
... Procrustes: rmse 0.0002492233  max resid 0.0005931013 
... Similar to previous best
Run 34 stress 0.1448504 
... Procrustes: rmse 0.0002816716  max resid 0.000671023 
... Similar to previous best
Run 35 stress 0.1609319 
Run 36 stress 0.1552838 
Run 37 stress 0.1537681 
Run 38 stress 0.1554601 
Run 39 stress 0.1459783 
Run 40 stress 0.1448503 
... Procrustes: rmse 0.0002199271  max resid 0.0005219861 
... Similar to previous best
*** Best solution repeated 5 times

stressplot(park1s.mds, main="Shepard plot")

park1s.mds

Call:
metaMDS(comm = park1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     park1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1448502 
Stress type 1, weak ties
Best solution was repeated 5 times in 40 tries
The best solution was from try 22 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

ordiplot(park1s.mds$points,type="n",display="sites",choices=c(1,3))
points(park1s.mds$points,pch=as.numeric(park$symb))

note that ordihull might be useful to draw boundaries around each group??

Check homogeneity of dispersions

Check first with file of abundances, then proportional abundances

park1.disp <- betadisper(park1.bc,park$treat)
anova(park1.disp)

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups     1 0.01698 0.016977  0.4559 0.5081
Residuals 18 0.67021 0.037234               

park1s.disp <- betadisper(park1s.bc,park$treat)
anova(park1s.disp)

Analysis of Variance Table

Response: Distances
          Df   Sum Sq   Mean Sq F value  Pr(>F)  
Groups     1 0.022832 0.0228325  3.5176 0.07704 .
Residuals 18 0.116837 0.0064909                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Do anosim on abundances and prop. abundances

park1.ano <- anosim(park1.bc, park$treat,permutations=999)
summary(park1.ano)

Call:
anosim(x = park1.bc, grouping = park$treat, permutations = 999) 
Dissimilarity: bray 

ANOSIM statistic R: 0.678 
      Significance: 0.001 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.111 0.152 0.228 0.358 

Dissimilarity ranks between and within classes:
        0%   25%   50%    75% 100%   N
Between  8 95.75 135.5 165.25  190 100
ALAN     1 15.00  36.0 113.00  158  45
no ALAN 10 45.00  64.0  81.00  160  45

park1s.ano <- anosim(park1s.bc, park$treat,permutations=999)
summary(park1s.ano)

Call:
anosim(x = park1s.bc, grouping = park$treat, permutations = 999) 
Dissimilarity: bray 

ANOSIM statistic R: 0.2989 
      Significance: 0.001 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0884 0.1145 0.1374 0.1649 

Dissimilarity ranks between and within classes:
        0%  25% 50%   75% 100%   N
Between  5 72.5 113 151.5  190 100
ALAN     1 25.0  51  77.0  177  45
no ALAN  8 59.0 111 143.5  187  45

# do permanova on abundances
park1.ado <- adonis2(park1.bc~treat,data=park,permutations=999)
print(park1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = park1.bc ~ treat, data = park, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
treat     1   1.6916 0.46132 15.415  0.001 ***
Residual 18   1.9753 0.53868                  
Total    19   3.6669 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# do permanova on prop abundances
park1s.ado <- adonis2(park1s.bc~treat,data=park,permutations=999)
print(park1s.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = park1s.bc ~ treat, data = park, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
treat     1   0.7251 0.12406 2.5493  0.001 ***
Residual 18   5.1195 0.87594                  
Total    19   5.8446 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mvabund

Run using negative binomial for raw abundances and double standardized abundances from earlier

park1mv <- mvabund(park1)
park1mv.mv <- manyglm(park1mv~park$treat,family="negative.binomial")
plot(park1mv.mv)

anova(park1mv.mv)

Time elapsed: 0 hr 0 min 6 sec
Analysis of Deviance Table

Model: park1mv ~ park$treat

Multivariate test:
            Res.Df Df.diff   Dev Pr(>Dev)    
(Intercept)     19                           
park$treat      18       1 106.7    0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

park1smv <- mvabund(park1s)
park1smv.mv <- manyglm(park1smv~park$treat,family="negative.binomial")

Warning: Non-integer data are fitted to the negative.binomial model.

plot(park1smv.mv)

anova(park1smv.mv)

Time elapsed: 0 hr 0 min 1 sec
Analysis of Deviance Table

Model: park1smv ~ park$treat

Multivariate test:
            Res.Df Df.diff   Dev Pr(>Dev)    
(Intercept)     19                           
park$treat      18       1 10.32    0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

SIMPER analysis for group differences

Again, abundances then proportional abundances

park1.sim <- simper(park1, park$treat, permutations=1000)
summary(park1.sim)

Contrast: ALAN_no ALAN 

          average        sd     ratio       ava       avb cumsum      p    
caen      0.53070   0.20750   2.55780 134.20000  16.40000  0.751 0.0010 ***
chiron    0.08170   0.04078   2.00370  22.50000   5.00000  0.867 0.0040 ** 
dolich    0.01470   0.01466   1.00550   2.20000   3.00000  0.888 1.0000    
cerat     0.00790   0.01279   0.61550   1.00000   0.20000  0.899 0.0519 .  
aphid     0.00750   0.00928   0.81220   1.20000   0.30000  0.909 0.0919 .  
ephyd     0.00660   0.00852   0.78060   1.10000   0.40000  0.919 0.9940    
baet      0.00620   0.00540   1.15350   1.20000   0.80000  0.928 1.0000    
chryso    0.00520   0.00625   0.82480   1.00000   0.40000  0.935 1.0000    
culicid   0.00450   0.00600   0.74470   0.90000   0.20000  0.941 0.8601    
tetrag    0.00440   0.00853   0.51620   0.50000   0.10000  0.948 0.1429    
hydro     0.00410   0.00580   0.71390   0.50000   0.30000  0.954 0.9920    
tipulid   0.00300   0.00474   0.62510   0.10000   0.40000  0.958 1.0000    
muscid    0.00280   0.00628   0.43820   0.60000   0.00000  0.962 0.0010 ***
empid     0.00190   0.00480   0.40060   0.20000   0.00000  0.964 0.0010 ***
formic    0.00190   0.00360   0.52220   0.00000   0.30000  0.967 1.0000    
polycen   0.00180   0.00356   0.51820   0.00000   0.30000  0.970 1.0000    
sciomyz   0.00180   0.00456   0.39200   0.20000   0.00000  0.972 0.0010 ***
helod     0.00170   0.00277   0.62610   0.30000   0.10000  0.975 0.9940    
taban     0.00150   0.00457   0.32730   0.10000   0.00000  0.977 0.0010 ***
cosmop    0.00150   0.00457   0.32730   0.10000   0.00000  0.979 0.0010 ***
coen      0.00140   0.00289   0.49910   0.20000   0.10000  0.981 0.9980    
delph     0.00140   0.00324   0.43420   0.30000   0.00000  0.983 0.0010 ***
notonec   0.00130   0.00249   0.53090   0.20000   0.10000  0.985 0.9950    
leptoc    0.00110   0.00315   0.36010   0.10000   0.10000  0.986 1.0000    
noctuid   0.00110   0.00267   0.42410   0.00000   0.20000  0.988 1.0000    
hydroph   0.00100   0.00282   0.34170   0.10000   0.10000  0.989 1.0000    
hydropt   0.00090   0.00271   0.34770   0.10000   0.10000  0.991 1.0000    
aran      0.00090   0.00271   0.34770   0.10000   0.10000  0.992 1.0000    
corix     0.00090   0.00182   0.48520   0.20000   0.00000  0.993 0.0010 ***
limneph   0.00080   0.00167   0.49110   0.20000   0.00000  0.994 0.0010 ***
phryg     0.00070   0.00269   0.26640   0.00000   0.10000  0.995 1.0000    
cicad     0.00070   0.00269   0.26640   0.00000   0.10000  0.996 1.0000    
agrom     0.00070   0.00251   0.27120   0.00000   0.10000  0.997 1.0000    
mymar     0.00070   0.00251   0.27120   0.00000   0.10000  0.998 1.0000    
halip     0.00060   0.00222   0.27940   0.00000   0.10000  0.999 1.0000    
acrid     0.00060   0.00175   0.33090   0.10000   0.00000  1.000 0.0010 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 1000

park1s.sim <- simper(park1s, park$treat, permutations=10000)
summary(park1s.sim)

Contrast: ALAN_no ALAN 

        average      sd   ratio     ava     avb cumsum      p    
dolich  0.05568 0.05036 1.10550 0.06115 0.14459  0.069 0.0451 *  
tipulid 0.05435 0.06561 0.82830 0.03055 0.09695  0.137 0.1938    
formic  0.04838 0.07860 0.61550 0.00000 0.09676  0.197 0.0001 ***
caen    0.03616 0.01885 1.91820 0.08926 0.01838  0.242 0.0005 ***
baet    0.03615 0.02236 1.61730 0.07487 0.06517  0.287 0.2027    
helod   0.03376 0.04641 0.72760 0.05221 0.02576  0.329 0.3398    
polycen 0.03299 0.05306 0.62170 0.00000 0.06598  0.370 0.0001 ***
cerat   0.03128 0.03457 0.90500 0.05570 0.02632  0.409 0.3629    
hydroph 0.03085 0.07579 0.40700 0.01319 0.05115  0.448 0.4751    
coen    0.02717 0.04624 0.58750 0.03912 0.02293  0.481 0.4763    
ephyd   0.02609 0.03149 0.82850 0.04708 0.01636  0.514 0.3082    
aphid   0.02539 0.02526 1.00530 0.05701 0.01516  0.545 0.0281 *  
chiron  0.02316 0.01568 1.47660 0.06636 0.02707  0.574 0.0099 ** 
tetrag  0.02269 0.02822 0.80400 0.04104 0.01191  0.602 0.1924    
chryso  0.02090 0.01881 1.11120 0.03550 0.02412  0.628 0.3283    
empid   0.02003 0.04668 0.42900 0.04005 0.00000  0.653 0.0001 ***
noctuid 0.01971 0.04172 0.47230 0.00000 0.03941  0.678 0.0001 ***
notonec 0.01917 0.03336 0.57470 0.02914 0.01387  0.702 0.4788    
leptoc  0.01892 0.04045 0.46780 0.01752 0.02383  0.725 0.4649    
hydro   0.01827 0.02060 0.88690 0.02471 0.02474  0.748 0.7008    
hydropt 0.01755 0.03856 0.45520 0.01319 0.02456  0.770 0.4700    
aran    0.01755 0.03856 0.45520 0.01319 0.02456  0.792 0.4700    
halip   0.01674 0.05049 0.33170 0.00000 0.03349  0.812 0.0001 ***
culicid 0.01439 0.01842 0.78120 0.02561 0.00793  0.830 0.1964    
sciomyz 0.01412 0.02846 0.49610 0.02824 0.00000  0.848 0.0001 ***
phryg   0.01328 0.04004 0.33170 0.00000 0.02656  0.865 0.0001 ***
cicad   0.01328 0.04004 0.33170 0.00000 0.02656  0.881 0.0001 ***
corix   0.01316 0.02770 0.47510 0.02632 0.00000  0.897 0.0001 ***
limneph 0.01259 0.02624 0.47980 0.02518 0.00000  0.913 0.0001 ***
agrom   0.01147 0.03457 0.33170 0.00000 0.02293  0.927 0.0001 ***
mymar   0.01147 0.03457 0.33170 0.00000 0.02293  0.942 0.0001 ***
muscid  0.01124 0.02506 0.44870 0.02249 0.00000  0.956 0.0001 ***
delph   0.01078 0.02556 0.42180 0.02157 0.00000  0.969 0.0001 ***
acrid   0.00982 0.02960 0.33170 0.01964 0.00000  0.981 0.0001 ***
taban   0.00753 0.02270 0.33170 0.01505 0.00000  0.991 0.0001 ***
cosmop  0.00753 0.02270 0.33170 0.01505 0.00000  1.000 0.0001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 10000

Generate graphs

Note: graphs may not run at the moment because of a problem in vegan 2.6.2 relating to the scores function when dealing with mds objects. Expecting a fix soon (as of 7 June 2022). Problem did not exist prior to this version of vegan

br=c("ALAN","no ALAN")
la=c("ALAN", "No ALAN")
a<-as.data.frame(scores(park1.mds))
a<-cbind(park[c(1:2)],a)   #Add mesocosms & treatments from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, shape=treat) )+
  geom_point()+
  labs(y="MDS2", x="MDS1", title="Raw")+
  scale_shape_manual(values=sym2,
                     name="Treat",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )+
  theme_qk()
as<-as.data.frame(scores(park1s.mds))
as<-cbind(park[c(1:2)],as)   #Add mesocosms & treatments from original data file

p2<-ggplot(data=as, aes(x=NMDS1, y=NMDS2, shape=treat) )+
  geom_point()+
  labs(y=NULL, x="MDS1", title="Standardized")+
  scale_shape_manual(values=sym2,
                     name="Treat",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()

p4<-p1+p2+plot_layout(guides='collect')&theme(axis.text=element_blank())

p4

# with ellipses
p1a<-p1+stat_ellipse(geom="polygon", alpha= 0.2)
p2a<-p2+stat_ellipse(geom="polygon", alpha= 0.2)
p4a<-p1a+p2a+plot_layout(guides='collect')

p4a

Color versions

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=treat, fill=treat) )+
  geom_point()+
  labs(y="MDS2", x="MDS1", title="Raw")+
  scale_color_uchicago(name="Treat",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )+
  scale_fill_uchicago(name="Treat",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL))+
  theme_qk()

p2<-ggplot(data=as, aes(x=NMDS1, y=NMDS2, color=treat, fill=treat) )+
  geom_point()+
  labs(y=NULL, x="MDS1", title="Standardized")+
  scale_color_uchicago(
                     name="Treat",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  scale_fill_uchicago(name="Treat",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL))+
  theme_qk()

p4<-p1+p2+plot_layout(guides='collect')&theme(axis.text=element_blank())

p4

# with ellipses
p1a<-p1+stat_ellipse(geom="polygon", alpha= 0.2)
p2a<-p2+stat_ellipse(geom="polygon", alpha= 0.2)
p4a<-p1a+p2a+plot_layout(guides='collect')&theme(axis.text=element_blank())

p4a