# R/driver-utils.R In fido: Bayesian Multinomial Logistic Normal Regression

#### Documented in alralr_arrayalrInvalrInv_arrayclr_arraycreate_default_ilr_basegather_arraysummarise_posterior

```###Logratio functions

#' Create a default ILR base
#' @param D	the number of parts (e.g., number of columns in untransformed data)
#' @return A matrix
#' @export
create_default_ilr_base <- function(D){
qr.Q(qr(create_alr_base(D, D)))
}

#' Compute the ALR of a matrix
#'
#' @param x A matrix where the rows are the samples
#' @param d Index of column used as a reference. Defaults to last column
#'
#' @return matrix
#' @export
alr <- function(x, d=NULL){
x <- vec_to_mat(x)
if (is.null(d)) d <- ncol(x)
B <- create_alr_base(ncol(x), d, inv=FALSE)
glr(x, B)
}

#' Compute the inverse ALR of a matrix
#'
#' @param y An ALR transformed matrix
#' @param d Index of column used as a reference. Defaults to last column
#'
#' @return matrix
#' @export
alrInv <- function(y, d=NULL){
y <- vec_to_mat(y)
if (is.null(d)) d <- ncol(y)+1
B <- create_alr_base(ncol(y)+1, d, inv=TRUE)
glrInv(y, B)
}

create_alr_base <- function(D,d, inv=FALSE){
if (d < 1 | d > D) stop("invalid d given D")
B <- diag(D)
if (!inv){
B[d,] <- rep(-1, D)
} else {
B[d,] <- rep(0, D)
}
B[,-d]
}

glr <- function(x, V){
log(x) %*% V
}

glrInv <- function(y, V){
tmp <- exp(y %*% t(V))
miniclo(tmp)
}

glr_array <- function(x, V, parts, dimname = colnames(V)){
if (parts == 1){
dn <- dimnames(x)
d <- dim(x)
x <- matrix(x, d[1], prod(d[-1]))
x <- t(V) %*% log(x)
d[1] <- ncol(V)
dim(x) <- d
if (!is.null(dn)){
if (!is.null(dimname)){
dn[[1]] <- dimname
} else {
dn[1] <- list(NULL)
}
dimnames(x) <- dn
}
return(x)
} else {
f <- function(x) glr(x, V)
return(array_apply_1D_function(x, parts, f, dimname))
}
}

#' Compute the CLR of an array
#'
#' @param x multidimensional array in index
#' @param parts index of dimension of `x` that represents parts
#'
#' @return array
#' @export
clr_array <- function(x, parts){
n.parts <- dim(x)[parts]
V <- create_clr_base(n.parts)
glr_array(x, V, parts)
}

clr <- function(x){
x <- vec_to_mat(x)
glr(x, create_clr_base(ncol(x)))
}

create_clr_base <- function(D, inv=FALSE){
if(!inv){
M <- matrix(-1, D, D) + D*diag(D)
return(M/D)
} else {
return(diag(D))
}
}

alrInv_array <- function(y, d=dim(y)[coords]+1, coords){
B <- create_alr_base(dim(y)[coords]+1, d, inv=TRUE)
glrInv_array(y, B, coords)
}

glrInv_array <- function(y, V, coords, dimname = rownames(V)){
if (coords==1){
dn <- dimnames(y)
d <- dim(y)
y <- matrix(y, d[1], prod(d[-1]))
y <- exp(V %*% y)
y <- t(miniclo(t(y)))
d[1] <- nrow(V)
dim(y) <- d
if (!is.null(dn)){
if (!is.null(dimname)) { dn[[1]] <- dimname } else {dn[1] <- list(NULL)}
dimnames(y) <- dn
}
return(y)
} else {
f <- function(y) glrInv(y, V)
return(array_apply_1D_function(y, coords, f, dimname))
}
}

clrInv_array <- function(y, coords){
n.coords <- dim(y)[coords]
V <- diag(n.coords) # Not efficient but reuses code...
glrInv_array(y, V, coords)
}

ilr_array <- function(x, V=NULL, parts){
n.parts <- dim(x)[parts]
if (is.null(V)) V <- create_default_ilr_base(n.parts)
glr_array(x, V, parts)
}

ilrInv_array <- function(y, V=NULL, coords){
n.coords <- dim(y)[coords]
if (is.null(V)) V <- qr.Q(qr(create_alr_base(n.coords+1, n.coords+1)))
glrInv_array(y, V, coords)
}

#' Compute the ALR of an array
#'
#' @param x multidimensional array in simplex
#' @param d Index of column used as a reference. Defaults to last column
#' @param parts index of dimension of `x` that represents parts
#'
#' @return array
#' @export
alr_array <- function(x, d=dim(x)[parts], parts){
B <- create_alr_base(dim(x)[parts], d, inv=FALSE)
glr_array(x, B, parts)
}

#' Compute the ALR of an array
#'
#' @param y multidimensional ALR transformed array
#' @param d Index of column used as a reference. Defaults to last column
#' @param coords index of dimension of `x` that represents coordinates
#'
#' @return array
#' @export
alrInv_array <- function(y, d=dim(y)[coords]+1, coords){
B <- create_alr_base(dim(y)[coords]+1, d, inv=TRUE)
glrInv_array(y, B, coords)
}

clrvar2alrvar <- function(Sigma, d2){
D <- nrow(Sigma)
V <- create_alr_base(D, d2, inv=FALSE)
return(t(V)%*%Sigma%*%V)
}

alrvar2clrvar <- function(Sigma, d1){
D <- nrow(Sigma)+1
G1 <- create_alr_base(D, d1, inv=TRUE) - 1/D
return(G1%*%Sigma%*%t(G1))
}

alrvar2alrvar <- function(Sigma, d1, d2){
S <- alrvar2clrvar(Sigma, d1)
clrvar2alrvar(S, d2)
}

alrvar2ilrvar <- function(Sigma, d1, V2){
S <- alrvar2clrvar(Sigma, d1)
clrvar2ilrvar(S, V2)
}

ilrvar2alrvar <- function(Sigma, V1, d2){
S <- ilrvar2clrvar(Sigma, V1)
clrvar2alrvar(S, d2)
}

alrvar2varmat <- function(Sigma, d1){
S <- alrvar2clrvar(Sigma, d1)
clrvar2varmat(S)
}

ilrvar2varmat <- function(Sigma, V){
Sigma <- ilrvar2clrvar(Sigma, V)
clrvar2varmat(Sigma)
}

clrvar2varmat <- function(Sigma){
varmat <- matrix(0, nrow(Sigma), ncol(Sigma))
for (i in 1:dim(Sigma)[1]){
for (j in 1:dim(Sigma)[2]){
varmat[i,j] <- Sigma[i,i] + Sigma[j,j] - 2*Sigma[i,j]
}
}
}

clrvar2ilrvar <- function(Sigma, V){
t(V) %*% Sigma %*% V
}

ilrvar2clrvar <- function(Sigma, V){
V %*% Sigma %*% t(V)
}

ilrvar2ilrvar <- function(Sigma, V1, V2){
t(V2) %*% V1 %*% Sigma %*% t(V1) %*% V2
}

ilrvar2phi <- function(Sigma, V){
Sigma.clr <- ilrvar2clrvar(Sigma,V)
phi <- Sigma.clr
phi[] <- 0
for(i in 1:dim(Sigma.clr)[1]){
for (j in 1:dim(Sigma.clr)[2]){
phi[i,j] <- Sigma.clr[i,i] + Sigma.clr[j,j] - 2*Sigma.clr[i,j]
}
}
return(phi)
}

##Helper functions

#' Gather Multidimensional Array to Tidy Tibble
#'
#' @param a multidimensional array
#' @param value unquoted name of column with values (defaults to "var")
#' @param ... unquoted dimension names (defaults to "dim_1", "dim_2", etc...)
#' @param .id if specified, name for column created with name of a captured
#'
#' @return data.frame
#' @export
#' @import dplyr purrr tidyr
#' @importFrom rlang quos enquo quo_name sym syms
#'
#' @examples
#' a <- array(1:100, dim =c(10, 5, 2))
#' gather_array(a, sequence, A, B, C)
gather_array <- function(a, value, ..., .id=NULL){
qs <- rlang::quos(...)
if (missing(value)) {
evalue <- rlang::sym("var")}
else {
evalue <- rlang::enquo(value)
}
len <- length(qs)
d <- dim(a)

# Default Values
if (len > 0) {
dimnames <- purrr::map(qs, rlang::quo_name) %>%
purrr::as_vector()
} else {
dimnames <- paste0("dim_", 1:length(d))
}

l <- list()
for (i in 1:length(d)){
l[[i]] <- 1:d[i]
}
names(l) <- dimnames
tidy <- expand.grid(l)
tidy[[rlang::quo_name(evalue)]] <- a[as.matrix(tidy)]
if (!is.null(.id)) tidy[[.id]] <- rlang::expr_name(a)
return(tidy)
}

#' Shortcut for summarize variable with quantiles and mean
#'
#' @param data tidy data frame
#' @param var variable name (unquoted) to be summarised
#' @param ... other expressions to pass to summarise
#'
#' @return data.frame
#' @export
#' @details Notation: \code{pX} refers to the \code{X}\% quantile
#' @import dplyr
#' @importFrom stats quantile
#' @importFrom rlang quos quo UQ
#' @examples
#' d <- data.frame("a"=sample(1:10, 50, TRUE),
#'                 "b"=rnorm(50))
#'
#' # Summarize posterior for b over grouping of a and also calcuate
#' # minmum of b (in addition to normal statistics returned)
#' d <- dplyr::group_by(d, a)
#' summarise_posterior(d, b, mean.b = mean(b), min=min(b))
summarise_posterior <- function(data, var, ...){
qvar <- enquo(var)
qs <- quos(...)

data %>%
summarise(p2.5 = quantile(!!qvar, prob=0.025),
p25 = quantile(!!qvar, prob=0.25),
p50 = quantile(!!qvar, prob=0.5),
mean = mean(!!qvar),
p75 = quantile(!!qvar, prob=0.75),
p97.5 = quantile(!!qvar, prob=0.975),
!!!qs)
}

vec_to_mat <- function(x){
if (is.vector(x)) {
n <- names(x)
x <- matrix(x, nrow = 1)
colnames(x) <- n
}
x
}

#' @importFrom rlang :=
array_apply_1D_function <- function(a, dimno, f, dimname=NULL){

d <- dim(a)
ndim <- length(d)
sdim <- sym(paste0("dim_", dimno))
sdim_other <- syms(paste0("dim_", (1:ndim)[(1:ndim) != dimno]))

# Store Dimnames
dn <- dimnames(a)

# Actual Computation
var <- "var"
ga <- a %>%
gather_array(!!var) %>%

indicies <- ga %>%
select(!!!sdim_other)

b <- ga %>%
select(-contains("dim")) %>%
as.matrix() %>%
f %>%
`colnames<-`(., 1:ncol(.)) %>%
as.data.frame() %>%
bind_cols(indicies, .) %>%
gather(!!sdim, var, -contains("dim")) %>%
mutate(!!quo_name(sdim)  := as.integer(!!sdim)) %>%

# Update Dimnames
if (!is.null(dn) || !is.null(dimname)){
if (!is.null(dn)){
dn[dimno] <- list(NULL)
} else {
dn <- list()
for (i in 1:length(d)) dn[[i]] <- NULL
}
if (!is.null(dimname)){
dn[[dimno]] <- dimname
}
dimnames(b) <- dn
} else {
names(dim(b)) <- NULL
}
return(b)
}

dd <- dim(a)
if (d > length(dd)+1) stop("d must be <= length(dim(a))+1")
passed=FALSE
for (i in 1:(length(dd)+1)) {
if (i==d) { ad[i] <- 1; passed <- TRUE }
else if (passed) { ad[i] <- dd[i-1] }
}
}

evalue <- rlang::enquo(value)
qs <- rlang::quos(...)
l <- length(qs)

# Default Values
if (l == 0) {
cn <- colnames(data)
cn <- cn[grepl("dim_", cn)]
# Validation of Defaults
consecutive1 <- strsplit(cn, "_") %>%
purrr::map(~.x[2]) %>%
purrr::map(as.integer) %>%
as_vector()
if (!setequal(consecutive1, 1:length(consecutive1))) {
stop("default dimnetion names must have consecutive integer suffixes")
}
##
cn <- cn[match(consecutive1, 1:length(consecutive1))]
qs <- rlang::syms(cn)
}
if (missing(value)) evalue <- rlang::sym("var")

tidy_dim <- data %>%
select(!!!qs)
unique_dim <- tidy_dim %>%
as.list() %>%
purrr::map(unique)
length_dim <- unique_dim %>%
purrr::map(length)

# Input validation - Must be sequential integers
class_dim <- data %>%
select(!!!qs) %>%
sapply(class)
if (!all(class_dim=="integer")) stop("Dimension indexes must be integers")

consecutive2 <- unique_dim %>%
map2(map(length_dim, ~1:.x), setequal) %>%
as_vector() %>%
all()
if (!consecutive2) stop("Dimension indexes must be consecutive")
#####

a <- array(NA, dim = length_dim)
a[as.matrix(tidy_dim)] <- pull(data, rlang::quo_name(evalue))

return(a)
}

##Closure functions

#' Closure operator
#'
#' @param x vector or matrix (rows are samples, parts are columns) of data in simplex
#'
#' @return x with row entries divided by sum of row (converts vectors to row matricies)
#' @export
#'
#' @examples
#' x <- matrix(runif(30), 10, 3)
#' x <- miniclo(x)
miniclo <- function(x){
if (is.vector(x)) {
n <- names(x)
x <- matrix(x, nrow = 1)
colnames(x) <- n
}
(x/rowSums(x))
}

#' Closure Operation applied to array on margin
#'
#'
#' @param x multidimensional array
#' @param parts index of dimension of \code{x} that represents parts (e.g., compositional variables)
#'
#' @return array
#' @export
#'
#' @examples
#' x <- array(1:100, dim=c(10, 5, 2))
#' miniclo_array(x, parts=2)
miniclo_array <- function(x, parts){
array_apply_1D_function(x, parts, miniclo)
}
```

## Try the fido package in your browser

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

fido documentation built on June 22, 2024, 9:36 a.m.