Function to estimate temporal changes in isotopic values
estimateIntervals.Rd
Given renewal rates for different isotopic probes over time, a model
estimates a normal distribution of the isotopic values for each time interval. Out of
the given rates of renewal, first the influence of each time interval on the final
isotopic value is calculated. This proportion is used as a regressor in a fully Bayesian
linear model. The variance is estimated as rbf-Kernel-Matrix, where the hyperparameters are
chosen such that the correlation between time intervals close to each other is about 0.5.
The estimation is implemented using link[rstan]{rstan}
.
Usage
estimateIntervals(
renewalRates,
timeVars,
boneVars = NULL,
indVar = "",
isoMean,
isoSigma,
renewalRatesSD = NULL,
mu_df = 3,
mu_mean = 0,
mu_sd = 20,
rho_mean = 1,
rho_sd = 0.25,
alpha_mean = 2,
alpha_sd = 0.5,
mc = TRUE,
adapt_delta = 0.98,
max_treedepth = 15,
chains = 4,
iter = 2000,
burnin = 500,
...
)
Arguments
- renewalRates
A dataframe specifying the renewal rates of different probes for each time interval. The renewalRates should be between 0 and 100 (percentages). The dataframe should include a column specifying a time-index (e.g. 1, 2, 3, ...) as well as columns for the different bones. The renewal rates should contain the times of origin, containing 100.
- timeVars
A character string specifying the name of the column indicating the time.
- boneVars
A vector of character strings indicating the relevant variables containing the renewal rates of bones and teeth. If not given, the variables with bones are all variables of the dataframe apart from the time variables.
- indVar
A character string specifying the name of the column indicating the individual.
- isoMean
A numeric number indicating the mean of the isotopic values measured.
- isoSigma
A numeric, positive number indicating the standard deviation of the isotopic values measured.
- renewalRatesSD
An optional dataframe specifying the renewal rates uncertainties in standard deviations of different probes for each time interval.
- mu_df
Hyperparameter for the mean of the interval estimates: degrees of freedom of a non-standardized t-Student distribution. Defaults to 1.
- mu_mean
Hyperparameter for the mean of the interval estimates: mean of a non-standardized t-Student distribution. Defaults to 0.
- mu_sd
Hyperparameter for the mean of the interval estimates: standard deviation of a non-standardized t-Student distribution. Defaults to 10.
- rho_mean
Hyperparameter for the length scale of the rbf-kernel: mean of a normal. Defaults to 1.
- rho_sd
Hyperparameter for the length scale of the rbf-kernel: standard deviation of a normal. Defaults to 0.25.
- alpha_mean
Hyperparameter for the signal variance of the rbf-kernel: mean of a normal. Defaults to 2.
- alpha_sd
Hyperparameter for the signal variance of the rbf-kernel: standard deviation of a normal. Defaults to 0.5.
- mc
A boolean indicating if multiple cores should be used. If
TRUE
, which is the default, 4 cores are used.- adapt_delta
A numeric value between 0 and 1 controlling the samplers behavior. Defaults to 0.9999. See
stan
for more details.- max_treedepth
A numeric, positive value controling the NUTS sampler. Defaults to 25. See
stan
for more details.- chains
Number of chains for mcmc. Defaults to 4
- iter
Number of iterations per chain for mcmc. Defaults to 2000
- burnin
Number of burnin iterations per chain for mcmc. Defaults to 500
- ...
Arguments passed to rstand sampling
Value
A fitted object of class TemporalIso.
Examples
if (FALSE) {
#dataset with renewal rates
data <- data.frame(
intStart = 0:5,
intEnd = 1:6,
bone1 = c(100, 50, 20, 10, 5, 2),
bone2 = c(100, 10, 5, 1, 1, 1),
tooth1 = c(0, 100, 0, 0, 0, 0),
tooth2 = c(0, 0, 100, 0, 0, 0)
)
# isotope-values and sd
y_mean <- c(-10, -7, -12, -9)
y_sigma <- c(2, 1.5, 2.5, 2.5)
###############################
#estimate values
###############################
fit <- estimateIntervals(renewalRates = data,
timeVars = "intStart",
boneVars = c("bone1", "bone2", "tooth1", "tooth2"),
isoMean = y_mean,
isoSigma = y_sigma)
print(fit)
# implemented Methods
plotTime(fit)
#get estimates for specific time points
estimateTimePoint(fit, time = seq(0,5, by = 0.5))
###############################
#shift point estimation
###############################
plotTime(fit, plotShifts = TRUE, threshold = 0.5)
plotTime(fit, plotShifts = TRUE, threshold = 0.5,
type = TRUE, probability = 0.5)
getShiftTime(fit, threshold = 0.5)
###############################
#Staying time estimation
###############################
estimatedStayTimes <- getSiteStayTimes(object = fit,
siteMeans = c(-8, -10),
siteSigma = c(1, 1.5))
###############################
#compute resulting isotopic values out of historic ones
###############################
testDat <- data.frame(
t = 1:3,
bone1 = c(100, 50, 10, 10, 10, 10),
bone2 = c(0, 20, 20, 20, 20, 20),
teeth1 = c(0, 0, 0, 10, 50, 60),
mean = c(1, 3, 3,0, -1, - 5),
sd = c(1, 3, 1, 1.5, 2, 3)
)
computeResult(
data = testDat,
timeVars = "t",
boneVars = c("bone1","bone2","teeth1"),
meanVar = "mean",
sdVar = "sd"
)
shiny::runApp(paste0(system.file(package = "OsteoBioR"),"/app"))
}