R与MATLAB中高维数据结构的方法
bri*_*tar 20 arrays matlab methodology r multidimensional-array
题
在R中进行探索性分析的重复试验中积累的分类标签构建多变量数据的正确方法是什么？我不想回到MATLAB.
说明
我喜欢R的分析函数和语法(以及令人惊叹的绘图)比MATLAB更好,并且一直在努力重构我的东西.但是,我一直对工作中数据的组织方式感到困惑.
MATLAB
我通常使用多个时间序列来重复多次试验,这些试验存储在SERIESxSAMPLESxTRIALS 的大矩阵 秩-3张量多维数组中.这偶尔适用于一些不错的线性代数东西,但是当涉及另一个变量,即CLASS时,它是笨拙的.通常,类标签存储在尺寸为1x的另一个矢量中TRIALS.
在分析方面,我基本上尽可能少地绘制,因为需要做很多工作来组合一个非常好的情节,教你很多关于MATLAB中的数据.(我不是唯一一个有这种感觉的人).
[R
在R中我一直坚持尽可能接近MATLAB结构,但是当试图保持类标签分离时,事情变得非常复杂; 即使我只使用它们的属性,我也必须继续将标签传递给函数.所以我所做的就是通过CLASS将数组分成一个数组列表.这增加了我所有apply()功能的复杂性,但在保持一致性(和错误)方面似乎是值得的.
另一方面,R对于张量/多维数组似乎并不友好.只是为了与他们合作,你需要抓住abind图书馆.关于多变量分析的文档,比如这个例子似乎是假设你有一个巨大的2-D数据点表,比如一些长期的中世纪滚动数据框,并且没有提到如何从我所在的位置获得'那里' .
一旦我对绘制的数据进行绘图和分类,它就不是一个大问题了,因为到那时我一直在使用像TRIALSxFEATURES这样的形状的数据框架友好结构(melt对此有很大帮助).另一方面,如果我想快速生成探索阶段的散点图矩阵或latticist直方图集(即统计矩,分离,类/方差,直方图等),我必须停下来弄清楚如何我将apply()这些巨大的多维数组转换为这些库所理解的东西.
如果我继续在丛林中捣乱,为此提出临时解决方案,我要么永远不会变得更好,要么我最终会以自己奇怪的巫术方式结束,这对任何人都没有意义.
那么,对于R中的探索性分析,在重复试验中积累的分类标签构建多变量数据的正确方法是什么？拜托,我不想回到MATLAB.
额外奖励:我倾向于在多个主题的相同数据结构上重复这些分析.有没有比将代码块包装到for循环更好的通用方法？
Bro*_*ieG 19
正如已经指出的,许多更强大的分析和可视化工具依赖于长格式的数据.当然,对于受益于矩阵代数的转换,您应该将数据保存在数组中,但只要您想要对数据的子集进行并行分析,或者根据数据中的因子绘制内容,您真的想要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
现在我们有一个三维数组.让我们melt把它变成一个data.table(注意melt操作data.frames,这基本上data.table是s sans bells&whistles,所以我们必须首先融化,然后转换为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
有几点需要了解:
Jdata.table从矢量创建一个- 尝试将一个行的子集
data.table与另一个数据表进行子集(即使用adata.table作为大括号后的第一个参数[.data.table)导致data.table将外连接(dt在本例中)连接到内部表(在此情况下)创建的内容表J在这种情况下,飞行).连接key在外部的列上完成data.table,您可能已经注意到我们在之前的melt/data.table转换步骤中定义了.
你必须阅读文档才能完全理解正在发生的事情,但想到J(trial.labs, class=trial.class)有效地创建一个data.tablewith data.table(trial.labs, class=trial.class),除了J只在内部使用时才有效[.data.table.
所以现在,通过一个简单的步骤,我们将类数据附加到值.同样,如果您需要矩阵代数,请先对您的阵列进行操作,然后在两到三个简单命令中切换回长格式.正如评论中所指出的,除非你有充分的理由这样做,否则你可能不希望从长格式到格式格式来回转换.
一旦进入data.table,您可以非常轻松地对数据进行分组/聚合(类似于split-apply-combine样式的概念).假设我们要获得每个汇总统计数据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 ).如果您的表/键的顺序不同(如果您注意到的话,就像这里的情况一样),您需要重新排序列或确保两个表在所需的列上都有相同顺序的列(我们在这里做了什么).列不正确顺序的原因是因为我们在第一次连接时添加了类,该类加入tr并导致该列成为第一个data.table.
- 所有这些重塑产生了相当多的副本.如果数据量很大,那么性能会更高.因此,与代数的性能增益存在交易.数据越大,您需要关注的副本就越多. (4认同)
Tro*_*roy 11
也许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().有用的,您可以与%.%操作员一起链接这些.
对于其他示例,我将使用内置的nasatbl_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
nd和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
尺寸 - 或不......
我尝试动态添加全新维度(有效地扩展现有多维数据集以及其他维度并使用您的立方体克隆或修改原始数据$dims[$...]),但发现行为有点不一致.无论如何最好避免这种情况,并在操作它之前先构建你的立方体.如果我到达任何地方,会让你发布.
持久性
显然,让解释器访问多维数据库的主要问题之一是有可能在不合时宜的击键时误操作它.所以我想早点经常坚持:
tempfilename<-gsub("[ :-]","",paste0("DBX",(Sys.time()),".cub"))
# save:
save(arr.cube,file=tempfilename)
# load:
load(file=tempfilename)
希望有所帮助!
- @Troy,你有可能分享一些你发现`data.table`较慢的基准吗？(也许是一个要点并分享链接)？很高兴看到是否可以有任何改进. (4认同)
- 这是一个很好的答案和超级信息!如果我和R在另外六个月后,我觉得我甚至可以充分利用它:(你对这种结构的比较速度/内存使用情况有什么信息吗？`data.table`它的参考电话非常诱人. (3认同)