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
devtools::source_url("https://raw.githubusercontent.com/mjkeough/mjkeough.github.io/refs/heads/main/R/libraries.R")   #This loads a set of packages used often

devtools::source_url("https://raw.githubusercontent.com/mjkeough/mjkeough.github.io/refs/heads/main/R/appearance.R")  #Graphics tweaks

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.0314781086 
Run 1 stress 0.03147916732 
... Procrustes: rmse 0.0009162213019  max resid 0.002573579273 
... Similar to previous best
Run 2 stress 0.03147824936 
... Procrustes: rmse 0.0006511648466  max resid 0.001823701035 
... Similar to previous best
Run 3 stress 0.03147826474 
... Procrustes: rmse 8.269200682e-05  max resid 0.0002333509795 
... Similar to previous best
Run 4 stress 0.03147804377 
... New best solution
... Procrustes: rmse 0.0005376042751  max resid 0.001505825555 
... Similar to previous best
Run 5 stress 0.04150746776 
Run 6 stress 0.04150779055 
Run 7 stress 0.04150747985 
Run 8 stress 0.04025041497 
Run 9 stress 0.04150747706 
Run 10 stress 0.04150755613 
Run 11 stress 0.03148315131 
... Procrustes: rmse 0.0008968505207  max resid 0.002517014572 
... Similar to previous best
Run 12 stress 0.0402503331 
Run 13 stress 0.03148097289 
... Procrustes: rmse 0.0007796402469  max resid 0.002183605502 
... Similar to previous best
Run 14 stress 0.03154474849 
... Procrustes: rmse 0.004097556535  max resid 0.01195077784 
Run 15 stress 0.03368096676 
Run 16 stress 0.03244801114 
Run 17 stress 0.0314781642 
... Procrustes: rmse 6.604590811e-05  max resid 0.0001860086218 
... Similar to previous best
Run 18 stress 0.03197404232 
... Procrustes: rmse 0.0120005642  max resid 0.03267575851 
Run 19 stress 0.04150829536 
Run 20 stress 0.03570233717 
*** Best solution repeated 4 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.03147804377 
Stress type 1, weak ties
Best solution was repeated 4 times in 20 tries
The best solution was from try 4 (random start)
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.1448504913 
Run 1 stress 0.1459783974 
Run 2 stress 0.1537681638 
Run 3 stress 0.1766276124 
Run 4 stress 0.1548081397 
Run 5 stress 0.1448504997 
... Procrustes: rmse 0.0007745907072  max resid 0.001822508025 
... Similar to previous best
Run 6 stress 0.1537680528 
Run 7 stress 0.1552734079 
Run 8 stress 0.1550377245 
Run 9 stress 0.145978259 
Run 10 stress 0.1448502487 
... New best solution
... Procrustes: rmse 0.0004035091964  max resid 0.0009831655366 
... Similar to previous best
Run 11 stress 0.1547491124 
Run 12 stress 0.1459782229 
Run 13 stress 0.1459785746 
Run 14 stress 0.1448503331 
... Procrustes: rmse 0.0002484488014  max resid 0.0006681234278 
... Similar to previous best
Run 15 stress 0.1459781267 
Run 16 stress 0.145978326 
Run 17 stress 0.1551463838 
Run 18 stress 0.1553990628 
Run 19 stress 0.1459783702 
Run 20 stress 0.14485031 
... Procrustes: rmse 0.0002054943791  max resid 0.0005974124936 
... Similar to previous best
Run 21 stress 0.1459782639 
Run 22 stress 0.1448503685 
... Procrustes: rmse 0.0002864081346  max resid 0.0007413505381 
... Similar to previous best
Run 23 stress 0.1459784582 
Run 24 stress 0.1547491536 
Run 25 stress 0.145978345 
Run 26 stress 0.1537684861 
Run 27 stress 0.145978299 
Run 28 stress 0.1459787331 
Run 29 stress 0.1459784516 
Run 30 stress 0.1459783164 
Run 31 stress 0.1459784266 
Run 32 stress 0.1552819078 
Run 33 stress 0.1548081105 
Run 34 stress 0.1459783393 
Run 35 stress 0.1459785045 
Run 36 stress 0.1658671891 
Run 37 stress 0.1459784109 
Run 38 stress 0.1763639878 
Run 39 stress 0.1547493436 
Run 40 stress 0.145978561 
*** Best solution repeated 4 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.1448502487 
Stress type 1, weak ties
Best solution was repeated 4 times in 40 tries
The best solution was from try 10 (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.01697682 0.016976819 0.45595 0.50811
Residuals 18 0.67021403 0.037234113                

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.022832467 0.0228324670 3.51759 0.077037 .
Residuals 18 0.116837023 0.0064909457                   
---
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.112 0.154 0.223 0.344 

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.2988889 
      Significance: 0.001 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0891 0.1150 0.1364 0.1568 

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
Permutation: free
Number of permutations: 999

adonis2(formula = park1.bc ~ treat, data = park, permutations = 999)
         Df  SumOfSqs         R2        F Pr(>F)    
Model     1 1.6916059 0.46132073 15.41506  0.001 ***
Residual 18 1.9752701 0.53867927                    
Total    19 3.6668759 1.00000000                    
---
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
Permutation: free
Number of permutations: 999

adonis2(formula = park1s.bc ~ treat, data = park, permutations = 999)
         Df  SumOfSqs         R2       F Pr(>F)    
Model     1 0.7250769 0.12405974 2.54935  0.001 ***
Residual 18 5.1195016 0.87594026                   
Total    19 5.8445785 1.00000000                   
---
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 3 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.6619    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.31607    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
caen      0.53073710   0.20750121   2.55775440 134.20000000  16.40000000 0.75123
chiron    0.08169960   0.04077522   2.00365870  22.50000000   5.00000000 0.86687
dolich    0.01473630   0.01465528   1.00552580   2.20000000   3.00000000 0.88773
cerat     0.00787050   0.01278722   0.61550080   1.00000000   0.20000000 0.89887
aphid     0.00753970   0.00928324   0.81218300   1.20000000   0.30000000 0.90954
ephyd     0.00664750   0.00851621   0.78057120   1.10000000   0.40000000 0.91895
baet      0.00622850   0.00539966   1.15350300   1.20000000   0.80000000 0.92777
chryso    0.00515730   0.00625262   0.82481640   1.00000000   0.40000000 0.93507
culicid   0.00447040   0.00600293   0.74471020   0.90000000   0.20000000 0.94140
tetrag    0.00440360   0.00853100   0.51618280   0.50000000   0.10000000 0.94763
hydro     0.00413960   0.00579859   0.71388990   0.50000000   0.30000000 0.95349
tipulid   0.00296020   0.00473538   0.62512810   0.10000000   0.40000000 0.95768
muscid    0.00275280   0.00628151   0.43823190   0.60000000   0.00000000 0.96158
empid     0.00192420   0.00480348   0.40058360   0.20000000   0.00000000 0.96430
formic    0.00187970   0.00359958   0.52220550   0.00000000   0.30000000 0.96696
polycen   0.00184710   0.00356411   0.51824890   0.00000000   0.30000000 0.96958
sciomyz   0.00178880   0.00456287   0.39203540   0.20000000   0.00000000 0.97211
helod     0.00173320   0.00276821   0.62610920   0.30000000   0.10000000 0.97456
taban     0.00149740   0.00457476   0.32731400   0.10000000   0.00000000 0.97668
cosmop    0.00149740   0.00457476   0.32731400   0.10000000   0.00000000 0.97880
coen      0.00144050   0.00288607   0.49911020   0.20000000   0.10000000 0.98084
delph     0.00140790   0.00324262   0.43419600   0.30000000   0.00000000 0.98283
notonec   0.00132110   0.00248828   0.53093330   0.20000000   0.10000000 0.98470
leptoc    0.00113500   0.00315211   0.36007570   0.10000000   0.10000000 0.98631
noctuid   0.00113060   0.00266564   0.42412650   0.00000000   0.20000000 0.98791
hydroph   0.00096370   0.00282052   0.34167350   0.10000000   0.10000000 0.98927
hydropt   0.00094090   0.00270596   0.34769980   0.10000000   0.10000000 0.99060
aran      0.00094090   0.00270596   0.34769980   0.10000000   0.10000000 0.99194
corix     0.00088120   0.00181616   0.48521560   0.20000000   0.00000000 0.99318
limneph   0.00082200   0.00167388   0.49105700   0.20000000   0.00000000 0.99435
phryg     0.00071650   0.00268957   0.26641010   0.00000000   0.10000000 0.99536
cicad     0.00071650   0.00268957   0.26641010   0.00000000   0.10000000 0.99637
agrom     0.00068070   0.00250987   0.27120330   0.00000000   0.10000000 0.99734
mymar     0.00068070   0.00250987   0.27120330   0.00000000   0.10000000 0.99830
halip     0.00062130   0.00222367   0.27938610   0.00000000   0.10000000 0.99918
acrid     0.00057870   0.00174866   0.33094120   0.10000000   0.00000000 1.00000
               p    
caen    0.000999 ***
chiron  0.003996 ** 
dolich  1.000000    
cerat   0.035964 *  
aphid   0.092907 .  
ephyd   0.995005    
baet    1.000000    
chryso  1.000000    
culicid 0.854146    
tetrag  0.157842    
hydro   0.994006    
tipulid 1.000000    
muscid  0.000999 ***
empid   0.000999 ***
formic  1.000000    
polycen 1.000000    
sciomyz 0.000999 ***
helod   0.990010    
taban   0.000999 ***
cosmop  0.000999 ***
coen    1.000000    
delph   0.000999 ***
notonec 0.997003    
leptoc  1.000000    
noctuid 1.000000    
hydroph 1.000000    
hydropt 1.000000    
aran    1.000000    
corix   0.000999 ***
limneph 0.000999 ***
phryg   1.000000    
cicad   1.000000    
agrom   1.000000    
mymar   1.000000    
halip   1.000000    
acrid   0.000999 ***
---
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.05567753 0.05036353 1.10551280 0.06114571 0.14459229 0.06928 0.0443956
tipulid 0.05434679 0.06560897 0.82834390 0.03054523 0.09694927 0.13691 0.1920808
formic  0.04838040 0.07859915 0.61553340 0.00000000 0.09676081 0.19712 0.0001000
caen    0.03615812 0.01885013 1.91818890 0.08925553 0.01838462 0.24211 0.0001000
baet    0.03615482 0.02235550 1.61726740 0.07487352 0.06517397 0.28710 0.2178782
helod   0.03376458 0.04640781 0.72756250 0.05221490 0.02575724 0.32912 0.3311669
polycen 0.03298760 0.05306306 0.62166790 0.00000000 0.06597519 0.37017 0.0001000
cerat   0.03128276 0.03456680 0.90499440 0.05570346 0.02631800 0.40910 0.3570643
hydroph 0.03085007 0.07579031 0.40704510 0.01318681 0.05115069 0.44749 0.4725527
coen    0.02716705 0.04624199 0.58749760 0.03912138 0.02293446 0.48129 0.4756524
ephyd   0.02608681 0.03148699 0.82849480 0.04707879 0.01636238 0.51375 0.3050695
aphid   0.02539183 0.02525769 1.00531100 0.05701335 0.01515914 0.54535 0.0308969
chiron  0.02315855 0.01568403 1.47656860 0.06636337 0.02706614 0.57417 0.0082992
tetrag  0.02268828 0.02822100 0.80395020 0.04104293 0.01191421 0.60240 0.1962804
chryso  0.02090415 0.01881287 1.11116230 0.03550069 0.02411942 0.62841 0.3288671
empid   0.02002624 0.04668232 0.42898980 0.04005248 0.00000000 0.65334 0.0001000
noctuid 0.01970658 0.04172432 0.47230440 0.00000000 0.03941316 0.67786 0.0001000
notonec 0.01917122 0.03335689 0.57473050 0.02914493 0.01387370 0.70171 0.4716528
leptoc  0.01892133 0.04045124 0.46775650 0.01751780 0.02382842 0.72526 0.4765523
hydro   0.01827123 0.02060112 0.88690470 0.02470708 0.02474372 0.74800 0.7012299
hydropt 0.01755325 0.03856413 0.45517030 0.01318681 0.02455704 0.76984 0.4646535
aran    0.01755325 0.03856413 0.45517030 0.01318681 0.02455704 0.79168 0.4646535
halip   0.01674495 0.05048793 0.33166250 0.00000000 0.03348990 0.81252 0.0001000
culicid 0.01439120 0.01842130 0.78122570 0.02561191 0.00792619 0.83043 0.1963804
sciomyz 0.01412053 0.02846080 0.49613970 0.02824107 0.00000000 0.84800 0.0001000
phryg   0.01328101 0.04004377 0.33166250 0.00000000 0.02656203 0.86453 0.0001000
cicad   0.01328101 0.04004377 0.33166250 0.00000000 0.02656203 0.88105 0.0001000
corix   0.01316086 0.02769928 0.47513370 0.02632172 0.00000000 0.89743 0.0001000
limneph 0.01259072 0.02623983 0.47983220 0.02518144 0.00000000 0.91310 0.0001000
agrom   0.01146723 0.03457499 0.33166250 0.00000000 0.02293446 0.92737 0.0001000
mymar   0.01146723 0.03457499 0.33166250 0.00000000 0.02293446 0.94164 0.0001000
muscid  0.01124399 0.02505917 0.44869770 0.02248799 0.00000000 0.95563 0.0001000
delph   0.01078405 0.02556459 0.42183530 0.02156810 0.00000000 0.96905 0.0001000
acrid   0.00981884 0.02960492 0.33166250 0.01963769 0.00000000 0.98127 0.0001000
taban   0.00752713 0.02269514 0.33166250 0.01505425 0.00000000 0.99063 0.0001000
cosmop  0.00752713 0.02269514 0.33166250 0.01505425 0.00000000 1.00000 0.0001000
           
dolich  *  
tipulid    
formic  ***
caen    ***
baet       
helod      
polycen ***
cerat      
hydroph    
coen       
ephyd      
aphid   *  
chiron  ** 
tetrag     
chryso     
empid   ***
noctuid ***
notonec    
leptoc     
hydro      
hydropt    
aran       
halip   ***
culicid    
sciomyz ***
phryg   ***
cicad   ***
corix   ***
limneph ***
agrom   ***
mymar   ***
muscid  ***
delph   ***
acrid   ***
taban   ***
cosmop  ***
---
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