covls.R 12.1 KB
Newer Older
1
#' @include procmod.R
2
#' @include procmod_frame.R
3
#' @include multivariate.R
4
#' @import foreach
5
#' @import stats
6 7 8 9 10
#'
#' @author Christelle Gonindard-Melodelima
#' @author Eric Coissac
NULL

11 12
.has_doParallel <- is.element("doParallel",installed.packages())
if (.has_doParallel) require(doParallel)
13

14 15
#' Compute the trace of a square matrix.
#'
16 17 18
#' The trace of a square matrix is defined as the sum
#' of its diagonal elements.
#'
19 20 21 22
#' @param X a square matrix
#' @return the trace of X
#'
#' @examples
23 24
#' m <- matrix(1:16, nrow = 4)
#' ProcMod:::.Trace(m)
25 26
#' @note Internal function do not use.
#'
27
#' @rdname internal.Trace
28 29 30
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#'
31
.Trace <- function(X) sum(diag(X))
32

33 34 35 36 37 38 39 40 41 42 43
#' Convert a covariance matrix to a correlation matric
#'
#' @param c a square covariance/variance matrix
#' @return correlation matrix corresponding to \code{c}
#'
#' @note Internal function do not use.
#'
#' @rdname internal.var2cor
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#'
44 45 46 47 48 49
.var2cor <- function(c) {
  v <- sqrt(diag(c))
  vv <- v %o% v
  c / vv
}

