1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220 | #' Assess variable importance
#'
#' Assess importance of variables on a given PC component by computing the
#' Crone-Crosby distance between original sample positions and sample positions
#' obtained by a random permutation of the variables.
#'
#' @details
#' \code{plotVar.kernel.pca} produces a barplot for each block. The variables
#' for which the importance has been computed with
#' \code{\link{kernel.pca.permute}} are displayed. The representation is limited
#' to the \code{ndisplay} most important variables.
#'
#' @param kpca.result a kernel.pca object returned by the
#' \code{\link{kernel.pca}} function.
#' @param ncomp integer. Number of KPCA components used to compute the
#' importance. Default: \code{1}.
#' @param ... list of character vectors. The parameter name must be the kernel
#' name to be considered for permutation of variables. Provided vectors length
#' has to be equal to the number of variables of the input dataset. A kernel is
#' performed on each unique variables values. Crone-Crosby distances are
#' computed on each KPCA performed on resulted kernels or meta-kernels and can
#' be displayed using the \code{\link{plotVar.kernel.pca}}.
#' @param directory character. To limit computational burden, this argument
#' allows to store / read temporary computed kernels.
#'
#' @return \code{kernel.pca.permute} returns a copy of the input
#' \code{kpca.result} results and add values in the three entries:
#' \code{cc.distances}, \code{cc.variables} and \code{cc.blocks}.
#'
#' @author Jerome Mariette <jerome.mariette@@inrae.fr>
#' Nathalie Vialaneix <nathalie.vialaneix@@inrae.fr>
#' @references Mariette J. and Villa-Vialaneix N. (2018). Unsupervised multiple
#' kernel learning for heterogeneous data integration. \emph{Bioinformatics},
#' \bold{34}(6), 1009-1015. DOI: \doi{10.1093/bioinformatics/btx682}
#'
#' Crone L. and Crosby D. (1995). Statistical applications of a metric on
#' subspaces to satellite meteorology. \emph{Technometrics}, \bold{37}(3),
#' 324-328.
#' @seealso \code{\link{compute.kernel}}, \code{\link{kernel.pca}}
#' @export
#' @examples
#' data(TARAoceans)
#'
#' # compute one kernel for the psychem dataset
#' phychem.kernel <- compute.kernel(TARAoceans$phychem, kernel.func = "linear")
#' # perform a KPCA
#' kernel.pca.result <- kernel.pca(phychem.kernel)
#'
#' # compute importance for all variables in this kernel
#' kernel.pca.result <- kernel.pca.permute(kernel.pca.result,
#' phychem = colnames(TARAoceans$phychem))
#'
kernel.pca.permute = function(kpca.result, ncomp = 1, ..., directory = NULL) {
#-- checking general input parameters --#
permutations <- list(...)
blocks <- names(permutations)
if (!inherits(kpca.result, "kernel.pca")) {
stop(paste0("'kpca.result' should be an instance of 'kernel.pca' object",
" returned by the kernel.pca function."), call. = FALSE)
}
if (length(blocks) == 0) {
stop("At least permutation variables of 1 block should be provided.",
call. = FALSE)
}
# it the kpca has been performed on a meta kernel
is.metakernel <- inherits(kpca.result$kernel, "metaKernel")
if (is.metakernel) {
# test if there is no identity kernels, i.e. kernel with X = NULL
sapply(blocks, function(b.label) {
if (is.null(kpca.result$kernel$X[[b.label]]$X)) {
stop(paste0("No permutation can be done on block '", b.label,
"' as there is no dataset available for this block in the",
" 'kpca.result' object."), call. = FALSE)
}
})
sapply(blocks, function(b.label) {
if (!b.label %in% names(kpca.result$kernel$X)) {
stop(paste0("Block '", b.label, "' does not exists as a block in the ",
"'kpca.result' object."), call. = FALSE)
}
})
sapply(blocks, function(b.label) {
if (!ncol(kpca.result$kernel$X[[b.label]]$X) == length(permutations[[b.label]])) {
stop(paste0("'permutations' vector length for block '", b.label,
"' should be equal to the number of variables in '",
b.label, "' dataset: ",
ncol(kpca.result$kernel$X[[b.label]]$X), "."), call. = FALSE)
}
})
} else {
# test if there its a kidentity kernel, i.e. kernel with X = NULL
if (is.null(kpca.result$kernel$X)) {
stop(paste0("No permutation can be done as no dataset is available in ",
"'kpca.result' object."), call. = FALSE)
}
if (length(blocks) > 1) {
stop(paste0("Only 1 'permutations' vector should be provided as the ",
"kernel.pca has been performed on a single kernel."),
call. = FALSE)
}
if (!ncol(kpca.result$kernel$X) == length(permutations[[1]])) {
stop(paste0("'permutations' vector length should be equal to the number ",
"of variables in the input dataset: ",
ncol(kpca.result$kernel$X), "."), call. = FALSE)
}
}
if (!is.null(directory) && !file.exists(file.path(directory))) {
stop(paste0("directory '", directory, "' does not exist."), call. = FALSE)
}
if (ncomp > kpca.result$ncomp) {
ncomp <- kpca.result$ncomp
}
#-- make permutations ------------------#
all.pdist <- all.variable <- all.block <- c()
for (block in blocks) {
if (is.metakernel) {
K <- kpca.result$kernel$X[[block]]
} else {
K <- kpca.result$kernel
}
for (variable in as.vector(unique(permutations[[block]]))) {
if (!is.na(variable)) {
if (!variable %in% kpca.result$cc.variables) {
# permutes values of variable
data.permuted <- K$X
for (col in which(permutations[[block]] == variable)) {
data.permuted[, col] <- sample(data.permuted[, col], replace = FALSE)
}
# if the directory parameter is set, try to load the permuted kernel
f.path <- file.path(directory, paste0(variable, ".rds"))
if (!is.null(directory) && file.exists(f.path)) {
kernel.permuted <- readRDS(f.path)
} else {
# recompute the kernel
compute.args <- K$kernel.args
compute.args$X <- data.permuted
compute.args$kernel.func <- K$kernel.func
kernel.permuted <- do.call("compute.kernel", compute.args)
if (!is.null(directory)) {
saveRDS(kernel.permuted, f.path)
}
}
if (is.metakernel) {
# change the considered kernel
all.kernels <- kpca.result$kernel$X
all.kernels[[block]] <- kernel.permuted
all.kernels.scaled <- lapply (all.kernels, function(x) {
x.cosinus <- sweep(sweep(x$kernel, 2, sqrt(diag(x$kernel)), "/"), 1,
sqrt(diag(x$kernel)), "/")
t(t(x.cosinus - colSums(x.cosinus) / nrow(x.cosinus)) -
rowSums(x.cosinus) / nrow(x.cosinus)) + sum(x.cosinus) /
nrow(x.cosinus)^2
})
# compute the new meta.kernel with permuted variables
meta.kernel.permuted <- lapply(as.list(1:length(all.kernels.scaled)),
function(x) {
all.kernels.scaled[[x]] * kpca.result$kernel$weights[x]
})
meta.kernel.permuted <- Reduce("+", meta.kernel.permuted)
kernel.tmp <- list(kernel = meta.kernel.permuted,
weights = kpca.result$kernel$weights)
class(kernel.tmp) <- c("kernel", "metaKernel")
} else {
kernel.tmp <- kernel.permuted
}
# re performe the KPCA
kpca.permuted <- kernel.pca(kernel.tmp, ncomp = kpca.result$ncomp)
# compute the crone and crosby distance
pdist <- Pdist(list(O2P(kpca.result$x[,1:ncomp]),
O2P(kpca.permuted$x[,1:ncomp])),
weights = "constant")
all.pdist <- c(all.pdist, pdist)
all.variable <- c(all.variable, variable)
all.block <- c(all.block, block)
} else {
stop(paste0("'", variable, "'", "variable has already been permuted",
" for this kernel.pca object!"))
}
}
}
}
kpca.result$cc.distances <- c(kpca.result$cc.distances, all.pdist)
kpca.result$cc.variables <- c(kpca.result$cc.variables, all.variable)
kpca.result$cc.blocks <- c(kpca.result$cc.blocks, all.block)
return(invisible(kpca.result))
}
|