Nothing
bhnoneq <- function(year = NULL, mlen = NULL, ss = NULL, K = NULL, Linf = NULL,
Lc = NULL, nbreaks = NULL, styrs = NULL, stZ = NULL, graph = TRUE) {
if (is.null(mlen)) stop("mean length vector does not exist")
if (is.null(year)) stop("year vector does not exist")
if (is.null(ss)) stop("numbers vector does not exist")
if (!is.numeric(mlen)) stop("vector is not numeric")
if (is.null(stZ)) stop("Initial Z vector does not exist")
if (is.null(K)) stop("K not specified")
if (is.null(Linf)) stop("Linf not specified")
if (is.null(Lc)) stop("Lc not specified")
if (length(mlen) != length(year)) stop("vectors have different lengths")
if (is.null(nbreaks)) stop("Number of mortality breaks not specified")
if (nbreaks > 0 & is.null(styrs)) stop("Starting guesses for years of mortality breaks not specified")
if (nbreaks > 0 & length(styrs) != nbreaks) stop("Number of starting guesses for years of mortality breaks not equal to number of mortality breaks")
if (nbreaks > 0 & length(stZ) != nbreaks+1) stop("Number of starting guesses for Z not equal to number of mortality breaks plus one")
dat <- data.frame(year = year, mlen = mlen, ss = ss)
sigma<-NULL; Lpred<-NULL;results<-NULL
return_LL <- function(x, dat, nbreaks, K, Linf, Lc) {
year <- dat$year
Lbar <- dat$mlen
ss <- dat$ss
Z <- x[1:(nbreaks+1)]
ggyr <- x[(nbreaks + 2):length(x)]
nyr <- sum(!is.na(Lbar))
count <- length(year)
dy <- array(0, dim = c(nbreaks, count))
for(i in 1:nbreaks) {
for(j in 1:count) {
dy[i, j] <- ifelse(ggyr[i] >= j, 0, j - ggyr[i])
}
}
#Quang correction 1/12/2016
if(nbreaks > 1) {
for(i in 1:(nbreaks - 1)) {
for(j in 1:count) dy[i, j] <- dy[i, j] - dy[i + 1, j]
}
}
a <- array(0, dim = c(nbreaks + 1, count))
s <- array(0, dim = c(nbreaks + 1, count))
r <- array(0, dim = c(nbreaks + 1, count))
denom <- rep(0, count)
numersum <- rep(0, count)
numerator <- rep(0, count)
for(m in 1:count) {
for(i in 1:(nbreaks + 1)) {
a[i, m] <- 1
r[i, m] <- 1
if (i < (nbreaks + 1)) s[i, m] <- 1 - exp(-(Z[nbreaks + 2 - i] + K) * dy[nbreaks + 1 - i, m])
if (i == (nbreaks + 1)) s[i, m] <- 1
for (j in 1:(i - 1)) {
if (i > 1) {
a[i, m] <- a[i, m] * exp(-Z[nbreaks + 2 - j] * dy[nbreaks + 1 - j, m])
r[i, m] <- r[i, m] * exp(-(Z[nbreaks + 2 - j] + K) * dy[nbreaks + 1 - j, m])
}
}
if (i <= nbreaks) {
denom[m] <- denom[m] + a[i, m] *
((1 - exp(-Z[nbreaks + 2 - i] * dy[nbreaks + 1 - i, m]))/Z[nbreaks + 2 - i])
}
if (i == (nbreaks + 1)) {
denom[m] <- denom[m] + a[i, m]/Z[nbreaks + 2 - i]
}
numersum[m] <- numersum[m] + r[i, m] * s[i, m]/(Z[nbreaks + 2 - i] + K)
}
}
numerator <- Linf * (denom - (1 - Lc/Linf) * numersum)
Lpred <<- numerator/denom
sigma <<- sqrt(sum(ss * (Lbar - Lpred)^2, na.rm = TRUE)/nyr)
LL <- -nyr * log(sigma) - sum(ss * (Lbar - Lpred)^2, na.rm = TRUE)/(2 * sigma^2)
return(-1 * LL)
}
return_LL_eq <- function(x, dat, K, Linf, Lc) {
year <- dat$year
Lbar <- dat$mlen
ss <- dat$ss
Z <- x
nyr <- sum(!is.na(Lbar))
count <- length(year)
Lpred <<- rep(Linf * (1 - (Z/(Z+K))*(1 - Lc/Linf)), count)
sigma <<- sqrt(sum(ss * (Lbar - Lpred)^2, na.rm = TRUE)/nyr)
LL <- -nyr * log(sigma) - sum(ss * (Lbar - Lpred)^2, na.rm = TRUE)/(2 * sigma^2)
return(-1 * LL)
}
if(nbreaks == 0) {
results <- optim(stZ, return_LL_eq, method = "BFGS", dat = dat, K = K, Linf = Linf, Lc = Lc,
control = list(maxit = 1e+06, abstol = 1e-07), hessian = TRUE)
}
else {
g_year <- styrs - year[1] + 1
stpar <- c(stZ, g_year)
results <- optim(stpar, return_LL, method = "BFGS", dat = dat, nbreaks = nbreaks, K = K, Linf = Linf, Lc = Lc,
control = list(maxit = 1e+06, abstol = 1e-07), hessian = TRUE)
}
npar <- length(results$par) + 1
AIC <- 2 * results$value + 2 * npar
var <- solve(results$hessian)
SE <- sqrt(diag(var))
if(nbreaks == 0) {
tab <- data.frame(Parameter = c("Z", "sigma", "AIC"),
Estimate = c(results$par, sigma, AIC),
StdErr = c(SE, NA, NA))
}
else {
vname <- c(paste("Z", seq(1, nbreaks + 1, 1), sep = ""),
paste("Y", seq(1, nbreaks, 1), sep = ""), "sigma", "AIC")
est <- c(results$par)
est[(nbreaks +2):length(est)] <- est[(nbreaks +2):length(est)] + year[1] - 1
tab <- data.frame(Parameter = vname, Estimate = c(est, sigma, AIC), StdErr = c(SE, NA, NA))
}
if(graph) {
par(mfrow = c(1, 2))
plot(year, mlen, xlab = "Year", ylab = "Mean Length", typ = "b")
lines(year, Lpred, col = "red")
plot(year, mlen - Lpred, xlab = "Year", ylab = "Residual", typ = "b")
abline(h = 0, col = "red")
}
if(results$convergence == 0) status <- paste("Optimizer converged after",results$counts[1],"function calls.")
else status <- "Optimizer did not converge."
if(all(eigen(results$hessian)$values > 0)) hess.status <- "Hessian matrix is positive definite."
else hess.status <- "Hessian matrix is not positive definite. Check parameter estimates."
return(list(summary = tab, convergence = status, hessian = hess.status,
results = data.frame(year = year, observed = mlen, predicted = Lpred, residual = mlen - Lpred)))
}
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.