Nothing
#' Likelihood function for Cormack-Jolly-Seber model
#'
#' For a given set of parameters and data, it computes -log Likelihood value.
#'
#' This function uses a FORTRAN subroutine (cjs.f) to speed up computation of the likelihood but the
#' result can also be obtained wholly in R with a small loss in precision. See R code below.
#' The R and FORTRAN code uses the likelihood formulation of Pledger et al.(2003).
#' \preformatted{
#' get.p=function(beta,dm,nocc,Fplus)
#' {
#' # compute p matrix from parameters (beta) and list of design matrices (dm)
#' # created by function create.dm
#' ps=cbind(rep(1,nrow(dm)/(nocc-1)),
#' matrix(dm%*%beta,ncol=nocc-1,nrow=nrow(dm)/(nocc-1),byrow=TRUE))
#' ps[Fplus==1]=plogis(ps[Fplus==1])
#' return(ps)
#' }
#' get.Phi=function(beta,dm,nocc,Fplus)
#' {
#' # compute Phi matrix from parameters (beta) and list of design matrices (dm)
#' # created by function create.dm
#' Phis=cbind(rep(1,nrow(dm)/(nocc-1)),
#' matrix(dm%*%beta,ncol=nocc-1,nrow=nrow(dm)/(nocc-1),byrow=TRUE))
#' Phis[Fplus==1]=plogis(Phis[Fplus==1])
#' return(Phis)
#' }
#'#################################################################################
#'# cjs.lnl - computes likelihood for CJS using Pledger et al (2003)
#'# formulation for the likelihood. This code does not cope with fixed parameters or
#'# loss on capture but could be modified to do so. Also, to work directly with cjs.r and
#'# cjs.accumulate call to process.ch would have to have all=TRUE to get Fplus and L.
#'# Arguments:
#'# par - vector of beta parameters
#'# imat - list of freq, indicator vector and matrices for ch data created by process.ch
#'# Phi.dm - list of design matrices; a dm for each capture history
#'# p.dm - list of design matrices; a dm for each capture history
#'# debug - if TRUE show iterations with par and -2lnl
#'# time.intervals - intervals of time between occasions
#'# Value: -LnL using
#'#################################################################################
#'cjs.lnl=function(par,model_data,Phi.links=NULL,p.links=NULL,debug=FALSE,all=FALSE) {
#' if(debug)cat("\npar = ",par)
#' #extract Phi and p parameters from par vector
#' nphi=ncol(model_data$Phi.dm)
#' np=ncol(model_data$p.dm)
#' beta.phi=par[1:nphi]
#' beta.p=par[(nphi+1):(nphi+np)]
#' #construct parameter matrices (1 row for each capture history and a column
#' #for each occasion)
#' Phis=get.Phi(beta.phi,model_data$Phi.dm,nocc=ncol(model_data$imat$chmat),
#' model_data$imat$Fplus)
#' if(!is.null(model_data$time.intervals))
#' {
#' exponent=cbind(rep(1,nrow(Phis)),model_data$time.intervals)
#' Phis=Phis^exponent
#' }
#' ps=get.p(beta.p,model_data$p.dm,nocc=ncol(model_data$imat$chmat),
#' model_data$imat$Fplus)
#' if(debug)cat("\npar = ",par)
#' # Compute probability of dying in interval from Phis
#' M=cbind((1-Phis)[,-1],rep(1,nrow(Phis)))
#' # compute cummulative survival from release across each subsequent time
#' # and the cummulative probability for detection (capture) across each time
#' Phi.cumprod=1-model_data$imat$Fplus + Phis*model_data$imat$Fplus
#' cump=(1-model_data$imat$Fplus)+model_data$imat$Fplus*
#' (model_data$imat$chmat*ps+(1-model_data$imat$chmat)*(1-ps))
#' for (i in 2:ncol(cump))
#' {
#' Phi.cumprod[,i]=Phi.cumprod[,i-1]*Phi.cumprod[,i]
#' cump[,i]=cump[,i-1]*cump[,i]
#' }
#' # compute prob of capture-history
#' pch=rowSums(model_data$imat$L*M*Phi.cumprod*cump)
#' lnl=-sum(model_data$imat$freq*log(pch))
#' if(debug)cat("\n-2lnl = ",2*lnl)
#' return(lnl)
#'}
#'}
#'
#' @param par vector of parameter values
#' @param model_data a list that contains: 1)imat-list of vectors and matrices constructed by
#' \code{\link{process.ch}} from the capture history data, 2)Phi.dm design matrix for Phi constructed by \code{\link{create.dm}},
#' 3)p.dm design matrix for p constructed by \code{\link{create.dm}},
#' 4)Phi.fixed matrix with 3 columns: ch number(i), occasion number(j),
#' fixed value(f) to fix phi(i,j)=f, 5)p.fixed matrix with 3 columns: ch number(i), occasion number(j), and
#' 6) time.intervals intervals of time between occasions if not all 1
#' fixed value(f) to fix p(i,j)=f
#' @param Phi.links vector of links for each parameter
#' @param p.links vector of links for each parameter
#' @param debug if TRUE will printout values of \code{par} and function value
#' @param all if TRUE, returns entire list rather than just lnl; can be used to
#' extract reals
#' @param cjsenv environment for cjs to update iteration counter
#' @return either -log likelihood value if \code{all=FALSE} or the entire
#' list contents of the call to the FORTRAN subroutine if \code{all=TRUE}. The
#' latter is used from \code{\link{cjs_admb}} after optimization to extract the real
#' parameter estimates at the final beta values.
#' @author Jeff Laake
#' @references Pledger, S., K. H. Pollock, et al. (2003). Open
#' capture-recapture models with heterogeneity: I. Cormack-Jolly-Seber model.
#' Biometrics 59(4):786-794.
#'
# cjs.lnl=function(par,model_data,Phi.links=NULL,p.links=NULL,debug=FALSE,all=FALSE,cjsenv,nc,cl)
cjs.lnl=function(par,model_data,Phi.links=NULL,p.links=NULL,debug=FALSE,all=FALSE,cjsenv)
{
f_eval=get("markedfunc_eval",envir=cjsenv)+1
assign("markedfunc_eval", f_eval, envir = cjsenv)
if(debug)cat("par = ",par,"\n")
nocc=model_data$imat$nocc
nphi=ncol(model_data$Phi.dm)
np=ncol(model_data$p.dm)
# allow for possibility that all parameter values are fixed for Phi or p
if(nphi!=0)
{
beta.phi=par[1:nphi]
Phibeta=as.vector(model_data$Phi.dm%*%beta.phi)
}
else
Phibeta=0
if(np!=0){
beta.p=par[(nphi+1):(nphi+np)]
pbeta=as.vector(model_data$p.dm%*%beta.p)
}
else
pbeta=0
Phibeta=matrix(Phibeta,ncol=nocc-1,nrow=nrow(model_data$Phi.dm)/(nocc-1),byrow=TRUE)
pbeta=matrix(pbeta,ncol=nocc-1,nrow=nrow(model_data$p.dm)/(nocc-1),byrow=TRUE)
lnl=.Fortran("cjs",as.double(model_data$imat$chmat),as.double(Phibeta),as.double(pbeta),
as.double(model_data$imat$first),as.double(model_data$imat$last),as.double(model_data$imat$freq),
as.integer(model_data$imat$loc),as.double(model_data$Phi.fixed),as.double(model_data$p.fixed),
as.double(model_data$time.intervals),as.integer(nrow(model_data$imat$chmat)),
as.integer(ncol(model_data$imat$chmat)),as.integer(nrow(model_data$Phi.fixed)),
as.integer(nrow(model_data$p.fixed)),lnl=double(1),p0=double(nrow(model_data$imat$chmat)),PACKAGE="marked")
if(debug)
{
cat("-2lnl = ",2*lnl$lnl,"\n")
} else
if((f_eval-100*floor(f_eval/100))==0)
{
cat("\r Number of evaluations: ",f_eval," -2lnl:",formatC(2*lnl$lnl,digits=10))
flush.console()
}
if(all)
return(lnl)
else
return(lnl$lnl)
}
# attempt to parallelize the code
#cjs.parallel=function(i,chmat,Phibeta,pbeta,first,last,freq,loc,Phi.fixed,p.fixed,time.intervals,vec)
#{
#ir=vec[i,1]:vec[i,2]
#lnl=.Fortran("cjs",as.double(chmat[ir,]),as.double(Phibeta[ir,]),as.double(pbeta[ir,]),
# as.double(first[ir]),as.double(last[ir]),as.double(freq[ir]),
# as.integer(loc[ir]),as.double(Phi.fixed),as.double(p.fixed),
# as.double(time.intervals[ir,]),as.integer(nrow(chmat[ir,])),
# as.integer(ncol(chmat[ir,])),as.integer(nrow(Phi.fixed)),
# as.integer(nrow(p.fixed)),lnl=double(1),p0=double(nrow(chmat[ir,])),PACKAGE="marked")
#return(lnl)
#}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.