simDSpat: Simulate a distance sample from a specified spatial point...
In DSpat: Spatial Modelling for Distance Sampling Data

Description Usage Arguments Value Author(s) See Also Examples

This is a wrapper function that calls all of the functions needed to simulate and sample a point process over a defined study.area with a specified covariates on a grid. In sequence it calls create.lines, lines_to_strips, simPts, and sample.points.

simDSpat(study.area=owin(xrange=c(0,100),yrange=c(0,100)),covariates,
         angle=90,nlines=10,spacing=10,width=1,int.formula=~1,
         int.par=1,model="exp",cor.par=NULL,EN=1000,detfct=hndetfct,
         det.formula=~1,det.par=log(width/3),showplot=FALSE,showlines=FALSE,
         showpts=FALSE,pts=NULL,...)

`study.area`	owin class defining area
`covariates`	a matrix with columns x,y and any number of covariates x and y are the mid points of the grid cells; the order of the rows must match the formulation for function im
`angle`	angle of rotation in degrees anticlockwise from x-axis
`nlines`	number of lines
`spacing`	spacing distance between centerlines
`width`	full transect width
`int.formula`	formula for deriving expected intensity from covariates
`int.par`	parameters for intensity formula
`model`	either "exp" or "gauss" for exponential or Gaussian correlation
`cor.par`	parameters controlling clustering of points cor.par[1] sigma^2 cor.par[2]=alpha where cov(y1,y2)=sigma^2exp(-d^p/alpha)* and d is the distance between y1 and y2 and p=1 for exp and p=2 for gauss; if it is not specified then no additional clustering is included.
`EN`	expected number of points
`detfct`	detection function name
`det.formula`	formula of covariates to use for scale of distance if det.formula=~-1, uses a strip transect
`det.par`	parameters for the detection function
`showplot`	if TRUE show plot of the simulated points
`showlines`	if TRUE show lines and transects on the plot
`showpts`	if TRUE show points on the plot
`pts`	if not NULL use these points rather than generating new ones; this allows generation of a single set of points and evaluation of different sampling designs or intensity
`...`	parameters, if any, passed to plot

a list with elements

`lines`	lines dataframe with label,x0,y0,x1,y1,width where x0,y0 is beginning and x1,y1 is end of the line
`observations`	a dataframe of the coordinates of the observed points

Jeff Laake

simCovariates,simPts

