IPFit: Fit a model of Internesting Period for marine turtles.

In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

Fit a model of Internesting Period for marine turtles.

Description

This function fits a model of internesting period using maximum likelihood or using Metropolis-Hastings algorithm with Bayesian model.
The fit using maximum likelihood is not the best strategy because the objective function is based on a stochastic model (and then a single set of parameters does not produce exactly the same output each time). The use of Metropolis-Hastings algorithm (a Markov chain Monte Carlo method) should be prefered.

Usage

IPFit(
  x = NULL,
  fixed.parameters = NULL,
  data = stop("Formated data must be provided"),
  method = c("Nelder-Mead", "BFGS"),
  control = list(trace = 1, REPORT = 100, maxit = 500),
  itnmax = c(500, 100),
  hessian = TRUE,
  verbose = TRUE,
  parallel = TRUE,
  model = c("MH", "ML"),
  parametersMH,
  n.iter = 10000,
  n.chains = 1,
  n.adapt = 100,
  thin = 30,
  trace = TRUE,
  adaptive = TRUE,
  adaptive.lag = 500,
  adaptive.fun = function(x) {
     ifelse(x > 0.234, 1.3, 0.7)
 },
  intermediate = NULL,
  filename = "intermediate.Rdata"
)

Arguments

x	Initial parameters to be fitted
fixed.parameters	Parameters that are fixed.
data	Data as a vector
method	Method to be used by optimx()
control	List of controls for optimx()
itnmax	A vector with maximum iterations for each method.
hessian	Logical to estimate SE of parameters
verbose	If TRUE, show the parameters for each tested model
parallel	If TRUE, will use parallel computing
model	Can be ML for Maximum likelihood or MH for Metropolis Hastings
parametersMH	The priors. See MHalgoGen
n.iter	See MHalgoGen
n.chains	See MHalgoGen
n.adapt	See MHalgoGen
thin	See MHalgoGen
trace	See MHalgoGen
adaptive	See MHalgoGen
adaptive.lag	See MHalgoGen
adaptive.fun	See MHalgoGen
intermediate	See MHalgoGen
filename	See MHalgoGen

Details

IPFit fit a model of Internesting Period for marine turtles.

Value

Return a list of class IP with the fit informations and the fitted model.

Author(s)

Marc Girondot [e-mail]

See Also

Other Model of Internesting Period: IPModel(), IPPredict(), plot.IP(), summary.IP()

Examples

## Not run: 
library(phenology)
# Example
data <- structure(c(`0` = 0, `1` = 47, `2` = 15, `3` = 6, `4` = 5, `5` = 4, 
                    `6` = 2, `7` = 5, `8` = 57, `9` = 203, `10` = 205, `11` = 103, 
                    `12` = 35, `13` = 24, `14` = 12, `15` = 10, `16` = 13, `17` = 49, 
                    `18` = 86, `19` = 107, `20` = 111, `21` = 73, `22` = 47, `23` = 30, 
                    `24` = 19, `25` = 17, `26` = 33, `27` = 48, `28` = 77, `29` = 83, 
                    `30` = 65, `31` = 37, `32` = 27, `33` = 23, `34` = 24, `35` = 22, 
                    `36` = 41, `37` = 42, `38` = 44, `39` = 33, `40` = 39, `41` = 24, 
                    `42` = 18, `43` = 18, `44` = 22, `45` = 22, `46` = 19, `47` = 24, 
                    `48` = 28, `49` = 17, `50` = 18, `51` = 19, `52` = 17, `53` = 4, 
                    `54` = 12, `55` = 9, `56` = 6, `57` = 11, `58` = 7, `59` = 11, 
                    `60` = 12, `61` = 5, `62` = 4, `63` = 6, `64` = 11, `65` = 5, 
                    `66` = 6, `67` = 7, `68` = 3, `69` = 2, `70` = 1, `71` = 3, `72` = 2, 
                    `73` = 1, `74` = 2, `75` = 0, `76` = 0, `77` = 3, `78` = 1, `79` = 0, 
                    `80` = 2, `81` = 0, `82` = 0, `83` = 1), Year = "1994", 
                    Species = "Dermochelys coriacea", 
                    location = "Yalimapo beach, French Guiana", 
                    totalnumber = 2526L, class = "IP")
  par(mar=c(4, 4, 1, 1)+0.4)
  plot(data, xlim=c(0,100))
  text(100, 190, labels=bquote(italic(.(attributes(data)$Species))), pos=2)
  text(100, 150, labels=attributes(data)$location, pos=2, cex=0.8)
  text(100, 110, labels=paste0(as.character(attributes(data)$totalnumber), " females"), pos=2)
  
######### Fit using Maximum-Likelihood

