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

#### Documented in Tagloss_model

```#' Tagloss_model returns the daily rate of tag loss.
#' @title Return the daily rate of tag loss.
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Return the daily rate of tag loss if hessian is null or a data.frame with distribution of daily 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 model The model of parameter to be used, can be 1, 2, L1, L2, R1 or R2
#' @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 method Can be NULL, "delta", "SE", "Hessian", "MCMC", or "PseudoHessianFromMCMC"
#' @param replicates Number of replicates to estimate se of output
#' @description This function compute a model of daily tag loss rate for days t
#' based on a set of parameters, par or a fitted tag loss model in x.\cr
#' Parameters are described in \code{\link{Tagloss_fit}}.
#' @family Model of Tag-loss
#' @examples
#' \dontrun{
#' library(phenology)
#'
#' # Example
#' t <- 1:1000
#' par <- c(D1=200, D2D1=100, D3D2=200,
#'          A=-logit(0.02), B=-logit(0.05), C=-logit(0.07))
#' y <- Tagloss_model(t, par, model="1")
#' plot(x=t, y, type="l")
#' par <- c(D1_1=200, D2D1_1=100, D3D2_1=200,
#'          A_1=-logit(0.02), B_1=-logit(0.05), C_1=-logit(0.07))
#' y <- Tagloss_model(t, par, model="1")
#' phenology:::plot.Tagloss(x=list(), t=1:1000, fitted.parameters=par, model="1")
#'
#' # Fig1A in Rivalan et al. 2005 (note an error for a0; a0 must be negative)
#' par <- c(a0=-1E5, a1=-2000, a2=0, a3=2*max(t), a4=0.1)
#' y <- Tagloss_model(t, par)
#' plot(x=t, y, type="l")
#'
#' # Fig1B in Rivalan et al. 2005
#' par <- c(a0=-0.5, a1=-2000, a2=-0.001, a3=0, a4=0.1)
#' y <- Tagloss_model(t, par)
#' plot(x=t, y, type="l")
#'
#' # Fig1C in Rivalan et al. 2005
#' par <- c(a0=-1, a1=-6, a2=0, a3=0, a4=0)
#' y <- Tagloss_model(t, par)
#' plot(x=t, y, type="l")
#'
#' # Fig1D in Rivalan et al. 2005
#' par <- c(a0=-1, a1=-6, a2=0, a3=0, a4=0.1)
#' y <- Tagloss_model(t, par)
#' plot(x=t, y, type="l")
#'
#' # Fig1E in Rivalan et al. 2005
#' par <- c(a0=-0.1, a1=-10, a2=-0.2, a3=60, a4=0.1)
#' y <- Tagloss_model(t, par)
#' plot(x=t, y, type="l")
#'
#' # Fig1F in Rivalan et al. 2005
#' par <- c(a0=-0.1, a1=-10, a2=0.2, a3=60, a4=0.1)
#' y <- Tagloss_model(t, par)
#' plot(x=t, y, type="l")
#'
#' # Example with fitted data
#' data_f_21 <- Tagloss_format(outLR, model="21")
#' # Without the N20 the computing is much faster
#' data_f_21_fast <- subset(data_f_21, subset=(is.na(data_f_21\$N20)))
#' par <- c('D1_2' = 49.086835072129126,
#'          'D2D1_2' = 1065.0992647723231,
#'          'D3D2_2' = 6.15531475922079,
#'          'A_2' = 5.2179675647973758,
#'          'B_2' = 8.0045560376751386,
#'          'C_2' = 8.4082505219581876,
#'          'D1_1' = 177.23337287498103,
#'          'D2D1_1' = 615.42690323741033,
#'          'D3D2_1' = 2829.0806609455867,
#'          'A_1' = 28.500118091731551,
#'          'B_1' = 10.175426055942701,
#'          'C_1' = 6.9616630417169398)
#' o <- Tagloss_fit(data=data_f_21_fast, fitted.parameters=par)
#' t <- 1:10
#' y <- Tagloss_model(t, o\$par, model="1")
#' y <- Tagloss_model(t, x=o, method=NULL, model="1")
#' y <- Tagloss_model(t, x=o, method="Hessian", model="1", replicates=1000)
#' }
#' @export

Tagloss_model <- function(t=NULL                                                      ,
par=NULL                                                    ,
Hessian=NULL                                                ,
mcmc = NULL                                                 ,
model_before = NULL                                         ,
model_after = NULL                                          ,
model=stop("You must specify which tag loss rate you want."),
method=NULL                                                 ,
replicates=NULL                                             ,
x=NULL                                                      ) {

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) & is.null(par)) {
stop("Both par and x cannot be null at the same time")
}

if (is.null(t)) {
if (is.null(x)) {
stop("Both t and x cannot be null at the same time")
} else {
t <- 1:x\$days.maximum
}
}

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
}

if (!is.null(model)) {
model <- toupper(model)
model <- match.arg(model, choices = c("1", "2", "L1", "L2", "R1", "R2"))
} else {
stop("Model cannot be NULL.")
}

if (method=="hessian" & is.null(Hessian)) {
stop("Hessian method cannot be used with an empty Hessian matrix.")
}

if ((method=="mcmc" | method=="pseudohessianfrommcmc") & is.null(mcmc)) {
stop("MCMC methods cannot be used without MCMC parameter.")
}

days.maximum <- max(t)

if (method=="null") {

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

p1 <- NULL
p2 <- NULL
pL1 <- NULL
pR1 <- NULL
pL2 <- NULL
pR2 <- NULL

# Casale 2017 Model 1
if (any(names(par) == "CasaleModelIc_1")) p1 <- rep(par["CasaleModelIc_1"], length(t))
if (any(names(par) == "CasaleModelIc_2")) p2 <- rep(par["CasaleModelIc_2"], length(t))
if (any(names(par) == "CasaleModelIc")) p2 <- p1 <- rep(par["CasaleModelIc"], length(t))

# Casale 2017 Model II
if (any(names(par) == "CasaleModelIIa0_1")) p1 <- (1-par["CasaleModelIIa4_1"])/(1+exp(par["CasaleModelIIa0_1"]*(par["CasaleModelIIa1_1"]-t))) + par["CasaleModelIIa4_1"]
if (any(names(par) == "CasaleModelIIa0_2")) p2 <- (1-par["CasaleModelIIa4_2"])/(1+exp(par["CasaleModelIIa0_2"]*(par["CasaleModelIIa1_2"]-t))) + par["CasaleModelIIa4_2"]
if (any(names(par) == "CasaleModelIIa0")) p1 <- p2 <- (1-par["CasaleModelIIa4"])/(1+exp(par["CasaleModelIIa0"]*(par["CasaleModelIIa1"]-t))) + par["CasaleModelIIa4"]

# Casale 2017 Model III
if (any(names(par) == "CasaleModelIIIa0_1")) p1 <- par["CasaleModelIIIa4_1"]/(1+exp(par["CasaleModelIIIa0_1"]*(par["CasaleModelIIIa1_1"]-t)))
if (any(names(par) == "CasaleModelIIIa0_2")) p2 <- par["CasaleModelIIIa4_2"]/(1+exp(par["CasaleModelIIIa0_2"]*(par["CasaleModelIIIa1_2"]-t)))
if (any(names(par) == "CasaleModelIIIa0")) p1 <- p2 <- par["CasaleModelIIIa4"]/(1+exp(par["CasaleModelIIIa0"]*(par["CasaleModelIIIa1"]-t)))

# Casale 2017 Model IV
if (any(names(par) == "CasaleModelIVa0_1")) {
mint <- par["CasaleModelIVa4_1"]/(1+exp(par["CasaleModelIVa2_1"]*(par["CasaleModelIVa3_1"]-t)))
p1 <- (0.01-mint)/(1+exp(par["CasaleModelIVa0_1"]*(par["CasaleModelIVa1_1"]-t))) + mint
}
if (any(names(par) == "CasaleModelIVa0_2")) {
mint <- par["CasaleModelIVa4_2"]/(1+exp(par["CasaleModelIVa2_2"]*(par["CasaleModelIVa3_2"]-t)))
p2 <- (0.01-mint)/(1+exp(par["CasaleModelIVa0_2"]*(par["CasaleModelIVa1_2"]-t))) + mint
}
if (any(names(par) == "CasaleModelIVa0")) {
mint <- par["CasaleModelIVa4"]/(1+exp(par["CasaleModelIVa2"]*(par["CasaleModelIVa3"]-t)))
p1 <- p2 <- (0.01-mint)/(1+exp(par["CasaleModelIVa0"]*(par["CasaleModelIVa1"]-t))) + mint
}
# Casale 2017 Model 5
if (any(names(par) == "CasaleModelVa0_1")) {
maxt <- (0.01-par["CasaleModelVa4_1"])/(1+exp(par["CasaleModelVa2_1"]*(par["CasaleModelVa3_1"]-t))) + par["CasaleModelVa4_1"]
p1 <- maxt/(1+exp(par["CasaleModelVa0_1"]*(par["CasaleModelVa1_1"]-t)))
}
if (any(names(par) == "CasaleModelVa0_2")) {
maxt <- (0.01-par["CasaleModelVa4_2"])/(1+exp(par["CasaleModelVa2_2"]*(par["CasaleModelVa3_2"]-t))) + par["CasaleModelVa4_2"]
p2 <- maxt/(1+exp(par["CasaleModelVa0_2"]*(par["CasaleModelVa1_2"]-t)))
}
if (any(names(par) == "CasaleModelVa0")) {
maxt <- (0.01-par["CasaleModelVa4"])/(1+exp(par["CasaleModelVa2"]*(par["CasaleModelVa3"]-t))) + par["CasaleModelVa4"]
p1 <- p2 <- maxt/(1+exp(par["CasaleModelVa0"]*(par["CasaleModelVa1"]-t)))
}

# modÃ¨le Rivalan 2005
if (any(names(par) == "a0_1")) {
delta <- par["delta_1"]
if (is.na(delta)) delta <- 0
mint <- par["a4_1"]/(1+exp(par["a2_1"]*(par["a3_1"]-t)))
p1 <- (1-mint)/(1+exp(par["a0_1"]*(par["a1_1"]-t))) + mint + delta
}
if (any(names(par) == "a0_2")) {
delta <- par["delta_2"]
if (is.na(delta)) delta <- 0
mint <- par["a4_2"]/(1+exp(par["a2_2"]*(par["a3_2"]-t)))
p2 <- (1-mint)/(1+exp(par["a0_2"]*(par["a1_2"]-t))) + mint + delta
}
if (any(names(par) == "a0")) {
delta <- par["delta"]
if (is.na(delta)) delta <- 0
mint <- par["a4"]/(1+exp(par["a2"]*(par["a3"]-t)))
p1 <- p2 <- (1-mint)/(1+exp(par["a0"]*(par["a1"]-t))) + mint + delta
}

suffixe <- c("", "_1", "_2", "_L1", "_R1", "_L2", "_R2")

for (suf in suffixe) {

if (any(names(par) == paste0("D1", suf))) {
if (!is.na(par[paste0("D1", suf)])) D1 <- abs(par[paste0("D1", suf)]) else D1 <- 0
if (!is.na(par[paste0("D2D1", suf)])) D2D1 <- abs(par[paste0("D2D1", suf)]) else D2D1 <- max(t)+1
D2 <- D1 + D2D1
if (!is.na(par[paste0("D3D2", suf)])) D3D2 <- abs(par[paste0("D3D2", suf)]) else D3D2 <- max(t)+1
D3 <- D2 + D3D2
if (!is.na(par[paste0("A", suf)])) A <- invlogit(-par[paste0("A", suf)]) else A <- 0
if (A == 0) A <- 1E-9
if (A == 1) A <- 1-1E-9
if (!is.na(par[paste0("B", suf)])) B <- invlogit(-par[paste0("B", suf)]) else B <- 0
if (B == 0) B <- 1E-9
if (B == 1) B <- 1-1E-9
if (!is.na(par[paste0("C", suf)])) C <- invlogit(-par[paste0("C", suf)]) else C <- 0
if (C == 0) C <- 1E-9
if (C == 1) C <- 1-1E-9

if (is.na(par[paste0("delta", suf)])) delta <- 0 else delta <- par[paste0("delta", suf)]
if (D1==0) D1_2 <- 0.1 else D1_2 <- 2*D1
if (D3==D2) D3D2 <- (D3-0.1)*D2 else D3D2 <-D3-D2

suffi <- gsub("_", "", suf)
if (suf=="") {
p1 <- p2 <- ifelse(t<=D1,
(1+cos(pi*((t+D1)/(D1_2))))*(A-B)+B,
ifelse(t<=D2, B,
ifelse(t<=D3, (1+cos(pi*((t-D2)/(D3D2))))*0.5*(B-C)+C, C))) + delta
} else {
assign(paste0("p", suffi),  ifelse(t<=D1,
(1+cos(pi*((t+D1)/(D1_2))))*(A-B)+B,
ifelse(t<=D2, B,
ifelse(t<=D3, (1+cos(pi*((t-D2)/(D3D2))))*0.5*(B-C)+C, C))) + delta)
}
}
}

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

if (is.null(model)) {
z <- p1
} else {
z <- get(paste0("p", model))
}

return(as.numeric(ifelse(z<0, 1E-9, ifelse(z>1, 1-1E-9, z))))

}

# Si je suis ici c'est que je dois calculer une erreur standard

if (method == "delta") {
if (!(is.element('nlWaldTest', installed.packages()[,1]))) stop("nlWaldTest package must be installed to use the delata method.")

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

VCov <- solve(Hessian)
par_hess <- par
par_hess <- par_hess[colnames(VCov)]

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

if (any(diag(VCov)<0)) stop("Hessian matrix is not correct; probably you are not at maximum likelihood.")

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

par <- c(...)

# print(d(par))

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

return(tgm)
}

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, ")")))),
theta=namess)
} else {
ez = numericDeriv(quote(eval(parse(text = ltext[i]))),
theta=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, "\"")
}

if (!is.null(model)) {
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_g2(", 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)
}

# Ce n'est pas null et pas delta; je peux utiliser RandomFromHessianOrMCMC

dfr <- RandomFromHessianOrMCMC(method = method, mcmc=mcmc, Hessian = Hessian,
fitted.parameters = par, replicates = replicates,
fn=Tagloss_model, 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.