Menu
I've plotted centroids for a Gower similarity index using the vegan{} package in R and would like to color code the ellipse fills based on factors in my original dataset. The centroids are for Month-Site combinations (and these are the only centroids I want to plot; see data sample below), but I want to color code the ellipses by Month and then, in a second plot, by site, to show differences between factors. At the moment the code looks like this:
dat.fac<-dat[,1:7] ##factors
dat.mgower <- vegdist(decostand(dat[,8:34], "log"), "altGower")
ord<-capscale(dat.mgower~1) #Does the PCO
plot(ord, display = "sites", type = "n")
#for i in
stats<-with(dat.fac,
ordiellipse(ord, Site_TP, kind="se", conf=0.95,
lwd=2, draw = "polygon",col="#0000ff22",border = "blue"))
The MDS plot with the centroids looks like this:
I would like to be able to color code the fills of the centroids based on a) Site and b) season (May vs. September). I have tried using col=c("#0000ff22","#CAFF7022",...) but this doesn't work because the coordinates for the ellipses are in ord. I've also tried
with(dat.fac,ordiellipse(ord, Site_TP[1:54], kind="se", conf=0.95,
lwd=2, draw = "polygon",col="#0000ff22",border = "blue"))
to get the first half of the ceentroids to be one color (the first half are all from the same month), but this results in the error "Error in cov.wt(X, W) : 'x' must contain finite values only." Help! Thanks.
The data looks like this:
structure(list(Core_num = structure(c(39L, 50L, 61L, 70L, 71L,
72L, 73L, 74L, 75L, 76L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L,
93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L,
105L, 106L, 107L, 108L, 112L, 113L, 115L, 116L, 117L, 118L, 119L,
120L, 121L, 122L, 123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L,
131L, 132L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 17L, 18L,
19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L,
32L, 33L, 34L, 35L, 36L, 37L, 38L, 40L, 44L, 45L, 46L, 47L, 48L,
49L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 62L, 63L,
64L, 65L, 66L), .Label = c("152", "153", "154", "155", "156",
"157", "158", "159", "160", "161", "162", "163", "164", "165",
"166", "167", "168", "169", "170", "171", "172", "173", "174",
"175", "176", "177", "178", "179", "180", "181", "182", "183",
"184", "185", "186", "187", "188", "189", "19", "190", "191",
"192", "193", "194", "195", "196", "197", "198", "199", "20",
"200", "201", "202", "203", "204", "205", "206", "207", "208",
"209", "21", "210", "211", "212", "213", "214", "215", "216",
"217", "22", "23", "24", "25", "26", "27", "28", "29", "30",
"31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41",
"42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52",
"53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63",
"64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74",
"75", "76", "77", "78", "79", "80", "81", "82", "83", "84", "L1",
"L10", "L11", "L12", "L13", "L14", "L15", "L16", "L17", "L18",
"L2", "L20", "L21", "L22", "L23", "L24", "L25", "L26", "L27",
"L28", "L3", "L4", "L5", "L6", "L7", "L9"), class = "factor"),
Cage_num = c(47L, 46L, 45L, 2L, 24L, 1L, 23L, 25L, 3L, 61L,
5L, 6L, 27L, 49L, 4L, 7L, 48L, 29L, 28L, 26L, 8L, 10L, 51L,
50L, 9L, 11L, 30L, 33L, 31L, 32L, 52L, 53L, 54L, 55L, 18L,
60L, 41L, 38L, 63L, 19L, 16L, 17L, 40L, 44L, 43L, 21L, 20L,
22L, 39L, 62L, 65L, 66L, 42L, 64L, 9L, 54L, 31L, 52L, 10L,
53L, 33L, 32L, 8L, 11L, 55L, 30L, 66L, 64L, 21L, 20L, 43L,
42L, 44L, 65L, 63L, 18L, 16L, 22L, 62L, 61L, 38L, 19L, 60L,
17L, 39L, 40L, 41L, 3L, 24L, 25L, 2L, 23L, 1L, 46L, 47L,
45L, 5L, 28L, 48L, 6L, 51L, 50L, 4L, 7L, 29L, 49L, 26L, 27L
), Treatment = structure(c(3L, 3L, 3L, 1L, 2L, 1L, 2L, 2L,
1L, 3L, 1L, 1L, 2L, 3L, 1L, 1L, 3L, 2L, 2L, 2L, 1L, 1L, 3L,
3L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 3L, 2L, 2L,
3L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 3L, 3L, 3L, 2L,
3L, 1L, 3L, 2L, 3L, 1L, 3L, 2L, 2L, 1L, 1L, 3L, 2L, 3L, 3L,
1L, 1L, 2L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 3L, 3L, 2L, 1L, 3L,
1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 3L, 3L, 3L, 1L, 2L,
3L, 1L, 3L, 3L, 1L, 1L, 2L, 3L, 2L, 2L), .Label = c("C",
"CC", "NC"), class = "factor"), Site = structure(c(2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 7L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 8L, 8L, 1L, 1L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 6L, 6L, 6L, 6L,
6L, 7L, 7L, 6L, 6L, 6L, 6L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L,
6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L), .Label = c("M1", "M2", "R07", "R08_Feb", "R08_Jul",
"R09", "R10", "R11", "R11_5m"), class = "factor"), Days_Since_Removal = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 291, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 291, 291, 291,
291, 291, 291, 291, 291, 291, 669, 669, 669, 669, 669, 291,
291, 669, 669, 669, 669, 86, 86, 86, 86, 86, 86, 86, 86,
86, 86, 86, 86, 806, 806, 806, 806, 806, 806, 806, 806, 428,
428, 428, 806, 428, 428, 428, 428, 428, 428, 428, 428, 428,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0), Time_point = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("May",
"September"), class = "factor"), Site_TP = structure(c(3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 10L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 14L, 14L, 1L, 1L, 14L, 14L, 14L, 14L,
14L, 14L, 14L, 14L, 14L, 14L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 8L, 8L, 8L, 8L, 8L, 10L, 10L, 8L, 8L, 8L,
8L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
15L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 11L, 11L, 11L, 9L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("M1_May", "M1_September", "M2_May", "M2_September",
"R07_May", "R08_Feb_May", "R08_Jul_May", "R09_May", "R09_September",
"R10_May", "R10_September", "R11_5m_May", "R11_5m_September",
"R11_May", "R11_September"), class = "factor"), Spionidae = c(0,
2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 19, 7, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Nereididae = c(0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Syllidae = c(0, 0, 0, 0,
0, 0, 0, 0, 0, 2, 8, 0, 0, 0, 3, 1, 0, 0, 0, 4, 4, 0, 0,
0, 1, 0, 1, 4, 1, 0, 0, 1, 1, 2, 2, 0, 1, 0, 3, 0, 0, 0,
3, 1, 22, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 5, 3, 0, 0, 0, 0, 9, 5, 0, 0,
0, 0, 29, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0), Opheliidae = c(0, 2, 1, 1, 3,
0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 14, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 3, 0, 14, 1, 0, 1, 0, 0, 1, 1, 0, 2,
2, 3, 5, 0, 1, 0, 0, 10, 0, 4, 7, 0, 2, 1, 5, 0, 0, 3, 3,
0, 1, 5, 0, 8, 11, 0, 0, 13, 9, 5, 4, 6, 1, 19, 3, 2, 5,
1, 1, 0, 0, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1, 1), Cossuridae = c(6, 29, 8, 7, 5,
0, 10, 31, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 4, 7, 6, 12, 0, 0, 19, 3, 4, 2, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0), Sternaspidae = c(0, 0, 0, 0, 0,
0, 0, 0, 6, 0, 0, 11, 5, 0, 1, 9, 1, 0, 7, 0, 0, 0, 13, 13,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 10, 0, 0, 0, 24, 1, 27,
5, 31, 5, 7, 4, 0, 2, 17), Sabellariidae = c(0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0), Sabellidae = c(1, 2, 1, 2, 0,
0, 3, 0, 2, 3, 2, 0, 0, 0, 0, 0, 2, 1, 0, 1, 0, 6, 0, 0,
0, 0, 0, 0, 0, 0, 0, 2, 0, 3, 0, 2, 0, 0, 0, 2, 0, 0, 0,
4, 16, 29, 0, 1, 0, 3, 2, 1, 7, 3, 24, 1, 1, 11, 1, 2, 1,
50, 0, 18, 27, 0, 0, 6, 1, 1, 43, 0, 4, 2, 1, 4, 153, 119,
2, 2, 2, 0, 0, 0, 0, 3, 0, 0, 4, 0, 0, 0, 0, 1, 0, 0, 0,
0, 4, 0, 0, 0, 2, 2, 1, 0, 0, 0), Capitellidae = c(1, 3,
3, 10, 7, 1, 10, 6, 17, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 3,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0, 0, 4, 0, 0, 0, 0,
0, 1, 0, 0, 0, 4, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 3, 1, 3,
1, 3, 14, 1, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0,
0, 0, 6, 6, 2, 1, 0, 4, 2, 2, 1, 0, 1, 0, 0, 4, 0, 0, 1,
1, 1, 0, 0, 0, 0, 0, 0, 0, 4, 1, 0), Dorvillidae = c(0, 0,
0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0,
2, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,
0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Cirratulidae = c(0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Oligochaeta = c(19, 35,
14, 41, 31, 12, 42, 38, 27, 16, 28, 2, 3, 2, 32, 0, 12, 12,
2, 7, 2, 4, 1, 20, 9, 14, 35, 5, 1, 0, 1, 4, 0, 6, 264, 48,
27, 8, 22, 0, 3, 133, 42, 5, 40, 21, 0, 2, 101, 4, 2, 0,
0, 11, 20, 4, 21, 22, 1, 9, 64, 5, 0, 38, 8, 20, 6, 5, 5,
35, 8, 7, 13, 9, 93, 20, 25, 12, 28, 63, 21, 71, 271, 4,
204, 8, 51, 14, 11, 10, 10, 1, 5, 4, 6, 4, 2, 3, 1, 0, 0,
1, 0, 1, 2, 9, 24, 1), Nemertea = c(0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0), Nematoda = c(2, 0, 0, 4, 5, 5, 0, 1, 13,
0, 0, 0, 1, 0, 0, 2, 0, 8, 3, 3, 2, 1, 4, 4, 0, 0, 18, 0,
1, 1, 1, 8, 2, 12, 1, 0, 4, 0, 1, 0, 0, 6, 0, 1, 4, 0, 0,
0, 1, 0, 0, 2, 6, 0, 0, 0, 1, 15, 0, 0, 20, 1, 0, 1, 2, 0,
0, 0, 0, 1, 16, 5, 1, 1, 0, 0, 35, 3, 1, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 1, 0, 1, 16, 2, 1, 0, 18, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0), Platyhelminthes = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0), Sipuncula = c(0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 25, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0,
0, 0, 0, 0), Ostracoda = c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, 5,
7, 0, 2, 11, 5, 5, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0), Decapoda = c(0, 0, 0, 4, 0, 0, 2, 0, 0, 0, 0, 0,
1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 1, 1, 0, 0, 4, 0, 1, 1, 0, 0, 1, 0, 2, 1, 0, 0, 1, 0,
0, 0, 17, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 2,
0), Amphipoda = c(1, 1, 5, 17, 8, 75, 22, 0, 34, 4, 38, 0,
27, 24, 89, 2, 1, 5, 3, 42, 88, 3, 1, 21, 6, 1, 6, 60, 6,
11, 0, 4, 11, 17, 15, 0, 13, 5, 44, 9, 0, 0, 0, 11, 20, 28,
8, 45, 0, 0, 18, 2, 14, 8, 3, 0, 2, 6, 6, 0, 2, 14, 2, 2,
6, 3, 37, 37, 5, 20, 63, 12, 10, 49, 6, 3, 150, 83, 7, 3,
5, 20, 15, 4, 19, 1, 6, 0, 3, 0, 0, 0, 14, 0, 1, 2, 3, 0,
0, 0, 1, 16, 3, 3, 0, 3, 1, 1), Isopoda = c(0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0), Copepoda = c(0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 2, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 7, 0, 0, 1, 0, 0, 0,
1, 0, 8, 1, 0, 2, 5, 1, 5, 2, 2, 0, 2, 4, 1, 0, 0, 0, 0,
0, 1, 1, 0, 0, 1, 0, 0, 0, 5, 3, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0), Tanaidacea = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 6, 0, 1, 0, 0, 0, 0, 3,
0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 2, 0, 1, 1, 6, 2, 0, 1, 0, 2, 1, 3, 1, 0, 0, 0, 0, 2,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0), HYRDOZOA = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0), ANTHOZOA = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0), Cnidarian_larvae = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0), Tunicata = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
Mollusca = c(0L, 0L, 0L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 2L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 14L, 2L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L,
0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("Core_num",
"Cage_num", "Treatment", "Site", "Days_Since_Removal", "Time_point",
"Site_TP", "Spionidae", "Nereididae", "Syllidae", "Opheliidae",
"Cossuridae", "Sternaspidae", "Sabellariidae", "Sabellidae",
"Capitellidae", "Dorvillidae", "Cirratulidae", "Oligochaeta",
"Nemertea", "Nematoda", "Platyhelminthes", "Sipuncula", "Ostracoda",
"Decapoda", "Amphipoda", "Isopoda", "Copepoda", "Tanaidacea",
"HYRDOZOA", "ANTHOZOA", "Cnidarian_larvae", "Tunicata", "Mollusca"
), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 19L,
20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L,
33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 46L, 47L, 49L,
50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L,
63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L,
76L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L,
95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 103L, 104L, 105L,
109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 118L, 119L,
120L, 121L, 122L, 123L, 124L, 125L, 126L, 127L, 128L, 129L), class = "data.frame")
973
To create a 95% confidence ellipse, we instead place . 95 in cell H9 and the formula =CHISQ. INV(H9,2) in cell H8 (resulting in a chi-square value of 5.99).
The following figure shows a 95% confidence ellipse for a set of 2D normally distributed data samples. This confidence ellipse defines the region that contains 95% of all samples that can be drawn from the underlying Gaussian distribution.
Ideally, the function could return the x and y coordinates used in the ellipses. Here's a simple modification of ordiellipse() that adds xy to each list element of the function output:
ordiellipse2 <-
function (ord, groups, display = "sites", kind = c("sd", "se"),
conf, draw = c("lines", "polygon", "none"),
w = weights(ord, display), col = NULL, alpha = 127,
show.groups, label = FALSE, ...)
{
weights.default <- function(object, ...) NULL
kind <- match.arg(kind)
draw <- match.arg(draw)
pts <- scores(ord, display = display, ...)
## ordiellipse only works with 2D data (2 columns)
pts <- as.matrix(pts)
if (ncol(pts) > 2)
pts <- pts[ , 1:2, drop = FALSE]
if (ncol(pts) < 2)
stop("ordiellipse needs two dimensions")
w <- eval(w)
if (length(w) == 1)
w <- rep(1, nrow(pts))
if (is.null(w))
w <- rep(1, nrow(pts))
## make semitransparent fill
if (draw == "polygon" && !is.null(col))
col <- rgb(t(col2rgb(col)), alpha = alpha, maxColorValue = 255)
if (!missing(show.groups)) {
take <- groups %in% show.groups
pts <- pts[take, , drop = FALSE]
groups <- groups[take]
w <- w[take]
}
out <- seq(along = groups)
inds <- names(table(groups))
res <- list()
if (label)
cntrs <- names <- NULL
## Remove NA scores
kk <- complete.cases(pts)
for (is in inds) {
gr <- out[groups == is & kk]
if (length(gr) > 2) {
X <- pts[gr, ]
W <- w[gr]
mat <- cov.wt(X, W)
if (kind == "se")
mat$cov <- mat$cov/mat$n.obs
if (missing(conf))
t <- 1
else t <- sqrt(qchisq(conf, 2))
xy <- vegan:::veganCovEllipse(mat$cov, mat$center, t)
if (draw == "lines")
vegan:::ordiArgAbsorber(xy, FUN = lines,
col = if(is.null(col)) par("fg") else col,
...)
else if (draw == "polygon")
vegan:::ordiArgAbsorber(xy[, 1], xy[, 2], col = col, FUN = polygon,
...)
if (label && draw != "none") {
cntrs <- rbind(cntrs, mat$center)
names <- c(names, is)
}
mat$scale <- t
mat$xy <- xy # this is the only substantive change I made
res[[is]] <- mat
}
}
if (label && draw != "none") {
if (draw == "lines")
vegan:::ordiArgAbsorber(cntrs[,1], cntrs[,2], labels=names, col = col,
FUN = text, ...)
else
vegan:::ordiArgAbsorber(cntrs, labels = names, col = NULL,
FUN = ordilabel, ...)
}
class(res) <- "ordiellipse"
invisible(res)
}
Now you can create your stats object, per your example, but specify draw = "none", because you'll want to play withstats` on your own to have full control of the ellipse plotting (no need to specify graphical arguments):
stats<-with(dat.fac, ordiellipse2(ord, Site_TP, kind = "se", conf = 0.95, draw = "none"))
notice now that stats has the same form, but there is an element xy for each ellipse:
names(stats[[1]])
[1] "cov" "center" "n.obs" "wt" "scale" "xy"
You can plot these in a loop using the polygon() function, but better would be to place all the xy coordinates in a single matrix or data.frame, so that you can just make a single call to polygon(), but with vectorized arguments, such as col, lwd. I'll extract the xy matrices using lapply(), then collapse them together as a data.frame. Key is the separate the individual ellipses with NAs, so the polygon drawing knows where one stops and the next starts:
xy.mat <- t(as.data.frame(lapply(stats, function(x){
t(rbind(x$xy, NA))
})))
Plotting becomes easy:
plot(ord, display = "sites", type = "n")
polygon(xy.mat[, 1], xy.mat[, 2], col = c(rep("#0000ff20", 5), rep("#ff000030", 5)), border = c(rep("blue", 5), rep("red", 5)))
Now it's just a matter of ordering the graphical arguments to polygon() correctly.
98
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
