\documentclass{article}
\newcommand{\R}{{\sf R}}
\newcommand{\codef}[1]{{\tt #1}}
\title{\R\ code for mid-domain models}
\begin{document}

\section{Spatial mid-domain model}

Detailed explanation of the model is provided in the text. Some of the
terms used in the code are: 
\codef{rangepar}=nesting range parameters (gamma distribution shape and scale in km);
\codef{beachsize}= beach length (km); 
\codef{segsize}=the size of each beach segment (km);
\codef{rangectrs}=center point of nesting ranges; 
\codef{nclutch}=number of clutches per turtle;
\codef{nturtles}=number of turtles; 
\codef{loc}=location;
\codef{nsim}=number of simulations;
\codef{nsimvec}=vector of simulation results 
\codef{rtolvec}=vector for tolerance values for the
optimization; 
\codef{par}=parameter.

Observed mean number of clutches at each spatial location on the beach:
<<>>=
X<-c(1130,950,842,670,583,338,238,1080,1339,1447,2066,2052,2592,2858,
     2887,2923,3319,3074,3730,3650,3434,2844,2808,2765,2635,2693,
     2390,2160,2210,1994,1843,1649,1404,1253,1202,943)
@ 

Choose random range size(s) from the specified
gamma distribution:
<<>>=
rrange <- function(rangepar,n=1) {
   rgamma(n,shape=rangepar[1],scale=rangepar[2])
}
@ 

Simulating nesting range sizes, and random placement of nesting ranges and 
clutches: 
<<>>=
simdist2S <- function(rangepar,nturtles=24000,nclutch=3,
                      beachsize=18,segsize=0.5) {
  segvec <- numeric(beachsize/segsize)
  rangesize <- rrange(rangepar,n=nturtles)  
  rangectrs <- runif(nturtles,min=rangesize/2,max=beachsize-rangesize/2)
  minvec <- rep(rangectrs-rangesize/2,rep(nclutch,nturtles))
  maxvec <- rep(rangectrs+rangesize/2,rep(nclutch,nturtles))
  loc <- runif(nturtles*nclutch,min=minvec,max=maxvec)
  as.vector(table(cut(loc,seq(0,beachsize,by=segsize))))
}
@ 

Function to determine mean distribution of nests and 95\% confidence intervals
for specified range parameters:
<<>>=
simdist2.env <- function(rangepars,
                         nsim=1000,beachsize=18,segsize=0.5,rpt=0,...) {
  nseg <- round(beachsize/segsize)
  simres <- matrix(nrow=nseg,ncol=nsim)
  for (i in 1:nsim) {
    simres[,i] <- simdist2S(rangepars,beachsize=beachsize,segsize=segsize,...)
    if (rpt>0) {
      if (i %% rpt == 0) cat(i,"\n")
    }
  }
  meanval <- apply(simres,1,mean)
  conf <- apply(simres,1,quantile,c(0.025,0.975))
  res <- list(x=1:nseg,mean=meanval,lower=conf[1,],upper=conf[2,],n=nsim)
  class(res) <- "simdist.env"
  res
}
@ 

Calculate the likelihood of the observed numbers of clutches
based on the simulated (mean) nest 
distribution in segments along the beach: 
<<>>=
simdist2.lik <- function(rangepars,data=X,
                        nsim=500,beachsize=18,segsize=0.5) {
  mvals <- simdist2.env(rangepars,nsim=nsim,
                        beachsize=beachsize,
                        segsize=segsize,nturtles=24000,
                        nclutch=3)
  -sum(dpois(X,mvals$mean,log=TRUE))
}
@

Using the Nelder-Mead optimization method (default
for \R\ \codef{optim()} function) with a series of tolerance values and 
numbers of simulations to determine the gamma parameter values for the best-
fitting mean curve for the observed data :

<<>>=
optfun <- function(start,fn=simdist2.lik,data=X,
                   nsimvec=c(10000,15000,20000),
                   rtolvec=c(0.001,0.001,0.001),trace=0) {
  resmat <- matrix(nrow=length(nsimvec),ncol=7)
  O1 <- list(par=start)
  for (i in 1:length(nsimvec)) {
    O1 <- optim(par=O1$par,fn=fn,data=data,nsim=nsimvec[i],
                control=list(reltol=rtolvec[i],trace=trace))
    resmat[i,] <- c(O1$par,O1$value,rtolvec[i],nsimvec[i],O1$convergence,
                    O1$counts[1])
  }
  dimnames(resmat) <- list(NULL,
                           c("shape","scale","lik",
                             "tol","nsim","conv","evals"))
  list(opt=O1,trace=resmat)
}
@ 

Determining the best gamma parameter estimates by starting with shape and scale 
parameters equivalent to 1:

<<eval=FALSE>>=
O220  <- optfun(c(1,1))
O220t <- O220$trace
s.020 <- simdist2.env(c(1,1),nsim=20000,nturtles=24000)
s.120 <- simdist2.env(O220$trace[1,1:2],nsim=20000,nturtles=24000)
s.220 <- simdist2.env(O220$trace[2,1:2],nsim=20000,nturtles=24000)
s.320 <- simdist2.env(O220$trace[3,1:2],nsim=20000,nturtles=24000)
@ 
                        
Plotting simulation output:
<<>>=
plot.simdist.env <- function(x,add=FALSE,col=par("fg"),lty=c(2,1,2),...) {
  if (!add) { matplot(x$x,cbind(x$lower,x$mean,x$upper),lty=lty,
        type="l",xlab="Distance",ylab="Density",col=col)
  } else {
    matlines(x$x,cbind(x$lower,x$mean,x$upper),lty=lty,col=col)
  }
}
@ 

<<eval=FALSE>>=
plot(X,ylim=c(0,4000),pch=16)
plot(s.0,add=TRUE,col=2)
plot(s.1,add=TRUE,col=3)
plot(s.2,add=TRUE,col=4)
plot(s.3,add=TRUE,col=5)
@ 

\section{Temporal mid-domain model}

Parameters:
\codef{seasonlt}=length of season (days);
\codef{nclutches}=number of clutches; 
\codef{lambda}=the mean number of clutches laid by a turtle; 
\codef{loc}=date each clutch was laid; 
\codef{nestint}=nesting interval; 
\codef{rangectr}=range center; 
\codef{nsim}=number of simulations. 

Simulating nesting range size and random placement of nesting ranges and clutches 
for a total of 72,000 clutches in the beach:
<<>>=
simdist <- function(clutches=72000,seasonlt=123, lambda=3) {
  loc <- numeric(clutches) 
  e <- 1 
  i=0
  while(i <clutches) {
    nclutches=rpois(1,lambda)
    if(nclutches>0 & nclutches<8){
      nestint =rbeta(1,shape1=4.2,shape2=5.5)
      nestint=round((nestint*12)+7)
      rangesize<-nestint*(nclutches-1)
      rangectr=runif(1,min=rangesize/2,max=seasonlt-rangesize/2)
      cat(rangesize," ",nestint,"\n")
      for (j in seq(from=-rangesize/2, to=rangesize/2, by=nestint))  {
        loc[e] <- rangectr+j
        e <- e+1
      }
      i<-i+nclutches
    }
  }
  sort(loc)
}
@ 

Determining mean distribution of nests and 95\% confidence intervals:
<<>>=
simdist.env <- function(nsim=10000,seasonlt=123,npts=seasonlt,rpt=10,...) {
  simres <- matrix(nrow=npts,ncol=nsim)
  for (i in 1:nsim) {
simres[,i] <-as.numeric(table(factor(round(simdist(seasonlt=seasonlt,...))
,levels=1:seasonlt)))
    if (rpt>0) {
      if (i %% rpt == 0) cat(i,"\n")
    }
  }
  x <- seq(0,seasonlt,length=npts)
  meanval <- apply(simres,1,mean)
  conf <- apply(simres,1,quantile,c(0.025,0.975))
  res <- list(x=x,mean=meanval,lower=conf[1,],upper=conf[2,],n=nsim)
  class(res) <- "simdist.env"
  res
}
@ 

Running the model for 1000 simulations:
<<eval=FALSE>>=
L1e <- simdist.env(nsim=1000)
@ 


\end{document}