# R/length2.R In DoE.base: Full Factorials, Orthogonal Arrays and Base Utilities for DoE Packages

#### Documented in length2

```length2 <- function(design, with.blocks=FALSE, J=FALSE){
## function to calculate number of generalized words of length 3
## according to Xu and Wu 2001 Annals

## it might be helpful to locate non-zeros
## this is so far not done
n <- nrow(design)
if (!(is.data.frame(design) | is.matrix(design))) stop("design must be a data frame or a matrix")
if (is.matrix(design)) design <- as.data.frame(design)
if (!"design" %in% class(design)){
for (i in 1:ncol(design)){
design[,i] <- factor(design[,i])
contrasts(design[,i]) <- "contr.XuWu"
}
nlevels <- sapply(as.list(design), function(obj) nlevels(obj))
##names(nlevels) <- colnames(design)
fo <- formula("~.", data=design)
}
else {
di <- design.info(design)
nlevels <- di\$nlevels
if (is.null(nlevels)){
if (length(grep("FrF2",di\$type))>0 | length(grep("pb",di\$type))>0 )
nlevels <- rep(2,length(di\$factor.name))
}
## orthogonal contrasts
design <- change.contr(design, "contr.XuWu")

## if blocked and requested, accomodate blocks
if (with.blocks & !is.null(di\$block.name)){
if (!is.factor(design[[di\$block.name]])) design[[di\$block.name]] <- factor(design[[di\$block.name]])
contrasts(design[[di\$block.name]]) <- "contr.XuWu"
fo <- formula(paste("~",paste(c(di\$block.name,names(di\$factor.names)),collapse="+")), data=design)
nlevels <- c(di\$nblocks,nlevels)
}
else
fo <- formula(paste("~",paste(names(di\$factor.names),collapse="+")), data=design)
}

## create model matrix

mm <- model.matrix(fo,design)
## store column allocation to factors
zuord <- attr(mm, "assign")
## no longer necessary, already done by contr.XuWu
## --> works for imbalanced designs as well
## normalize columns to Euclidean norm sqrt(n)
## mm <- sqrt(n)*mm %*% diag(1/sqrt(colSums(mm^2)))
## remove intercept column
mm <- mm[,-1]
zuord <- zuord[-1]
##colnames(mm) <- zuord   ## takes too long
nfac <- max(zuord)
if (nfac < 2) return(0)

## 2fi columns
## store all in one matrix, faster but requires large matrix
paare <- nchoosek(nfac, 2)
anz <- 0
for (i in 1:ncol(paare))
anz <- anz+prod(nlevels[paare[,i,drop=FALSE]]-1)
vec2 <- rep(NA, anz)
## names(vec2) <- rep("sp",length(vec2))   ## takes too long

zaehl <- 1   ## column accessor

for (i in 1:(nfac - 1)){
icols <- which(zuord==i)
for (j in (i+1):(nfac)){
jcols <- which(zuord==j)
for (a in icols){
for (b in jcols){
vec2[zaehl] <- sum(mm[,a]*mm[,b])
##names(vec2)[zaehl] <- paste(i,j,k,sep=":")
zaehl <- zaehl+1
}}
}
}
if (J) {
names(vec2) <- unlist(apply(paare,2,function(obj) rep(paste((1:nfac)[obj],collapse=":"), prod(nlevels[obj]-1))) )
return(abs(vec2)/n)
}
else
sum(vec2^2)/(n^2)
}
```

## Try the DoE.base package in your browser

Any scripts or data that you put into this service are public.

DoE.base documentation built on Aug. 23, 2018, 1:03 a.m.