[Rd] cmdscale problem

From: Christian Kleiber <kleiber_at_statistik.uni-dortmund.de>
Date: Sun 09 Jan 2005 - 02:00:28 EST


Dear R developers,

there appears to be a small problem with function cmdscale: for non-Euclidean distance matrices, using option add=FALSE (the default), cmdscale misses the smallest eigenvalue. This affects GOF statistic g.1 (See Mardia, Kent + Bibby (1979): Multivariate Analysis, eq. (14.4.7). The corresponding formula in Cox + Cox (2001): Multidimensional Scaling, 2nd ed., p 38, would seem to contain a misprint, it should be n instead of n-1.)

Example:

R> cmdscale(eurodist, eig=TRUE)$GOF
[1] 0.7968344 0.8679134

The full set of eigenvalues can be obtained via

R> H <- function(n) { diag(n) - matrix(1, n, n)/n } R> eigen( -.5 * H(21) %*% (as.matrix(eurodist))^2 %*% H(21) , + only.values=TRUE)$values

  [1]  1.953838e+07  1.185656e+07  1.528844e+06
  [4]  1.118742e+06  7.893472e+05  5.816552e+05
  [7]  2.623192e+05  1.925976e+05  1.450845e+05

[10] 1.079673e+05 5.139484e+04 3.538963e-10
[13] -9.496124e+03 -5.305820e+04 -1.322166e+05
[16] -2.573360e+05 -3.326719e+05 -5.162523e+05
[19] -9.191491e+05 -1.006504e+06 -2.251844e+06

whereas cmdscale allows access of the first 20 of these via

R> cmdscale(eurodist, k=20, eig=TRUE)$eig

  [1]  1.953838e+07  1.185656e+07  1.528844e+06
  [4]  1.118742e+06  7.893472e+05  5.816552e+05
  [7]  2.623192e+05  1.925976e+05  1.450845e+05

[10] 1.079673e+05 5.139484e+04 4.656613e-10
[13] -9.496124e+03 -5.305820e+04 -1.322166e+05
[16] -2.573360e+05 -3.326719e+05 -5.162523e+05
[19] -9.191491e+05 -1.006504e+06

Warning messages:
...

Suggestion: replace

     evalus <- e$values[-n]

with

     evalus <- e$values

This yields

R> cmdscale(eurodist, eig=TRUE)$GOF
[1] 0.7537543 0.8679134

Minor points:

May I suggest to allow access of all eigenvalues, instead of just the first k as delivered by cmdscale(... , eig=TRUE)$eig, in order to permit inspection of the severity of the problem? (I am aware that a comparison of g.1 and g.2 indicates the presence of negative eigenvalues.) It might also be useful to deliver a warning whenever some eigenvalues are < 0, using something like

if (any(e$values < 0))

       warning("some eigenvalues are < 0")

Also, the documentation says that cmdscale(...)$eig returns "the n-1 eigenvalues computed during the scaling process if 'eig' is true", whereas in fact it returns the first k.

Incidentally, there is an undocumented feature: depending on the value of if(eig || x.ret || add) the value is a list including a component 'ac' (add. constant) in addition to those mentioned in the documentation.

Best regards,
Christian Kleiber

Version: R 2.0.1
OS: Win XP Pro

-- 

*****************************************

Dr. Christian Kleiber
FB Statistik
Universitaet Dortmund
Vogelpothsweg 78
D-44227 DORTMUND
Germany

Tel: ++49-(0)231-755 5419
Fax: ++49-(0)231-755 5284
E-mail: kleiber@statistik.uni-dortmund.de

______________________________________________
R-devel@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Sun Jan 09 01:11:45 2005

This archive was generated by hypermail 2.1.8 : Sun 09 Jan 2005 - 01:17:20 EST