sumscal {lazy.mds}

R Documentation

SUMSCAL

Description

Given symmetric matrices, B_k, k=1,2, ... , nG, this function minimizes the LS criterion

RSS = sum_k ssq(B_k - Xc diag(W_k)^2 Xc')^2

with respect to Xc and diagonal matrices W_k, k=1,2, ... , nG using two-stage least squares method.

Usage

sumscal(
  B,
  ndim = 2,
  maxiter = 50,
  eps = 1e-07,
  print = 2,
  estm = 0,
  Xall = NULL
)

Arguments

B	Input scalar product array
ndim	# of dimensions
maxiter	max # of iterations
eps	convergence criterion for B
print	= 1 to print result
estm	= 1 to subtract column mean from X when Xall is given.
Xall	An Array of X_k's, or NULL

Details

Let B_k, k=1,2, ... , nG, be the scalar product matrices calculated from the Euclidean distance matrices, D_k, as

B_k = -2 J D_k^2 J

where J is the column centering operator.
This program find those Xc and diagonal W_k matrices that approximately minimize

RSS = sum_k ssq( B_k - Xc diag(W_k)^2 Xc' )^2 .

Note that, if the D_k matrices above were calculated from the X_k matrices which can be expressed as

X_k = Xc diag(W_k) T_k + Error

we have

B_k == X_k X_k' == Xc diag(W_k)^2 Xc'

Therefore, given nG X_k matrices, we may find Xc, W_k and T_k matrices by SUMSCAL by first calculating the Euclidean distance matrices from X_k and then double centering them to get B_k matrices. and finally applying SUMSCAL to those B_k matrices to obtain Xc and W_k matrices.
Given Xc and W_k matrices, T_k matrices can be found by orthogonal Procrustes rotation.


When X_k are present as array Xall, they have priority over B.
That is, B_k matries will be calculated from Xall array/

Xc matrix will be normalized so that diag(t(Xc)%*%Xc) = n I and the columns will be sorted according to the magnitude of diag(t(W)%*%W).

Value

A list of
X=Xc, W, T, rmseB, rmseX

References

De Leeuw, Jan, and Pruzansky, Sandra (1978) A New Computational Method to Fit the Weighted Euclidean Distance Model. Psychometrika, v43 n4 p479-9

Examples

set.seed(1701)
# generate O from X and W
X <- matrix( c( -1,1, 0,1, 1,1,  -1,0, 0,0, 1,0, -1,-1, 0,-1, 1,-1  )
, ,2, byrow=1 )
dimnames(X) <- list(paste("s",1:nrow(X),sep="")
, paste("d",1:ncol(X),sep=""))
X <- standard( X )
W <- matrix( c( 1,1, 2,1, 1,1.5 ),,2, byrow=1 )
rownames(W) <- paste("sbj",1:nrow(W),sep=""); colnames(W) <- colnames(X)
resgen <- gendatamds( X=X, W=W^2, stderror=0.5 )
O <- vechinv(resgen$vO, nodiag=1, array=1)  # This is n x n x nG.
# double center O
n <- nrow(X); nG <- nrow(W)
J <- diag(n)-(1/n)*matrix(1,n,1)%*%matrix(1,1,n)
B <- O
for( k in 1:nG ){
 B[,,k] <- -0.5*J%*%O[,,k]^2%*%J
}
Print(O, fmt=".2")

# res1 <- indscal( B, ndim=2, print=1 )
res2 <- sumscal( B, ndim=2, print=1 )

---