par <- c(meanIP = 9.8229005713237623, 
         sdIP = 0.079176011861863474, 
         minIP = 6.8128364577569309, 
         pAbort = 1.5441529841959203, 
         meanAbort = 2.7958742380756121, 
         sdAbort = 0.99370406770777175, 
         pCapture = -0.80294884905867658, 
         meanECF = 4.5253772889275758, 
         sdECF = 0.20334743335612529)

fML <- IPFit(x=par, 
             fixed.parameters=c(N=20000),
             data=data, 
             verbose=FALSE, 
             model="ML")

# Plot the fitted ECF
plot(fML, result="ECF")

# Plot the Internesting Period distribution
plot(fML, result="IP")

# Plot the distribution of days between tentatives
plot(fML, result="Abort", xlim=c(0, 15))
#' 
######### Fit using ML and non parametric ECF

par <- c(ECF.2 = 0.044151921569961131, 
         ECF.3 = 2.0020778325280748, 
         ECF.4 = 2.6128345101617083, 
         ECF.5 = 2.6450582416622375, 
         ECF.6 = 2.715198206774927, 
         ECF.7 = 2.0288031327239904, 
         ECF.8 = 1.0028041546528881, 
         ECF.9 = 0.70977432157689235, 
         ECF.10 = 0.086052204035003091, 
         ECF.11 = 0.011400419961702518, 
         ECF.12 = 0.001825219438794076, 
         ECF.13 = 0.00029398731859899116, 
         ECF.14 = 0.002784886479846703, 
         meanIP = 9.9887100433529721, 
         sdIP = 0.10580250625108811, 
         minIP = 6.5159124624132048, 
         pAbort = 2.5702251748938956, 
         meanAbort = 2.2721679285648841, 
         sdAbort = 0.52006431730489933, 
         pCapture = 0.079471782729506113)
         
fML_NP <- IPFit(x=par, 
             fixed.parameters=c(N=20000),
             data=data, 
             verbose=FALSE, 
             model="ML")
             
par <- fML_NP$ML$par

fML_NP <- IPFit(x=par, 
             fixed.parameters=c(N=1000000),
             data=data, 
             verbose=FALSE, 
             model="ML")
             
par <- c(ECF.2 = 0.016195025683080871, 
         ECF.3 = 2.0858089267994315, 
         ECF.4 = 3.1307578727979348, 
         ECF.5 = 2.7495760827322622, 
         ECF.6 = 2.8770821670450939, 
         ECF.7 = 2.1592708144943145, 
         ECF.8 = 1.0016227335391867, 
         ECF.9 = 0.80990178270345259, 
         ECF.10 = 0.081051214954249967, 
         ECF.11 = 0.039757901443389344, 
         ECF.12 = 6.3324056808464527e-05, 
         ECF.13 = 0.00037500864146146936, 
         ECF.14 = 0.0010383506745475582, 
         meanIP = 10.004121090603523, 
         sdIP = 0.10229422354470977, 
         minIP = 6.5051758088487883, 
         pAbort = 2.5335985958484839, 
         meanAbort = 2.3145895392189173, 
         sdAbort = 0.51192514362374153, 
         pCapture = 0.055440514236842105, 
         DeltameanIP = -0.046478049165483697)

fML_NP_Delta <- IPFit(x=par, 
             fixed.parameters=c(N=20000),
             data=data, 
             verbose=FALSE, 
             model="ML")
             
par <- fML_NP_Delta$ML$par
             
fML_NP_Delta <- IPFit(x=par, 
             fixed.parameters=c(N=1000000),
             data=data, 
             verbose=FALSE, 
             model="ML")
             
# Test for stability of -Ln L value according to N
grandL.mean <- NULL
grandL.sd <- NULL
N <- c(10000, 20000, 30000, 40000, 50000,
            60000, 70000, 80000, 90000,  
            100000, 200000, 300000, 400000, 500000, 
            600000, 700000, 800000, 900000,  
            1000000)
for (Ni in N) {
    print(Ni)
    smallL <- NULL
    for (replicate in 1:100) {
         smallL <- c(smallL, 
         getFromNamespace(".IPlnL", ns="phenology")
                 (x=par, fixed.parameters=c(N=Ni), data=data))
    }
    grandL.mean <- c(grandL.mean, mean(smallL))
    grandL.sd <- c(grandL.sd, sd(smallL))
}

grandL.mean <- c(242.619750064524, 239.596145944548, 238.640010536147, 237.965573853263, 
237.727506424543, 237.240740566494, 237.527948232993, 237.297225856515, 
237.17073080938, 237.103397800143, 236.855939567838, 
236.704861853456, 236.82264801458, 236.606065021519, 236.685930841831, 
236.697562908131, 236.568003663293, 236.58097471402, 236.594282543024
)
grandL.sd <- c(6.54334049298099, 3.04916614991682, 2.57932397492509, 2.15990307710982, 
1.59826856034413, 1.54505295915354, 1.59734964880484, 1.41845032728396, 
1.43096821211286, 1.20048923027244, 0.912467350448495, 
0.75814052890774, 0.668841336554019, 0.539505594152166, 0.554662419326559, 
0.501551009304687, 0.415199780254872, 0.472274287714195, 0.386237047201706
)

