# R/Tagloss_cumul.R In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

#### Documented in Tagloss_cumul

```#' Tagloss_cumul returns the cumulative rate of tag loss.
#' @title Return the cumulative rate of tag loss.
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Return the cumulative rate of tag loss if hessian is null or a data.frame with distribution of cumulative rate of tag loss if hessian is not null.
#' @param t Time for which values of model must be estimated
#' @param x A Tagloss fitted model
#' @param par Parameters
#' @param Hessian Hessian matrix of parameters
#' @param mcmc A mcmc result
#' @param method Can be NULL, "delta", "SE", "Hessian", "MCMC", or "PseudoHessianFromMCMC"
#' @param model_before Function to be used before estimation of daily tagloss rate
#' @param model_after Function to be used after estimation of daily tagloss rate
#' @param model The model of parameter to use, can be N2, N1 or N0; or NLR, NL0, N0R, or N00 or NULL if hessian is NULL.
#' @param replicates Number of replicates to estimate se of output for resampling method
#' @description This function compute a model of cumulative tag loss rate for days t
#' based on a set of parameters, par.\cr
#' If hessian is not null, it will estimate standard error of the output using numerical delta method is replicates
#' is null or using resampling if replicates is not null.\cr
#' Parameters are described in \code{\link{Tagloss_fit}}.
#' @family Model of Tag-loss
#' @examples
#' \dontrun{
#' # Example
#' library(phenology)
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 1 - Fig 1D
#' par <- c(a0_2=-5.43E-2, a1_2=-103.52, a2_2=0, a3_2=0, a4_2=5.62E-4)
#' (y <- Tagloss_cumul(t=(1:6)*365, par=par))
#' plot(y[, "time"], y[, "N2"], type="l", bty="n",
#'      xlab="Days after tagging", ylab="N2 proportion")
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 2 - Fig 1E
#' par <- c(a0_2=-6.80E-2, a1_2=-81.15, a2_2=-2.20E-4, a3_2=6348.01, a4_2=1.65E-3)
#' (y <- Tagloss_cumul(t=(1:6)*365, par=par))
#' plot(y[, "time"], y[, "N2"], type="l",
#'      xlab="Days after tagging", ylab="N2 proportion")
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 3 - Fig 1F
#' par <- c(a0_2=-6.93E-2, a1_2=-78.92, a2_2=8.45E-4, a3_2=-16272.76, a4_2=2.87E-4)
#' (y <- Tagloss_cumul(t=(1:6)*365, par=par))
#' plot(y[, "time"], y[, "N2"], type="l",
#'      xlab="Days after tagging", ylab="N2 proportion")
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 4 - Fig 1C
#' par <- c(a0_2=-1.68E-3, a1_2=-4141.68, a2_2=0, a3_2=0, a4_2=0)
#' (y <- Tagloss_cumul(t=(1:6)*365, par=par))
#' plot(y[, "time"], y[, "N2"], type="l",
#'      xlab="Days after tagging", ylab="N2 proportion")
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 5 - Fig 1B
#' par <- c(a0_2=-3.77E-4, a1_2=-2000, a2_2=-0.001, a3_2=0, a4_2=5.00E-8)
#' (y <- Tagloss_cumul(t=(1:6)*365, par=par))
#' plot(y[, "time"], y[, "N2"], type="l",
#'      xlab="Days after tagging", ylab="N2 proportion")
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 6 - Fig 1A
#' par <- c(a0_2=-1E5, a1_2=-2000, a2_2=0, a3_2=4000, a4_2=8.34E-4)
#' (y <- Tagloss_cumul(t=(1:6)*365, par=par))
#' plot(y[, "time"], y[, "N2"], type="l",
#'      xlab="Days after tagging", ylab="N2 proportion")
#'
#' # Data from Rivalan et al. 2005 - Table 2, line 1 - Fig 1D
#' # With tagloss rate dependency on tage number
#' par <- c(a0_2=-5.43E-2, a1_2=-103.52, a2_2=0, a3_2=0, a4_2=5.62E-4,
#'          a0_1=-5.43E-2, a1_1=-103.52, a2_1=0, a3_1=0, a4_1=5.62E-4, delta_1=3.2E-4)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#'
#' p2 <- Tagloss_model(t=1:(6*365), par=par, model="2")
#' p1 <- Tagloss_model(t=1:(6*365), par=par, model="1")
#' par(mar=c(4, 5, 2, 1))
#' plot(x=1:(6*365), y=p2, bty="n", type="l", las=1, ylim=c(0,0.003), ylab="")
#' mtext("Daily tag loss", side=2, line=4)
#' lines(x=1:(6*365), y=p1, col="red")
#' legend("topright", legend=c("2>1", "1>0"), lty=1, col=c("black", "red"))
#'
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' # Without tagloss rate dependency on tag number
#' par <- c(a0_2=-5.43E-2, a1_2=-103.52, a2_2=0, a3_2=0, a4_2=5.62E-4,
#'          a0_1=-5.43E-2, a1_1=-103.52, a2_1=0, a3_1=0, a4_1=5.62E-4)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' #### Data from Casale et al. 2017
#' # Table 1 - Model II
#' par <- c(CasaleModelIIa0_2=-0.0511, CasaleModelIIa1_2=-100, CasaleModelIIa4_2=0.00014)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' # Table 1 - Model IV
#' par <- c(CasaleModelIVa0_2=-0.0132, CasaleModelIVa1_2=-100,
#'          CasaleModelIVa2_2=0.0327, CasaleModelIVa3_2=109.98,
#'          CasaleModelIVa4_2=0.00011)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' # Table 1 - Model I
#' par <- c(CasaleModelIc_2=0.00027)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' # Table 1 - Model III
#' par <- c(CasaleModelIIIa0_2=1.14E-10, CasaleModelIIIa1_2=-110.04,
#'          CasaleModelIIIa4_2=0.00055)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' # Table 1 - Model V
#' par <- c(CasaleModelVa0_2=4.04E-10, CasaleModelVa1_2=-90,
#'          CasaleModelVa2_2=-0.0326, CasaleModelVa3_2=100.31,
#'          CasaleModelVa4_2=0.00006)
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par,
#'                          model="Cumul",
#'                          decoration = TRUE)
#' Tagloss_cumul(t=(1:6)*365, par=par)
#'
#' }
#' @export

Tagloss_cumul <- function(t                                                     ,
par=NULL                                              ,
Hessian=NULL                                          ,
mcmc=NULL                                             ,
method=NULL                                           ,
model_before = NULL                                   ,
model_after=NULL                                      ,
model=NULL                                            ,
replicates=NULL                                       ,
x=NULL                                                ) {

# par=NULL; Hessian=NULL; model_before = NULL; model_after=NULL; model="cumul"; replicates=NULL; x=NULL

if (!is.null(model)) model <- toupper(model)
# D'abord je fais sans hessian

if (is.null(x) & is.null(par)) {
stop("Both par and x cannot be null at the same time")
}

if (is.null(replicates)) replicates <- 0
if (is.null(method) | (replicates == 0)) method <- "null"
method <- tolower(method)
method <- match.arg(method, choices = c("null",
"delta",
"se",
"hessian",
"mcmc",
"pseudohessianfrommcmc"))

if (!is.null(x)) {
if (is.null(par)) par <- c(x\$par, x\$fixed.par)
if (is.null(Hessian)) Hessian <- x\$hessian
if (is.null(model_before)) model_before <- x\$model_before
if (is.null(model_after)) model_after <- x\$model_after
}

days.maximum <- max(t)

if (method == "null") {

if (!is.null(model_before)) eval(parse(text=model_before), envir= environment())

if (any(grepl("_R2", names(par)))) pR2 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model="R2") else pR2 <- NA
if (any(grepl("_R1", names(par)))) pR1 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model="R1") else pR1 <- NA
if (any(grepl("_L2", names(par)))) pL2 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model="L2") else pL2 <- NA
if (any(grepl("_L1", names(par)))) pL1 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model="L1") else pL1 <- NA

if (any(grepl("_2", names(par)))) p2 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model="2") else p2 <- NA
if (any(grepl("_1", names(par)))) p1 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model="1") else p1 <- NA

if (!any(grepl("_", names(par)))) p1 <- p2 <- Tagloss_model(1:days.maximum, par, model_before = model_before, model_after = model_after, model=NULL)

if (!is.null(model_after)) eval(parse(text=model_after), envir= environment())

if (is.na(p1[1]) & !is.na(p2[1])) p1 <- p2
if (is.na(p2[1]) & !is.na(p1[1])) p2 <- p1

if (is.na(pR1[1]) & !is.na(pR2[1])) pR1 <- pR2
if (is.na(pR2[1]) & !is.na(pR1[1])) pR2 <- pR1

if (is.na(pL1[1]) & !is.na(pL2[1])) pL1 <- pL2
if (is.na(pL2[1]) & !is.na(pL1[1])) pL2 <- pL1

if ((!is.na(p1[1])) & (!is.na(p2[1]))) {
################
# modÃ¨le sans RL
################

d <- data.frame(time=1:days.maximum, N2=rep(x = 1, times=days.maximum), N1=rep(x=0, times=days.maximum), N0=rep(x=0, times=days.maximum))
if (days.maximum >1) {
for (i in 2:days.maximum) {
d[i, "N2"] <- d[i-1, "N2"] * (1-p2[i-1])
d[i, "N1"] <- d[i-1, "N1"] * (1-p1[i-1]) + d[i-1, "N2"] * p2[i-1] * (1-p1[i-1])
d[i, "N0"] <- d[i-1, "N0"] + d[i-1, "N1"] * p1[i-1] + d[i-1, "N2"] * p2[i-1] * p1[i-1]
}
}
d[, -1] <- d[, -1]/rowSums(d[, -1])

} else {
################
# modÃ¨le avec RL
################

d <- data.frame(time=1:days.maximum, NLR=rep(x = 1, times=days.maximum), NL0=rep(x=0, times=days.maximum), N0R=rep(x=0, times=days.maximum), N00=rep(x=0, times=days.maximum))
if (days.maximum >1) {
for (i in 2:days.maximum) {
d[i, "NLR"] <- d[i-1, "NLR"] * (1-pR2[i-1]) * (1-pL2[i-1])
d[i, "NL0"] <- d[i-1, "NL0"] * (1-pR1[i-1]) + d[i-1, "NLR"] * pL2[i-1] * (1-pR1[i-1])
d[i, "N0R"] <- d[i-1, "N0R"] * (1-pL1[i-1]) + d[i-1, "NLR"] * pR2[i-1] * (1-pL1[i-1])
d[i, "N00"] <- d[i-1, "N00"] + d[i-1, "NL0"] * pL1[i-1] + d[i-1, "N0R"] * pR1[i-1] + d[i-1, "NLR"] * pL2[i-1] * pR1[i-1] + d[i-1, "NLR"] * pR2[i-1] * pL1[i-1]
}
}
d[, -1] <- d[, -1]/rowSums(d[, -1])
}

d <- d[match(t, d[, "time"]), ]
return(d)

}

if (method == "delta") {

# Bon, j'ai une matrice... mais maintenant je dois l'utliser !

VCov <- solve(Hessian)

par_hess <- par[colnames(VCov)]

gh <- data.frame(time=numeric(), value=numeric(),
mean=numeric(), se=numeric(),
'X2.5'=numeric(), 'X97.5'=numeric())

suppressPackageStartupMessages(  requireNamespace("nlWaldTest"))
message("Estimation of distribution using delta method")

try_N <- function(..., Time, model_before=NULL, model_after=NULL, pfixed=NULL, model) {

par <- c(...)

tgm <- Tagloss_cumul(Time, par=c(par, pfixed),
model_before=model_before,
model_after=model_after,
Hessian = NULL, model = NULL)

return(tgm[1, model])
}

nlConfint2 <- function (obj = NULL, texts, level = 0.95, coeff = NULL, Vcov = NULL,
df2 = NULL, x = NULL, parameters=NULL)
{
if (!is.null(obj)) {
co = try(coef(obj), silent = T)
cond = attr(co, "condition")
if (is.null(coeff) && (is.null(cond)))
coeff = co
vc = try(vcov(obj), silent = T)
cond2 = attr(vc, "condition")
if (is.null(Vcov) && (is.null(cond2)))
Vcov = vc
}
if (is.null(coeff)) {
if (is.null(obj))
mess = "Both  'obj' and 'coeff' are missing"
else {
clm <- class(obj)
part1 <- "There are no coef() methods for model objects of class \""
mess <- paste0(part1, clm, "\".\nInput the 'coeff' parameter.")
}
stop(mess)
}
if (is.null(Vcov)) {
if (is.null(obj))
mess = "Both  'obj' and 'Vcov' are missing"
else {
clm <- class(obj)
part1 <- "There are no vcov() methods for model objects of class \""
mess <- paste0(part1, clm, "\".\nInput the 'Vcov' parameter.")
}
stop(mess)
}
if (length(texts) > 1)
{ kkk = texts[1]
} else kkk = strsplit(texts[1], ";")[[1]]
kkkfl = as.formula(paste("~", kkk[1]))
vvss = setdiff(all.vars(kkkfl), "x")
texts = getFromNamespace(".smartsub", ns="nlWaldTest")(vvss, "b", texts)
if (length(texts) > 1)
{ ltext0 = texts
} else ltext0 = strsplit(texts, ";")[[1]]
texts1 = gsub("[", "", texts, fixed = T)
texts1 = gsub("]", "", texts1, fixed = T)
if (length(texts1) > 1)
{ ltext = texts1
} else ltext = strsplit(texts1, ";")[[1]]
r = length(ltext)
n = length(coeff)
namess = paste0("b", 1:n)
for (j in 1L:n) assign(namess[j], coeff[j])
if (!is.null(x)) {
nx = length(x)
namesx = paste0("x", 1:nx)
for (j in 1L:nx) assign(namesx[j], x[j])
}
hess = c()
for (i in 1L:r) {
if (!is.null(parameters)) {
fli <- as.formula(paste("~", paste0(gsub(")", "", ltext[i]), ", ",
parameters, ")")))
} else {
fli <- as.formula(paste("~", ltext[i]))
}
z = try(deriv(as.formula(fli), namess), silent = T)
if (inherits(z, "try-error")) {
tei = as.character(i)
tri2 = ", numerical derivatives were used in delta-method"
wate = paste0("Note: For function ", i, tri2)
message(wate)
if (!is.null(parameters)) {
ez = numericDeriv(quote(eval(parse(text = paste0(gsub(")", "", ltext[i]), ", ",
parameters, ")")))),
namess)
} else {
ez = numericDeriv(quote(eval(parse(text = ltext[i]))),
namess)
}
} else ez = eval(z)
hess = rbind(hess, hessj)
}

ddd = hess %*% Vcov %*% t(hess)
# On retire ce test
# matr = chol2inv(chol(ddd))
ses = sqrt(diag(ddd))

trydf = identical(df2, T)
if (trydf) {
isdf = try(df.residual(obj), silent = T)
df2 = isdf
if (is.null(df2)) {
wn = "Note: Failed to extract df for denomenator; z-intervals applied"
message(wn)
}
}
getFromNamespace(".getint", ns="nlWaldTest")(as.numeric(Rb), ltext0, ses, level, df = df2)
}

for (ti in t) {

if (is.null(model_before)) {
mb <- "model_before=NULL"
} else {
mb <- paste0("model_before=\"", model_before, "\"")
}

if (is.null(model_after)) {
mb <- paste0(mb, ", model_after=NULL")
} else {
mb <- paste0(mb, ", model_after=\"", model_after, "\"")
}

mb <- paste0(mb, ", model=\"", model, "\"")

mb <- paste0(mb, ", pfixed=")
mb <- paste0(mb, paste0("c(", paste(names(par_add), "=", unname(par_add), collapse = ", "), ")"))
}

ghi <- suppressMessages ( nlConfint2(texts=paste0(c("try_N(", paste0("b[", seq_along(par_hess), "], ", collapse=""),"Time=", ti, ")"),
collapse = ""),
parameters = mb,
level = 0.95, coeff = par_hess,
Vcov = VCov, df2 = TRUE)
)

gh <- rbind(gh, data.frame(time=ti, value=ghi[, 1], mean=NA, sd=NA,
'X2.5'=ghi[, 2], 'X97.5'=ghi[, 3]))

}

return(gh)

}

# Je suis en SE, hessian, mcmc ou peudohessian

dfr <- RandomFromHessianOrMCMC(method = method, mcmc=mcmc, Hessian = Hessian,
fitted.parameters = par, replicates = replicates,
fn=Tagloss_cumul, ParTofn = "par",  model_before = model_before,
model_after = model_after,
model=model, x=x, t=t)
return(dfr\$quantiles)

}
```

## Try the phenology package in your browser

Any scripts or data that you put into this service are public.

phenology documentation built on Sept. 11, 2024, 6:07 p.m.