Nothing
F.cjs.estim <- function(capture, survival, histories, cap.init, sur.init,
group, nhat.v.meth=1,
c.hat=-1.0, df=NA, intervals=rep(1,ncol(histories)-1),
conf=0.95, link="logit",
control=mra.control() ){
start.sec <- proc.time()[3]
run.date <- Sys.time()
# ---- Verify some of the inputs
if( missing(histories) ){
stop("Capture histories matrix must be specified")
}
if( missing(capture) ){
stop("Capture covariates must be specified")
}
if( missing(survival) ){
stop("Survival covariates must be specified")
}
if( length(union( unique(histories), c(0,1,2))) > 3 ) stop("Capture histories must consist of 0's, 1's, and 2's only.")
if( !is.numeric(nhat.v.meth) || !(nhat.v.meth %in% c(1,2,3)) ){
stop("value of 'nhat.v.meth' must be 1, 2, or 3")
}
# Validity of cov.meth and algorithm and other inputs is determined in mra.control()
# ---- Initialize some variables
hist.name <- deparse(substitute(histories))
cr.call <- match.call()
nan <- nrow( histories )
ns <- ncol( histories )
if( length(intervals) < (ns-1)){
stop(paste("Too few time intervals specified. INTERVALS vector should have length", ns-1))
} else if(length(intervals) >= (ns-1)){
intervals <- intervals[1:(ns-1)]
intervals <- c(intervals, 0) # Make this vector length ns, but we never use the last element in estimation.
}
SVD_ZERO = 0.5E-6 # This is the value that determines when a singular value is zero.
# I.e., when a singular value is less than this, it is considered
# zero. This is key when computing rank of the
# variance-covariance matrix. Mark uses 0.5e-6.
# ---- Get the X and Y matricies. After this, the x and Y matricies are in huge NAN by (ns*nx) matricies.
covars <- F.cr.model.matrix( capture, survival, nan, ns )
nx <- covars$n.cap.covars #nx and ny include intercept. This is total number of parameters
ny <- covars$n.sur.covars
# Set initial values if missing or short
if( missing(cap.init) ){
cap.init <- rep(0,nx)
} else if(length(cap.init) < (nx) ){
cap.init <- c(cap.init, rep(0, nx-length(cap.init)))
}
if( missing(sur.init) ){
sur.init <- rep(0,ny)
} else if(length(sur.init) < (ny) ){
sur.init <- c(sur.init, rep(0, ny-length(sur.init)))
}
# Set up the tolerance vector, if not specified, or if not long enough
if( length(control$tol) < (nx+ny) ){
control$tol <- rep(control$tol, trunc((nx+ny) / length(control$tol))+1)[1:(nx+ny)]
} else if( length(control$tol > (nx+ny)) ){
control$tol <- control$tol[1:(nx+ny)]
}
# Fix up the group variable
if( missing( group )){
group <- rep(1, nan)
ng <- 1
} else {
ng <- length(unique(group))
}
# Transfer over c.hat
vif <- c.hat
# Do the estimation, but first allocate room for answers
loglik <- chisq.vif <- df.vif <- 0
parameters <- se.param <- rep(0, nx + ny )
covariance <- matrix( 0, nx+ny, nx+ny )
p.hat <- se.p.hat <- s.hat <- se.s.hat <- matrix( 0, nan, ns )
exit.code <- cov.code <- 0
n.hat <- se.n.hat <- rep(0, ns)
# on entry to .Fortran, maxfn is maximum number of function evals. On return, maxfn is actual number of evaluations.
maxfn <- control$maxfn
# Re-code the link specification to integers
if( link=="logit" ){
link.code <- 1
} else if( link == "sine" ){
link.code <- 2
} else if( link == "hazard" ){
link.code <- 3
} else {
stop("Unknown link function specified.")
}
if(control$trace) cat( "Calling MRA DLL to maximize likelihood. Please wait...\n")
ans <- .Fortran( "cjsmod",
nan = as.integer(nan),
ns = as.integer(ns),
nx = as.integer(nx),
ny = as.integer(ny),
ng = as.integer(ng),
histories = as.integer(histories),
group = as.integer(group),
algorithm = as.integer(control$algorithm),
cov.meth = as.integer(control$cov.meth),
link = as.integer(link.code),
nhat.v.meth = as.integer(nhat.v.meth),
capX = as.double(covars$capX),
survX = as.double(covars$survX),
cap.init = as.double(cap.init),
sur.init = as.double(sur.init),
maxfn = as.integer(maxfn),
beta.tol.vec= as.double(control$tol),
loglik = as.double(loglik),
vif = as.double(vif),
chisq.vif = as.double(chisq.vif),
df.vif = as.double(df.vif),
parameters = as.double(parameters),
se.param = as.double(se.param),
covariance = as.double(covariance),
p.hat = as.double(p.hat),
se.p.hat = as.double(se.p.hat),
s.hat = as.double(s.hat),
se.s.hat = as.double(se.s.hat),
n.hat = as.double(n.hat),
se.n.hat = as.double(se.n.hat),
exit.code = as.integer(exit.code),
cov.code = as.integer(cov.code),
intervals = as.double(intervals),
PACKAGE="mra"
)
if( control$trace ) cat(paste("Returned from MRAWIN. Fitting details written to MRA.LOG.\n", sep=""))
## ----- Reset missing standard errors to NA
ans$se.param[ ans$se.param < 0 ] <- NA
# ----- R does not preserve the matrix structures in .Fortran call. Put matricies,
# which are now vectors, back to matricies.
covariance <- matrix( ans$covariance, nrow=nx+ny )
ans$p.hat <- matrix( ans$p.hat, nrow=nan )
ans$se.p.hat <- matrix( ans$se.p.hat, nrow=nan )
ans$s.hat <- matrix( ans$s.hat, nrow=nan )
ans$se.s.hat <- matrix( ans$se.s.hat, nrow=nan )
# ----- Fix up number of parameters.
if( is.na(df) ){
# Singular values of hessian are more stable than of covariance
hess <- tryCatch(solve(covariance),
error=function(x){1},
warning=function(x){2})
if(is.matrix(hess)){
# likely everything ok with var-covar, but could have NA in hess
svals <- tryCatch(svd(hess)$d,
error=function(x){NA},
warning=function(x){NA})
if(!all(is.na(svals))){
svals = svals / max(svals, na.rm=TRUE)
df <- sum(svals >= SVD_ZERO, na.rm=TRUE)
} else {
df <- NA
ans$cov.code <- 1
}
} else {
df <- NA
ans$cov.code <- 1
}
} else if( df <= 0 ){
df <- nx+ny # assume full rank
} # else use the values supplied by user (unchanged from input)
# ----- Work out exit codes
if( ans$exit.code==0 ){
exit.mess = "FAILURE: Initial Hessian not positive definite"
} else if( ans$exit.code == 1 ){
exit.mess = "SUCCESS: Convergence criterion met"
} else if( ans$exit.code == 2 ){
exit.mess = "FAILURE: G'dX > 0, rounding error"
} else if( ans$exit.code == 3 ){
exit.mess = "FAILURE: Likelihood evaluated too many times"
} else if( ans$exit.code == -1 ){
exit.mess = "FAILURE: Algorithm 2 not implimented."
} else {
exit.mess = "Unknown exit code"
}
if(control$algorithm == 1){
alg.mess <- "Optimization by line search."
} else {
alg.mess <- "Unknown optimization routine."
}
if(control$cov.meth == 1){
cov.mess = "Covariance from numeric derivatives."
} else if (control$cov.meth == 2){
cov.mess = "Covariance from optimization Hessian."
} else {
cov.mess = "Unkown covariance method."
}
if(ans$cov.code == 0){
cov.mess = paste( cov.mess, "SUCCESS: Non-singular covariance matrix.")
} else if( ans$cov.code == 1) {
cov.mess = paste( cov.mess, "ERROR: COVARIANCE MATRIX IS SINGULAR.")
} else {
cov.mess = paste( cov.mess, "ERROR: NON-NEGATIVE DEFINITE COVARIANCE MATRIX.")
}
# do some printouts
if( ans$exit.code != 1 ){
cat(exit.mess)
cat("\n")
} else if( ans$cov.code != 0 ){
cat(cov.mess)
cat("\n")
} else if( control$trace ){
cat(paste( alg.mess, "\n", sep="" ))
cat(paste( exit.mess, "\n", sep="" ))
cat(paste( cov.mess, "\n", sep="" ))
}
# ----- Fix up capture and survival estimates
# Wipe out the first capture probability and the last survival probability. These are computed
# by mrawin5 if there are covariates out there, but they are not used in the likelihood, and
# we can't really estimate them.
ans$p.hat[,1] <- NA
ans$s.hat[,ns] <- NA
ans$se.s.hat[,ns] <- NA
ans$se.p.hat[,1] <- NA
capcoef <- ans$parameters[1:nx]
se.capcoef <- ans$se.param[1:nx]
surcoef <- ans$parameters[(nx+1):(nx+ny)]
se.surcoef <- ans$se.param[(nx+1):(nx+ny)]
names(capcoef) <- covars$cap.vars
names(se.capcoef) <- covars$cap.vars
names(surcoef) <- covars$sur.vars
names(se.surcoef) <- covars$sur.vars
dimnames(covariance) <- list( c( paste("cap.",names( capcoef ),sep=""), paste("sur.",names( surcoef ), sep="")),
c( paste("cap.",names( capcoef ),sep=""), paste("sur.",names( surcoef ), sep="")))
# ----- Now that df is computed, recompute fit statistics
dev <- -2*ans$loglik
aic <- -2*ans$loglik + 2*df
qaic <- -2*(ans$loglik/ans$vif) + 2*df
aicc <- aic + (2*df*(df+1))/(nan - df - 1)
qaicc <- qaic + (2*df*(df+1))/(nan - df - 1)
# ----- LEAVE THIS CODE HERE FOR FUTURE REFERENCE
# ----- Compute EBC of Peterson (1986), Stat & Prob letters, p227
#qf <- t(parameters) %*% covariance %*% parameters
#q1 <- qf/((df)*nan)
#q2 <- qf/(df)
#ebc <- -2*loglik + (df)*log(nan) +
# (df)*log( max(c( 1/nan, q1 ))) +
# (df)*min( c( 1, q2 ) )
# ----- Put all the "auxillary" info together
aux <- list( call=cr.call,
nan=nan,
ns=ns,
nx=nx,
ny=ny,
link=link,
cov.name=c(names(capcoef), names(surcoef)),
ic.name=hist.name,
mra.version=packageDescription("mra")$Version,
R.version=R.version.string,
run.date=run.date )
# ----- Fix up the estimates of N.
names(ans$n.hat) <- dimnames(histories)[[2]]
num.observed <- c( t( rep(1, nrow(histories)) ) %*% (histories >= 1) )
crit.value <- qnorm( 1-((1-conf)/2) )
lower.ci <- ans$n.hat - crit.value * ans$se.n.hat
lower.ci <- ifelse(lower.ci < num.observed, num.observed, lower.ci)
upper.ci <- ans$n.hat + crit.value * ans$se.n.hat
ans$n.hat[1] <- NA # can't estimate the first N
ans$se.n.hat[1] <- NA
lower.ci[1] <- NA
upper.ci[1] <- NA
# ----- Compute execution time here so can store in output.
ex.time <- (proc.time()[3] - start.sec) / 60
# ----- Done. Put into a 'CR' object.
ans <- list( histories=histories,
aux=aux,
intervals=intervals[1:(ns-1)],
loglike=ans$loglik,
deviance=dev,
aic=aic,
qaic=qaic,
aicc=aicc,
qaicc=qaicc,
vif=ans$vif,
chisq.vif=ans$chisq.vif,
vif.df=ans$df.vif,
parameters=ans$parameters,
se.param=ans$se.param,
capcoef=capcoef,
se.capcoef=se.capcoef,
surcoef=surcoef,
se.surcoef=se.surcoef,
covariance=covariance,
p.hat=ans$p.hat,
se.p.hat=ans$se.p.hat,
s.hat=ans$s.hat,
se.s.hat=ans$se.s.hat,
df=df,
df.estimated=ans$df.estimated,
control=control,
message=c(alg.mess,exit.mess,cov.mess),
exit.code=ans$exit.code,
cov.code=ans$cov.code,
fn.evals=ans$maxfn,
ex.time = ex.time,
n.hat =ans$n.hat,
se.n.hat=ans$se.n.hat,
n.hat.lower = lower.ci,
n.hat.upper = ceiling(upper.ci),
n.hat.conf = conf,
nhat.v.meth = ans$nhat.v.meth,
num.caught=num.observed)
class(ans) <- c("cjs", "cr")
# Compute fitted and residual components, if converged
if( ans$exit.code == 1 ){
ans$fitted <- predict( ans )
ans$residuals <- residuals( ans, type="pearson" )
ans$resid.type <- "pearson"
} else {
ans$fitted <- NA
ans$residuals <- NA
ans$resid.type <- NA
}
if( control$trace ){
if( ex.time < 0.01666667 ){
cat("\t(Execution time < 1 second)\n")
} else {
cat(paste("\t(Execution time =", round(ex.time,2), "minutes)\n"))
}
}
ans
}
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.