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