show.groups argument in ordiellipse() to use different colours and line types for groups

Question

I am trying to produce a plot like the one shown at the bottom of this question: Adding ellipses to a principal component analysis (PCA) plot. However, I am neither familiar with ggplot nor with modifying the function ordiellipse() itself like shown here: Color-coding 95% confidence ellipses for centroids.

Therefore, I tried using the argument show.groups in ordiellipse() but I cannot figure out how to tell this argument what I want.

I am analyzing a table with 83 sites that contains abundance data for about 418 species. I run a DCA with the function decorana() and downweighted rare species on the data. Afterwards I am fitting mean Ellenberg indicator values onto the ordination.

vali <-read.table(file = "VegOrdi2012.txt", header=T)
attach(vali)

calcul.env<-read.table(file = "indicator-values-2012.txt",header=T)
attach(calcul.env)
names(calcul.env)

#Downweighting of rare species
calcul.dca1<-decorana(vali,iweigh=1)
calcul.dca1

#Fitting of environmental values
habi<-envfit(calcul.dca1~Habitat+Light+Temp+Continent+Moisture+pH+Nutrients+Salt,calcul.env,permu=999)
habi

It is the factor Habitat that I need to visualize in a diagram by means of ordiellipse() together with all the 83 plots shown as points and their plot number as text. Habitat consists of nine categories which comprise:

with(calcul.env, levels(Habitat))
[1] "Brown"  "Emb"    "Eph"    "FBS"    "GreeS"  "RiBa"   "Tracks" "TraffA"
[9] "UrbanL"

I was able to color and label them without drawing ellipses around them:

#Assigning colors to the nine categories
colvec <- c("black", "darkviolet", "green3", "darkorange", "magenta", "deepskyblue", "red", "blue", "forestgreen")

plot(calcul.dca1, type = "n", xlim=c(-1.5,1.6),ylim=c(-3,3))
ordipointlabel(calcul.dca1, display="site", cex=0.7, col = colvec[Habitat],bg = colvec[Habitat],pch = 21)

# Creating a legend
with(calcul.env, legend("topright", legend = levels(Habitat), cex=0.65,bty = "n",col = colvec, pch = 21, pt.bg = colvec))

But how do I tell the show.groups argument how to color the different Habitat types as well as assigning a different line type other than the standard one to the nine ellipses?

So far, I came up with this plot by using:

plot(calcul.dca, type="n", xlim=c(-1.5,1.6),ylim=c(-1,1))

ordiellipse(calcul.dca, Habitat, display="site", kind = "sd",
conf = 0.95, label = T, col="black", cex=0.7)

But it is rather elusive as it only shows the ellipses without the plots even though I added display="site". Also, because they're heavily overlapping it is hard to see which one is which. Therefore, it would be great to color each habitat type differently.

In case you need an extract of my DCA ordination data, here it comes:

dput(calcul.dca1)

structure(list(rproj = structure(c(0.622196100508291, 0.0425187062298239, 
0.809118643435233, 1.42649346881288, 0.942469303820393, 2.14549756552088, 
2.41922446503012, 1.28982244065239, 2.38544988614361, 1.40637644094225, 
1.91277135864581, 0.605410167146248, 0.806888393988314, 0.683853698855826, 
0.487565297706596, 1.33517822183407, 1.0239859712969, 1.11996240234112, 
0.878649568698446, 1.53912253203697, 1.10967929396299, 1.23838793138041, 
0.423890604035373, 0.624635178433215, 0.625774941613906, 0.873372685928585, 
1.02496140772642, 1.44785313041614, 1.08257665268759, 0.349688854561329, 
0.391547839427374, 0.730346245879122, 0.286646737460428, 0.731484625760109, 
1.08559070145785, 1.53243971551768, 0.85664914035342, 1.81525177285522, 
1.5519555131377, 1.39064843441504, 0.760013397658093, 1.22801365648165, 
2.21178067417998, 1.79390604164473, 1.65361126768449, 1.73966835197937, 
1.9237320812287, 1.80114956024378, 2.23293893445506, 2.02322760613128, 
...
0.857428410615221, 1.48594779863814, 0.935954361867925, 1.0524685508641, 
0.898391402052793), .Dim = c(84L, 4L), .Dimnames = list(c("p1", 
"p2", "p3", "p4", "p5", "p6", "p7", "p8", "p9", "p10", "p11", 
"p12", "p13", "p14", "p15", "p16", "p17", "p18", "p19", "p20", 
"p21", "p22", "p23", "p24", "p25", "p26", "p27", "p28", "p29", 
"p30", "p31", "p32", "p33", "p34", "p35", "p36", "p37", "p38", 
"p39", "p40", "p41", "p42", "p43", "p44", "p45", "p46", "p47", 
"p48", "p49", "p50", "p51", "p52", "p53", "p54", "p55", "p56", 
"p57", "p58", "p59", "p61", "p62", "p63", "p64", "p65", "p66", 
"p67", "p68", "p69", "p70", "p71", "p73", "p74", "p75", "p76", 
"p77", "p78", "p79", "p80", "p81", "p82", "p83", "p85", "p86", 
"p87"), c("DCA1", "DCA2", "DCA3", "DCA4"))), cproj = structure(c(1.25651540418674, 
4.10086480672813, 0.985170303000948, 0.356752380731314, 0.0871005927968299, 
-0.65389620717045, 0.416500256857967, -0.0687562150601063, -0.448057138258203, 
-0.851379298478351, -4.22395990436862, 3.38456796429771, 3.50656010401288, 
3.0970679364606, 3.76068022295243, 4.15832492473722, 1.73508588573266, 
-0.0979648246992866, -1.54594775838685, 0.0871005366833882, 2.10389514743405, 
0.705016209734685, -1.54594765746973, -0.264737082168187, 1.28992362813467, 
...
2.37831113227017, 2.6918186855265, 3.96654774047496), .Dim = c(414L, 
4L), .Dimnames = list(c("Ace_cam", "Ace_pla", "Ace_pse", "Ach_mil", 
"Aeg_pod", "Aes_hip", "Agr_cap", "Agr_sto", "All_pet", "Alo_pra", 
"Alo_myo", "Ama_alb", "Ama_bli", "Ama_pow", "Ama_ret", "Amb_art", 
"Ame_lam", "Ang_arc", "Anc_arv", "Ant_cau", "Ant_syl", "Ant_tin", 
"Ant_vul", "Ape_spi", "Aqu_vul", "Ara_tha", "Arc_lap", "Arc_min", 
"Are_ser", "Arr_ela", "Art_bie", "Art_cam", "Art_vul", "Asp_rut", 
...

Answer 1

Try adding each ellipse separately and assign the colors from colvec to each Habitat group separately using show.groups. For example:

with(calcul.env, ordiellipse(calcul.dca, Habitat, kind="se", conf=0.95,
                       lwd=2, col="black", show.groups = "Brown"))
with(calcul.env, ordiellipse(calcul.dca, Habitat, kind="se", conf=0.95,
                       lwd=2, col="darkviolet", show.groups = "Emb"))
with(calcul.env, ordiellipse(calcul.dca, Habitat, kind="se", conf=0.95,
                       lwd=2, col="green3", show.groups = "Eph"))

and so forth. You can make this easier by creating a loop if you want. See http://www.mail-archive.com/r-help@r-project.org/msg149786.html