Eric Coissac committed
50 51
#' Procrustean Correlation, and Variance / Covariance Matrices.
#'
52
#' \code{varls}, \code{corls}, \code{corls_partial} compute the procrustean
Eric Coissac committed
53 54
#' variance / covariance, correlation, or partial correlation matrices
#' between a set of real matrices and \code{\link[stats]{dist}} objects.
55
#'
56
#' Procrustean covariance between two matrices X and Y, is defined as the sum
Eric Coissac committed
57 58 59 60 61 62 63 64
#' of the singular values of the X'Y matrix \insertCite{Gower:71:00,Lingoes:74:00}{ProcMod}.
#' Both the X and Y matrices must have the same number of rows.
#'
#' The variances and covariances and correlations are corrected
#' to avoid over fitting \insertCite{Coissac-Eric:19:00}{ProcMod}.
#'
#' Partial correlation coefficients are computed by inverting the correlation followed
#' by a normalisation by the diagonal of the inverted matrix.
65
#'
Eric Coissac committed
66 67
#' The inputs must be numeric matrices or \code{\link[stats]{dist}} object.
#' The set of input matrices can be aggregated un a
68
#' \code{\link[ProcMod]{procmod_frame}}.
69
#'
Eric Coissac committed
70
#' Before computing the coefficients, matrices are projected into an
71 72
#' orthogonal space using the \code{\link[ProcMod]{ortho}} function.
#'
Eric Coissac committed
73 74 75 76 77 78 79 80 81
#' The denominator n - 1 is used which gives an unbiased estimator of the
#' (co)variance for i.i.d. observations.
#'
#' Scaling a covariance matrix into a correlation one can be achieved in many ways,
#' mathematically most appealing by multiplication with a diagonal matrix from left
#' and right, or more efficiently by using sweep(.., FUN = "/") twice.
#' The \code{\link[stats]{cov2cor}} function is even a bit more efficient,
#' and provided mostly for didactical reasons.
#'
82 83 84 85 86 87 88 89
#' @references{
#'  \insertRef{Gower:71:00}{ProcMod}
#'
#'  \insertRef{Lingoes:74:00}{ProcMod}
#'
#'  \insertRef{Coissac-Eric:19:00}{ProcMod}
#' }
#'
90
#' @param ...   the set of matrices or a \code{\link[ProcMod]{procmod_frame}}
91
#'              object.
92 93
#' @param nrand number of randomisation used to estimate the mean
#'              covariance observed between two random matrix.
94 95 96
#'              If rand is \code{NULL} or equal to \code{0}, no correction
#'              is estimated and the raw procrustean covariances are
#'              estimated.
97 98
#' @param p_adjust_method the multiple test correction method used
#'              to adjust p values. \code{p_adjust_method}
99 100 101 102 103
#'              belongsone of the folowing values: \code{"holm"},
#'              \code{"hochberg"}, \code{"hommel"}, \code{"bonferroni"},
#'              \code{"BH"}, \code{"BY"}, \code{"fdr"},  \code{"none"}.
#'              The default is,set to \code{"holm"}.
#'
104
#' @return a \code{procmod_varls} object which corresponds to a numeric
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
#'              matrix annotated by several attributes.
#'
#'              The following attribute is always added:
#'
#'              - \code{nrand} an integer value indicating the number of
#'                randomisations used to estimate the mean of the random
#'                covariance.
#'
#'              When \code{nrand} is greater than 0 a couple of attributes
#'              is added:
#'
#'              - \code{rcovls} a numeric matrix containing the estimation
#'                of the mean of the random covariance.
#'
#'              - \code{p.value} a numeric matrix containing the estimations
#'                of the p.values of tests checking that the observed
#'                covariance is larger than the mean of the random covariance.
#'                p.values are corrected for multiple tests according to the
123
#'                method specified by the \code{p_adjust_method} parameter.
124 125
#'
#' @examples
126
#' # Build Three matrices of 3 rows.
127 128 129 130 131
#' A <- simulate_matrix(10,3)
#' B <- simulate_matrix(10,5)
#' C <- simulate_correlation(B,10,r2=0.6)
#'
#' # Computes the variance covariance matrix
132 133
#' varls(A, B, C)
#' varls(A = A, B = B, C = C)
134
#'
135
#' data = procmod_frame(A = A, B = B, C = C)
136
#' varls(data)
Eric Coissac committed
137
#'
138 139 140 141 142 143 144
#' # Computes the correlation matrix
#' corls(data, nrand = 500)
#'
#' # Computes the partial correlation matrix
#' corls_partial(data)
#' corls_partial(data, nrand = 0)
#'
145 146
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
Eric Coissac committed
147
#'
148 149
#' @seealso \code{\link[stats]{p.adjust}}
#'
Eric Coissac committed
150 151 152 153
#' @rdname varls
#' @name varls
#' @aliases varls
#' @aliases corls
154
#' @aliases corls_partial
155
#' @export
156 157
varls <- function(...,
                  nrand = 100,
158
                  p_adjust_method = "holm") {
159 160
  if (!is.null(nrand) && length(nrand) > 1) {
    stop("nrand must be a single numeric value or the NULL value")
161
  }
162

163 164 165
  if (nrand == 0) nrand <- NULL

  xs <- list(...)
166

167 168
  if (length(xs) == 1) {
    x <- xs[[1]]
169
    if (is_procmod_frame(x)) {
170
      xs <- x
171
    } else {
172
      xs <- procmod_frame(x)
173 174 175
    }
  }
  else {
176
    xs <- as_procmod_frame(xs)
177 178
  }

179
  x_names <- names(xs)
180

181
  xs <- ortho(xs)
182

183 184 185
  nx <- length(xs)
  n <- nrow(xs)
  n1 <- n - 1
186

187
  cov_xxs <- Matrix::Matrix(0, nrow = nx, ncol = nx)
188 189

  # Computes the standard Covls covariance matrix
190
  for (i in seq_len(nx))
191 192
    for (j in i:nx)
      cov_xxs[i, j] <- sum(svd(crossprod(xs[[i]],xs[[j]]))$d) / n1
193

194
  cov_xxs <- as.matrix(Matrix::forceSymmetric(cov_xxs, uplo = "U"))
195

196 197 198
  # Computes Covls under null hypothesis
  if (!is.null(nrand)) {
    n_r_greater <- array(0, dim = dim(cov_xxs))
199

200 201 202 203
    v_xs <- vector(mode = "list", nx)
    for (i in seq_len(nx))
      v_xs[[i]] <- var(xs[[i]])

204 205 206 207 208 209
    if (.has_doParallel && getDoParRegistered()) {
        `%dp%` <- `%dopar%`
      }
    else{
        `%dp%` <- `%do%`
      }
Eric Coissac committed
210

211
    s_cov_xxs <- foreach(k = seq_len(nrand),
212
                         .combine = cbind) %dp% {
213
      s1_cov_xxs <- matrix(0, nrow = nx, ncol = nx)
214
      r_xs <- vector(mode = "list", nx)
215 216 217 218 219 220 221 222 223 224 225
      r_ys <- vector(mode = "list", nx)
      for (i in seq_len(nx)) {
        r_xs[[i]] <- MASS::mvrnorm(n,
                                   rep(0,nrow(v_xs[[i]])),
                                   Sigma = v_xs[[i]],
                                   empirical = TRUE)
        r_ys[[i]] <- MASS::mvrnorm(n,
                                   rep(0,nrow(v_xs[[i]])),
                                   Sigma = v_xs[[i]],
                                   empirical = TRUE)
      }
226 227 228

      for (i in seq_len(nx))
        for (j in i:nx)
229
          s1_cov_xxs[i, j] <- sum(svd(crossprod(r_xs[[i]],r_ys[[j]]))$d)
230

231
      s1_cov_xxs  / n1
232

233
                         }
234

235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
    dim(s_cov_xxs) = c(nx, nx, nrand)

    s_cov_xxs <- s_cov_xxs

    r_cov_xxs <- apply(s_cov_xxs,
      MARGIN = c(1,2),
      FUN = mean
    )

    for (k in seq_len(nrand))
      n_r_greater <- n_r_greater + (s_cov_xxs[, , k] >= cov_xxs)

    r_cov_xxs <- as.matrix(Matrix::forceSymmetric(r_cov_xxs, uplo = "U"))
    cov_xxs <- pmax(cov_xxs - r_cov_xxs, 0)

    colnames(r_cov_xxs) <- x_names
    rownames(r_cov_xxs) <- x_names
252

253
    p_values <- n_r_greater / nrand
254

255
    c_p_values <- p.adjust(p_values[upper.tri(p_values, diag = FALSE)],
256
      method = p_adjust_method
257 258 259 260 261 262 263 264 265 266 267
    )

    p_values[upper.tri(p_values, diag = FALSE)] <- c_p_values
    p_values <- as.matrix(Matrix::forceSymmetric(p_values, uplo = "U"))

    colnames(p_values) <- x_names
    rownames(p_values) <- x_names

    attr(cov_xxs, "rcovls") <- as.matrix(r_cov_xxs)
    attr(cov_xxs, "p.value") <- p_values
    attr(cov_xxs, "nrand") <- nrand
268
  }
269

270 271
  colnames(cov_xxs) <- x_names
  rownames(cov_xxs) <- x_names
272

273
  make_subS3Class(cov_xxs, "procmod_varls")
274 275 276
}


