Search code examples
rloopspoint-cloudsrgl

How to build multiple 3D ellipses grouped by IDs?

How to build multiple 3D ellipses around 3D points grouped by their IDs?

Following example shows the wanted result:

# load data
data(iris)
head(iris)
# add column
iris$ID <- iris$Species
# change column values
within(iris, ID <- factor(ID, labels = 1:3))

# add variables
x <- sep.l <- iris$Sepal.Length
y <- pet.l <- iris$Petal.Length
z <- sep.w <- iris$Sepal.Width

x.set <- iris$Sepal.Length[(iris$Species == "setosa")]
y.set <- iris$Petal.Length[(iris$Species == "setosa")]
z.set <- iris$Sepal.Width[(iris$Species == "setosa")]

x.ver <- iris$Sepal.Length[(iris$Species == "versicolor")]
y.ver <- iris$Petal.Length[(iris$Species == "versicolor")]
z.ver <- iris$Sepal.Width[(iris$Species == "versicolor")]

x.vir <- iris$Sepal.Length[(iris$Species == "virginica")]
y.vir <- iris$Petal.Length[(iris$Species == "virginica")]
z.vir <- iris$Sepal.Width[(iris$Species == "virginica")]

# create 3D ellipses
ellips.set <- ellipse3d(cov(cbind(x.set, y.set, z.set)), centre = c(mean(x.set), mean(y.set), mean(z.set)), level = 0.9)
ellips.ver <- ellipse3d(cov(cbind(x.ver, y.ver, z.ver)), centre = c(mean(x.ver), mean(y.ver), mean(z.ver)), level = 0.9)
ellips.vir <- ellipse3d(cov(cbind(x.vir, y.vir, z.vir)), centre = c(mean(x.vir), mean(y.vir), mean(z.vir)), level = 0.9)

# plot
plot3d(x, y, z, col="black", box = FALSE, type ="s", radius = 0.10)

plot3d(ellips.set, col = "red", alpha = 0.5, add = TRUE, type = "wire")
plot3d(ellips.ver, col = "green", alpha = 0.5, add = TRUE, type = "wire")
plot3d(ellips.vir, col = "gold", alpha = 0.5, add = TRUE, type = "wire")

However, instead of creating ellipses manually they should be created in a loop. 

   for (i in 1:length(iris$ID)) {
      ellipse.iris <- ellipse3d(cov(cbind(
           x=iris$Sepal.Length[(iris$ID == i)], y=iris$Petal.Length[(iris$ID == i)], z=iris$Sepal.Width[(iris$ID == i)])), 
           centre=c(mean(iris$Sepal.Length), mean(iris$Petal.Length), mean(levels(iris$Sepal.Width))), level = 0.95)
    }

Doesn´t work. It´s throwing following error message:

Error in chol.default(cov) : the leading minor of order 1 is not positive definite

Solution

  • You have a few problems.

    First, you are calculating the covariance of single rows of the matrix. It would be better to have

    for (id in unique(iris$ID))

    The center calculation then needs to be changed to

    centre=c(mean(iris$Sepal.Length[iris$ID == id]), 
       mean(iris$Petal.Length[iris$ID == id]), 
       mean(iris$Sepal.Width[iris$ID == id]))

    And then you need to plot it. So the loop should be

    for (id in unique(iris$ID)) {
      ellipse.iris <- ellipse3d(cov(cbind(
           x=iris$Sepal.Length[(iris$ID == id)], y=iris$Petal.Length[(iris$ID == id)], z=iris$Sepal.Width[(iris$ID == id)])), 
           centre=c(mean(iris$Sepal.Length[iris$ID == id]), 
             mean(iris$Petal.Length[iris$ID == id]), 
             mean(iris$Sepal.Width[iris$ID == id])), level = 0.95)
      shade3d(ellipse.iris, alpha = 0.3)
}