pvermees / IsoplotR

An alternative to Isoplot written in R
64 stars 22 forks source link

IsoPlotR package gives different 2-sigma error to IsoplotR online #212

Closed MichelineCampbell closed 6 months ago

MichelineCampbell commented 6 months ago

Hi,

Many thanks for writing IsoPlotR!

In trying to figure out how to correct some U-Th dates, I have been attempting to replicate the corrections in Treble, P.C., Baker, A., Abram, N.J., Hellstrom, J.C., Crawford, J., Gagan, M.K., Borsato, A., Griffiths, A.D., Bajo, P., Markowska, M., Priestley, S.C., Hankin, S., Paterson, D., 2022. Ubiquitous karst hydrological control on speleothem oxygen isotope variability in a global study. Commun Earth Environ 3, 1–10. https://doi.org/10.1038/s43247-022-00347-3 (Supplementary Table 5).

The first two samples were analysed in 2017, giving corrected dates of 0.263 +- 0.07 and 0.406 +- 0.067 ka BP 1950.

IsoPlotR online gives comparable results to those presented in the paper (attempted first two dates only). However, in IsoPlotR (v 6.0) in R 4.3.1 the 2sigma error (oerr = 2) is double to that which both the online version and Treble et al. 2022 report.

The first two samples were analysed in 2017, giving corrected dates of 0.263 +- 0.07 and 0.406 +- 0.067 ka BP 1950.

The input data are:

t <- structure(list(x = structure(c(0.001205, 0.0010015, 1e-05, 7.7e-06, 
1.0604, 1.0564, 0.004, 0.0034, 0.00361, 0.00488, 0.00065, 0.00063, 
0, 0, 0, 0, 0, 0), dim = c(2L, 9L), dimnames = list(NULL, c("Th232U238", 
"errTh232U238", "U234U238", "errU234U238", "Th230U238", "errTh230U238", 
"rXY", "rXZ", "rYZ"))), format = 2, U8Th2 = 0, Th02i = c(0.33, 
0.25), Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0)), class = "ThU")

and the code I used to run the corrections:

age(
  t,
  isochron = FALSE,
  Th0i = 3, #Th0i = 3 says to use assumed initial 
  exterr = FALSE, # don't propagate the external error
  oerr = 2, #caculate output uncertainty as 2sigma absolute uncertainties
  sigdig = 2 # number of significant digits fir the uncertainty estimate
  )

results:

       t err[t]   48_0 err[48_0] cov[t,48_0]
1 0.3310 0.1341 1.0605    0.0080    -1.1e-09
2 0.4708 0.1307 1.0565    0.0068    -1.0e-09

331 - (2017-1950) = 264 = .264 ka BP 1950
471 - (2017-1950) = 404 = .404 ka BP 1950

Input to isoplotR online:

IsoPlotRhelp1

Reults from isoplotR online IsoPlotRhelp2

That the error from the offline version is double that of the online version (and of the Treble et al. results), suggests there is a bug in the calculation of the 2-sigma error in the offline version. I have had a look through the source files, but haven't been able to find it, sorry!

Please let me know if there is anything else you need. I will also keep looking for the bug.

Many thanks, Micha

pvermees commented 6 months ago

Dear Micha,

Thank you for the feedback! Fortunately, the discrepancy between the online version and the command line code is not caused by a bug. IsoplotR internally stores all data at 1-sigma. Your code bypasses the read.data() function and thereby omits the conversion from 2-sigma to 1-sigma. Here is an alternative solution that fixes the issue:

tab <- cbind(
 Th232U238=c(0.001205, 0.0010015),
 errTh232U238=c(1e-05, 7.7e-06),
 U234U238=c(1.0604, 1.0564),
 errU234U238=c(0.004, 0.0034),
 Th230U238=c(0.00361,0.00488),
 errTh230U238=c(0.00065, 0.00063),
 rXY=c(0,0), rXZ=c(0,0), rYZ=c(0,0)
)

t <- as.ThU(
 x=tab,
 format=2,
 ierr=2,
 U8Th2 = 0, 
 Th02i = c(0.33,0.25), 
 Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0)
)

age(
  t,
  isochron = FALSE,
  Th0i = 3, #Th0i = 3 says to use assumed initial 
  exterr = FALSE, # don't propagate the external error
  oerr = 2, #calculate output uncertainty as 2sigma absolute uncertainties
  sigdig = 2 # number of significant digits fir the uncertainty estimate
 )

I hope this helps,

Pieter

MichelineCampbell commented 6 months ago

Dear Pieter,

This is so helpful, thanks so much! I didn't read deeply enough into read.data to see that the data were stored at 1-sigma, I feel pretty silly now.

Thanks again for your help, and the quick reply :)

Cheers, Micha

On Mon, Feb 5, 2024 at 7:37 PM Pieter Vermeesch @.***> wrote:

Dear Micha,

Thank you for the feedback! Fortunately, the discrepancy between the online version and the command line code is not caused by a bug. IsoplotR internally stores all data at 1-sigma. Your code bypasses the read.data() function and thereby omits the conversion from 2-sigma to 1-sigma. Here is an alternative solution that fixes the issue:

tab <- cbind( Th232U238=c(0.001205, 0.0010015), errTh232U238=c(1e-05, 7.7e-06), U234U238=c(1.0604, 1.0564), errU234U238=c(0.004, 0.0034), Th230U238=c(0.00361,0.00488), errTh230U238=c(0.00065, 0.00063), rXY=c(0,0), rXZ=c(0,0), rYZ=c(0,0) )

t <- as.ThU( x=tab, format=2, ierr=2, U8Th2 = 0, Th02i = c(0.33,0.25), Th02U48 = c(0, 0, 1e+06, 0, 0, 0, 0, 0, 0) )

age( t, isochron = FALSE, Th0i = 3, #Th0i = 3 says to use assumed initial exterr = FALSE, # don't propagate the external error oerr = 2, #calculate output uncertainty as 2sigma absolute uncertainties sigdig = 2 # number of significant digits fir the uncertainty estimate )

I hope this helps,

Pieter

— Reply to this email directly, view it on GitHub https://github.com/pvermees/IsoplotR/issues/212#issuecomment-1926464564, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADNQYFEDX4GXEWRIHTCZTNDYSCK35AVCNFSM6AAAAABCZOCL2OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMRWGQ3DINJWGQ . You are receiving this because you authored the thread.Message ID: @.***>

pvermees commented 6 months ago

No worries. Thanks for using the command line :-)

pvermees commented 6 months ago

Issue closed.