"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 }