2

在下面编辑以显示一个非常简洁的解决方案——感谢哈德利·威克姆。

我有一个非常具体的查询，但它也与我想纠正的 R 知识中的一些普遍缺陷有关。我还想（如果可能的话）不仅解决我的问题，而且以一种优雅而有效的方式来解决问题——也许我将目光投向了高处。任何人都可以回答我的具体问题，但也可以推荐一个好的来源来了解更多信息吗？非常感谢任何帮助。似乎 Hadley Wickham 在这里遇到了类似的问题 - http://www.slideshare.net/hadley/plyr-one-data-analytic-strategy - 但这些是演示文稿中的幻灯片，我很难理解这些幻灯片他们自己。

我正在尝试操作存储在 R 列表中的 MCMC 输出。数据分为五年，每年我有四组。目标是绘制这些。为了使问题易于处理，这里是十次迭代的输出，就像这样。

iterations      [,1]       [,2]      [,3]       [,4]
      [1,] 49.184181  4.3515983 16.051958 -14.896019
      [2,] 45.910362  2.1738066 17.161775 -29.880989
      [3,] 14.575248  7.9476606  8.385455 -34.753004
      [4,] 55.029604  2.3422748 16.366960 -66.182627
      [5,] 25.338546  8.3039173 16.937638 -26.697235
      [6,] 48.633115  0.4698014 16.130142 -65.659992
      [7,]  1.356642  3.0249349  2.388576  -1.700559
      [8,] 49.831352 -2.0644832 15.403726 -23.378055
      [9,] 13.057886 -2.8856576 11.481152 -36.697754
     [10,] 50.889166  2.6846852 15.763382 -23.049868

, , 2


 iterations       [,1]      [,2]      [,3]       [,4]
      [1,] 51.6134663 15.659392 17.218244 -47.864892
      [2,] 46.0545981 17.067779 18.158151 -38.336587
      [3,] 16.5690775 10.386358 10.991029 -30.225820
      [4,] 55.5724832 14.840466 15.556193 -54.432882
      [5,] 26.1064404  5.656579 15.063810  -5.085942
      [6,] 57.3084200 12.551751 16.212203 -52.459935
      [7,]  0.9825892  6.651478  1.999976  -5.350995
      [8,] 56.1117252  3.204124 16.011812 -21.179722
      [9,] 15.4204854  5.761157 12.594028 -43.691113
     [10,] 50.1407397 16.404882 15.990908 -26.019990

, , 3


iterations      [,1]      [,2]      [,3]       [,4]
      [1,] 53.521436 24.340327 16.073063 -20.939950
      [2,] 46.040969 21.025351 16.535917 -47.611395
      [3,] 19.276578 16.575285 14.824175 -18.432136
      [4,] 58.050774 20.886686 15.944355 -37.646286
      [5,] 26.008007 11.449253 13.027001 -56.572886
      [6,] 61.474771 18.270354 15.879238 -31.316868
      [7,]  1.515227  1.434234  3.568761  -1.328706
      [8,] 61.725967 19.212081 16.717331 -18.993349
      [9,] 15.303739  6.939953 11.940742 -54.261739
     [10,] 47.968838 20.070758 17.168400 -48.598802

, , 4


 iterations      [,1]      [,2]      [,3]       [,4]
      [1,] 51.952695 24.267668 17.867717 -28.129743
      [2,] 49.680524 22.914727 16.001512 -44.434294
      [3,] 18.519755 17.961953 15.831455 -57.110802
      [4,] 59.652211 21.655724 16.876315 -24.965724
      [5,] 29.091609 20.831196 15.546565 -59.272164
      [6,] 62.190041 21.112490 15.759867 -19.910655
      [7,]  3.116584  1.187595  1.050807  -7.721749
      [8,] 61.384355 27.331487 16.646250 -17.793893
      [9,] 16.320224 14.321294 13.726538 -47.748184
     [10,] 47.676867 27.325987 17.056364 -31.032911

, , 5


iterations      [,1]      [,2]      [,3]      [,4]
      [1,] 55.326522 33.737691 19.698060 -46.34804
      [2,] 51.122038 31.055026 19.668949 -64.52942
      [3,] 22.036674 17.577561 13.546166 -85.24881
      [4,] 60.481009 34.300432 16.903054 -25.19277
      [5,] 29.168884 26.811356 16.066908 -37.56252
      [6,] 54.221450 28.760434 16.480317 -36.42441
      [7,]  3.672456  1.571084  2.397663 -10.91522
      [8,] 56.223306 30.730421 18.185858 -28.30282
      [9,] 16.955258 16.699139 18.101711 -36.85851
     [10,] 48.220404 29.749342 17.557532 -38.22831

