2

Background

I have a dataframe of probability distributions that I would like to calculate statistical summaries for:

priors <- structure(list(name = c("theta1", "theta2", "theta3", "theta4", 
  "theta5"), distn = c("gamma", "beta", "lnorm", "weibull", "gamma"), 
   parama = c(2.68, 4, 1.35, 1.7, 2.3), paramb = c(0.084, 7.2, 0.69, 0.66, 3.9),
   another_col = structure(c(3L, 4L, 5L, 1L, 2L
   ), .Label = c("1", "2", "a", "b", "c"), class = "factor")), 
   .Names = c("name", "distn", "parama", "paramb", "another_col"), row.names = c("1",
   "2", "3", "4", "5"), class = "data.frame")

Approach

Step 1: I wrote a function to calculate the summaries and returning mean(lcl, ucl)

 summary.stats <- function(distn, A, B) {
  if (distn == 'gamma'  ) ans <- c(A*B,                       qgamma(c(0.05, 0.95), A[ ], B))
  if (distn == 'lnorm'  ) ans <- c(exp(A + 1/2 * B^2),        qlnorm(c(0.05, 0.95), A, B))
  if (distn == 'beta'   ) ans <- c(A/(A+B),                   qbeta( c(0.05, 0.95), A, B))
  if (distn == 'weibull') ans <- c(mean(rweibull(10000,A,B)), qweibull(c(0.05, 0.95), A, B))
  if (distn == 'norm'   ) ans <- c(A,                         qnorm( c(0.05, 0.95), A, B))
  ans <- (signif(ans, 2))
  return(paste(ans[1], ' (', ans[2], ', ', ans[3],')', sep = ''))
}

Step 2: I would like to add a new column to my dataframe called stats

priors$stats <- ddply(priors, 
                     .(name, distn, parama, paramb), 
                     function(x)  summary.stats(x$distn, x$parama, x$paramb))$V1

Question 1:

what is the proper way to do this? I get an error when I try

                ddply(priors, 
                     .(name, distn, parama, paramb),
                     transform, 
                     stats = function(x)  summary.stats(x$distn, x$parama, x$paramb))

Question 2: (extra credit)

Is there a more efficient way to code the summary.stats function, i.e., with less 'if's'?

update

Thanks to Shane and Joshua for clearing this up for me.

I also found a question that should be helpful for others trying to do a plyr operation on every row of a dataframe

Community
  • 1
  • 1
David LeBauer
  • 31,011
  • 31
  • 115
  • 189

2 Answers2

4

Here's a cleaned-up version of your summary.stats that uses switch instead. I also added the name "stats" to the output, since that seems to be the thing tripping you up.

summaryStats <- function(distn, A, B) {
  CI <- c(0.05, 0.95)
  FUN <- get(paste("q",distn,sep=""))
  ans <- switch(distn,
    gamma   = A*B,
    lnorm   = exp(A + 1/2 * B^2),
    beta    = A/(A+B),
    weibull = mean(rweibull(10000,A,B)),
    norm    = A)
  ans <- c(ans, FUN(CI, A, B))
  ans <- (signif(ans, 2))
  out <- c(stats=paste(ans[1], ' (', ans[2], ', ', ans[3],')', sep=''))
  return(out)
}

I'm not sure how to do this with plyr, but you can do it with boring ol' sapply like this:

priors$stats <- sapply(1:nrow(priors),
  function(i) with(priors[i,], summaryStats(distn, parama, paramb) ))
Joshua Ulrich
  • 173,410
  • 32
  • 338
  • 418
4

I could be missing something, but using Josh's function and your data, this works fine.

priors <- ddply(priors, 
  .(name, distn, parama, paramb), 
  function(x)  summaryStats(x$distn, x$parama, x$paramb))
colnames(priors)[5] <- "stats"

What do you want your output to look like?

> priors
    name   distn parama paramb            stats
1 theta1   gamma   2.68  0.084   0.23 (7.8, 69)
2 theta2    beta   4.00  7.200 0.36 (0.15, 0.6)
3 theta3   lnorm   1.35  0.690    4.9 (1.2, 12)
4 theta4 weibull   1.70  0.660 0.59 (0.12, 1.3)
5 theta5   gamma   2.30  3.900    9 (0.12, 1.3)

Edit

Sorry, didn't read your whole comment. Then this should work (in my example here, I leave out one column):

ddply(priors, .(distn, parama, paramb), function(x) 
   data.frame(x, stats=summaryStats(x$distn, x$parama, x$paramb)))
Shane
  • 98,550
  • 35
  • 224
  • 217
  • I think it's that he wants the new column name to be "stats". I've added that to my version of his function. – Joshua Ulrich Dec 09 '10 at 22:34
  • @Shane, thanks for your answer. What I am confused about is how to get `ddply` to output a single column that I can assign to a new column in priors, e.g. `priors$stats <- ddply(....)$V1`; I suspect that using `ddply()$V1` isn't correct use. The 'actual' `priors` dataframe has other columns that I would like to retain without having to specify them all in `ddply()`. – David LeBauer Dec 09 '10 at 22:36
  • 1
    You can avoid the name reassignment at the end by doing it within the function: `data.frame(stats=summaryStats(...))` – Shane Dec 09 '10 at 22:44
  • 1
    @Shane add some other columns to `priors` before running `ddply`. They won't be included in the output without specifying them in `.variables`. – Joshua Ulrich Dec 09 '10 at 22:45
  • I'll put an extra column in the `priors` example data frame – David LeBauer Dec 09 '10 at 22:47
  • @David: Sorry...should read your whole comment. yes, in that case, just `merge(priors, ddply(....))` should work (it will match on the same columns). You would only need to do this *if there are multiple rows per the grouping columns*. – Shane Dec 09 '10 at 22:49
  • @David Correction, look at my updated solution. Just add in the original into your return value. – Shane Dec 09 '10 at 22:59