title: “Q & K Box 16.1” output: html_notebook —

Worked example of PCoA: invertebrates in artificial ponds This example continues with the pond invertebrate example used in Boxes 15.2 and 15.3 and uses the data file lemminvert.csv.

Preliminaries

Load core set of packages and vegan, ape

source("../R/libraries.R")
library(vegan)
library(ape)

Read data

lemminvert <- read.csv("../data/lemminvert.csv")

Do PCoA on lemmens invert data

Use B-C matrix on raw abundances

# Remove pond classification variables
lemminvert.bc <- vegdist(lemminvert[,4:32],method="bray")
# do PCoA
lemminvert.pco <- pcoa(lemminvert.bc)
print(lemminvert.pco)
$correction
[1] "none" "1"   

$note
[1] "No correction was applied to the negative eigenvalues"

$values
       Eigenvalues    Relative_eig   Rel_corr_eig   Broken_stick Cum_corr_eig Cumul_br_stick
1   1.425605593977  0.263894086652 0.216855943873 0.173588509751 0.2168559439   0.1735885098
2   1.266281010728  0.234401486766 0.193714195654 0.125969462132 0.4105701395   0.2995579719
3   0.686166627796  0.127016417653 0.109453118588 0.102159938322 0.5200232581   0.4017179102
4   0.454073168665  0.084053559150 0.075741758515 0.086286922449 0.5957650166   0.4880048327
5   0.358117768704  0.066291239236 0.061804325503 0.074382160544 0.6575693421   0.5623869932
6   0.297913465130  0.055146810670 0.053059706144 0.064858351020 0.7106290483   0.6272453442
7   0.249986262908  0.046274998362 0.046098324391 0.056921843084 0.7567273727   0.6841671873
8   0.200085412339  0.037037843682 0.038850272043 0.050119121996 0.7955776447   0.7342863093
9   0.157519751844  0.029158507246 0.032667649262 0.044166741043 0.8282452940   0.7784530503
10  0.127511876510  0.023603744494 0.028309033150 0.038875735752 0.8565543271   0.8173287861
11  0.102437725035  0.018962264178 0.024667035862 0.034113830990 0.8812213630   0.8514426171
12  0.086648234474  0.016039468976 0.022373626975 0.029784826661 0.9035949900   0.8812274437
13  0.052157452222  0.009654874584 0.017363872784 0.025816572693 0.9209588627   0.9070440164
14  0.043769562909  0.008102191010 0.016145539626 0.022153569030 0.9371044024   0.9291975855
15  0.033994869408  0.006292795883 0.014725774444 0.018752208486 0.9518301768   0.9479497940
16  0.024405460446  0.004517698809 0.013332921670 0.015577605311 0.9651630985   0.9635273993
17  0.018274390907  0.003382775515 0.012442389496 0.012601414835 0.9776054880   0.9761288141
18  0.000000000000  0.000000000000 0.008610318709 0.009800294387 0.9862158067   0.9859291085
19 -0.008108364421 -0.001500940675 0.006387395720 0.007154791741 0.9926032024   0.9930839002
20 -0.023412578925 -0.004333906347 0.005202865160 0.004648526077 0.9978060676   0.9977324263
21 -0.031567746639 -0.005843510787 0.002193932429 0.002267573696 1.0000000000   1.0000000000
22 -0.052283422430 -0.009678192950 0.000000000000 0.000000000000 1.0000000000   1.0000000000
23 -0.067388044891 -0.012474212105 0.000000000000 0.000000000000 1.0000000000   1.0000000000

$vectors
           Axis.1         Axis.2          Axis.3           Axis.4          Axis.5         Axis.6          Axis.7          Axis.8          Axis.9         Axis.10          Axis.11
