IPFit: 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"
)

Value

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

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

Author

Marc Girondot marc.girondot@gmail.com

Details

IPFit fit a model of Internesting Period for marine turtles.

Examples

Run this code

if (FALSE) {
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")

}

Run the code above in your browser using DataLab