unmarked (version 0.9-1)

parboot: Parametric bootstrap method for fitted models inheriting class.

Description

Simulate datasets from a fitted model, refit the model, and generate a sampling distribution for a user-specified fit-statistic.

Arguments

object

a fitted model inheriting class "unmarkedFit"

statistic

a function returning a vector of fit-statistics. First argument must be the fitted model. Default is sum of squared residuals.

nsim

number of bootstrap replicates

report

print fit statistic every 'report' iterations during resampling

...

Additional arguments to be passed to statistic

Value

• An object of class parboot with three slots:
• call
  parboot call
• t0
  Numeric vector of statistics for original fitted model.
• t.star
  nsim by length(t0) matrix of statistics for each simulation fit.

Details

This function simulates datasets based upon a fitted model, refits the model, and evaluates a user-specified fit-statistic for each simulation. Comparing this sampling distribution to the observed statistic provides a means of evaluating goodness-of-fit or assessing uncertainty in a quantity of interest.

Examples

data(linetran)
(dbreaksLine <- c(0, 5, 10, 15, 20)) 
lengths <- linetran$Length

ltUMF <- with(linetran, {
	unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4), 
	siteCovs = data.frame(Length, area, habitat), dist.breaks = dbreaksLine,
	tlength = lengths*1000, survey = "line", unitsIn = "m")
    })

# Fit a model
(fm <- distsamp(~area ~habitat, ltUMF))

# Function returning three fit-statistics.
fitstats <- function(fm) {
    observed <- getY(fm@data)
    expected <- fitted(fm)
    resids <- residuals(fm)
    sse <- sum(resids^2)
    chisq <- sum((observed - expected)^2 / expected)
    freeTuke <- sum((sqrt(observed) - sqrt(expected))^2)
    out <- c(SSE=sse, Chisq=chisq, freemanTukey=freeTuke)
    return(out)
    }

(pb <- parboot(fm, fitstats, nsim=25, report=1))
plot(pb, main="")


# Local population size (derived parameter)
Nhat <- function(fm) {
    lengths <- fm@data@tlength
    strip.width <- max(fm@[e-mail])
    plot.area.ha <- lengths * (2*strip.width) / 1e4
    N <- sum(predict(fm, type="state")$Predicted*plot.area.ha, na.rm=TRUE)
    }
    
(pb.N <- parboot(fm, Nhat, nsim=25, report=5))
plot(pb.N)