"screeplot.faccr"<-
function(x, min.factor = 1, max.factor = 0, method = "wl1fit", drawplot = T, 
	tol = 0.01, iter.max = 10)
{
#
# P. Filzmoser (Vienna, Austria) and C. Croux (Brussels, Belgium)
#
# July 29th, 1999
#
# Scree plot (i.e. plot of the explained cumulative variance) for robust factor analysis
# using criss-cross regressions
#
#
# x ............ data matrix
# min.factor ... minimum number of factors for the scree plot
# max.factor ... maximum number of factors for the scree plot (by default p/2)
# method ....... regression method: "lsfit", "l1fit", "wl1fit","lmsreg","ltsreg", "rreg", ...
# drawplot ..... if TRUE, the scree plot is shown
# tol .......... tolerance for convergence
# iter.max ..... maximum number of iteration steps
#
#
# 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
#       k .... 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
# 
#    rexpl .... cumulative variance explained by "min.factor" to "max.factor" factors
#
#
# Begin of the procedure:
#
	n <- dim(x)[1]
	p <- dim(x)[2]
	if(max.factor == 0)
		kmax <- min(n, p) - 1
	else kmax <- max.factor
	if(min.factor > kmax)
		stop("\n min.factor is larger than max.factor!\n")
	kmin <- min.factor
	if(method == "wl1fit") {
		if((n <= 2 * kmax) | (p <= 2 * kmax))
			stop("\nThe number of rows and columns should at least be > 2*kmax for WL1fit !\n"
				)
	}
	rexpl <- NULL
	for(k in kmin:kmax) {
		if(method == "lsfit") {
			result <- factanal.ccr(x, factors = k, method = "lsfit",
				iter.max = iter.max, tol = tol, orth = F)
			xs <- scale(x, center = result$center, scale = result$
				scale)
			resid <- xs - result$scores %*% t(result$loadings)
			objopt <- sqrt(sum(resid * resid))
			objzero <- sqrt(sum(xs * xs))
			rexpl <- c(rexpl, 1 - (objopt/objzero)^2)
		}
		if(method == "l1fit") {
			result <- factanal.ccr(x, factors = k, method = "l1fit",
				iter.max = iter.max, tol = tol, orth = F)
			xs <- scale(x, center = result$center, scale = result$
				scale)
			resid <- xs - result$scores %*% t(result$loadings)
			objopt <- sum(abs(resid))
			objzero <- sum(abs(xs))
			rexpl <- c(rexpl, 1 - (objopt/objzero)^2)
		}
		if(method == "wl1fit") {
			result <- factanal.ccr(x, factors = k, method = 
				"wl1fit", iter.max = iter.max, tol = tol, orth
				 = F)
			xs <- scale(x, center = result$center, scale = result$
				scale)
			resid <- xs - result$scores %*% t(result$loadings)
			objopt <- t(result$rw) %*% abs(resid) %*% result$cw
			objzero <- t(result$rw) %*% abs(xs) %*% result$cw
			rexpl <- c(rexpl, 1 - (objopt/objzero)^2)
		}
	}
	if(kmax == (min(n, p) - 1)) {
		kmax <- kmax + 1
		rexpl <- c(rexpl, 1)
	}
	if(drawplot == T) {
		if(kmin == 1)
			plot(0:kmax, c(0, rexpl), type = "b", lty = 2, ylim = c(
				0, 1), xlab = "Number of Factors", ylab = 
				"Fitting measure")
		else barplot(rexpl)
	}
	rexpl
}