1  -0.27508949873 -0.05212328360  0.014015261009  0.1089673486792  0.035340804174  0.04824095155 -0.014069004387  0.292721869963  0.093090742756 -0.016760875159 -0.0771717188901
2  -0.09225316967 -0.27722834929  0.061415789541  0.0473687348225 -0.014826832565 -0.10030236281  0.039336944245 -0.050133809209 -0.013650266348  0.030936268285  0.0841374212898
3   0.04497220159 -0.31710565088 -0.084121681319 -0.0774249988121  0.001257019948 -0.00189993911 -0.001147952734  0.013503447194  0.012011736860  0.006199503153 -0.0379881711385
4  -0.07478040177  0.24794033551  0.319389075209 -0.1006848698509  0.031128948998 -0.23107137032  0.017170290126  0.022282248093 -0.005392612792  0.042208335824  0.0680942813806
5   0.31007361349  0.08830565568  0.021372311118 -0.0003292858567 -0.129666829527  0.08401290834 -0.197277645404 -0.208081746479  0.142381544796  0.009046744912 -0.0191102746926
6  -0.32990486973 -0.06914779567 -0.199921267176 -0.1273830032493 -0.232877701405  0.05954403021  0.127489871590 -0.073283125535 -0.025272096839  0.010578503228 -0.0166339515141
7  -0.44598071727  0.29070151329 -0.274800125698  0.0895717959311  0.192311510708  0.12642290946 -0.052232268283 -0.057558285243  0.035020854632  0.151701029465  0.0679684067238
8  -0.43179013888  0.42182479954 -0.039661088651 -0.0239355519179  0.054820950193  0.10144271665 -0.095990226546 -0.012643732944 -0.034131211524 -0.085917379936 -0.0409222323760
9  -0.28813448977 -0.13905574476 -0.003969375199  0.1489252052391  0.056658377557 -0.13799278633  0.166888623105 -0.072845742333  0.059205424439 -0.023137704395  0.0167152262695
10 -0.10720543595  0.17517697996  0.125995564097 -0.2939766412196  0.171589947791 -0.20075787163 -0.110118066006 -0.058258007961 -0.076564633497 -0.025259139660 -0.0824635294075
11 -0.15999745080  0.18903652281  0.094658999259 -0.2286412221754 -0.331423881925  0.08502821640  0.056741781296  0.069252398508  0.026809553047 -0.072983418776  0.0731475047807
12  0.31619846568  0.34300840034 -0.379288368019  0.2345505531504 -0.205494740240 -0.26142606460 -0.081116519509  0.068074632278 -0.068960843274  0.014359483696 -0.0304303075437
13  0.16302805891 -0.12061021666  0.157297592675  0.0801902280596 -0.036058548523  0.10929101893 -0.113442354380  0.010290075316 -0.087205322781  0.127911781954  0.0720795336929
14  0.02245818427 -0.25678113517 -0.013446781631  0.0039316964956  0.024786725307  0.05418826419  0.056792833546 -0.101069363236 -0.066314475601 -0.009021090850 -0.1611182860690
15  0.05266626691  0.19087431823  0.340790209886  0.2939955826411 -0.001252034785  0.01002821687  0.111603394481 -0.080981931773  0.073957974889 -0.037058692027 -0.0297391084593
16  0.36841006387  0.21297423012 -0.001289891340 -0.0120565163727  0.094633825400  0.17163352729  0.136319300543  0.015080513545 -0.234460895500 -0.034937950432  0.0428994073724
17  0.41642243507  0.11965964413 -0.140501002168 -0.0672173936658  0.173650106425  0.03719521799  0.040264493211  0.008953128593  0.149334015750 -0.148583650503  0.0940098057162
18 -0.08282681422 -0.17334241330  0.150216458267  0.2376349553081 -0.023016186857  0.06416433425 -0.096221881798  0.006200815476 -0.081833642670 -0.119123141298 -0.0216478876730
19  0.03112215703 -0.23821347506  0.098020719035  0.0044359546117  0.017430630881 -0.01197600800 -0.072103659890  0.099638238257 -0.021250395252  0.037799817469  0.0791786453877
20  0.38304152033  0.17693621408 -0.037443612040 -0.0955406276020  0.064253824025  0.03419850918  0.197679390536  0.069185519034  0.073377090669  0.118659773979 -0.0880331530917
21 -0.05236358691 -0.37639442872 -0.230788679681 -0.0266194377851  0.041032275117 -0.08766014130  0.070288487425 -0.059778356269 -0.018101052674 -0.031542237915  0.0574310241043
22  0.12617652088 -0.30377158448 -0.128618034598 -0.1041302858626  0.085723592360  0.01250784500 -0.149846702455  0.065171266229  0.019202232509 -0.055732270775 -0.0001278653959
23  0.10575708569 -0.13266453607  0.150677927422 -0.0916322205682 -0.070001783055  0.03518787777 -0.037009128711  0.034279948497  0.048746278406  0.110656309760 -0.0502747704666
           Axis.12        Axis.13          Axis.14          Axis.15          Axis.16         Axis.17