plot_errbar(x=N, y=grandL.mean, errbar.y = 2*grandL.sd, 
            xlab="N", ylab="-Ln L (2 SD)", bty="n", las=1)
             
# Plot the fitted ECF
plot(fML_NP_Delta, result="ECF")

# Plot the Internesting Period distribution
plot(fML_NP_Delta, result="IP")

# Plot the distribution of days between tentatives
plot(fML_NP_Delta, result="Abort", xlim=c(0, 15))

print(paste("Probability of capture", invlogit(-fML_NP_Delta$ML$par["pCapture"])))
# Confidence interval at 95%
print(paste(invlogit(-fML_NP_Delta$ML$par["pCapture"]-1.96*fML_NP_Delta$ML$SE["pCapture"]), "-", 
invlogit(-fML_NP_Delta$ML$par["pCapture"]+1.96*fML_NP_Delta$ML$SE["pCapture"])))

print(paste("Probability of abort", invlogit(-fML_NP_Delta$ML$par["pAbort"])))
# Confidence interval at 95%
print(paste(invlogit(-fML_NP_Delta$ML$par["pAbort"]-1.96*fML_NP_Delta$ML$SE["pAbort"]), "-", 
invlogit(-fML_NP_Delta$ML$par["pAbort"]+1.96*fML_NP_Delta$ML$SE["pAbort"])))
             
compare_AIC(parametric=fML$ML, 
            nonparameteric=fML_NP$ML, 
            nonparametricDelta=fML_NP_Delta$ML)

######### Fit using Metropolis-Hastings algorithm
# ECF.1 = 1 is fixed
par <- c(ECF.2 = 0.044151921569961131, 
         ECF.3 = 2.0020778325280748, 
         ECF.4 = 2.6128345101617083, 
         ECF.5 = 2.6450582416622375, 
         ECF.6 = 2.715198206774927, 
         ECF.7 = 2.0288031327239904, 
         ECF.8 = 1.0028041546528881, 
         ECF.9 = 0.70977432157689235, 
         ECF.10 = 0.086052204035003091, 
         ECF.11 = 0.011400419961702518, 
         ECF.12 = 0.001825219438794076, 
         ECF.13 = 0.00029398731859899116, 
         ECF.14 = 0.002784886479846703, 
         meanIP = 9.9887100433529721, 
         sdIP = 0.10580250625108811, 
         minIP = 6.5159124624132048, 
         pAbort = 2.5702251748938956, 
         meanAbort = 2.2721679285648841, 
         sdAbort = 0.52006431730489933, 
         pCapture = 0.079471782729506113)

df <- data.frame(Density=rep("dunif", length(par)), 
Prior1=c(rep(0, 13), 8, 0.001, 0, -8, 0, 0.001, -8), 
Prior2=c(rep(10, 13), 12, 1, 10, 8, 2, 1, 8), 
SDProp=unname(c(ECF.2 = 6.366805760909012e-05, 
                ECF.3 = 6.366805760909012e-05, 
                ECF.4 = 6.366805760909012e-05, 
                ECF.5 = 6.366805760909012e-05, 
                ECF.6 = 6.366805760909012e-05, 
                ECF.7 = 6.366805760909012e-05, 
                ECF.8 = 6.366805760909012e-05, 
                ECF.9 = 6.366805760909012e-05, 
                ECF.10 = 6.366805760909012e-05, 
                ECF.11 = 6.366805760909012e-05, 
                ECF.12 = 6.366805760909012e-05, 
                ECF.13 = 6.366805760909012e-05, 
                ECF.14 = 6.366805760909012e-05, 
                meanIP = 6.366805760909012e-05, 
                sdIP = 6.366805760909012e-05, 
                minIP = 6.366805760909012e-05, 
                pAbort = 6.366805760909012e-05, 
                meanAbort = 6.366805760909012e-05, 
                sdAbort = 6.366805760909012e-05, 
                pCapture = 6.366805760909012e-05)),               
Min=c(rep(0, 13), 8, 0.001, 0, -8, 0, 0.001, -8), 
Max=c(rep(10, 13), 12, 1, 10, 8, 2, 1, 8),
Init=par, stringsAsFactors = FALSE)
rownames(df)<- names(par)

fMH <- IPFit(parametersMH=df, 
fixed.parameters=c(N=10000),
data=data, 
verbose=FALSE, 
n.iter = 10000,
n.chains = 1, n.adapt = 100, thin = 1, trace = TRUE,
adaptive = TRUE, 
model="MH")

# Plot the fitted ECF
plot(fMH, result="ECF")


## End(Not run)

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