2

First off, I'm new to both R and Stack Overflow, so I apologize in advance for the inevitable faux pas or silly mistakes I make in this question and will do what I can to remedy them. In the meantime, please excuse my naïveté.

I'm using the R package Bchron to calibrate radiocarbon dates (i.e., turning 14C ages with error into probability distributions from a calibration curve, and especially, resulting in a point estimate (e.g., a mean) for that probability distribution). My goal is to get a vector of the weighted average ages resulting from calibrating multiple dates simultaneously. In this case, the weighted average is the sum of the products of the calibration age grid and the probability densities.

Let's say we are calibrating just three dates as follows:

library(Bchron)

ages1 <- BchronCalibrate(ages=c(3445,11553,7456), 
                    ageSds=c(50,230,110),
                    calCurves=c('intcal13','intcal13','shcal13'),
                    ids=c('sample1','sample2','sample3'))

Bchron outputs calibrations (i.e., from BchronCalibrate()) as a list of vectors. It's easy enough to get the mean one at a time:

print (sum(ages1$'sample1'$ageGrid * ages1$'sample1'$densities))

But when we calibrate multiple dates we get a list of dates (designated by "ids" above), each comprised of the calibration list for that date. Given the nested lists, I'm having a very hard time figuring out how to get multiple means in a vector from my list(s). Snooping around Stack Exchange, I've found lots of examples for iterating over a list, or even doing nested apply functions, but unfortunately, given my relative inexperience, I can't quite seem to adapt them to my problem. The following is what I've got so far.

Given that the ageGrid and densities vectors are always the 4th and 5th list items for each date in the parent list, I think it should be simple... something like

 for(sample.age in ages1) {
 sapply(ages1[[sample.age]],sum([4]*[5]))
 }

Given all the errors I've gotten in my various attempts, though, it's clear my syntax or (more probably) my entire line of thinking on this is mixed up.

And even if I can get the basic for loop to work, how do I ask it to return a vector?

I appreciate any help on this subject, and will especially appreciate answers that are simple enough for a novice like me to understand.

Cheers, Matt

m.s.bolton
  • 47
  • 5
  • From a statistical and physics point of view, I would have guessed that it would be necessary to know the scintillator counts that formed the basis of the various estimates. – IRTFM Mar 25 '18 at 22:23

1 Answers1

3

The R sapply and lapply functions are disguised loops. (And you only need one loop if I understand the question.) I'm eschewing your effort to use the numeric index and am instead suggesting you stick to using character indices as they are less prone to being misinterpreted.

sapply(ages1, function(x) sum(x[['ageGrid']]*x[['densities']]))

#  sample1   sample2   sample3 
# 3713.214 13372.815  8217.632 

The sapply function delivered those variously labeled sample-items one at a time to the function that extracted and multiplied them. Upon completion of hte multiplication it gave them back the names of the lists from which they were calculated. It's often useful at the beginning stages of R practice to look at the structure of the object you are dealing with:

str(ages1)
List of 3
 $ sample1:List of 5
  ..$ ages     : num 3445
  ..$ ageSds   : num 50
  ..$ calCurves: chr "intcal13"
  ..$ ageGrid  : num [1:473] 3470 3471 3472 3473 3474 ...
  ..$ densities: num [1:473] 1.01e-05 1.01e-05 1.01e-05 1.01e-05 1.01e-05 ...
 $ sample2:List of 5
  ..$ ages     : num 11553
  ..$ ageSds   : num 230
  ..$ calCurves: chr "intcal13"
  ..$ ageGrid  : num [1:1489] 12711 12712 12713 12714 12715 ...
  ..$ densities: num [1:1489] 1.05e-05 1.11e-05 1.17e-05 1.23e-05 1.30e-05 ...
 $ sample3:List of 5
  ..$ ages     : num 7456
  ..$ ageSds   : num 110
  ..$ calCurves: chr "shcal13"
  ..$ ageGrid  : num [1:714] 7860 7861 7862 7863 7864 ...
  ..$ densities: num [1:714] 1.00e-05 1.04e-05 1.08e-05 1.13e-05 1.17e-05 ...
 - attr(*, "class")= chr "BchronCalibratedDates"
thelatemail
  • 91,185
  • 12
  • 128
  • 188
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Thank you so much for the answer. This is great. I knew it would be "simple", but was really struggling with using a function in sapply. And let's just say I'm only barely grasping data indexes in lists anyhow (e.g., the need for double square brackets to actually get the data inside). – m.s.bolton Mar 26 '18 at 01:50
  • Yeah. That double bracket thisng is confusing. It's essential to know whether you're dealing with a list or a vector. In this case, you had the additional potential for confusion because it was a list of lists rather than just a list. – IRTFM Mar 26 '18 at 01:53