/
calc_CommonDose.R
206 lines (188 loc) · 7.15 KB
/
calc_CommonDose.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
#' Apply the (un-)logged common age model after Galbraith et al. (1999) to a
#' given De distribution
#'
#' Function to calculate the common dose of a De distribution.
#'
#' **(Un-)logged model**
#'
#' When `log = TRUE` this function
#' calculates the weighted mean of logarithmic De values. Each of the estimates
#' is weighted by the inverse square of its relative standard error. The
#' weighted mean is then transformed back to the dose scale (Galbraith &
#' Roberts 2012, p. 14).
#'
#' The log transformation is not applicable if the
#' De estimates are close to zero or negative. In this case the un-logged model
#' can be applied instead (`log = FALSE`). The weighted mean is then
#' calculated using the un-logged estimates of De and their absolute standard
#' error (Galbraith & Roberts 2012, p. 14).
#'
#' @param data [RLum.Results-class] or [data.frame] (**required**):
#' for [data.frame]: two columns with De `(data[,1])` and De error `(data[,2])`
#'
#' @param sigmab [numeric] (*with default*):
#' additional spread in De values.
#' This value represents the expected overdispersion in the data should the sample be
#' well-bleached (Cunningham & Walling 2012, p. 100).
#' **NOTE**: For the logged model (`log = TRUE`) this value must be
#' a fraction, e.g. 0.2 (= 20 \%). If the un-logged model is used (`log = FALSE`),
#' sigmab must be provided in the same absolute units of the De values (seconds or Gray).
#'
#' @param log [logical] (*with default*):
#' fit the (un-)logged central age model to De data
#'
#' @param ... currently not used.
#'
#' @return
#' Returns a terminal output. In addition an
#' [RLum.Results-class] object is returned containing the
#' following element:
#'
#' \item{.$summary}{[data.frame] summary of all relevant model results.}
#' \item{.$data}{[data.frame] original input data}
#' \item{.$args}{[list] used arguments}
#' \item{.$call}{[call] the function call}
#'
#' The output should be accessed using the function [get_RLum]
#'
#' @section Function version: 0.1.1
#'
#' @author
#' Christoph Burow, University of Cologne (Germany)
#'
#' @seealso [calc_CentralDose], [calc_FiniteMixture],
#' [calc_FuchsLang2001], [calc_MinDose]
#'
#' @references
#' Galbraith, R.F. & Laslett, G.M., 1993. Statistical models for
#' mixed fission track ages. Nuclear Tracks Radiation Measurements 4, 459-470.
#'
#' Galbraith, R.F., Roberts, R.G., Laslett, G.M., Yoshida, H. & Olley,
#' J.M., 1999. Optical dating of single grains of quartz from Jinmium rock
#' shelter, northern Australia. Part I: experimental design and statistical
#' models. Archaeometry 41, 339-364.
#'
#' Galbraith, R.F. & Roberts, R.G., 2012. Statistical aspects of equivalent dose and error calculation and
#' display in OSL dating: An overview and some recommendations. Quaternary
#' Geochronology 11, 1-27.
#'
#' **Further reading**
#'
#' Arnold, L.J. & Roberts, R.G., 2009. Stochastic modelling of multi-grain equivalent dose
#' (De) distributions: Implications for OSL dating of sediment mixtures.
#' Quaternary Geochronology 4, 204-230.
#'
#' Bailey, R.M. & Arnold, L.J., 2006. Statistical modelling of single grain quartz De distributions and an
#' assessment of procedures for estimating burial dose. Quaternary Science
#' Reviews 25, 2475-2502.
#'
#' Cunningham, A.C. & Wallinga, J., 2012. Realizing the potential of fluvial archives using robust OSL chronologies.
#' Quaternary Geochronology 12, 98-106.
#'
#' Rodnight, H., Duller, G.A.T., Wintle, A.G. & Tooth, S., 2006. Assessing the reproducibility and accuracy
#' of optical dating of fluvial deposits. Quaternary Geochronology, 1 109-120.
#'
#' Rodnight, H., 2008. How many equivalent dose values are needed to
#' obtain a reproducible distribution?. Ancient TL 26, 3-10.
#'
#' @examples
#'
#' ## load example data
#' data(ExampleData.DeValues, envir = environment())
#'
#' ## apply the common dose model
#' calc_CommonDose(ExampleData.DeValues$CA1)
#'
#' @md
#' @export
calc_CommonDose <- function(
data,
sigmab,
log=TRUE,
...
) {
##============================================================================##
## CONSISTENCY CHECK OF INPUT DATA
##============================================================================##
if(missing(data)==FALSE){
if(is(data, "data.frame") == FALSE & is(data,"RLum.Results") == FALSE){
stop("[calc_CentralDose] Error: 'data' object has to be of type
'data.frame' or 'RLum.Results'!")
}else{
if(is(data, "RLum.Results") == TRUE){
data <- get_RLum(data, "data")
}
}
}
try(colnames(data)<- c("ED","ED_Error"), silent = TRUE)
if(colnames(data[1])!="ED"||colnames(data[2])!="ED_Error") {
cat(paste("Columns must be named 'ED' and 'ED_Error'"), fill = FALSE)
stop(domain=NA)
}
if(!missing(sigmab)) {
if(sigmab <0 | sigmab >1) {
cat(paste("sigmab needs to be given as a fraction between",
"0 and 1 (e.g. 0.2)"), fill = FALSE)
stop(domain=NA)
}
}
##============================================================================##
## ADDITIONAL ARGUMENTS
##============================================================================##
settings <- list(verbose = TRUE)
settings <- modifyList(settings, list(...))
##============================================================================##
## CALCULATIONS
##============================================================================##
# set default value of sigmab
if (missing(sigmab)) sigmab<- 0
# calculate yu = log(ED) and su = se(logED)
if (log) {
yu<- log(data$ED)
su<- sqrt( (data$ED_Error/data$ED)^2 + sigmab^2 )
}
else {
yu<- data$ED
su<- sqrt((data$ED_Error)^2 + sigmab^2)
}
# calculate weights
wu<- 1/su^2
delta<- sum(wu*yu)/sum(wu)
n<- length(yu)
#standard error
sedelta<- 1/sqrt(sum(wu))
if (!log) {
sedelta<- sedelta/delta
}
if (log){
delta<- exp(delta)
}
##============================================================================##
## TERMINAL OUTPUT
##============================================================================##
if (settings$verbose) {
cat("\n [calc_CommonDose]")
cat(paste("\n\n----------- meta data --------------"))
cat(paste("\n n: ",n))
cat(paste("\n log: ",if(log==TRUE){"TRUE"}else{"FALSE"}))
cat(paste("\n----------- dose estimate ----------"))
cat(paste("\n common dose: ", round(delta,2)))
cat(paste("\n SE: ", round(delta*sedelta, 2)))
cat(paste("\n rel. SE [%]: ", round(sedelta*100,2)))
cat(paste("\n------------------------------------\n\n"))
}
##============================================================================##
## RETURN VALUES
##============================================================================##
summary<- data.frame(de=delta,
de_err=delta*sedelta)
call<- sys.call()
args<- list(log=log, sigmab=sigmab)
newRLumResults.calc_CommonDose<- set_RLum(
class = "RLum.Results",
data = list(summary = summary,
data = data,
args = args,
call = call))
invisible(newRLumResults.calc_CommonDose)
}