CODEWITHSUNDEEP
rplyr
Lindsay
Lindsay

Reputation: 153

Manipulating list objects in R - processing MCMC output

EDITED BELOW TO SHOW A REALLY NEAT SOLUTION -- THANKS TO HADLEY WICKHAM.

I have a very specific query, but it also relates to some general shortcomings in my R knowledge which I would like to rectify. I'd like also (if possible) not just solve my problem but do so in an elegant and efficient way---maybe I am setting my sights to high. Can anyone both answer my specific queries, but also recommend a good source to find out more? Any help greatly appreciated. It seems Hadley Wickham has wrestled with a similar problem here - http://www.slideshare.net/hadley/plyr-one-data-analytic-strategy - but these are slides from a presentation, and I struggle to understand the slides by themselves.

I am trying to manipulate MCMC output stored in a list in R. The data are grouped into five years, and for each year I have four groups. The goal is to plot these. To make the problem tractable, here is the output for just ten iterations, like so. 

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

Some further information:

> 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

What I am looking for (for the immediate problem) is a way of naming the dimensions (i.e. the groups and the years) of this (with the dimnames() command), and second, taking some summary values from each column (group) in each of the five years. Something that will apply the following to each of the four columns for each of the five years:

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

Any help greatly appreciated. Also, as I said, if anyone can recommend a good source on problems like this, I might not have to ask questions like this so often in future.

Note added: Based on the helpful answer below, I have figured out part of my problem. I can name the dimensions as follows:

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

This makes the following (just showing year 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

So that works. A further question: how can I just name the groups and years---I have not much interest in naming the iterations, and indeed I want to be able to output different numbers of iterations without changing my code. In other words is there a logical way to skip over naming the iterations. If I do...

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

...then I get an error message...

> 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

On the other thing, applying a function. I can do the following, but that gives me I think the mean and quantiles across all years:

> 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
 >

On the other hand, I can do the following, and apply my function to just one year at a time:

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

Now obviously, that would solve my problem -- I would just have to type five lines of code. But to figure out the deeper problem, is there a way of getting means and quantiles a year at a time?

Thanks.

Note added 17 March 2013. Thanks to Hadley Wickham's marvellous plyr package, I seem to have a solution---and thanks Zach for turning me onto it. 

library(plyr)

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

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

This gives the following output. 

> 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

All that is left now is to plot this (well, and several other versions of the same model). I am having a play with ggplot.

Upvotes: 2

Views: 902

Answers (2)

Thruston
Thruston

Reputation: 1557

You want to get your data into a data frame instead of a matrix, and then use the formula interface to aggregate.

Ideally you want to get your MCMC output in a form that you can read directly into a data frame, but if you are stuck with the matrix, then use melt or reshape + as.data.frame or just do something like this (assuming you have a matrix called M with the three dimensions discussed above):

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))

in order to get a data frame that looks like this:

  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
...

Now you can use aggregate to apply your function, like this:

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

to get

   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
...

and now you can slice and dice as you wish

aggregate(spend~year,d,myfunc)
aggregate(spend~agency,d,myfunc)
etc...

Upvotes: 1

agstudy
agstudy

Reputation: 121578

I don't know the dimensions of your array , but here an example:

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

Using dimnames here is a good idea since you have many dimensions, this will help you to understand the output of your aggregation function. You need just to spply a list of names with the good dimensions.

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

Then using apply to get means by iteration

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

To get means by group over years

> 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

Upvotes: 1

Related Questions