# R/meanspv.R In vdg: Variance Dispersion Graphs and Fraction of Design Space Plots

#### Documented in expmatmeanspv

#' Compute Mean Spherical SPV
#'
#' Computes the matrix of spherical region moments for a given model formula and a vector of radii, and uses this to
#' calculate the mean spherical SPV for each of the radii. The function \code{expmat} calculates
#' the matrix containing the exponents of each model factor within each model term as columns.
#' Only simple formulae are allowed for. Only products of terms should be included in
#' of \code{\link{sqrt}}. Models should contain only monomial terms.
#'
#' @param formula model formula
#' @param radii numeric vector or radii at which to calculate the matrix of spherical region moments
#' @param FtF.inv inverse of F'F, where F is the design matrix
#' @param n integer giving the number of design runs
# ' @return a matrix of variables by terms containing the exponents of the varaibles in every term
#' @author Pieter C. Schoonees
#' @references
#' Pieter C. Schoonees, Niel J. le Roux, Roelof L.J. Coetzer (2016). Flexible Graphical Assessment of
#' Experimental Designs in R: The vdg Package. \emph{Journal of Statistical Software}, 74(3), 1-22.
#'  \doi{10.18637/jss.v074.i03}.
#' @examples
#' f1 <- formula(~ x1*x2)
#' expmat(f1)
#' f2 <- update(f1, ~ . + I(x1^2) + I(x2^2))
#' expmat(f2)
#' f3 <- update(f2, ~ . + I(x2^0.4))
#' expmat(f3)
#' f4 <- update(f3, ~ . + I(x1^2):I(x2^2))
#' expmat(f4)
#' f5 <- update(f4, ~ . + I(x1^3*x2^0.5))
#' expmat(f5)
meanspv <- function(formula, radii, FtF.inv, n){
expmat <- expmat(formula = formula)
nterms <- ncol(expmat)

sigfun <- function(delta, r){
if(any(delta %% 2 != 0)) return(rep(0, length(r)))
#     if(all(delta == 0)) return(rep(1, length(r)))
m <- length(delta)
sdelta <- sum(delta)
out <- r^sdelta * gamma(m/2) * prod(gamma((delta + 1)/2)) / (pi^(m/2) *
gamma((sdelta + m)/2))
return(out)
}

lowers <- t(apply(combn(nterms, 2), 2, function(x) sigfun(,x[1]] +
diags <- t(apply(2*expmat, 2, sigfun, r = radii))
smom <- mapply(as.data.frame(lowers), as.data.frame(diags),
FUN = function(x, y){tmp <- matrix(0, length(y), length(y));
tmp[lower.tri(tmp)] <- x;
tmp <- tmp + t(tmp);
diag(tmp) <- y;
dimnames(tmp) <- rep(list(colnames(expmat)), 2);
return(tmp)
}, SIMPLIFY = FALSE, USE.NAMES = FALSE)
mspv <- mapply(FUN = function(x, y) n * sum(diag(x %*% y)), smom,
list(FtF.inv))
class = c("meanspv", "list"))
return(out)
}
#' @rdname meanspv
#' @export
expmat <- function(formula){
vars <- all.vars(formula)
terms <- terms(formula)
factors <- attr(terms, "factors")
test <- tnames <- colnames(factors)

for(i in vars) test <- gsub(i, replacement = "", x = test)
test <- gsub("[[:digit:]^()I.:* -]", "", test)
if(!all(test == ""))
stop("Characters other that variable names and the regexp class '[[:digit:]^()I.:* -]' found in formula")

varrows <- na.omit(match(vars, rownames(factors)))
Irows <- seq_len(nrow(factors))[-varrows]
Inames <- rownames(factors)[Irows]
expmat <- factors[varrows, , drop = FALSE]
expmat[expmat == 2] <- 1

getExp <- function(term, var){
if(!grepl(var, term)) return(0L)
sterm <- gsub("^I\\(", "", term)
sterm <- gsub(")\$", "", sterm)
sterm <- unlist(strsplit(sterm, "*", fixed = TRUE))
sterm <- sterm[grepl(var, sterm)]
if(length(sterm) != 1) stop("A variable should only occur once in each term")
tmatch <- match(var, sterm)
if(!is.na(tmatch)) if(tmatch == 1) return(1L)
sterm <- gsub("[()]", "", sterm)
sterm <- unlist(strsplit(sterm, "^", fixed = TRUE))
if(length(sterm) != 2) stop("Unkown use of '^' operator")
return(as.numeric(sterm[2]))
}

whichI <- match(Inames, tnames)
for(i in whichI) for(j in seq_along(vars[varrows]))
expmat[j, i] <- getExp(tnames[i], vars[varrows][j])
whichInt <- setdiff(seq_along(tnames), match(rownames(factors), tnames))
whichIntI <- whichInt[grep("I\\(", tnames[whichInt])]
nrInt <- length(whichIntI)
if(nrInt){
sint <- strsplit(tnames[whichIntI], ":", fixed = TRUE)
intPos <- lapply(sint, match, Inames)
for(i in seq_along(intPos)){
index <- intPos[[i]][!is.na(intPos[[i]])]
, whichIntI[i]] <- , whichIntI[i]] +
rowSums(, Inames[index], drop = FALSE])
}
}
if(attr(terms, "intercept")) expmat <- cbind("(Intercept)" = 0, expmat)
return(expmat)
}

## Try the vdg package in your browser

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

vdg documentation built on July 8, 2022, 1:08 a.m.