一些进一步的信息：

> str(a.type)
List of 1
 $ a_type: num [1:10, 1:4, 1:5] 49.2 45.9 14.6 55 25.3 ...
  ..- attr(*, "dimnames")=List of 3
  .. ..$ iterations: NULL
  .. ..$           : NULL
  .. ..$           : NULL

我正在寻找（针对当前问题）是一种命名维度（即组和年份）的方法（使用dimnames()命令），其次，从每个列（组）中获取一些汇总值五年。将以下内容应用于五年中每一年的四列中的每一列：

 myfunc <- function(x)c(mean(x),
                   quantile(x,c(.025,.975)))

非常感谢任何帮助。另外，正如我所说，如果有人可以推荐此类问题的良好来源，那么我将来可能不必经常问这样的问题。

补充说明：根据下面的有用答案，我已经解决了部分问题。我可以将尺寸命名如下：

dimnames(a.type[[1]])=list(paste('iter',1:10,sep=''),                       ## 10 iterations
               paste(c("Delivery", "Other", "Regulatory", "Transfer")),     ## 4 groups
               paste('Year',1:5,sep=''))                                    ## 5 Years

这使得以下内容（仅显示第 1 年）：

> a.type
$a_type
, , Year1
        Delivery      Other Regulatory   Transfer 
iter1  49.184181  4.3515983  16.051958 -14.896019
iter2  45.910362  2.1738066  17.161775 -29.880989
iter3  14.575248  7.9476606   8.385455 -34.753004
iter4  55.029604  2.3422748  16.366960 -66.182627
iter5  25.338546  8.3039173  16.937638 -26.697235
iter6  48.633115  0.4698014  16.130142 -65.659992
iter7   1.356642  3.0249349   2.388576  -1.700559
iter8  49.831352 -2.0644832  15.403726 -23.378055
iter9  13.057886 -2.8856576  11.481152 -36.697754
iter10 50.889166  2.6846852  15.763382 -23.049868

所以这行得通。另一个问题：我如何才能命名组和年份——我对命名迭代没有太大兴趣，实际上我希望能够在不更改代码的情况下输出不同数量的迭代。换句话说，是否有一种逻辑方法可以跳过命名迭代。如果我做...

dimnames(a.type[[1]])=list(,                       ## 
               paste(c("Delivery", "Other", "Regulatory", "Transfer")), ## 4 groups
               paste('Year',1:5,sep=''))                                ## 5 Years

...然后我收到一条错误消息...

