############################################################# # STAT 718A Computer Practical 2 ############################################################# library(shapes) # we'll do some more analysis of the digit 3 data plotshapes(digit3.dat, joinline=c(1:13) ) # Carry out Procrustes analysis (calculated using complex eigenvectors) out <- procGPA( digit3.dat , eigen2d = TRUE ) out\$rotated # contains the Procrustes registered data plotshapes ( out\$rotated , joinline=c(1:13) ) # procrustes registered data out\$mshape # contains the Procrustes mean shape lines( out\$mshape, col=2 ) # adds the mean shape to the plot (in red) plot( out\$size , out\$rho ) # NOTE `rho' is called the Riemannian distance # and it is given by sin rho = d_F , # where d_F is the full Procrustes distance # given in our notes. # here out\$rho is a vector of shape distances to the mean text( out\$size , out\$rho , c(1:30) ) # is there a possible outlier ? plotshapes( digit3.dat , digit3.dat[,,1] , joinline= 1:13) # To see all the calculations given by the command procGPA() # type summary( out ) #and look at help( procGPA ) # Shape variability print( out\$rmsd1 ) #full Procrustes distances print( out\$rho ) #Riemannian distances # We'll now look at the first few principal components out\$percent # plot first two PC scores plot(out\$rawscores[,1],out\$rawscores[,2], xlab="pc1",ylab="pc2") # pairwise plots - do you see any interesting structure ? pairs( cbind( out\$size, out\$rho, out\$rawscores[,1:3] ) , label=c("size","rho","pc1","pc2","pc3") ) # Let us interpret the structure of shape variability in the # first 3 PCs shapepca ( out , 1:3 , type="r" , joinline=1:13 ) # row plot shapepca ( out , 1:3 , type="v" , joinline=1:13 ) #vector plot shapepca ( out , 1:3 , type="s" , joinline=1:13 ) #superposition plot shapepca ( out , 1:3 , type="m" , joinline=1:13 ) # movie plot ################################################################## # Now let us remove the 1st digit as it looks a bit strange. newdigit3.dat <- digit3.dat[,,-1] # check dimension dim(newdigit3.dat) #### carry out the above analysis on the digit3 data but with the #### outlier removed. #Is there a big difference in mean shape and shape variability #with or without the outlier ? #################################################################### #Let's now compare the Procrustes mean shape and Bookstein mean shapes. #Consider the female Chimpanzee data help(panf.dat) dim(panf.dat) proc <- procGPA( panf.dat ) plotshapes( panf.dat , proc\$mshape ) book <- bookstein2d( panf.dat ) plotshapes( book\$bshpv , book\$mshape ) final <- procOPA( proc\$mshape , book\$mshape ) plotshapes( proc\$mshape , final\$Bhat) print( riemdist( proc\$mshape, book\$mshape ) ) # is there much difference ?? #####################################################################