1SCFTs <- function(design, digits=3, all=TRUE, resk.only=TRUE, kmin=NULL, kmax=ncol(design), 2 regcheck=FALSE, arft=TRUE, cancors=FALSE, with.blocks=FALSE){ 3 ## returns list with SCFTs and ARFTs 4 ## either for all full resolution subsets or for all subsets 5 if ("design" %in% class(design)) { 6 fn <- names(factor.names(design)) 7 if (with.blocks) fn <- c(fn, design.info(design)$block.name) 8 design <- design[,fn] 9 nfac <- length(fn) 10 } 11 else{ 12 nfac <- ncol(design) 13 fn <- 1:nfac 14 } 15 16 nlev <- levels.no(design) 17 dfs <- nlev-1 18 ks <- which(round(GWLP(design, kmax=kmax)[-1],8)>0) 19 if (length(ks)==0) { 20 hilf <- list(list(SCFT=cbind(SC=0, frequency=sum(nlev) - kmax), 21 ARFT=cbind(aveR2=0, frequency=kmax), cancor1=0)) 22 names(hilf) <- kmax 23 return(hilf) 24 } 25 k <- min(ks) ## here, k is resolution 26 if (k < 2) stop("resolution of design must be at least 2") 27 kminset <- FALSE ## kminset TRUE: user has chosen kmin 28 if (is.null(kmin)) kmin <- k 29 else { 30 kminset <- TRUE 31 redu <- ks[ks >= kmin] ## reduced set of levels 32 message(paste("Check sets of sizes ", paste(redu, collapse=","))) 33 if (length(redu)==0) return() 34 kmin <- min(redu) 35 k <- kmin 36 ks <- redu ## added March 2016, after commenting out the interim 0 calculations 37 } 38 if (k >= kmin) 39 { 40 ## should always be reached 41 if (!"design" %in% class(design)){ 42 ## make sure these are factors 43 if (!is.data.frame(design)) design <- as.data.frame(design) 44 faktoren <- sapply(design, is.factor) 45 keinfaktor <- which(!faktoren) 46 if (length(keinfaktor)>0) for (i in keinfaktor) design[[i]] <- as.factor(design[[i]]) 47 } 48 49 ### calculations for k=kmin 50 k <- kmin 51 nGRindi <- choose(nfac-1,k-1) 52 sel <- nchoosek(nfac-1,k-1) 53 auswahl <- 1:nGRindi 54 erg3 <- array(NA, c(nfac, nGRindi, max(nlev)-1)) 55 selproj <- nchoosek(nfac,k) 56 ## added March 2016 57 ## determine resolution k projections 58 GWLPs <- round(apply(selproj, 2, function(obj) GWLP(design[,obj])[-1]),4) 59 ## column for each projection 60 selproj <- apply(selproj, 2, function(obj) paste(obj, collapse=":")) 61 ergproj <- rep(NA, length(selproj)) 62 dimnames(erg3) <- list(factor=fn, others=apply(sel, 2, function(obj) paste(obj, collapse=":")), 1:(max(nlev)-1)) 63 for (i in 1:nfac){ 64 ## added March 2016 65 ## determine resolution k projections 66 if (resk.only){ 67 reskproj <- apply(GWLPs, 2, function(obj) all(obj[-k]==0)) 68 if (all(!reskproj)){ 69 message("no projections with resolution ", k, " or higher") 70 return() 71 } 72 } 73 74 seli <- sel 75 seli[seli>=i] <- seli[seli>=i]+1 76 proji <- apply(seli, 2, function(obj) paste(sort(c(i,obj)),collapse=":")) 77 indexes <- sapply(proji, function(obj) which(selproj==obj)) 78 berechn <- t(sapply(auswahl, function(obj){ 79 spalten <- c(i,seli[,obj]) 80 hilf2 <- design[,spalten] 81 mmX <- model.matrix(~., data=hilf2[,1,drop=FALSE])[,-1] ## first factor 82 83 if (k>2) 84 mmOther <- model.matrix(formula(substitute(~.^km1, list(km1 = k-1))), 85 data=hilf2[,-1,drop=FALSE])[,-1] 86 else 87 mmOther <- model.matrix(~., data=hilf2[,-1,drop=FALSE])[,-1] 88 hilfc <- cancor(mmOther, mmX)$cor 89 if (length(hilfc) < dfs[i]) 90 hilfc <- c(hilfc, rep(0,dfs[i]-length(hilfc))) 91 hilfc 92 })) 93 if (nrow(berechn)==1) berechn <- t(berechn) 94 erg3[i,1:nGRindi,1:dfs[i]] <- berechn 95 } 96 rund <- round(erg3, digits) 97 98 SCFT <- table(round(erg3^2, digits)) 99 SCFT <- cbind(SC=as.numeric(names(SCFT)), frequency=SCFT) 100 rownames(SCFT) <- rep("",nrow(SCFT)) 101 102 aus <- list(SCFT=SCFT) 103 if (arft) { 104 aveR2s <- apply(erg3, c(1, 2), function(obj) mean(obj^2, 105 na.rm = TRUE)) 106 ARFT <- table(round(aveR2s, digits)) 107 ARFT <- cbind(aveR2 = as.numeric(names(ARFT)), frequency = ARFT) 108 rownames(ARFT) <- rep("", nrow(ARFT)) 109 aus <- c(aus, list(ARFT=ARFT)) 110 } 111 if (cancors) aus <- c(aus, list(cancors=round(erg3^2, digits))) 112 113 if (regcheck) if (length(setdiff(aus$SCFT[,1], c(0,1)))>0){ 114 message("The design is not regular") 115 return(aus) 116 } 117 else message("The check for k=", k, " did not exclude regularity") 118 } 119 else aus <- list() 120 121 ## calculations for larger k, if requested 122 if (all){ 123 aus <- list(aus) 124 names(aus) <- k 125 ## March 2016: added condition 126 ## is this really needed, even for not only resk.only ? 127 ## it will provide interim 0 tables, which are unnecessary 128 ## given one knows the GWLP 129 ## therefore, commented out entirely 130 ## if (!resk.only) ks <- kmin:kmax ### entered in order to also obtain SCFTs for 131 ### numbers of factors with GWLP 0 132 for (k in ks[ks>kmin & ks<=kmax]){ 133 nGRindi <- choose(nfac-1,k-1) ## number of subsets excluding LHS 134 sel <- nchoosek(nfac-1,k-1) ## subsets excluding LHS 135 auswahl <- 1:nGRindi ## index vector of subsets exluding LHS 136 erg3 <- array(NA, c(nfac, nGRindi, max(nlev)-1)) 137 138 selproj <- nchoosek(nfac,k) ## k factor projections 139 140 ## determine resolution k projections 141 GWLPs <- round(apply(selproj, 2, function(obj) GWLP(design[,obj])[-1]),4) 142 ## column for each projection 143 if (resk.only){ 144 reskproj <- apply(GWLPs, 2, function(obj) all(obj[-k]==0)) 145 if (all(!reskproj)){ 146 message("no projections with resolution ", k, " or higher") 147 ks <- ks[ks<k] 148 break ## seemingly unnecessary break removed March 1 2018 149 ## and reentered as necessary June 2018 150 ## the break breaks the for loop 151 } 152 } 153 154 ## make into character descriptor 155 selproj <- apply(selproj, 2, function(obj) paste(obj, collapse=":")) 156 ergproj <- rep(NA, length(selproj)) 157 ## prepare result 158 dimnames(erg3) <- list(factor=fn, 159 others=apply(sel, 2, function(obj) 160 paste(obj, collapse=":")), 1:(max(nlev)-1)) 161 for (i in 1:nfac){ 162 seli <- sel 163 seli[seli>=i] <- seli[seli>=i]+1 164 proji <- apply(seli, 2, function(obj) paste(sort(c(i,obj)),collapse=":")) 165 indexes <- sapply(proji, function(obj) which(selproj==obj)) 166 if (resk.only) needi <- sapply(proji, function(obj) reskproj[which(selproj==obj)]) 167 else needi <- rep(TRUE, nGRindi) 168 if (any(needi)){ 169 berechn <- t(sapply(auswahl[needi], function(obj){ 170 spalten <- c(i,seli[,obj]) 171 hilf2 <- design[,spalten] 172 mmX <- model.matrix(~., data=hilf2[,1,drop=FALSE])[,-1] ## first factor 173 if (k>2) 174 mmOther <- model.matrix(formula(substitute(~.^km1, list(km1 = k-1))), 175 data=hilf2[,-1,drop=FALSE])[,-1] 176 else 177 mmOther <- model.matrix(~., data=hilf2[,-1,drop=FALSE])[,-1] 178 hilfc <- cancor(mmOther, mmX)$cor 179 if (length(hilfc) < dfs[i]) 180 hilfc <- c(hilfc, rep(0,dfs[i] - length(hilfc))) 181 hilfc 182 })) 183 if (nrow(berechn)==1) berechn <- t(berechn) 184 erg3[i,(1:nGRindi)[needi],1:dfs[i]] <- berechn 185 } 186 } 187 rund <- round(erg3, digits) 188 SCFT <- table(round(erg3^2, digits)) 189 SCFT <- cbind(SC=as.numeric(names(SCFT)), frequency=SCFT) 190 rownames(SCFT) <- rep("",nrow(SCFT)) 191 ausadd <- list(SCFT=SCFT) 192 193 if (arft) { 194 aveR2s <- apply(erg3, c(1, 2), function(obj) mean(obj^2, 195 na.rm = TRUE)) 196 ARFT <- table(round(aveR2s, digits)) 197 ARFT <- cbind(aveR2 = as.numeric(names(ARFT)), frequency = ARFT) 198 rownames(ARFT) <- rep("", nrow(ARFT)) 199 ausadd <- c(ausadd, list(ARFT=ARFT)) 200 } 201 202 if (cancors) ausadd <- c(ausadd, list(cancors=round(erg3^2, digits))) 203 aus <- c(aus, list(ausadd)) 204 205 if (regcheck) 206 if (length(setdiff(ausadd$SCFT[,1], c(0,1)))>0){ 207 message("The design is not regular") 208 names(aus) <- ks[1:which(ks==k)] 209 return(aus) 210 } 211 else message("The check for k=", k, " did not exclude regularity") 212 } 213 names(aus) <- ks 214 } 215 aus 216 217 } 218 219