1GRind <- function(design, digits=3, arft=TRUE, scft=TRUE, cancors=FALSE, with.blocks=FALSE){ 2 ## returns list of four or five 3 ## with first element GRs a vector of GR and GRind, 4 ## second element GR.i a matrix 2 x nfac elements 5 ## with rows for GRtot.i and GRind.i 6 ## third element the average R2 frequency table 7 ## fourth element the squared canonical correlations table 8 ## fifth element the canonical correlations 9 ## fifth element is an nfac x choose(nfac-1, R-1) x max(nlev)-1 array 10 ## returned only if cancors is set to TRUE 11 ## the canonical correlations are supplemented with 0es, 12 ## if there are fewer than the respective dfi 13 ## version 0.27: make sure only factors are included in calculations 14 ## with.blocks only used for class design 15 if ("design" %in% class(design)) { 16 fn <- names(factor.names(design)) 17 if (with.blocks) fn <- c(fn, design.info(design)$block.name) 18 design <- design[,fn] 19 nfac <- length(fn) 20 } 21 else{ 22 nfac <- ncol(design) 23 fn <- 1:nfac 24 } 25 26 nlev <- levels.no(design) 27 dfs <- nlev-1 28 ks <- which(round(GWLP(design)[-1],8)>0) 29 if (length(ks)==0) return(list(GRtot=Inf, GRind=Inf, GRind.i=rep(Inf,nfac), cancor1=0)) 30 k <- min(ks) 31 ## 10/08/15: resolution requirement lowered to 2, avoided warning in case of Inf 32 if (k<2) stop("resolution of design must be at least 2") 33 34 ## now there is a finite resolution of at least 2 35 if (!"design" %in% class(design)){ 36 ## make sure these are factors 37 if (!is.data.frame(design)) design <- as.data.frame(design) 38 faktoren <- sapply(design, is.factor) 39 keinfaktor <- which(!faktoren) 40 if (length(keinfaktor)>0) for (i in keinfaktor) design[[i]] <- as.factor(design[[i]]) 41 } 42 nGRindi <- choose(nfac-1,k-1) 43 sel <- nchoosek(nfac-1,k-1) 44 auswahl <- 1:nGRindi 45 erg3 <- array(NA, c(nfac, nGRindi, max(nlev)-1)) 46 selproj <- nchoosek(nfac,k) 47 selproj <- apply(selproj, 2, function(obj) paste(obj, collapse=":")) 48 ergproj <- rep(NA, length(selproj)) 49 dimnames(erg3) <- list(factor=fn, others=apply(sel, 2, function(obj) paste(obj, collapse=":")), 1:(max(nlev)-1)) 50 for (i in 1:nfac){ 51 seli <- sel 52 seli[seli>=i] <- seli[seli>=i]+1 53 proji <- apply(seli, 2, function(obj) paste(sort(c(i,obj)),collapse=":")) 54 indexes <- sapply(proji, function(obj) which(selproj==obj)) 55 berechn <- t(sapply(auswahl, function(obj){ 56 spalten <- c(i,seli[,obj]) 57 hilf2 <- design[,spalten] 58 mmX <- model.matrix(~., data=hilf2[,1,drop=FALSE])[,-1] ## first factor 59 if (k>2) 60 mmOther <- model.matrix(formula(substitute(~.^km1, list(km1 = k-1))), 61 data=hilf2[,-1,drop=FALSE])[,-1] 62 else 63 mmOther <- model.matrix(~., data=hilf2[,-1,drop=FALSE])[,-1] 64 hilfc <- cancor(mmOther, mmX)$cor 65 if (length(hilfc)<dfs[i]) 66 hilfc <- c(hilfc, rep(0,dfs[i]-length(hilfc))) 67 hilfc 68 })) 69 if (nrow(berechn)==1) berechn <- t(berechn) 70 erg3[i,1:nGRindi,1:dfs[i]] <- berechn 71 } 72 rund <- round(erg3, digits) 73 aveR2s <- apply(erg3,c(1,2),function(obj) mean(obj^2, na.rm=TRUE)) 74 if (arft){ 75 ARFT <- table(round(aveR2s,digits)) 76 ARFT <- cbind(aveR2=as.numeric(names(ARFT)), frequency=ARFT) 77 rownames(ARFT) <- rep("",nrow(ARFT)) 78 } 79 if (scft){ 80 SCFT <- table(round(erg3^2, digits)) 81 SCFT <- cbind(SC=as.numeric(names(SCFT)), frequency=SCFT) 82 rownames(SCFT) <- rep("",nrow(SCFT)) 83 } 84 aus <- list(GRs=c(GR=k+1-round(sqrt(max(apply(erg3,c(1,2), 85 function(obj) mean(obj^2, na.rm=TRUE)))),digits), 86 GRind=k+1-max(rund, na.rm=TRUE)), 87 GR.i=rbind(GRtot.i=k+1-round(sqrt(apply(aveR2s,1, max)), digits), 88 GRind.i=k+1-apply(rund, 1, max, na.rm=TRUE))) 89 if (arft) aus <- c(aus, list(ARFT=ARFT)) 90 if (scft) aus <- c(aus, list(SCFT=SCFT)) 91 if (cancors) aus <- c(aus, list(cancors=round(erg3^2, digits))) 92 class(aus) <- c("GRind", class(aus)) 93 aus 94 } 95 96