added test for continuous. #5

R/generalizedPCA.R

@@ -36,10 +36,10 @@
 #'    the null model estimates 0 for all natural parameters.}
 #' @export
 generalizedPCA <- function(x, k = 2, M = 4, family = c("gaussian", "binomial", "poisson"),
-                        majorizer = c("row", "all"), weights,
-                        quiet = TRUE, use_irlba = FALSE,
-                        max_iters = 1000, conv_criteria = 1e-5, random_start = FALSE,
-                        start_U, start_mu, main_effects = TRUE) {
+                           majorizer = c("row", "all"), weights,
+                           quiet = TRUE, use_irlba = FALSE,
+                           max_iters = 1000, conv_criteria = 1e-5, random_start = FALSE,
+                           start_U, start_mu, main_effects = TRUE) {
   use_irlba = use_irlba && requireNamespace("irlba", quietly = TRUE)
 
   family = match.arg(family)
@@ -111,7 +111,7 @@ generalizedPCA <- function(x, k = 2, M = 4, family = c("gaussian", "binomial", "
     cat("0 hours elapsed\n")
   }
 
-  for (m in 1:max_iters) {
+  for (m in seq_len(max_iters)) {
     last_U = U
     last_M = M
     last_mu = mu

README.md

@@ -1,6 +1,6 @@
 # Generalized PCA
 
-`generalizedPCA` is an R package which extends principal component analysis to other types of data, much like [generalized linear models](http://en.wikipedia.org/wiki/Generalized_linear_model) extends linear regression. The package, [logisticPCA](https://github.com/andland/logisticPCA), contains the extension to binary data, among other methods, and this package intends to generalize it to all exponential family distributions. Please note that it is still in the very early stages of development and the conventions will possibly change in the future.
+`generalizedPCA` is an R package which extends principal component analysis to other types of data, much like [generalized linear models](http://en.wikipedia.org/wiki/Generalized_linear_model) extends linear regression. The package [logisticPCA](https://github.com/andland/logisticPCA) contains the extension to binary data, among other methods, and this package intends to generalize it to all exponential family distributions. Please note that it is still in the very early stages of development and the conventions will possibly change in the future.
 
 ## Installation
 

tests/testthat/test-continuous-basics.R

@@ -0,0 +1,80 @@
+context("Continuous")
+
+# construct a low rank matrix in the natural parameter scale
+rows = 100
+cols = 10
+k = 1
+set.seed(1)
+mat_np = outer(rnorm(rows), rnorm(cols))
+
+mat = matrix(rnorm(rows * cols, c(mat_np)), rows, cols)
+
+gpca = generalizedPCA(mat, k = k, M = 4, family = "gaussian", main_effects = TRUE)
+
+pca_mu = colMeans(mat)
+pca = svd(scale(mat, center = pca_mu, scale = FALSE))
+gpca2 = gpca
+gpca$mu = pca_mu
+gpca$U = matrix(pca$v[, 1:k], cols, k)
+# expect_equal(gpca$mu, pca_mu)
+# expect_equal(pca$v[, 1], gpca$U[, 1])
+
+pred1 = predict(gpca, mat)
+pred1l = predict(gpca, mat, type = "link")
+pred1r = predict(gpca, mat, type = "response")
+fit1l = fitted(gpca, type = "link")
+fit1r = fitted(gpca, type = "response")
+
+pred2l = predict(gpca2, mat, type = "link")
+pred2r = predict(gpca2, mat, type = "response")
+
+test_that("correct classes", {
+  expect_is(gpca, "gpca")
+
+  expect_is(pred1, "matrix")
+  expect_is(pred1l, "matrix")
+  expect_is(pred1r, "matrix")
+  expect_is(fit1l, "matrix")
+  expect_is(fit1r, "matrix")
+})
+
+test_that("k = 1 dimensions", {
+  expect_equal(dim(gpca$U), c(cols, 1))
+  expect_equal(dim(gpca$PCs), c(rows, 1))
+  expect_equal(length(gpca$mu), cols)
+
+  expect_equal(dim(pred1), c(rows, 1))
+  expect_equal(dim(pred1l), c(rows, cols))
+  expect_equal(dim(pred1r), c(rows, cols))
+  expect_equal(dim(fit1l), c(rows, cols))
+  expect_equal(dim(fit1r), c(rows, cols))
+})
+
+test_that("k = 1 same fits as standard PCA", {
+  expect_equal(pred1r, pred2r)
+  expect_equal(pred1l, pred2l)
+  # Don't compare fitted, because they rely on $PCs
+})
+
+rm(gpca, pred1, pred1l, pred1r, fit1l, fit1r)
+
+k = 2
+gpca = generalizedPCA(mat, k = k, M = 4, family = "gaussian", main_effects = TRUE)
+
+pred1 = predict(gpca, mat)
+pred1l = predict(gpca, mat, type = "link")
+pred1r = predict(gpca, mat, type = "response")
+fit1l = fitted(gpca, type = "link")
+fit1r = fitted(gpca, type = "response")
+
+test_that("k = 2 dimensions", {
+  expect_equal(dim(gpca$U), c(cols, 2))
+  expect_equal(dim(gpca$PCs), c(rows, 2))
+  expect_equal(length(gpca$mu), cols)
+
+  expect_equal(dim(pred1), c(rows, 2))
+  expect_equal(dim(pred1l), c(rows, cols))
+  expect_equal(dim(pred1r), c(rows, cols))
+  expect_equal(dim(fit1l), c(rows, cols))
+  expect_equal(dim(fit1r), c(rows, cols))
+})