Search code examples
r3dpcargl

3D PCA plotting error: the leading minor of order 3 is not positive definite

I have been trying to troubleshoot a problem for the past couple of hours with no luck.

I have found one error message on another post that is similar here : Plotting Ellipse3d in R Plotly with surface ellipse

The only response highlights that it is because of a lack of samples. But my sample size = 3; which is the minimum for creating an ellipse.

Another response on another post highlights that it could be because the matrix being entered into the ellipse3d() is not a symmetric square matrix, but in fact it is because I take the cov(data_frame_of_interest) while entering it.

More about the problem:

I am trying to plot a 3d PCA using rgl package in R and prcomp(); where top_noSR3 is a prcomp() object. I am interested in the first 3 PC's for the following example.

#Here is the exact PCA object obtained from prcomp(). It is 
#top_noSR3$x[,c(1,2,3)].    

top_noSR3 <- list(x = structure(c(-1.51867921165966, 4.7156146538954, -0.812795773332441, 
-2.38413966890339, -2.57305046487183, 8.4906222835565, 1.66680581870767, 
-7.5843776373923, 3.35199279752431, 2.4452307290032, -7.33666903970592, 
1.53944551317846, -1.72127982564219, 2.68563509110228, -2.75720515666432, 
1.79284989120424, -0.962414945057766, -1.44498205066367, -0.270435252418625, 
2.67783224814005, -0.780985238414212, 7.3028581587737, -1.3641698814598, 
-5.15770303889969, 0.989625278145346, 0.101301456051966, 1.96037141973899, 
-3.05129815393625, -2.29377625653868, -4.3601221831666, 5.9554046825257, 
0.698493757179566, 2.3724852157055, 0.0505080025883657, -1.97013906903266, 
-0.452854149261224), .Dim = c(12L, 3L), .Dimnames = list(c("SR1.4", 
"SR1.5", "SR1.6", "SR1.7", "SR2.2", "SR2.3", "SR2.4", "SR2.5", 
"CI1.4", "CI1.5", "CI1.6", "CI1.7"), c("PC1", "PC2", "PC3"))))

#col_object is an array of colours used for colouring the points in the PCA:
col_object <- c("blue","red","green","purple","blue","red",   
"green","purple","blue","red","green","purple")

plot3d(top_noSR3$x[,c(1,2,3)], type = "s", size = 1, col = col_object)

#Divide object into 3 components for 3d ellipse plotting function

for_plot3d <- function(casetime){

  return(top_noSR3$x[casetime, c(1,2,3)])
}

BASE <- c("SR1.4", "SR2.2", "CI1.4")
EA <- c("SR1.5", "SR2.3", "CI1.5")
LA <- c("SR1.6", "SR2.4", "CI1.6")
FU <- c("SR1.7", "SR2.5", "CI1.7")

BASE_for3d <- for_plot3d(BASE)
EA_for3d <- for_plot3d(EA)
LA_for3d <- for_plot3d(LA)
FU_for3d <- for_plot3d(FU)

#The following arguments used above (BASE,EA, etc.) ^ are 
#timepoints/rownames of interest. length(BASE) = ... = length(FU) = 3


center4plot <- function(x){
  k <- mean(x[,1])
  h <- mean(x[,2])
  z <- mean(x[,3])
  return(c(k,h,z))
}

center_base <- center4plot(BASE_for3d)
center_EA <- center4plot(EA_for3d)
center_LA <- center4plot(LA_for3d)
center_FU <- center4plot(FU_for3d)

ellipse_base <- ellipse3d(cov(BASE_for3d),level = 0.95, centre=center_base)
ellipse_EA <- ellipse3d(cov(EA_for3d), level = 0.95, centre=center_EA)

#ERROR appears here, in ellipse_LA
ellipse_LA <- ellipse3d(cov(LA_for3d), level = 0.95, centre=center_LA)
ellipse_FU <- ellipse3d(cov(FU_for3d), level = 0.95, centre=center_FU)

I get the following message when i run the line with ellipse_LA: ellipse3d() Error in chol.default(cov) : the leading minor of order 3 is not positive definite

Thank you all for your time and effort.

Solution

  • When I run the current version of the code, I get the error you saw here:

    > ellipse_base <- ellipse3d(cov(BASE_for3d),level = 0.95, centre=center_base)
Error in chol.default(cov) : 
  the leading minor of order 3 is not positive definite

    If I look at the eigenvalues of that matrix, I see that the message is correct:

    > eigen(cov(BASE_for3d))
eigen() decomposition
$values
[1] 1.454977e+01 1.433775e+00 8.521084e-16

$vectors
           [,1]       [,2]       [,3]
[1,] 0.81230516  0.5213444  0.2614581
[2,] 0.04362183  0.3927276 -0.9186197
[3,] 0.58159906 -0.7576048 -0.2962727

    Notice that the smallest eigenvalue is approximately zero.

    The BASE_for3d object has 3 points in it. They all lie in a plane (as any 3 points must). Thus the covariance matrix for them is rank 2, not rank 3. You need to increase BASE to include at least 4 rows not all in a plane, not just 3. For example,

    BASE <- c("SR1.4", "SR2.2", "CI1.4", "SR1.5")

    This results in the following:

    BASE_for3d <- for_plot3d(BASE)
center_base <- center4plot(BASE_for3d)
ellipse_base <- ellipse3d(cov(BASE_for3d),level = 0.95, centre=center_base)
shade3d(ellipse_base, alpha = 0.2)

    screenshot