R 与 MATLAB 中的高维数据结构化方法
Posted
技术标签:
【中文标题】R 与 MATLAB 中的高维数据结构化方法【英文标题】：Methodology of high-dimensional data structuring in R vs. MATLAB 【发布时间】：2014-01-31 04:37:12 【问题描述】：问题
在 R 中使用在重复试验中积累的分类标签来构建多元数据以进行探索性分析的正确方法是什么？我不想溜回 MATLAB。
说明
我比 MATLAB 更喜欢 R 的分析函数和语法（以及令人惊叹的图表），并且一直在努力重构我的东西。但是，我一直对数据在工作中的组织方式感到困惑。
MATLAB
对我来说，使用在多次试验中重复的多元时间序列是典型的工作，这些时间序列存储在一个大的 matrix rank-3 张量 SERIESxSAMPLESxTRIALS 的多维数组中。这偶尔适用于一些不错的线性代数，但是当涉及到另一个变量，即 CLASS 时，它就显得笨拙了。通常类标签存储在另一个维度为 1xTRIALS 的向量中。
在分析方面，我基本上尽可能少地绘制，因为需要大量的工作才能拼凑出一个非常好的绘图，它可以教会你很多关于 MATLAB 中的数据的知识。 (I'm not the only one who feels this way)。
R
在 R 中，我一直尽可能地接近 MATLAB 结构，但是当试图将类标签分开时，事情变得非常复杂；即使我只使用它们的属性，我也必须继续将标签传递给函数。所以我所做的就是将数组按类分成数组列表。这增加了我所有 apply() 函数的复杂性，但在保持一致（以及排除错误）方面似乎是值得的。
另一方面，R 似乎对张量/多维数组并不友好。为了与他们合作，您需要获取 abind 库。 multivariate analysis, like this example 上的文档似乎是在假设您有一个巨大的 2-D 数据点表（如 一些中世纪长卷轴 一个数据框）的假设下运行的，并且没有提到如何获得“那里”我在哪里。
一旦我开始对处理后的数据进行绘图和分类，这并不是什么大问题，因为到那时我已经努力开发出具有类似 TRIALSxFEATURES 形状的数据框友好结构（melt 帮助很大）这）。另一方面，如果我想为探索阶段快速生成散点图矩阵或格子直方图集（即统计矩、分离、类内/类间方差、直方图等），我必须停下来弄清楚如何我要把apply()这些巨大的多维数组变成那些库能理解的东西。
如果我一直在丛林中四处奔波，为此想出临时解决方案，我要么永远不会变得更好，要么我最终会以自己奇怪的巫术方式做这件事对任何人都有意义。
那么正确的方法是用在 R 中为探索性分析而重复试验积累的分类标签来构造多元数据的方法是什么？拜托，我不想溜回 MATLAB。
奖励：我倾向于对多个主题的相同数据结构重复这些分析。有没有比将代码块包装到for 循环中更好的通用方法？
【问题讨论】：
许多 R 函数期望您的数据是“长格式”，即当您melt 数据时您得到的数据。也有例外，但一般来说，您应该将数据存储在 data.frames 中（或者如果它在 data.tables 中很大），其中包含值列和因子列（分类器）。但您可能还想看看在 R 中提供比基本 ts 对象更复杂的时间序列对象的包（例如 xts 包）。
这是意识形态的。看一下包 data.table，这使得这非常有效。
我已将核心问题的副本放在最上面，尤其是 MATLAB 用户可能不会急于一直阅读到底部以发现他们无法贡献。
关于您的“奖励”：将所有内容放在一个大 data.table 中并使用包的 by 语法。
我还建议阅读 vita.had.co.nz/papers/tidy-data.html - 它阐述了我在 R 中处理数据的理念。
【参考方案1】：
如前所述，许多更强大的分析和可视化工具都依赖于长格式数据。当然，对于受益于矩阵代数的转换，您应该将内容保存在数组中，但是一旦您想要对数据子集进行并行分析，或者根据数据中的因素绘制内容，您真的想要melt。
下面是一个示例，可帮助您开始使用 data.table 和 ggplot。
数组 -> 数据表
首先，让我们按照您的格式制作一些数据：
series <- 3
samples <- 2
trials <- 4
trial.labs <- paste("tr", seq(len=trials))
trial.class <- sample(c("A", "B"), trials, rep=T)
arr <- array(
runif(series * samples * trials),
dim=c(series, samples, trials),
dimnames=list(
ser=paste("ser", seq(len=series)),
smp=paste("smp", seq(len=samples)),
tr=trial.labs
)
)
# , , tr = Trial 1
# smp
# ser smp 1 smp 2
# ser 1 0.9648542 0.4134501
# ser 2 0.7285704 0.1393077
# ser 3 0.3142587 0.1012979
#
# ... omitted 2 trials ...
#
# , , tr = Trial 4
# smp
# ser smp 1 smp 2
# ser 1 0.5867905 0.5160964
# ser 2 0.2432201 0.7702306
# ser 3 0.2671743 0.8568685
现在我们有了一个 3 维数组。让我们将melt 转换为data.table（注意melt 在data.frames 上运行，基本上是data.tables 没有花里胡哨，所以我们必须先融化，然后转换为data.table）：
library(reshape2)
library(data.table)
dt.raw <- data.table(melt(arr), key="tr") # we'll get to what the `key` arg is doing later
# ser smp tr value
# 1: ser 1 smp 1 tr 1 0.53178276
# 2: ser 2 smp 1 tr 1 0.28574271
# 3: ser 3 smp 1 tr 1 0.62991366
# 4: ser 1 smp 2 tr 1 0.31073376
# 5: ser 2 smp 2 tr 1 0.36098971
# ---
# 20: ser 2 smp 1 tr 4 0.38049334
# 21: ser 3 smp 1 tr 4 0.14170226
# 22: ser 1 smp 2 tr 4 0.63719962
# 23: ser 2 smp 2 tr 4 0.07100314
# 24: ser 3 smp 2 tr 4 0.11864134
请注意，这是多么容易，我们所有的维度标签都逐渐变成了长格式。 data.tables 的花里胡哨之一是能够在 data.tables 之间进行索引合并（很像 mysql 索引连接）。所以在这里，我们将这样做以将class 绑定到我们的数据：
dt <- dt.raw[J(trial.labs, class=trial.class)] # on the fly mapping of trials to class
# tr ser smp value class
# 1: Trial 1 ser 1 smp 1 0.9648542 A
# 2: Trial 1 ser 2 smp 1 0.7285704 A
# 3: Trial 1 ser 3 smp 1 0.3142587 A
# 4: Trial 1 ser 1 smp 2 0.4134501 A
# 5: Trial 1 ser 2 smp 2 0.1393077 A
# ---
# 20: Trial 4 ser 2 smp 1 0.2432201 A
# 21: Trial 4 ser 3 smp 1 0.2671743 A
# 22: Trial 4 ser 1 smp 2 0.5160964 A
# 23: Trial 4 ser 2 smp 2 0.7702306 A
# 24: Trial 4 ser 3 smp 2 0.8568685 A
需要了解的几点：
J 从向量创建 data.table
尝试将一个data.table 的行与另一个数据表进行子集化（即使用data.table 作为[.data.table 中大括号后的第一个参数）导致data.table 左连接（用MySQL 的说法）外部在这种情况下，将表（在这种情况下为dt）复制到内部表（由J 即时创建的表）。连接是在外部data.table 的key 列上完成的，您可能已经注意到我们在前面的melt/data.table 转换步骤中定义的。
您必须阅读文档才能完全理解发生了什么，但认为J(trial.labs, class=trial.class) 实际上等同于使用data.table(trial.labs, class=trial.class) 创建data.table，除了J 仅在[.data.table 内部使用时有效.
所以现在，通过一个简单的步骤，我们将类数据附加到值上。同样，如果您需要矩阵代数，请先对数组进行操作，然后用两三个简单的命令切换回长格式。正如 cmets 中所述，您可能不希望在长格式和数组格式之间来回切换，除非您有充分的理由这样做。
一旦内容在data.table 中，您就可以很容易地对数据进行分组/聚合（类似于拆分-应用-组合样式的概念）。假设我们要获取每个class-sample 组合的汇总统计信息：
dt[, as.list(summary(value)), by=list(class, smp)]
# class smp Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1: A smp 1 0.08324 0.2537 0.3143 0.4708 0.7286 0.9649
# 2: A smp 2 0.10130 0.1609 0.5161 0.4749 0.6894 0.8569
# 3: B smp 1 0.14050 0.3089 0.4773 0.5049 0.6872 0.8970
# 4: B smp 2 0.08294 0.1196 0.1562 0.3818 0.5313 0.9063
在这里，我们只给 data.table 一个表达式 (as.list(summary(value))) 来评估每个数据的 class、smp 子集（如 by 表达式中指定的那样）。我们需要as.list，以便data.table 将结果重新组合为列。
您可以轻松计算类/样本/试验/系列变量的任意组合的矩（例如list(mean(value), var(value), (value - mean(value))^3）。
如果您想对数据进行简单的转换，使用data.table 非常容易：
dt[, value:=value * 10] # modify in place with `:=`, very efficient
dt[1:2] # see, `value` now 10x
# tr ser smp value class
# 1: Trial 1 ser 1 smp 1 9.648542 A
# 2: Trial 1 ser 2 smp 1 7.285704 A
这是就地转换，因此没有内存副本，因此速度很快。通常data.table 会尝试尽可能高效地使用内存，因此这是进行此类分析的最快方法之一。
从长格式绘图
ggplot 非常适合以长格式绘制数据。我不会详细介绍正在发生的事情，但希望这些图片能让您了解您可以做什么
library(ggplot2)
ggplot(data=dt, aes(x=ser, y=smp, color=class, size=value)) +
geom_point() +
facet_wrap( ~ tr)
ggplot(data=dt, aes(x=tr, y=value, fill=class)) +
geom_bar(stat="identity") +
facet_grid(smp ~ ser)
ggplot(data=dt, aes(x=tr, y=paste(ser, smp))) +
geom_tile(aes(fill=value)) +
geom_point(aes(shape=class), size=5) +
scale_fill_gradient2(low="yellow", high="blue", midpoint=median(dt$value))
数据表 -> 数组 -> 数据表
首先我们需要将acast（来自包reshape2）我们的数据表返回到一个数组中：
arr.2 <- acast(dt, ser ~ smp ~ tr, value.var="value")
dimnames(arr.2) <- dimnames(arr) # unfortunately `acast` doesn't preserve dimnames properly
# , , tr = Trial 1
# smp
# ser smp 1 smp 2
# ser 1 9.648542 4.134501
# ser 2 7.285704 1.393077
# ser 3 3.142587 1.012979
# ... omitted 3 trials ...
此时，arr.2 看起来就像 arr 所做的一样，除了值乘以 10。请注意，我们必须删除 class 列。现在，让我们做一些简单的矩阵代数
shuff.mat <- matrix(c(0, 1, 1, 0), nrow=2) # re-order columns
for(i in 1:dim(arr.2)[3]) arr.2[, , i] <- arr.2[, , i] %*% shuff.mat
现在，让我们回到 melt 的长格式。注意key 参数：
dt.2 <- data.table(melt(arr.2, value.name="new.value"), key=c("tr", "ser", "smp"))
最后，让我们重新加入dt 和dt.2。在这里你需要小心。 data.table 的行为是，如果外表没有键，内表将根据内表的所有键连接到外表。如果内表有键，data.table 将键连接键。这是一个问题，因为我们预期的外部表 dt 已经在之前的 tr 上有一个键，所以我们的连接将只发生在该列上。因此，我们需要将键放在外表上，或者重置键（我们在这里选择了后者）：
setkey(dt, tr, ser, smp)
dt[dt.2]
# tr ser smp value class new.value
# 1: Trial 1 ser 1 smp 1 9.648542 A 4.134501
# 2: Trial 1 ser 1 smp 2 4.134501 A 9.648542
# 3: Trial 1 ser 2 smp 1 7.285704 A 1.393077
# 4: Trial 1 ser 2 smp 2 1.393077 A 7.285704
# 5: Trial 1 ser 3 smp 1 3.142587 A 1.012979
# ---
# 20: Trial 4 ser 1 smp 2 5.160964 A 5.867905
# 21: Trial 4 ser 2 smp 1 2.432201 A 7.702306
# 22: Trial 4 ser 2 smp 2 7.702306 A 2.432201
# 23: Trial 4 ser 3 smp 1 2.671743 A 8.568685
# 24: Trial 4 ser 3 smp 2 8.568685 A 2.671743
注意data.table通过匹配键列进行连接，即-通过将外部表的第一个键列匹配到内部表的第一列/键，第二个到第二个，依此类推，不考虑列名（有一个 FR here）。如果您的表/键的顺序不同（如果您注意到的话，这里就是这种情况），您需要重新排序您的列或确保两个表在您想要的列上具有相同顺序的键（我们在这里做了什么）。列顺序不正确的原因是因为我们添加类的第一个连接，它在tr 上连接并导致该列成为data.table 中的第一个。
【讨论】：
非常好。我只想指出，使用apply() 进行任何您想要的平均非常简单有效，可以适当地折叠尺寸，之前 melting ...
@TrevorAlexander 我明天会把这个添加到答案中，但简而言之：reshape2::acast 允许您从长格式重构数组。然后你再做你的代数，melt，如上所示，最后将你的结果加入表格。
人们应该仔细考虑（对于大数据），如果使用数组带来的效率提升值得从 data.table 重塑为数组所需的时间。留在 data.table 框架中可能更有效。
@Roland，同意，除非问题难以表述/长格式慢，而矩阵格式容易/快速，你应该坚持data.table。速度非常快。
所有这些重塑都会产生相当多的副本。而且这会降低性能，如果数据很大，则更是如此。因此，需要权衡代数的性能增益。您的数据越大，您就越需要关心副本。【参考方案2】：
也许是 dplyr::tbl_cube ？
根据@BrodieG 的出色回答，我认为您可能会发现查看dplyr::tbl_cube 提供的新功能很有用。这本质上是一个多维对象，您可以从数组列表（正如您当前正在使用的那样）轻松创建它，它具有一些非常好的子集、过滤和汇总功能（重要的是，我认为）在“数据的立方体”视图和“表格”视图。
require(dplyr)
几个注意事项：
这是一个早期版本：随之而来的所有问题建议此版本在加载 dplyr 时卸载 plyr
将数组加载到立方体中
这是一个使用 arr 的示例，如另一个答案中所定义：
# using arr from previous example
# we can convert it simply into a tbl_cube
arr.cube<-as.tbl_cube(arr)
arr.cube
#Source: local array [24 x 3]
#D: ser [chr, 3]
#D: smp [chr, 2]
#D: tr [chr, 4]
#M: arr [dbl[3,2,4]]
请注意，D 表示维度和 M 度量，您可以拥有任意多个。
从多维到平面的轻松转换
您可以通过将数据返回为 data.frame 轻松地制作数据表格（如果您以后需要功能和性能优势，您可以简单地将其转换为 data.table）
head(as.data.frame(arr.cube))
# ser smp tr arr
#1 ser 1 smp 1 tr 1 0.6656456
#2 ser 2 smp 1 tr 1 0.6181301
#3 ser 3 smp 1 tr 1 0.7335676
#4 ser 1 smp 2 tr 1 0.9444435
#5 ser 2 smp 2 tr 1 0.8977054
#6 ser 3 smp 2 tr 1 0.9361929
子集
您显然可以为每个操作展平所有数据，但这对性能和实用性有很多影响。我认为这个包的真正好处是你可以在将它转换成对 ggplot 友好的表格格式之前“预先挖掘”你需要的数据的立方体，例如简单过滤以仅返回系列 1：
arr.cube.filtered<-filter(arr.cube,ser=="ser 1")
as.data.frame(arr.cube.filtered)
# ser smp tr arr
#1 ser 1 smp 1 tr 1 0.6656456
#2 ser 1 smp 2 tr 1 0.9444435
#3 ser 1 smp 1 tr 2 0.4331116
#4 ser 1 smp 2 tr 2 0.3916376
#5 ser 1 smp 1 tr 3 0.4669228
#6 ser 1 smp 2 tr 3 0.8942300
#7 ser 1 smp 1 tr 4 0.2054326
#8 ser 1 smp 2 tr 4 0.1006973
tbl_cube 目前可与dplyr 函数summarise()、select()、group_by() 和filter() 一起使用。有用的是，您可以将它们与 %.% 运算符链接在一起。
对于其余的示例，我将使用内置的nasa tbl_cube 对象，它有一堆气象数据（并演示了多个维度和度量）：
分组和汇总措施
nasa
#Source: local array [41,472 x 4]
#D: lat [dbl, 24]
#D: long [dbl, 24]
#D: month [int, 12]
#D: year [int, 6]
#M: cloudhigh [dbl[24,24,12,6]]
#M: cloudlow [dbl[24,24,12,6]]
#M: cloudmid [dbl[24,24,12,6]]
#M: ozone [dbl[24,24,12,6]]
#M: pressure [dbl[24,24,12,6]]
#M: surftemp [dbl[24,24,12,6]]
#M: temperature [dbl[24,24,12,6]]
下面是一个示例，展示了从多维数据集中拉回 已修改 数据的子集是多么容易，然后然后 将其展平以使其适合绘图：
plot_data<-as.data.frame( # as.data.frame so we can see the data
filter(nasa,long<(-70)) %.% # filter long < (-70) (arbitrary!)
group_by(lat,long) %.% # group by lat/long combo
summarise(p.max=max(pressure), # create summary measures for each group
o.avg=mean(ozone),
c.all=(cloudhigh+cloudlow+cloudmid)/3)
)
head(plot_data)
# lat long p.max o.avg c.all
#1 36.20000 -113.8 975 310.7778 22.66667
#2 33.70435 -113.8 975 307.0833 21.33333
#3 31.20870 -113.8 990 300.3056 19.50000
#4 28.71304 -113.8 1000 290.3056 16.00000
#5 26.21739 -113.8 1000 282.4167 14.66667
#6 23.72174 -113.8 1000 275.6111 15.83333
n-d 和 2-d 数据结构的一致表示法
遗憾的是，mutate() 功能尚未为 tbl_cube 实现，但看起来这只是（不多）时间的问题。不过，您可以在表格结果中使用它（以及适用于多维数据集的所有其他函数） - 使用完全相同的符号。例如：
plot_data.mod<-filter(plot_data,lat>25) %.% # filter out lat <=25
mutate(arb.meas=o.avg/p.max) # make a new column
head(plot_data.mod)
# lat long p.max o.avg c.all arb.meas
#1 36.20000 -113.8000 975 310.7778 22.66667 0.3187464
#2 33.70435 -113.8000 975 307.0833 21.33333 0.3149573
#3 31.20870 -113.8000 990 300.3056 19.50000 0.3033389
#4 28.71304 -113.8000 1000 290.3056 16.00000 0.2903056
#5 26.21739 -113.8000 1000 282.4167 14.66667 0.2824167
#6 36.20000 -111.2957 930 313.9722 20.66667 0.3376045
绘图 - 作为“喜欢”平面数据的 R 功能示例
然后您可以利用扁平化数据的优势使用ggplot() 进行绘图：
# plot as you like:
ggplot(plot_data.mod) +
geom_point(aes(lat,long,size=c.all,color=c.all,shape=cut(p.max,6))) +
facet_grid( lat ~ long ) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
对生成的平面数据使用 data.table
我不打算在这里扩展data.table 的使用，因为它在上一个答案中做得很好。显然，使用data.table 有很多充分的理由 - 对于这里的任何情况，您都可以通过对 data.frame 的简单转换来返回：
data.table(as.data.frame(your_cube_name))
动态处理您的多维数据集
我认为很棒的另一件事是能够将度量（切片/场景/班次，无论您想如何称呼它们）添加到您的多维数据集。我认为这将非常适合问题中描述的分析方法。这是arr.cube 的一个简单示例 - 添加一个额外的度量，它本身就是前一个度量的（当然很简单）功能。您可以通过 yourcube$mets[$...]
head(as.data.frame(arr.cube))
# ser smp tr arr
#1 ser 1 smp 1 tr 1 0.6656456
#2 ser 2 smp 1 tr 1 0.6181301
#3 ser 3 smp 1 tr 1 0.7335676
#4 ser 1 smp 2 tr 1 0.9444435
#5 ser 2 smp 2 tr 1 0.8977054
#6 ser 3 smp 2 tr 1 0.9361929
arr.cube$mets$arr.bump<-arr.cube$mets$arr*1.1 #arb modification!
head(as.data.frame(arr.cube))
# ser smp tr arr arr.bump
#1 ser 1 smp 1 tr 1 0.6656456 0.7322102
#2 ser 2 smp 1 tr 1 0.6181301 0.6799431
#3 ser 3 smp 1 tr 1 0.7335676 0.8069244
#4 ser 1 smp 2 tr 1 0.9444435 1.0388878
#5 ser 2 smp 2 tr 1 0.8977054 0.9874759
#6 ser 3 smp 2 tr 1 0.9361929 1.0298122
尺寸 - 或不...
我尝试过动态添加全新维度（有效地扩展现有多维数据集并添加额外维度并使用 yourcube$dims[$...] 克隆或修改原始数据），但我发现行为有点不一致。无论如何，最好避免这种情况，并在操作它之前先构建你的立方体。如果我到达任何地方，会及时通知你。
持久性
显然，解释器访问多维数据库的主要问题之一是可能会因不合时宜的击键而意外地对其进行窃听。所以我想只要尽早并经常坚持：
tempfilename<-gsub("[ :-]","",paste0("DBX",(Sys.time()),".cub"))
# save:
save(arr.cube,file=tempfilename)
# load:
load(file=tempfilename)
希望有帮助！
【讨论】：
这是一个很好的答案，而且信息量很大！如果我在 R 的道路上再过六个月，我觉得我什至可以充分利用它:(你有关于这种结构的比较速度/内存使用情况的任何信息吗？data.table 真的引人入胜的调用。
@TrevorAlexander 看起来 data.table 和 dplyr 在性能方面大致相当：有些事情更快/更慢。 （这是来自文档和我自己完成的一些有限的基准测试。）我喜欢 dplyr 中的语法，但这是个人喜好！此外，如果您想使用它，您可以随时将 dplyr 的任何输出转换为 data.table。我很想将它存储在一个多维数据集中，并展开您想要动态分析和报告的位，但这有点取决于您的数据。应用程序是什么？
@Troy，您是否可以分享一些您发现data.table 速度较慢的基准？ （也许是要点并分享链接）？很高兴看看是否有任何改进。
我的主要应用是 epoch 多元时间序列。所以每个 epoch 就像一个方形切片，以系列作为列向量。这听起来很合适，但我必须尝试一下。以上是关于R 与 MATLAB 中的高维数据结构化方法的主要内容，如果未能解决你的问题，请参考以下文章
拓端tecdat|R语言高维数据惩罚回归方法：主成分回归PCR岭回归lasso弹性网络elastic net分析基因数据
R语言高维数据的主成分pca t-SNE算法降维与可视化分析案例报告