Eric Coissac committed
277
#' @rdname varls
278
#' @export
279
corls <- function(..., nrand = 100,
280
                  p_adjust_method = "holm") {
281 282
  cov <- varls(...,
    nrand = nrand,
283
    p_adjust_method = p_adjust_method
284
  )
285
  s <- sqrt(diag(cov))
286 287
  vv <- s %o% s
  rls <- cov / vv
Eric Coissac committed
288
  class(rls) <- "matrix"
289 290


291 292
  if (!is.null(attr(cov, "rcovls"))) {
    attr(rls, "nrand") <- attr(cov, "nrand")
293
  }
294

295
  make_subS3Class(rls, "procmod_corls")
296 297
}

Eric Coissac committed
298
#' @rdname varls
299
#' @export
300
corls_partial <- function(..., nrand = 100) {
301
  rls <- corls(..., nrand = nrand)
302
  C <- solve(rls)
303 304 305
  S <- sqrt(diag(C))
  D <- S %o% S
  rp <- C / D
306

307 308
  attr(rp, "nrand") <- attr(rls, "nrand")
  attr(rp, "rcovls") <- attr(rls, "rcovls")
309

310
  make_subS3Class(rp, "procmod_corls")
311 312
}

313
#' Print procrustean Variance / Covariance Matrix.
314
#'
315
#' @param x a \code{procmod_varls}
316 317 318
#'          object
#' @param ... other parameters passed to other functions
#'
319 320 321 322 323 324 325 326 327 328 329 330
#' @examples
#' # Build Three matrices of 3 rows.
#' A <- simulate_matrix(10,3)
#' B <- simulate_matrix(10,5)
#' C <- simulate_correlation(B,10,r2=0.6)
#'
#' # Computes the variance covariance matrix
#' data <- procmod_frame(A = A, B = B, C = C)
#' v <- varls(data, nrand = 1000)
#'
#' print(v)
#'
331 332
#' @seealso \code{\link[ProcMod]{varls}}
#'
333 334 335
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
336
print.procmod_varls <- function(x, ...) {
337
  class(x) <- "matrix"
338
  attr(x, "nrand") <- NULL
339
  attr(x, "rcovls") <- NULL
340
  attr(x, "p.value") <- NULL
341 342 343
  print(x)
}