> dimnames(a.type[[1]][2:3])=list(#paste('iter',1:10,sep=''),                       ## 10 years
+                    paste(c("Delivery", "Other", "Regulatory", "Transfer")), ## 4 groups
+                    paste('Year',1:5,sep=''))                                ## 5 Years
Error in dimnames(a.type[[1]][2:3]) = list(paste(c("Delivery", "Other",  : 
 'dimnames' applied to non-array

另一方面，应用一个功能。我可以执行以下操作，但这让我认为所有年份的平均值和分位数：

> myfunc <- function(x)c(mean(x),
+                        quantile(x,c(.025,.975)))
>                      
> 
>                  
>                      
> a.type.bar <- apply(a.type[[1]], 2, myfunc)
> a.type.bar


   Delivery     Other Regulatory  Transfer
  38.351706 14.892788  14.450314 -34.61954
  2.5%   1.392323 -1.494269   2.087411 -66.06503
  97.5% 61.669447 33.134091  19.335254  -2.46227
 >

另一方面，我可以执行以下操作，并将我的功能一次仅应用于一年：

a.type.bar <- apply(a.type[[1]][,,1], 2, myfunc)

现在显然，这将解决我的问题——我只需要输入五行代码。但要弄清楚更深层次的问题，有没有办法一次获得平均值和分位数？

谢谢。

注意添加于 2013 年 3 月 17 日。感谢 Hadley Wickham 的奇妙 plyr 包，我似乎有了一个解决方案——感谢 Zach 让我接受它。

library(plyr)

myfunc <- function(x)c(mean(x),
                   quantile(x,c(.025,.975)))

summaries <- adply(a.type[[1]], 2:3, myfunc)

这给出了以下输出。

> summaries
       X1   X2           V1        2.5%       97.5%
1    Delivery 1996   78.6691388   39.912455   109.61078
2       Other 1996    4.3485461   -4.584758    16.61764
3  Regulatory 1996   19.6444938   14.135322    24.00373
4    Transfer 1996   -0.7922307 -195.263744   203.95175
5    Delivery 1997   79.6291215   29.853200   109.26860
6       Other 1997   14.3462871    5.607952    22.68043
7  Regulatory 1997   22.4131984   16.861994    30.09017
8    Transfer 1997 4392.7699174  991.168626  8426.64365
9    Delivery 1998   85.9237011   52.100181   115.78991
10      Other 1998   21.4735955    9.790307    37.40546
11 Regulatory 1998   25.5654754   19.558132    30.58021
12   Transfer 1998 6166.7374268 2456.330035 10249.00350
13   Delivery 1999   90.1843678   52.574874   128.28546
14      Other 1999   27.2028622   14.373959    38.54636
15 Regulatory 1999   28.8851480   20.913437    34.59272
16   Transfer 1999 8116.6049650 4186.782183 12030.65517
17   Delivery 2000   91.0299168   47.211931   125.35626
18      Other 2000   31.5885924   16.087480    46.28089
19 Regulatory 2000   31.7628775   21.082236    40.29969
20   Transfer 2000 9203.9975199 2349.851364 14382.00472

现在剩下的就是绘制这个（好吧，以及同一模型的其他几个版本）。我正在玩 ggplot。

4

2 回答 2

1

我不知道你的数组的尺寸，但这里有一个例子：

dat <- array(sample(1:5,10*4*5,rep=TRUE),c(10,4,5))

在这里使用dimnames是一个好主意，因为您有很多维度，这将帮助您了解聚合函数的输出。您只需要提供具有良好尺寸的名称列表。

dimnames(dat)=list(paste('year',1:10,sep=''),          ## 10 years
                   paste('group',letters[1:4],sep=''), ## 4 groups
                   paste('iter',1:5,sep=''))           ## 5 iterations

然后使用apply通过迭代获取手段

apply(dat,3,rowMeans)
       iter1 iter2 iter3 iter4 iter5
year1   2.25  3.00  3.75  3.00  3.00
year2   3.00  3.00  3.00  2.25  3.25
year3   3.75  3.50  3.50  3.50  3.50
year4   2.00  2.25  3.50  1.50  3.50
year5   2.50  2.50  3.50  2.75  3.50
year6   2.75  3.75  2.00  4.00  2.50
year7   3.50  2.50  3.50  2.50  2.75
year8   3.25  2.75  4.50  2.50  3.75
year9   4.50  3.25  3.25  3.00  2.25
year10  1.75  4.25  3.25  1.50  2.00

多年来按组获得手段

> apply(dat,3,colMeans)
       iter1 iter2 iter3 iter4 iter5
groupa   3.1   3.0   3.3   2.8   2.9
groupb   2.7   3.6   3.0   2.8   2.7
groupc   3.6   3.3   3.4   2.1   3.3
groupd   2.3   2.4   3.8   2.9   3.1
于 2013-03-14T18:53:04.027 回答
1

您想将数据放入数据框而不是矩阵，然后使用公式接口来aggregate.

理想情况下，您希望以可以直接读入数据框的形式获得 MCMC 输出，但如果您被矩阵卡住，则使用meltor reshape+as.data.frame或执行类似的操作（假设您有一个M使用三个上面讨论的尺寸）：

d<-data.frame(year=rep(1991:1995,each=40),
              agency=rep(c("D","O","T","R"),50),
              iteration=rep(0:9,5,each=4),
              spend=as.vector(M))

为了得到一个看起来像这样的数据框：

  year agency iteration      spend
1 1996      D         0  49.184181
2 1996      O         0   4.351598
3 1996      R         0  16.051958
4 1996      T         0 -14.896019   
5 1996      D         1  45.910362
6 1996      O         1   2.173807
7 1996      R         1  17.161775
...

现在您可以使用aggregate来应用您的功能，如下所示：

aggregate(spend~agency+year,d,myfunc)

要得到

   agency year   spend.V1 spend.2.5% spend.97.5%
1       D 1996  35.380610   3.989422   54.098005
2       O 1996   2.634854  -2.700893    8.223760
3       R 1996  13.607076   3.737874   17.111344
4       T 1996 -32.289610 -66.065034   -4.669537
5       D 1997  37.588003   4.231116   57.039164
6       O 1997  10.818397   3.755926   16.918627
...

现在你可以随心所欲地切片和切块

aggregate(spend~year,d,myfunc)
aggregate(spend~agency,d,myfunc)
etc...
于 2013-03-17T19:42:28.373 回答