Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Color-coding 95% confidence ellipses for centroids

Tags:

r

vegan

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: Centroids with Site_season titles

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")
like image 973
Margaret Avatar asked Aug 17 '12 01:08

Margaret


People also ask

How do you find the 95 confidence interval for an ellipse?

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).

What does a 95% confidence ellipse mean?

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.


1 Answers

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)))

enter image description here

Now it's just a matter of ordering the graphical arguments to polygon() correctly.

like image 98
tim riffe Avatar answered Oct 19 '22 22:10

tim riffe