344 345 346
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
347
`$.procmod_varls` <- function(x, name) {
348 349 350
  attr(x,name)
}

351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371
#' The Names of the elements of a Variance / Covariance Matrix.
#'
#' Returns the names of the elements associated to a \code{procmod_varls}
#' object.
#'
#' @param x a \code{procmod_varls} object
#'
#' @examples
#' # Build Three matrices of 3 rows.
#' A <- simulate_matrix(10,3)
#' B <- simulate_matrix(10,5)
#' C <- simulate_correlation(B,10,r2=0.6)
#'
#' # Computes the variance covariance matrix
#' data <- procmod_frame(A = A, B = B, C = C)
#' v <- varls(data, nrand = 1000)
#'
#' names(v)
#'
#' @seealso \code{\link[ProcMod]{varls}}
#'
372 373 374
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
375
names.procmod_varls <- function(x) {
Eric Coissac committed
376 377 378 379
  n <- names(attributes(x))
  bn <- grep(pattern = "^(dim|class)",
             x = n)
  n[-bn]
380 381 382
}


383
#' Print a procrustean Correlation Matrix.
384
#'
385
#' @param x a \code{procmod_corls}
386 387 388
#'          object
#' @param ... other parameters passed to other functions
#'
389 390 391 392 393 394 395 396 397 398 399 400
#' @examples
#' # Build Three matrices of 3 rows.
#' A <- simulate_matrix(10,3)
#' B <- simulate_matrix(10,5)
#' C <- simulate_correlation(B,10,r2=0.6)
#'
#' # Computes the correlation matrix
#' data <- procmod_frame(A = A, B = B, C = C)
#' cls <- corls(data, nrand = 1000)
#'
#' print(cls)
#'
401 402
#' @seealso \code{\link[ProcMod]{corls}}
#'
403 404 405
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
406
print.procmod_corls <- function(x, ...) {
407
  class(x) <- "matrix"
408
  attr(x, "nrand") <- NULL
409
  attr(x, "rcovls") <- NULL
410 411
  attr(x, "rcorls") <- NULL
  attr(x, "p.value") <- NULL
412 413
  print(x)
}
414 415 416 417

#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
418
`$.procmod_corls` <- function(x, name) {
419 420 421
  attr(x,name)
}

Eric Coissac committed
422

423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443
#' The Names of the elements of a Correlation Matrix
#'
#' Returns the names of the elements associated to a \code{procmod_corls}
#' object.
#'
#' @param x a \code{procmod_corls} object
#'
#' @examples
#' # Build Three matrices of 3 rows.
#' A <- simulate_matrix(10,3)
#' B <- simulate_matrix(10,5)
#' C <- simulate_correlation(B,10,r2=0.6)
#'
#' # Computes the correlation matrix
#' data <- procmod_frame(A = A, B = B, C = C)
#' cls <- corls(data, nrand = 1000)
#'
#' names(cls)
#'
#' @seealso \code{\link[ProcMod]{corls}}
#'
Eric Coissac committed
444 445 446
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
447
names.procmod_corls <- function(x) {
Eric Coissac committed
448 449 450 451 452 453 454
  n <- names(attributes(x))
  bn <- grep(pattern = "^(dim|class)",
             x = n)
  n[-bn]
}