1ICFT <- function (design, digits = 3, with.blocks = FALSE, 2 conc = TRUE, recode=TRUE) 3{ 4 ### function for detailed inspection for 5 ### a single combination of factors in terms of their 6 ### ICs 7 8 ### for overall inspection of entire designs, consider function ICFTs 9 10 if ("design" %in% class(design)) { 11 fn <- names(factor.names(design)) 12 if (with.blocks) 13 fn <- c(fn, design.info(design)$block.name) 14 design <- design[, fn] 15 nfac <- length(fn) 16 } 17 else { 18 nfac <- ncol(design) 19 fn <- 1:nfac 20 } 21 22 nlev <- levels.no(design) 23 dfs <- nlev - 1 24 25 if (!is.data.frame(design)) design <- as.data.frame(design) 26 if (recode) 27 for (i in 1:nfac){ 28 design[[i]] <- factor(design[[i]]) 29 contrasts(design[[i]]) <- contr.XuWu(nlev[i]) 30 } 31 32 k <- nfac 33 N <- nrow(design) 34 35 ns <- 1 36 auswahl <- 1 37 selproj <- sel <- matrix(1:k, ncol=1) 38 selproj <- paste(selproj, collapse = ":") 39 names(auswahl) <- selproj 40 41 hilf2 <- design 42 mm <- model.matrix(formula(substitute(~.^km1, 43 list(km1 = k))), data = hilf2) 44 mm <- mm[,-(1:(ncol(mm) - prod(dfs)))] 45 hilfc <- svd(mm) 46 ## spans of constant sv2s: from to cumcount 47 48 sv2s <- hilfc$d^2 49 mean.u2s <- colMeans(hilfc$u)^2 50 ICs <- sv2s*mean.u2s 51 52 ## check for and cure ambiguities 53 hilf2 <- table(round(sv2s,6)) 54 cumcounts <- cumsum(rev(hilf2)) 55 from <- c(1, cumcounts[-length(hilf2)]+1) 56 multi <- which(cumcounts-from > 0) 57 58 if (length(multi)>0){ 59 for (i in multi){ 60 ## rectify ambiguities 61 bereich <- from[i]:cumcounts[i] 62 li <- length(bereich) 63 rot <- HouseholderRotToOne(colMeans(hilfc$u)[bereich])$rot 64 #hilf2 <- rep(0, li) ## initialize output vector 65 if (conc){ 66 #hilf2[1] <- sum((sv2s*mean.u2s)[bereich]) 67 Q <- t(rot) ## matrix concentrating on first component 68 } 69 else { 70 ## even case 71 #hilf2 <- rep(mean((sv2s*colMeans(hilfc$u)^2)[bereich]),li) 72 Q <- rect_simplex(li) 73 Q <- t(rot)%*%Q 74 } 75 hilfc$u[,bereich] <- hilfc$u[,bereich]%*%Q 76 hilfc$v[,bereich] <- hilfc$v[,bereich]%*%Q 77 } 78 79 ## make colmeans of u non-negative (added March 9 2017) 80 ## for unique vector c 81 cmu <- colMeans(hilfc$u) 82 vorz <- sign(cmu) 83 hilfc$u <- hilfc$u%*%diag(vorz) 84 hilfc$v <- hilfc$v%*%diag(vorz) 85 86 ## redo after rotations 87 mean.u2s <- cmu^2 88 ICs <- sv2s*mean.u2s 89 } 90 91 rund <- round(ICs, digits) 92 93 ICFT <- table(unlist(rund)) 94 ICFT <- cbind(IC = as.numeric(names(ICFT)), frequency = ICFT) 95 rownames(ICFT) <- rep("", nrow(ICFT)) 96 aus <- list(ICFT = ICFT, 97 ICs = ICs, 98 sv2s = sv2s, 99 mean.u2s = mean.u2s, 100 mm = mm, 101 u = hilfc$u, 102 v = hilfc$v, 103 c.worst = hilfc$d*colMeans(hilfc$u)/sqrt(sum(sv2s*mean.u2s)) 104 ) 105 aus 106} 107