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