var.Np<-function(pars,mix.terms,width,z,zdim,pt,n,hessian,N,pa.vec,pa,x){
# calculate the variance of N and the average p
# most of this is taken from the summary.ds method in mrds
# by Jeff Laake, http://github.com/jlaake/mrds
ret<-list()
# function to calculate N, to be passed to DeltaMethod()
Nfct<-function(pars,mix.terms,width,z,zdim,pt,n){
mu<-mu.calc(pars,mix.terms,width,z,zdim,pt)
if(length(mu)!=n){
mu<-rep(mu,n)
}
if(!pt){
return(width*sum(1/mu))
}else{
return((pi*width^2)*sum(1/mu))
}
}
# quick function to be passed to DeltaMethod()
# this actually just returns n, this is a bit pointless here, but is
# retained for mrds comparability
pafct<-function(pars,mix.terms,width,z,zdim,pt,n){
mu<-mu.calc(pars,mix.terms,width,z,zdim,pt)
if(length(mu)!=n){
mu<-rep(mu,n)
}
return(sum(mu/mu))
}
######################################################
# calculate the variance-covariance matrix
# Compute varmat using first partial approach (pg 62 of Buckland et al 2002)
# set machine epsilon for numeric differentiation
eps<-(.Machine$double.eps)^(1/3)
#eps<-.0001 # from Jeff
fpar<-pars
fpar1<-pars
parmat<-NULL
# Compute first partial (numerically) of log(f(y)) for each observation
# for each parameter and store in parmat (n by length(fpar))
for (i in 1:length(fpar)){
deltap <- eps*fpar1[i]
if(deltap==0){
deltap <- eps
}
fpar[i] <- fpar1[i]- deltap
# x1<- flt(fpar,x,width,mix.terms,0,"hn",z,zdim,pt)
x1 <- log(rowSums(eval.pdf(fpar,x,width,mix.terms,showit=0,"hn",z,zdim,pt)))
fpar[i] <- fpar1[i]+deltap
#x2 <- flt(fpar,x,width, mix.terms,0,"hn",z,zdim,pt)
x2 <- log(rowSums(eval.pdf(fpar,x,width,mix.terms,showit=0,"hn",z,zdim,pt)))
parmat=cbind(parmat,(x2-x1)/(2*deltap))
}
# Compute varmat using first partial approach (pg 62 of Buckland et al 2002)
varmat <- matrix(0,ncol=length(fpar1),nrow=length(fpar1))
for(i in 1:length(fpar1)){
for(j in 1:length(fpar1)){
varmat[i,j] <- sum(parmat[,i]*parmat[,j])
}
}
# invert
vcov<-solvecov(varmat)$inv
######################################################
### calculate the variance, se and CV of N in the covered region
if(length(pa.vec)==1){
pa.vec<-rep(pa.vec,n)
}
Nhatvar.list<-DeltaMethod(pars,Nfct,vcov,eps,mix.terms=mix.terms,
width=width,z=z,zdim=zdim,pt=pt,n=n)
Nhatvar<-Nhatvar.list$variance + sum((1-pa.vec)/pa.vec^2)
# put them in the results list
ret$N.var<-Nhatvar
ret$N.se<-sqrt(Nhatvar)
ret$N.CV<-sqrt(Nhatvar)/N
cvN<-sqrt(Nhatvar)/N
### calculate variance, se and CV for average detectability
vc1.list<-DeltaMethod(pars,pafct,vcov,eps,mix.terms=mix.terms,
width=width,z=z,zdim=zdim,pt=pt,n=n)
vc1<-vc1.list$variance
vc2<-sum(1 - pa.vec)
covar<-sum((1 - pa.vec)/pa.vec)
var.pbar.list<-list(var=vc1+vc2,partial=vc1.list$partial,covar=covar)
covar<-t(Nhatvar.list$partial)%*%vcov%*%var.pbar.list$partial+
var.pbar.list$covar
var.pbar<-pa^2*(cvN^2 + var.pbar.list$var/n^2-2*covar/(n*N))
# put the results in the return list...
ret$pa.var<-var.pbar
ret$pa.se<-sqrt(var.pbar)
ret$pa.CV<-sqrt(var.pbar)/pa
ret$vcov<-vcov
return(ret)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.