1   0.077812751475 -0.03187609406  0.0747993853439  0.0229458889079  3.191437526e-02 -0.023395733393
2  -0.055389020205 -0.02693409203  0.0319426406035  0.0320850008010  6.107752417e-02 -0.008850538502
3   0.056771671469 -0.02083562147 -0.0596874914358  0.0166059756231  2.617688718e-02  0.084685887761
4   0.007603277406 -0.09018362834  0.0008408010569 -0.0427540870683  1.366959325e-02  0.006507392414
5   0.113171638751 -0.02733883261  0.0543010296619  0.0272783561690 -5.902123959e-03 -0.001998963742
6  -0.069720937732 -0.10382147589  0.0462979719501 -0.0251310005164 -3.766980936e-02 -0.008545967425
7   0.031038310577 -0.01817899508 -0.0777996506194  0.0002150316207  2.171755581e-02 -0.015933656584
8  -0.085803776447  0.04581164872  0.0506505185248 -0.0011925694844 -2.516961373e-02  0.043161774728
9   0.086763162142  0.08664456446  0.0248674574887 -0.0203362959273 -5.266172427e-02 -0.009495218628
10  0.025159387676  0.02701204871  0.0061483979561  0.0276462124724  4.746237933e-03 -0.010641406097
11  0.029825065777  0.07634609000 -0.0641719338721  0.0248279835189  3.365931620e-02 -0.011706550393
12 -0.014805107052  0.01006553732 -0.0166592573395  0.0033160138076 -6.011788656e-05 -0.006035505526
13 -0.074968839510  0.06994219518  0.0752736672419 -0.0110965392231  2.562162053e-02  0.013117182227
14 -0.047986408422 -0.01215291090 -0.0300863187059  0.0469379915693  2.901001748e-02 -0.041566765294
15 -0.069422763715  0.01646501394 -0.0416294147041  0.0165421357061  5.752265878e-03 -0.011251749630
16  0.109991608182 -0.01765826852  0.0169448267595 -0.0033367662936 -8.056615254e-03 -0.021018230875
17 -0.058997154393 -0.03445082756  0.0083719005343 -0.0064022798770  4.719126563e-03 -0.005503134976
18  0.024205741116 -0.04575699834 -0.0435289098724 -0.0564733030486 -2.583076013e-03  0.034924191941
19 -0.025491404524 -0.02418808296 -0.0327928713180  0.1028824232478 -9.550496510e-02  0.001061124666
20 -0.046846550385  0.02203677020  0.0056706103703 -0.0009173536798 -7.400224395e-03  0.040838534494
21  0.033080221105  0.05027730005  0.0336015217004 -0.0115034884680  1.736720210e-02  0.015715475865
22 -0.073378980009  0.02858229553 -0.0354879874747 -0.0553896798084 -6.481545128e-03 -0.041454110064
23  0.027388106717  0.02019236364 -0.0278668938503 -0.0867496500490 -3.394190726e-02 -0.022614032964

$trace
[1] 5.402188477

attr(,"class")
[1] "pcoa"
biplot(lemminvert.pco,lab=lemminvert$manag)

Improve plot

Plots used for QK use the ggplot classic theme, with some tweaks. Tweaks are consolidated into theme_QK; use this theme for figures and tweak the theme to avoid repetitive code changes.

source("../R/appearance.R")

Get labels back on

