Added OLS test, FW partialling test

DESCRIPTION

@@ -18,7 +18,6 @@ Suggests:
     tinytest,
     lmtest,
     sandwich,
-    MASS,
     SingleCellExperiment,
     SeuratObject,
     anndataR,

R/robust.R

@@ -11,12 +11,15 @@ ols_beta <- function(X, Y) {
     return(Matrix::Diagonal(x=1/Matrix::diag(X)) %*% Y)
   }
   X <- as.matrix(X)
+  Y <- as.matrix(Y)
   if (!is.null(nrow(Y))) {
     stopifnot(nrow(X)==nrow(Y))
   } else {
     stopifnot(nrow(X)==length(Y))
   }
-  return(r_ols_beta(X, Y))
+  beta <- r_ols_beta(X, Y)
+  dimnames(beta) <- list(colnames(X), colnames(Y))
+  return(beta)
 }
 
 #' Calculate OLS residuals
@@ -30,7 +33,9 @@ ols_resid <- function(X, Y, beta) {
   stopifnot(nrow(X)==nrow(Y))
   stopifnot(ncol(X)==nrow(beta))
   stopifnot(ncol(beta)==ncol(Y))
-  return(r_ols_resid(X, Y, beta))
+  res <- r_ols_resid(X, Y, beta)
+  dimnames(res) <- dimnames(Y)
+  return(res)
 }
 
 #' Calculate HC0 SE per-row

inst/tinytest/test_ols_beta.R

@@ -0,0 +1,22 @@
+set.seed(100)
+
+mat_nr <- c(100, 250, 500, 1000, 2000, 5000)
+mat_nx <- sample(seq(2, length(LETTERS)), length(mat_nr))
+mat_ny <- sample(seq(2, length(LETTERS)), length(mat_nr))
+apply(cbind(mat_nr, mat_nx, mat_ny), 1, function(shape) {
+    X <- matrix(rnorm(shape[["mat_nr"]] * shape[["mat_nx"]]), nrow=shape[["mat_nr"]])
+    colnames(X) <- paste0("X", 1:ncol(X))
+    Y <- matrix(rnorm(shape[["mat_nr"]] * shape[["mat_ny"]]), nrow=shape[["mat_nr"]])
+    colnames(Y) <- paste0("Y", 1:ncol(Y))
+    rownames(Y) <- seq_len(nrow(Y))
+    df <- as.data.frame(cbind(X, Y))
+    models <- lapply(setNames(colnames(Y), colnames(Y)), function(target) {
+        lm(as.formula(paste0(target, "~ 0 + ", paste0(colnames(X), collapse="+"))), df)
+    })
+    beta_real <- sapply(models, coef)
+    beta <- scdemon::ols_beta(X=X, Y=Y)
+    expect_equal(beta, beta_real)
+    res_real <- sapply(models, resid)
+    res <- scdemon::ols_resid(X=X, Y=Y, beta=beta)
+    expect_equal(res, res_real)
+})

inst/tinytest/test_partial.R

@@ -0,0 +1,29 @@
+set.seed(0)
+
+mat_nr <- c(100, 250, 500, 1000, 2000, 5000)
+mat_nc <- c(10, 1000, 400, 250, 555, 100)
+mat_nm <- as.integer(ceiling(sqrt(pmin(mat_nr, mat_nc))))
+mat_nb <- c(0, 1, 5, 10, 3, 50)
+apply(cbind(mat_nr, mat_nc, mat_nm, mat_nb), 1, function(shape) {
+    U <- matrix(rnorm(shape[["mat_nr"]] * shape[["mat_nm"]]), nrow=shape[["mat_nr"]])
+    V <- matrix(rnorm(shape[["mat_nm"]] * shape[["mat_nc"]]), ncol=shape[["mat_nc"]])
+    if (shape[["mat_nb"]] > 0) {
+        ### Make B such that U has a substantial component made from B
+        B <- U %*% matrix(rnorm(shape[["mat_nm"]] * shape[["mat_nb"]]), ncol=shape[["mat_nb"]])
+        B <- B + matrix(rnorm(prod(dim(B))), ncol=ncol(B))
+    } else {
+        B <- matrix(1, nrow=shape[["mat_nr"]])
+    }
+    #U <- U + 5 *  rowMeans(B) ### to add effect
+    ### Test that B is correctly partialled out
+    uvb <- U %*% V - B %*% solve(t(B) %*% B) %*% t(B) %*% U %*% V
+    V1 <- scdemon:::.robust_prepare(U, V, B, return_U=TRUE)
+    U1 <- attr(V1, "U")
+    expect_equal(U1 %*% V1, uvb)
+    ## plot(atanh(cor(V1)), atanh(cor(uvb)))
+### test accuracy
+    ## I = sample(nrow(U), as.integer(ceiling(sqrt(nrow(U)))), replace=FALSE)
+    ## V2 <- scdemon:::.robust_prepare(U[I,], V, B[I,], return_U=TRUE)
+    ## U2 <- attr(V2, "U")
+    ## plot(atanh(cor(V2)), atanh(cor(uvb)))
+})