"factanal.ccr"<-
function(x, factors = 2, method = "wl1fit", drawplot = F, alpha = 1, tol = 
	1e-005, iter.max = 50, orth = T, print = T)
{
#
# P. Filzmoser (Vienna, Austria) and C. Croux (Brussels, Belgium)
#
# Final version for "Robust Factor Analysis based on Criss-cross Regressions"
# July 29th, 1999
#
# x ......... data matrix
# factors.... dimension of the model fit
# method .... regression method: "lsfit", "l1fit", "wl1fit","lmsreg","ltsreg", "rreg", ...
# drawplot .. if TRUE, plots will be drawn
# alpha ..... adjusting constant for biplot representation
# tol ....... tolerance for convergence
# iter.max .. maximum number of iteration steps
# orth ...... orthogonalize f and l at the end for biplot representation (alpha=1)
# print ..... if TRUE, the iteration number and the convergence criterion is printed on the screen
#
#
# This program uses the following functions (not implemented in Splus):
#    prcomp.rob
#    weight.wl1
#    ccrfit
#  
#
# DESCRIPTION: 
#
# A robust fit of the factor analysis model can be computed.
# The factor analysis model is defined by 
#
# (x_ij - b_j) / s_j = sum_{t=1}^{k} (f_it*l_jt) + e_ij
#
# for i=1,...,n and j=1,...,p, 
# with
#       x_ij .... (i,j)-th element of the data matrix
#       b_j  .... robust column mean of x 
#       s_j  .... robust scale of x
#       f_it .... factor scores
#       l_jt .... loadings
#       e_ij .... error term
#       factors.. number of factors
#
# There are different possibilities for estimating all the parameters in the
# model:
#    method="lsfit" minimizes the sum over the squared residuals, 
#    method="l1fit" minimizes the sum over the absolute residuals, 
#    method="wl1fit" minimizes the sum over the weighted absolute residuals, the
#                    weights for rows and columns are computed by the routine
#                    "weight.wl1", 
# additionally it is possible to select one of "ltsreg", "lmsreg", and "rreg".
# It is suggested to use the methods L1-Regression ("l1fit") or Weighted L1-Regression
# ("wl1fit").
#
#
# OUTPUT:
#
# The output of the procedure is
# 
#    loadings  .... estimation of the loadings (p x k)
#    scores    .... estimation of the factor scores (n x k)
#    uniquenesses . estimation of the residual variances (p x 1)
#    residuals .... final residuals (n x p)
#    center    .... estimation of the column mean (p x 1)
#    scale     .... estimation of the column scale (p x 1)
#    rweight   .... weight for outlying scores
#    cweight   .... weight for outlying loadings
#
# With the option `drawplot=T', a plot is created which includes 
#   a biplot representation of loadings and scores (for k=2), 
#   a residual plot, 
#   box plots of the parameters `b' and `scale'.
#  
#
# Begin of the procedure:
#
	n <- dim(x)[1]
	if(method == "lsfit") {
		rscale <- diag(var(x))
		rscale <- sqrt(rscale)
		rmean <- apply(x, 2, mean)
	}
	else {
		rscale <- apply(x, 2, mad)
		rmean <- apply(x, 2, median)
	}
	xs <- scale(x, rmean, rscale)
	res <- ccrfit(xs, factors = factors, model = "m", method = method, 
		drawplot = drawplot, alpha = alpha, tol = tol, iter.max = 
		iter.max, orth = orth, print = T)
	if(factors == 1) {
		ressav <- xs - res$sco %*% t(res$load)
	}
	else {
		ressav <- xs - res$sco %*% t(res$load)
	}
	unique <- apply(ressav, 2, mad)
	ans <- list(loadings = res$load, scores = res$sco, uniquenesses = 
		unique, residuals = ressav, center = rmean, scale = rscale, 
		rweight = res$rweight, cweight = res$cweight, objective = res$
		obj, iterations = res$iter, method = method, call = res$call)
	ans
}