# extract coordinates
a<-as.data.frame(lemminvert.pco$vectors)
# Need to get management types on this plot
a<-cbind(lemminvert[c(1:3)],a)   #Add site names & symbols from original data file
br=c("nm","li","nf","yf")     #Categories from original file; define as object to avoid retyping
la=c("None", "Light", "No fish", "Young fish")   #Labels; define here to avoid retyping
p2<-ggplot(data=a, aes(x=Axis.1, y=Axis.2, shape=manag, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym4,
                     name="Management",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()+
  theme(legend.position = "top")
p2

# colour version using viridis palette d
p2a<-ggplot(data=a, aes(x=Axis.1, y=Axis.2, color=manag, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_color_viridis_d(
                     name="Management",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()+
  theme(legend.position = "top")
p2a

LS0tCnRpdGxlOiAiUUsgQm94IDE2LjEiCm91dHB1dDoKICBodG1sX25vdGVib29rCi0tLQoKdGl0bGU6ICJRICYgSyBCb3ggMTYuMSIgb3V0cHV0OiBodG1sX25vdGVib29rIC0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKYGBgCgpXb3JrZWQgZXhhbXBsZSBvZiBQQ29BOiBpbnZlcnRlYnJhdGVzIGluIGFydGlmaWNpYWwgcG9uZHMgVGhpcyBleGFtcGxlIGNvbnRpbnVlcyB3aXRoIHRoZSBwb25kIGludmVydGVicmF0ZSBleGFtcGxlIHVzZWQgaW4gQm94ZXMgWzE1LjJdKGxlbW1lbnNDQS5uYi5odG1sKSBhbmQgWzE1LjNdKGxlbW1lbnNSREEubmIuaHRtbCkgYW5kIHVzZXMgdGhlIGRhdGEgZmlsZSBbbGVtbWludmVydC5jc3ZdKC4uL2RhdGEvbGVtbWludmVydC5jc3YpLgoKIyMjIFByZWxpbWluYXJpZXMKCkxvYWQgY29yZSBzZXQgb2YgcGFja2FnZXMgYW5kIHZlZ2FuLCBhcGUKCmBgYHtyIHJlc3VsdHM9J2hpZGUnfQpzb3VyY2UoIi4uL1IvbGlicmFyaWVzLlIiKQpsaWJyYXJ5KHZlZ2FuKQpsaWJyYXJ5KGFwZSkKYGBgCgpSZWFkIGRhdGEKCmBgYHtyfQpsZW1taW52ZXJ0IDwtIHJlYWQuY3N2KCIuLi9kYXRhL2xlbW1pbnZlcnQuY3N2IikKYGBgCgojIyBEbyBQQ29BIG9uIGxlbW1lbnMgaW52ZXJ0IGRhdGEKClVzZSBCLUMgbWF0cml4IG9uIHJhdyBhYnVuZGFuY2VzCgpgYGB7cn0KIyBSZW1vdmUgcG9uZCBjbGFzc2lmaWNhdGlvbiB2YXJpYWJsZXMKbGVtbWludmVydC5iYyA8LSB2ZWdkaXN0KGxlbW1pbnZlcnRbLDQ6MzJdLG1ldGhvZD0iYnJheSIpCiMgZG8gUENvQQpsZW1taW52ZXJ0LnBjbyA8LSBwY29hKGxlbW1pbnZlcnQuYmMpCnByaW50KGxlbW1pbnZlcnQucGNvKQpiaXBsb3QobGVtbWludmVydC5wY28sbGFiPWxlbW1pbnZlcnQkbWFuYWcpCmBgYAoKIyMjIEltcHJvdmUgcGxvdAoKUGxvdHMgdXNlZCBmb3IgUUsgdXNlIHRoZSBnZ3Bsb3QgY2xhc3NpYyB0aGVtZSwgd2l0aCBzb21lIHR3ZWFrcy4gVHdlYWtzIGFyZSBjb25zb2xpZGF0ZWQgaW50byB0aGVtZV9RSzsgdXNlIHRoaXMgdGhlbWUgZm9yIGZpZ3VyZXMgYW5kIHR3ZWFrIHRoZSB0aGVtZSB0byBhdm9pZCByZXBldGl0aXZlIGNvZGUgY2hhbmdlcy4KCmBgYHtyfQpzb3VyY2UoIi4uL1IvYXBwZWFyYW5jZS5SIikKYGBgCgpHZXQgbGFiZWxzIGJhY2sgb24KCmBgYHtyfQojIGV4dHJhY3QgY29vcmRpbmF0ZXMKYTwtYXMuZGF0YS5mcmFtZShsZW1taW52ZXJ0LnBjbyR2ZWN0b3JzKQojIE5lZWQgdG8gZ2V0IG1hbmFnZW1lbnQgdHlwZXMgb24gdGhpcyBwbG90CmE8LWNiaW5kKGxlbW1pbnZlcnRbYygxOjMpXSxhKSAgICNBZGQgc2l0ZSBuYW1lcyAmIHN5bWJvbHMgZnJvbSBvcmlnaW5hbCBkYXRhIGZpbGUKYnI9Yygibm0iLCJsaSIsIm5mIiwieWYiKSAgICAgI0NhdGVnb3JpZXMgZnJvbSBvcmlnaW5hbCBmaWxlOyBkZWZpbmUgYXMgb2JqZWN0IHRvIGF2b2lkIHJldHlwaW5nCmxhPWMoIk5vbmUiLCAiTGlnaHQiLCAiTm8gZmlzaCIsICJZb3VuZyBmaXNoIikgICAjTGFiZWxzOyBkZWZpbmUgaGVyZSB0byBhdm9pZCByZXR5cGluZwpwMjwtZ2dwbG90KGRhdGE9YSwgYWVzKHg9QXhpcy4xLCB5PUF4aXMuMiwgc2hhcGU9bWFuYWcsICkgKSsKICBnZW9tX3BvaW50KCkrCiAgbGFicyh5PSJNRFMyIiwgeD0iTURTMSIpKwogIHNjYWxlX3NoYXBlX21hbnVhbCh2YWx1ZXM9c3ltNCwKICAgICAgICAgICAgICAgICAgICAgbmFtZT0iTWFuYWdlbWVudCIsCiAgICAgICAgICAgICAgICAgICAgIGJyZWFrcz1iciwKICAgICAgICAgICAgICAgICAgICAgbGFiZWxzPWxhLAogICAgICAgICAgICAgICAgICAgICBndWlkZSA9CiAgICAgICAgICAgICAgICAgICAgICAgICBndWlkZV9sZWdlbmQobGFiZWwudGhlbWUgPSBlbGVtZW50X3RleHQoc2l6ZT02KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGl0bGU9TlVMTCkKICAgICAgICAgICAgICAgICAgICAgKSsKICB0aGVtZV9xaygpKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiKQpwMgpgYGAKCmBgYHtyfQojIGNvbG91ciB2ZXJzaW9uIHVzaW5nIHZpcmlkaXMgcGFsZXR0ZSBkCnAyYTwtZ2dwbG90KGRhdGE9YSwgYWVzKHg9QXhpcy4xLCB5PUF4aXMuMiwgY29sb3I9bWFuYWcsICkgKSsKICBnZW9tX3BvaW50KCkrCiAgbGFicyh5PSJNRFMyIiwgeD0iTURTMSIpKwogIHNjYWxlX2NvbG9yX3ZpcmlkaXNfZCgKICAgICAgICAgICAgICAgICAgICAgbmFtZT0iTWFuYWdlbWVudCIsCiAgICAgICAgICAgICAgICAgICAgIGJyZWFrcz1iciwKICAgICAgICAgICAgICAgICAgICAgbGFiZWxzPWxhLAogICAgICAgICAgICAgICAgICAgICBndWlkZSA9CiAgICAgICAgICAgICAgICAgICAgICAgICBndWlkZV9sZWdlbmQobGFiZWwudGhlbWUgPSBlbGVtZW50X3RleHQoc2l6ZT02KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGl0bGU9TlVMTCkKICAgICAgICAgICAgICAgICAgICAgKSsKICB0aGVtZV9xaygpKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiKQpwMmEKYGBgCg==