1ICFTs <- function (design, digits = 3, resk.only = TRUE, 2 kmin = NULL, kmax = ncol(design), detail = FALSE, with.blocks = FALSE, 3 conc = TRUE) 4{ 5 if ("design" %in% class(design)) { 6 fn <- names(factor.names(design)) 7 if (with.blocks) 8 fn <- c(fn, design.info(design)$block.name) 9 design <- design[, fn] 10 nfac <- length(fn) 11 } 12 else { 13 nfac <- ncol(design) 14 fn <- 1:nfac 15 } 16 17 nlev <- levels.no(design) 18 dfs <- nlev - 1 19 20 if (!is.data.frame(design)) design <- as.data.frame(design) 21 for (i in 1:nfac){ 22 design[[i]] <- factor(design[[i]]) 23 contrasts(design[[i]]) <- contr.XuWuPoly(nlev[i]) 24 } 25 26 ks <- which(round(GWLP(design, kmax = kmax)[-1], 8) > 0) 27 N <- nrow(design) 28 if (length(ks) == 0) { 29 hilf <- list(list(ICFT = cbind(IC = 0, frequency = sum(nlev) - 30 kmax), IC1 = 0)) 31 names(hilf) <- kmax 32 return(hilf) 33 } 34 k <- min(ks) 35 if (k < 2) 36 stop("resolution of design must be at least 2") 37 kminset <- FALSE 38 if (is.null(kmin)) 39 kmin <- k 40 else { 41 if (!kmin == min(ks)) kminset <- TRUE 42 redu <- ks[ks >= kmin] 43 message(paste("Check sets of sizes ", paste(redu, collapse = ","))) 44 if (length(redu) == 0) 45 return() 46 kmin <- min(redu) 47 k <- kmin 48 ks <- redu 49 } 50 if (k >= kmin) { 51# if (!"design" %in% class(design)) { 52# if (!is.data.frame(design)) 53# design <- as.data.frame(design) 54# faktoren <- sapply(design, is.factor) 55# keinfaktor <- which(!faktoren) 56# if (length(keinfaktor) > 0) 57# for (i in keinfaktor) design[[i]] <- as.factor(design[[i]]) 58# } 59 k <- kmin 60 ns <- choose(nfac, k) 61 auswahl <- 1:ns 62 selproj <- sel <- nchoosek(nfac, k) 63 GWLPs <- round(apply(selproj, 2, function(obj) GWLP(design[, 64 obj])[-1]), 4) 65 selproj <- apply(selproj, 2, function(obj) paste(obj, 66 collapse = ":")) 67 names(auswahl) <- selproj 68 # ergproj <- rep(NA, length(selproj)) 69 # dimnames(erg3) <- list(factor = fn, others = apply(sel, 70 # 2, function(obj) paste(obj, collapse = ":")), 1:(max(nlev) - 71 # 1)) 72 if (resk.only) { 73 reskproj <- apply(GWLPs, 2, function(obj) all(obj[-k] == 74 0)) 75 if (all(!reskproj)) { 76 message("no projections with resolution ", 77 k, " or higher") 78 return() 79 }} 80 81 berechn <- lapply(auswahl, function(obj) { 82 hilf2 <- design[, sel[, obj]] 83 mmX <- model.matrix(formula(substitute(~.^km1, 84 list(km1 = k))), data = hilf2) 85 mmX <- mmX[,-(1:(ncol(mmX) - prod(dfs[sel[, obj]])))] 86 hilfc <- svd(mmX) 87 hilf2 <- table(round(hilfc$d^2,6)) 88 cumcounts <- cumsum(rev(hilf2)) 89 from <- c(1, cumcounts[-length(hilf2)]+1) 90 hilf2 <- rep(0, prod(dfs[sel[, obj]])) ## initialize output vector 91 if (conc) 92 hilf2[from] <- sapply(1:length(from), 93 function(obj) sum((hilfc$d^2*colMeans(hilfc$u)^2)[from[obj]:cumcounts[obj]]) 94 ) 95 else { 96 for (i in 1:length(from)){ 97 bereich <- from[i]:cumcounts[i] 98 hilf2[bereich] <- mean((hilfc$d^2*colMeans(hilfc$u)^2)[bereich]) 99 }} 100 list(hilf2, hilfc$d^2, colMeans(hilfc$u)^2) 101 }) 102 103 ICs <- lapply(berechn, function(obj) obj[[1]]) 104 sv2s <- lapply(berechn, function(obj) obj[[2]]) 105 mean.u2s <- lapply(berechn, function(obj) obj[[3]]) 106 rund <- lapply(ICs, function(obj) round(obj,digits)) 107 108 ICFT <- table(unlist(rund)) 109 ICFT <- cbind(IC = as.numeric(names(ICFT)), frequency = ICFT) 110 rownames(ICFT) <- rep("", nrow(ICFT)) 111 aus <- list(ICFT = ICFT) 112 113 if (detail) 114 aus <- c(aus, list(ICs = rund, sv2s = sv2s, mean.u2s = mean.u2s )) 115 116 if (!resk.only || kminset){ 117 aus <- list(aus); names(aus) <- k 118 if (!resk.only){ 119 ks <- kmin:kmax 120 if (length(ks)>1){ 121 for (k in (kmin+1):kmax){ 122 ns <- choose(nfac, k) 123 auswahl <- 1:ns 124 selproj <- sel <- nchoosek(nfac, k) 125 GWLPs <- round(apply(selproj, 2, function(obj) GWLP(design[, 126 obj])[-1]), 4) 127 selproj <- apply(selproj, 2, function(obj) paste(obj, 128 collapse = ":")) 129 names(auswahl) <- selproj 130 131 berechn <- lapply(auswahl, function(obj) { 132 hilf2 <- design[, sel[, obj]] 133 mmX <- model.matrix(formula(substitute(~.^km1, 134 list(km1 = k))), data = hilf2) 135 mmX <- mmX[,-(1:(ncol(mmX) - prod(dfs[sel[, obj]])))] 136 hilfc <- svd(mmX) 137 hilf2 <- table(round(hilfc$d^2,6)) 138 cumcounts <- cumsum(rev(hilf2)) 139 from <- c(1, cumcounts[-length(hilf2)]+1) 140 hilf2 <- rep(0, prod(dfs[sel[, obj]])) ## initialize output vector 141 if (conc) 142 hilf2[from] <- sapply(1:length(from), 143 function(obj) sum((hilfc$d^2*colMeans(hilfc$u)^2)[from[obj]:cumcounts[obj]]) 144 ) 145 else { 146 for (i in 1:length(from)){ 147 bereich <- from[i]:cumcounts[i] 148 hilf2[bereich] <- mean((hilfc$d^2*colMeans(hilfc$u)^2)[bereich]) 149 }} 150 list(hilf2, hilfc$d^2, colMeans(hilfc$u)^2) 151 }) 152 153 ICs <- lapply(berechn, function(obj) obj[[1]]) 154 sv2s <- lapply(berechn, function(obj) obj[[2]]) 155 mean.u2s <- lapply(berechn, function(obj) obj[[3]]) 156 rund <- lapply(ICs, function(obj) round(obj,digits)) 157 158 ICFT <- table(unlist(rund)) 159 ICFT <- cbind(IC = as.numeric(names(ICFT)), frequency = ICFT) 160 rownames(ICFT) <- rep("", nrow(ICFT)) 161 ausn <- list(ICFT = ICFT) 162 if (detail) 163 ausn <- c(ausn, list(ICs = rund, sv2s = sv2s, mean.u2s = mean.u2s )) 164 ausn <- list(ausn) 165 names(ausn) <- k 166 aus <- c(aus, ausn) 167 } 168 } 169 }} 170 } 171 else aus <- list(ICFT = NULL) 172 aus 173} 174