# Code stored in a function to speed up package checking run by typing do.simDSpat
# Some portions of this code will not pass the packge check on Linux and Mac and will issue
# an error that the polygons intersect even though when run as an example, the
# error is not encountered; so polygon checking is turned off
do.simDSpat=function()
{
# Now that it is in a function shouldn't need following line
spatstat.options(checkpolygons=FALSE)
study.area=owin(poly=list(x=c(0,40,40,100,100,0),y=c(0,0,40,40,100,100)))
covariates = simCovariates(hab.range=30, probs=c(1/3,2/3))
simdata=simDSpat(study.area,covariates,int.formula=~factor(habitat),
                   int.par=c(0,1,2),angle=45,nlines=10,width=3,det.par=.1)
sim.dspat=dspat(int.formula=~factor(habitat),study.area=study.area,
               obs=simdata$observations,lines=simdata$lines,
               covariates=covariates,epsvu=c(1,.05))
summary(sim.dspat)
AIC(sim.dspat)
coef(sim.dspat)
mu.B <- integrate.intensity(sim.dspat,dimyx=100,se=TRUE)
cat('Abundance =       ', round(mu.B$abundance,0), "\n")
cat('Standard Error =  ', round(mu.B$precision$se,0), "\n",
    '95 Percent Conf. Int. =   (', round(mu.B$precision$lcl.95,0), ',',
           round(mu.B$precision$ucl.95,0), ')', '\n')

mu.B <- integrate.intensity(sim.dspat,dimyx=100,se=TRUE,od=TRUE,reps=50)
cat('Abundance =       ', round(mu.B$abundance,0), "\n")
cat('Standard Error (corrected) =  ', round(mu.B$precision.od$se,0), "\n",
    '95 Percent Conf. Int.(corrected)  =   (', round(mu.B$precision.od$lcl.95,0), ',',
           round(mu.B$precision.od$ucl.95,0), ')', '\n')
plot(mu.B$lambda, main='Estimated Intensity')
plot(sim.dspat$lines.psp,lty=2,add=TRUE)
plot(owin(poly=sim.dspat$transect),add=TRUE)
plot(sim.dspat$model$Q$data,add=TRUE)
# Now sample with same point process realization with a different sampling angle
dev.new()
simdata=simDSpat(study.area,covariates,int.formula=~factor(habitat),
                   int.par=c(0,1,2),angle=90,nlines=10,width=3,pts=simdata$pts)
sim.dspat=dspat(int.formula=~factor(habitat),study.area=study.area,
               obs=simdata$observations,lines=simdata$lines,
               covariates=covariates,epsvu=c(1,.05))
mu.B <- integrate.intensity(sim.dspat,dimyx=100,se=TRUE)
cat('Abundance =       ', round(mu.B$abundance,0), "\n")
cat('Standard Error =  ', round(mu.B$precision$se,0), "\n",
    '95 Percent Conf. Int. =   (', round(mu.B$precision$lcl.95,0), ',',
           round(mu.B$precision$ucl.95,0), ')', '\n')
mu.B <- integrate.intensity(sim.dspat,dimyx=100,se=TRUE,od=TRUE,reps=50)
cat('Abundance  =       ', round(mu.B$abundance,0), "\n")
cat('Standard Error (corrected)=  ', round(mu.B$precision.od$se,0), "\n",
    '95 Percent Conf. Int. (corrected)=   (', round(mu.B$precision.od$lcl.95,0), ',',
           round(mu.B$precision.od$ucl.95,0), ')', '\n')
plot(mu.B$lambda, main='Estimated Intensity')
plot(sim.dspat$lines.psp,lty=2,add=TRUE)
plot(owin(poly=sim.dspat$transect),add=TRUE)
spatstat.options(checkpolygons=TRUE)
plot(sim.dspat$model$Q$data,add=TRUE)
# Sample with detection as a function of habitat
dev.new()
study.area=owin(poly=list(x=c(0,40,40,100,100,0),y=c(0,0,40,40,100,100)))
simdata=simDSpat(study.area,covariates,int.formula=~factor(habitat),
                   int.par=c(0,1,2),angle=45,nlines=10,width=3,
                   det.par=c(.1,.5,-.2),det.formula=~factor(habitat))
sim.dspat=dspat(int.formula=~factor(habitat),det.formula=~factor(habitat),
                study.area=study.area, obs=simdata$observations,lines=simdata$lines,
                covariates=covariates,epsvu=c(1,.05))
summary(sim.dspat)
AIC(sim.dspat)
coef(sim.dspat)
mu.B <- integrate.intensity(sim.dspat,dimyx=100,se=TRUE)
cat('Abundance =       ', round(mu.B$abundance,0), "\n")
cat('Standard Error =  ', round(mu.B$precision$se,0), "\n",
    '95 Percent Conf. Int. =   (', round(mu.B$precision$lcl.95,0), ',',
           round(mu.B$precision$ucl.95,0), ')', '\n')
plot(mu.B$lambda, main='Estimated Intensity')
plot(sim.dspat$lines.psp,lty=2,add=TRUE)
plot(owin(poly=sim.dspat$transect),add=TRUE)
plot(sim.dspat$model$Q$data,add=TRUE)
############################################################
# Generate example like Figure used in paper for simulations
# Note: it required a patch to plot.im from spatstat to
# fix the ribbon bar on the side.
#
#
#             if(is.null(list(...)$zlim))
#            {
#               ribbonvalues <- seq(vrange[1], vrange[2], length = ribn)
#               ribbonrange <- vrange
#               ribbonticks <- clamp(pretty(ribbonvalues), vrange)
#            }
#            else
#            {
#               zlim=list(...)$zlim
#               ribbonvalues <- seq(zlim[1], zlim[2], length = ribn)
#               ribbonrange <- zlim
#               ribbonticks <- clamp(pretty(ribbonvalues), zlim)
#            }
#
#############################################################
study.area=owin(poly=list(x=c(0,100,100,0),y=c(0,0,100,100)))
covariates = simCovariates(hab.range=30, probs=c(1/3,2/3))
postscript("Figure1.eps",horizontal=FALSE)
par(mfrow=c(2,1),mar=c(3, 1, 3, 1) + 0.1)
k=10
width=0.04*100/k
En=75
p=0.25
EN=En/(.04*p)
simdata=simDSpat(study.area,covariates,int.formula=~factor(habitat)+river,EN=EN,
                   int.par=c(0,1,2,-1),angle=90,nlines=k,width=width,
                   det.par=log(width/5),showplot=TRUE,col=gray(1-c(1:100)/120),
                   breaks=(0:100)*2.5/100,
                   zlim=c(0,2.5))
                   
lines(c(50,50),c(0,100),lty=2)
sim.dspat=dspat(int.formula=~factor(habitat)+river,study.area=study.area,
               obs=simdata$observations,lines=simdata$lines,
               covariates=covariates,epsvu=c(1,width/100))
summary(sim.dspat)
AIC(sim.dspat)
coef(sim.dspat)
mu.B <- integrate.intensity(sim.dspat,dimyx=100)
plot(mu.B$lambda, col=gray(1-c(1:100)/120), main='Estimated Intensity',
  breaks=(0:100)*2.5/100,zlim=c(0,2.5))
plot(sim.dspat$lines.psp,lty=2,add=TRUE)
plot(owin(poly=sim.dspat$transect),add=TRUE)
plot(sim.dspat$model$Q$data,add=TRUE)
dev.off()
spatstat.options(checkpolygons=TRUE)
}