# Fit a model of Internesting Period for marine turtles. 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. ```r 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 ## Returns Return a list of class IP with the fit informations and the fitted model. ## Details IPFit fit a model of Internesting Period for marine turtles. ## Examples ```r ## 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) ``` ## See Also Other Model of Internesting Period: `IPModel()`, `IPPredict()`, `plot.IP()`, `summary.IP()` ## Author(s) Marc Girondot marc.girondot@gmail.com

IPFit function