Nothing
# Morris' OAT sub-routines. (See main file morris.R)
#
# Gilles Pujol 2006
# Modified by Frank Weber (2016): support model functions
# returning a matrix or a 3-dimensional array.
#
# Adapted by Marcelo C. Pereira (2017) to reuse just
# the ee.oat function and its print/plot methods
# as CRAN does not accept linking unexported code
#
ee.oat <- function(X, y) {
# compute the elementary effects for a OAT design
p <- ncol(X)
r <- nrow(X) / (p + 1)
# if(is(y,"numeric")){
if(inherits(y, "numeric")){
one_i_vector <- function(i){
j <- ind.rep(i, p)
j1 <- j[1 : p]
j2 <- j[2 : (p + 1)]
# return((y[j2] - y[j1]) / rowSums(X[j2,] - X[j1,]))
return(solve(X[j2,] - X[j1,], y[j2] - y[j1]))
}
ee <- vapply(1:r, one_i_vector, FUN.VALUE = numeric(p))
ee <- t(ee)
# "ee" is now a (r times p)-matrix.
# } else if(is(y,"matrix")){
} else if(inherits(y, "matrix")){
one_i_matrix <- function(i){
j <- ind.rep(i, p)
j1 <- j[1 : p]
j2 <- j[2 : (p + 1)]
return(solve(X[j2,] - X[j1,],
y[j2, , drop = FALSE] - y[j1, , drop = FALSE]))
}
ee <- vapply(1:r, one_i_matrix,
FUN.VALUE = matrix(0, nrow = p, ncol = dim(y)[2]))
# Special case handling for p == 1 and ncol(y) == 1 (in this case, "ee" is
# a vector of length "r"):
if(p == 1 && dim(y)[2] == 1){
ee <- array(ee, dim = c(r, 1, 1))
}
# Transpose "ee" (an array of dimensions c(p, ncol(y), r)) to an array of
# dimensions c(r, p, ncol(y)) (for better consistency with the standard
# case that "class(y) == "numeric""):
ee <- aperm(ee, perm = c(3, 1, 2))
# } else if(is(y,"array")){
} else if(inherits(y, "array")){
one_i_array <- function(i){
j <- ind.rep(i, p)
j1 <- j[1 : p]
j2 <- j[2 : (p + 1)]
ee_per_3rd_dim <- sapply(1:(dim(y)[3]), function(idx_3rd_dim){
y_j2_matrix <- y[j2, , idx_3rd_dim]
y_j1_matrix <- y[j1, , idx_3rd_dim]
# Here, the result of "solve(...)" is a (p times dim(y)[2])-matrix or
# a vector of length dim(y)[2] (if p == 1):
solve(X[j2,] - X[j1,], y_j2_matrix - y_j1_matrix)
}, simplify = "array")
if(dim(y)[2] == 1){
# Correction needed if dim(y)[2] == 1, so "y_j2_matrix" and
# "y_j1_matrix" have been dropped to matrices (or even vectors, if also
# p == 1):
ee_per_3rd_dim <- array(ee_per_3rd_dim,
dim = c(p, dim(y)[2], dim(y)[3]))
} else if(p == 1){
# Correction needed if p == 1 (and dim(y)[2] > 1), so "y_j2_matrix" and
# "y_j1_matrix" have been dropped to matrices:
ee_per_3rd_dim <- array(ee_per_3rd_dim,
dim = c(1, dim(y)[2], dim(y)[3]))
}
# "ee_per_3rd_dim" is now an array of dimensions
# c(p, dim(y)[2], dim(y)[3]). Assign the corresponding names for the
# third dimension:
if(is.null(dimnames(ee_per_3rd_dim))){
dimnames(ee_per_3rd_dim) <- dimnames(y)
} else{
dimnames(ee_per_3rd_dim)[[3]] <- dimnames(y)[[3]]
}
return(ee_per_3rd_dim)
}
ee <- sapply(1:r, one_i_array, simplify = "array")
# Special case handling if "ee" has been dropped to a vector:
# if(is(ee,"numeric")){
if (inherits(ee, "numeric")){
ee <- array(ee, dim = c(p, dim(y)[2], dim(y)[3], r))
dimnames(ee) <- list(NULL, dimnames(y)[[2]], dimnames(y)[[3]], NULL)
}
# "ee" is an array of dimensions c(p, dim(y)[2], dim(y)[3], r), so it is
# transposed to an array of dimensions c(r, p, dim(y)[2], dim(y)[3]):
ee <- aperm(ee, perm = c(4, 1, 2, 3))
}
return(ee)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.