Add p-value computation

diff --git a/pkg/R/plot.R b/pkg/R/plot.R
index 93e9e65..b8024f5 100644 (file)
--- a/pkg/R/plot.R
+++ b/pkg/R/plot.R
@@ -8,8 +8,7 @@
 #
 .extractParam <- function(mr, x=1, y=1)
 {
 #
 .extractParam <- function(mr, x=1, y=1)
 {

-  if (is.list(mr[[1]]))
-  {
+  if (is.list(mr[[1]])) {

     # Obtain L vectors where L = number of res lists in mr
     return ( lapply( mr, function(mr_list) {
       sapply(mr_list, function(m) m[x,y])
     # Obtain L vectors where L = number of res lists in mr
     return ( lapply( mr, function(mr_list) {
       sapply(mr_list, function(m) m[x,y])


@@ -47,8 +46,7 @@ plotHist <- function(mr, x, y, ...)
   # Plot histograms side by side
   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
   args <- list(...)
   # Plot histograms side by side
   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
   args <- list(...)

-  for (i in 1:L)
-  {
+  for (i in 1:L) {

     hist(params[[i]], breaks=40, freq=FALSE,
       xlab=ifelse("xlab" %in% names(args), args$xlab, "Parameter value"),
       ylab=ifelse("ylab" %in% names(args), args$ylab, "Density"))
     hist(params[[i]], breaks=40, freq=FALSE,
       xlab=ifelse("xlab" %in% names(args), args$xlab, "Parameter value"),
       ylab=ifelse("ylab" %in% names(args), args$ylab, "Density"))


@@ -85,11 +83,8 @@ plotBox <- function(mr, x, y, ...)
   L <- length(params)
   # Plot boxplots side by side
   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
   L <- length(params)
   # Plot boxplots side by side
   par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))

-  args <- list(...)

   for (i in 1:L)
   for (i in 1:L)

-  {

     boxplot(params[[i]], ...)
     boxplot(params[[i]], ...)

-  }

 }
 
 #' plotCoefs
 }
 
 #' plotCoefs


@@ -124,10 +119,8 @@ plotCoefs <- function(mr, params, ...)
 
   params_hat <- matrix(nrow=d, ncol=K)
   stdev <- matrix(nrow=d, ncol=K)
 
   params_hat <- matrix(nrow=d, ncol=K)
   stdev <- matrix(nrow=d, ncol=K)

-  for (x in 1:d)
-  {
-    for (y in 1:K)
-    {
+  for (x in 1:d) {
+    for (y in 1:K) {

       estims <- .extractParam(mr, x, y)
       params_hat[x,y] <- mean(estims)
       # Another way to compute stdev: using distances to true params
       estims <- .extractParam(mr, x, y)
       params_hat[x,y] <- mean(estims)
       # Another way to compute stdev: using distances to true params

diff --git a/pkg/R/utils.R b/pkg/R/utils.R
index 5b9d999..e39119e 100644 (file)
--- a/pkg/R/utils.R
+++ b/pkg/R/utils.R
@@ -18,6 +18,33 @@ normalize <- function(x)
   sweep(x, 2, norm2, '/')
 }
 
   sweep(x, 2, norm2, '/')
 }
 

+#' pvalue
+#'
+#' Compute the p-values of the tests beta[i,j] == 0 vs != 0
+#'
+#' @param mr A list of matrices as output by multiRun()
+#'
+#' @return The matrix of p-values (same size as mr[[1]])
+#'
+#' @examples
+#' mr <- multiRun(...) #cf ?multiRun
+#' p <- pvalue(mr[[1]])
+#' @export
+pvalue <- function(mr)
+{
+  n = nrow(mr[[1]])
+  m = ncol(mr[[1]])
+  pval <- matrix(nrow=n, ncol=m)
+  for (x in 1:n) {
+    for (y in 1:m) {
+      coefs <- sapply(mr, function(m) m[x,y])
+      stat_test <- mean(coefs^2) / var(coefs)
+      pval[x, y] <- 1 - pchisq(stat_test, 1)
+    }
+  }
+  pval
+}
+

 # Computes a tensor-vector product
 #
 # @param Te third-order tensor (size dxdxd)
 # Computes a tensor-vector product
 #
 # @param Te third-order tensor (size dxdxd)