Skip to content

Commit

Permalink
Passing tests and examples
Browse files Browse the repository at this point in the history
  • Loading branch information
antagomir committed Nov 17, 2024
1 parent 4a2b84d commit 78daefe
Show file tree
Hide file tree
Showing 6 changed files with 48 additions and 35 deletions.
42 changes: 20 additions & 22 deletions R/decostand.R
Original file line number Diff line number Diff line change
Expand Up @@ -217,16 +217,11 @@
xx <- clog - means
attr(xx, "parameters") <- list("means" = means)

# Add the matrix completion step if na.rm=TRUE
# Replace NAs with zeroes if na.rm=TRUE
# Otherwise return the transformation with NAs
# (this is needed in some downstream methods)
# (to be completed separately)
if (na.rm && any(is.na(xx))) {
dimnams <- dimnames(xx)
# xx <- filling::OptSpace(xx)$X
# run OptSpace and pick the completed matrix
xx <- OptSpace(xx, ropt=ropt, niter=niter, tol=tol, verbose=verbose)$X
# add back dim names
dimnames(xx) <- dimnams
xx[is.na(xx)] <- 0
}

return(xx)
Expand Down Expand Up @@ -309,29 +304,32 @@
}


OptSpace <- function(A, ropt=3, niter=5, tol=1e-5, verbose=FALSE){
OptSpace <- function(x, ropt=3, niter=5, tol=1e-5, verbose=FALSE){

## Preprocessing : A : partially revelaed matrix
if (!is.matrix(A)){
stop("* OptSpace : an input A should be a matrix.")
## Preprocessing : x : partially revealed matrix
if (is.data.frame(x)) {
x <- as.matrix(x)
}
if (!is.matrix(x)){
stop("* OptSpace : an input x should be a matrix.")
}
if (any(is.infinite(A))){
stop("* OptSpace : no infinite value in A is allowed.")
if (any(is.infinite(x))){
stop("* OptSpace : no infinite value in x is allowed.")
}
if (!any(is.na(A))){
if (!any(is.na(x))){
stop("* OptSpace : there is no unobserved values as NA.")
}
idxna <- (is.na(A))
M_E <- array(0, c(nrow(A), ncol(A)))
M_E[!idxna] <- A[!idxna]
idxna <- (is.na(x))
M_E <- array(0, c(nrow(x), ncol(x)))
M_E[!idxna] <- x[!idxna]

## Preprocessing : size information
n <- nrow(A)
m <- ncol(A)
n <- nrow(x)
m <- ncol(x)

## Preprocessing : other sparse-related concepts
nnZ.E <- sum(!idxna)
E <- array(0,c(nrow(A), ncol(A))); E[!idxna] <- 1
E <- array(0,c(nrow(x), ncol(x))); E[!idxna] <- 1
eps <- nnZ.E/sqrt(m*n)

## Preprocessing : ropt : implied rank
Expand All @@ -346,7 +344,7 @@ OptSpace <- function(A, ropt=3, niter=5, tol=1e-5, verbose=FALSE){
} else {
r <- round(ropt)
if ((!is.numeric(r)) || (r<1) || (r>m) || (r>n)){
stop("* OptSpace: ropt should be an integer in [1,min(nrow(A),ncol(A))].")
stop("* OptSpace: ropt should be an integer in [1,min(nrow(x),ncol(x))].")
}
}

Expand Down
8 changes: 8 additions & 0 deletions R/vegdist.R
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,14 @@
x <- decostand(x, "clr", ...) # dots to pass possible pseudocount
if (method == 22) # robust.aitchison
x <- decostand(x, "rclr", ...) # No pseudocount for rclr
# Matrix completion with OptSpace
if (anyNA(x)) {
if (na.rm) {
x <- OptSpace(x, ...)$X
} else {
stop("missing values or zeroes are not allowed with argument 'na.rm = FALSE'")
}
}
if (binary)
x <- decostand(x, "pa")
N <- nrow(x)
Expand Down
11 changes: 5 additions & 6 deletions man/OptSpace.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -16,16 +16,15 @@ relative root mean square error (RMSE) \deqn{RMSE \le
C(\alpha)\sqrt{nr/|E|}}.
}
\usage{
OptSpace(x, ropt=3, niter=5, tol=1e-5, verbose=FALSE, ...)
OptSpace(x, ropt=3, niter=5, tol=1e-5, verbose=FALSE)
}
\arguments{
\item{A}{An \eqn{(n\times m)} matrix whose missing entries should be flagged as NA.}
\item{rop}{\code{NA} to guess the rank, or a positive integer as a pre-defined rank.}
\item{x}{An \eqn{(n\times m)} matrix whose missing entries should be flagged as NA.}
\item{ropt}{\code{NA} to guess the rank, or a positive integer as a pre-defined rank.}
\item{niter}{Maximum number of iterations allowed.}
\item{tol}{Stopping criterion for reconstruction in Frobenius norm.}
\item{verbose}{a logical value; \code{TRUE} to show progress, \code{FALSE} otherwise.}
\item{\dots}{Other arguments to the function (ignored).}
\item{verbose}{a logical value; \code{TRUE} to show progress, \code{FALSE} otherwise.}
}
Expand Down Expand Up @@ -69,7 +68,7 @@ data(varespec)
# rclr transformation with no matrix completion for the 0/NA entries
x <- decostand(varespec, method="rclr", na.rm=FALSE)
# Matrix completion for the missing values
OptSpace(x, ropt=3, niter=5, tol=1e-5, verbose=FALSE)
res <- OptSpace(x, ropt=3, niter=5, tol=1e-5, verbose=FALSE)
}
\keyword{ multivariate}
\keyword{ manip }
2 changes: 1 addition & 1 deletion man/decostand.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ sptrans <- wisconsin(varespec)
varespec.clr <- decostand(varespec, "clr", pseudocount=1)
# Robust CLR transformation for rows, no pseudocount necessary
varespec.rclr <- decostand(varespec, "rclr", na.rm=TRUE)
varespec.rclr <- decostand(varespec, "rclr")
# ALR transformation for rows, with pseudocount and reference sample
varespec.alr <- decostand(varespec, "alr", pseudocount=1, reference=1)
Expand Down
7 changes: 4 additions & 3 deletions man/vegdist.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -419,8 +419,9 @@ vare.dist <- vegdist(decostand(varespec, "norm"), "euclidean")
vare.dist <- vegdist(decostand(varespec, "log"), "altGower")
# Range standardization with "altGower" (that excludes double-zeros)
vare.dist <- vegdist(decostand(varespec, "range"), "altGower")
# Robust Aitchison distance equals to Euclidean distance for rclr transformed data
vare.dist <- vegdist(decostand(varespec, "rclr"), "euclidean")
vare.dist <- vegdist(varespec, "robust.aitchison")
# Robust Aitchison distance equals to Euclidean distance for rclr transformed data
# when matrix completion is used (rm.na=TRUE)
vare.dist <- vegdist(decostand(varespec, "rclr"), method="euclidean", na.rm=TRUE)
vare.dist <- vegdist(varespec, "robust.aitchison", na.rm=TRUE)
}
\keyword{ multivariate }
13 changes: 10 additions & 3 deletions tests/aitchison-tests.R
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,19 @@ relative <- decostand(testdata, "total")
testdata.with.pseudo <- testdata + 1
relative.with.pseudo <- decostand(testdata+1, "total")

# NAs correcly assigned
a1 <- testdata
a2 <- decostand(a1, "rclr", na.rm=FALSE)
a3 <- decostand(a1, "rclr", na.rm=TRUE)
all(is.na(a2[a1==0]))
all(!is.na(a3[a1==0]))

# aitchison equals to CLR + Euclid (pseudocount is necessary with clr)
a1 <- vegdist(testdata+1, method = "aitchison", na.rm=TRUE)
a2 <- vegdist(decostand(testdata+1, "clr"), method = "euclidean", na.rm=TRUE)
max(abs(a1-a2)) < 1e-6 # Tolerance

# Robust aitchison equals to rCLR + Euclid
# Robust aitchison equals to rCLR + Euclid when na.rm=TRUE (matrix completion applied)
# and works without pseudocount
a1 <- vegdist(testdata, method = "robust.aitchison", na.rm=TRUE)
a2 <- vegdist(decostand(testdata, "rclr"), method = "euclidean", na.rm=TRUE)
Expand All @@ -44,8 +51,8 @@ d2 <- vegdist(testdata, "robust.aitchison", na.rm=TRUE)
max(abs(d1-d2))<1e-6 # Tolerance

# Robust rclr+Euclid should give same result than robust.aitchison
d1 <- vegdist(decostand(testdata, "rclr"), "euclidean", na.rm=FALSE)
d2 <- vegdist(testdata, "robust.aitchison", na.rm=FALSE)
d1 <- vegdist(decostand(testdata+1, "rclr"), "euclidean", na.rm=FALSE)
d2 <- vegdist(testdata+1, "robust.aitchison", na.rm=FALSE)
max(abs(d1-d2))<1e-6 # Tolerance

# Compare the outcomes with an external package that also provides compositional transformations
Expand Down

0 comments on commit 78daefe

Please sign in to comment.