"screeplot.fanova"<-
function(x, min.factor = 1, max.factor = 0, model = "b", 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 robustly fitting the FANOVA model
#
#
# 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", ...
# model ........ "a" (additive), "m" (multiplicative), and "b" (both additive and multiplicative)
# 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 FANOVA model can be computed.
# The FANOVA model is defined by 
#
# x_ij = mu + a_i + b_j + sum_{t=1}^{k} (f_it*l_jt) + e_ij
#
# for i=1,...,n and j=1,...,p, 
# with
#       x_ij .... (i,j) of the twoway table
#       mu   .... overall mean  \
#       a_i  .... row effect     >  additive parameters
#       b_j  .... column effect /
#       f_it    .... factor scores \
#                                   >  multiplicative parameters
#       l_jt    .... loadings      /
#       e_ij .... residuals
#       k .... number of factors.
# The additive model (model="a") includes just the additive parameters, the
# multiplicative model (model="m") includes just the multiplicative parameters, 
# and model="b" accounts for bith, additive and multiplicative terms.
#
# 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(model == "a") {
		if(method == "lsfit") {
			result <- ccrfit(x, factors = k, model = "a", method = 
				"lsfit", iter.max = iter.max, tol = tol, orth
				 = F)
			resid <- result$res
			resid0 <- scale(x, scale = F)
			objopt <- sqrt(sum(resid^2))
			objzero <- sqrt(sum(resid0^2))
			rexpl <- 1 - (objopt/objzero)^2
		}
		if((method == "l1fit") || (method == "wl1fit")) {
			result <- ccrfit(x, factors = k, model = "a", method = 
				"l1fit", iter.max = iter.max, tol = tol, orth
				 = F)
			resid <- result$res
			resid0 <- scale(x, center = apply(x, 2, median), scale
				 = F)
			objopt <- sum(abs(resid))
			objzero <- sum(abs(resid0))
			rexpl <- 1 - (objopt/objzero)^2
		}
	}
	else {
		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
		if(method == "lsfit") {
			if(model == "b") {
				result0 <- ccrfit(x, factors = k, model = "a", 
				  method = "lsfit", iter.max = iter.max, tol = 
				  tol, orth = F)
				resid0 <- result0$res
			}
			else {
				resid0 <- x
			}
		}
		if(method == "l1fit") {
			if(model == "b") {
				result0 <- ccrfit(x, factors = k, model = "a", 
				  method = "l1fit", iter.max = iter.max, tol = 
				  tol, orth = F)
				resid0 <- result0$res
			}
			else {
				resid0 <- x
			}
		}
		if(method == "wl1fit") {
			if(model == "b") {
				result0 <- ccrfit(x, factors = k, model = "a", 
				  method = "wl1fit", iter.max = iter.max, tol
				   = tol, orth = F)
				resid0 <- result0$res
			}
			else {
				resid0 <- x
			}
		}
		for(k in kmin:kmax) {
			if(method == "lsfit") {
				result <- ccrfit(x, factors = k, model = model, 
				  method = "lsfit", iter.max = iter.max, tol = 
				  tol, orth = F)
				resid <- result$res
				objopt <- sqrt(sum(resid^2))
				objzero <- sqrt(sum(resid0^2))
				rexpl <- c(rexpl, 1 - (objopt/objzero)^2)
			}
			if(method == "l1fit") {
				result <- ccrfit(x, factors = k, model = model, 
				  method = "l1fit", iter.max = iter.max, tol = 
				  tol, orth = F)
				resid <- result$res
				objopt <- sum(abs(resid))
				objzero <- sum(abs(resid0))
				rexpl <- c(rexpl, 1 - (objopt/objzero)^2)
			}
			if(method == "wl1fit") {
				result <- ccrfit(x, factors = k, model = model, 
				  method = "wl1fit", iter.max = iter.max, tol
				   = tol, orth = F)
				resid <- result$res
				objopt <- t(result$rw) %*% abs(resid) %*% 
				  result$cw
				objzero <- t(result$rw) %*% abs(resid0) %*% 
				  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
}