Nothing
R/fit_optim_maxim_adtneh2bis.R
In curesurv: Mixture and Non Mixture Parametric Cure Models to Estimate Cure Indicators
Defines functions
fit.opt.maxim.adtneh2bis
fit.opt.maxim.adtneh2bis <- function(x = x,
d = d,
z_tau = z_tau,
z_alpha = z_alpha,
theta_init = theta_init,
theta_lower = theta_lower,
theta_upper = theta_upper,
pophaz = pophaz,
cumpophaz = cumpophaz,
method_opt = method_opt,
pophaz.alpha = pophaz.alpha,
maxit_opt = maxit_opt,
gradient = gradient,
hessian_varcov = hessian_varcov,
ncoor_des = ncoor_des,
optim_func = optim_func,
iter_eps = iter_eps,
trace = trace,
optim_fixed = optim_fixed,
optimizer = optimizer) {
fixed <- optim_fixed
theta <- theta2 <- theta_init
n_z_tau <- ncol(z_tau)
n_z_alpha <- ncol(z_alpha)
n_z_tau_ad <- n_z_tau - 1
n_z_alpha_ad <- n_z_alpha - 1
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if (pophaz.alpha) {
if (n_z_tau >= 1 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"),
'alpha')
} else if (n_z_tau >= 1 & n_z_alpha == 0) {
names(theta) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
names(theta2) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"),
'alpha')
} else if (n_z_tau == 0 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
'pophaz.alpha')
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha ,sep = "_"),
'beta',
'tau0',
'alpha')
}
else if (n_z_tau == 0 & n_z_alpha == 0) {
names(theta) <- c('alpha0','beta','tau0', 'pophaz.alpha')
names(theta2) <- c('alpha0','beta','tau0', 'alpha')
}
}else{
if (n_z_tau >= 1 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"))
}else if (n_z_tau == 0 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0')
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha,sep = "_"),
'beta',
'tau0')
} else if (n_z_tau >= 1 & n_z_alpha == 0) {
names(theta) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
names(theta2) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"))
} else if (n_z_tau == 0 & n_z_alpha == 0) {
names(theta) <- c('alpha0','beta','tau0')
names(theta2) <- c('alpha0','beta','tau0')
}
}
if (gradient == TRUE) {
if (n_z_alpha > 0 & n_z_tau > 0) {
adtneh_ll_fn_deriv_i <- Deriv::Deriv( ~ -sum(
d * log(c(exp(alpha_pophaz)) * pophaz +
ifelse((x <= (tau)),
(((1/(1+exp(-x / (tau)))))^ ((alpha) - 1) * (1 - ((1/(1+exp(-x / (tau)))))) ^ ((beta - 1))),
0)) -
ifelse((x <= (tau)),
(tau) * beta(alpha, beta) * stats::pbeta(u, alpha, beta),
(tau) * beta(alpha, beta) * stats::pbeta(1, alpha, beta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)),
c("alpha0","alpha_k","beta","tau0","tau_z"),
combine = "cbind",
function(alpha0, alpha_k, beta, tau0, tau_z, x, d, pophaz, cumpophaz,
alpha_pophaz, z_alpha, z_tau, n_z_alpha, n_z_tau)NULL)
adtneh_ll_fn_deriv_all <- function(theta
) {
tau0 <- theta[n_z_alpha + 2 + 1]
alpha0 <- theta[1]
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- theta[n_z_alpha + 2]
tau0 <- theta[n_z_alpha + 2 + 1]
tau_z <- theta[(n_z_alpha + 2 + 1 + 1):(n_z_alpha + 2 + n_z_tau + 1)]
alpha <- (alpha0 + z_alpha %*% alpha_k)
tau <- (tau0 + z_tau %*% tau_z)
u <- (1 / (1 + exp(-x / (tau))))
alpha_pophaz <- c(matrix(1, ncol = 1, nrow = length(pophaz)) %*% 0)
pophaz <- c(pophaz)
cumpophaz <- c(cumpophaz)
gradval <- colSums(adtneh_ll_fn_deriv_i(alpha0,
alpha_k,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz,
z_alpha,
n_z_alpha,
n_z_tau))
names(gradval) <- NULL
return(gradval)
}
}else if (n_z_alpha == 0 & n_z_tau > 0) {
} else if (n_z_alpha > 0 & n_z_tau == 0) {
out <- statmod::gauss.quad(20, "legendre")
adtneh_ll_fn_deriv_i1 <-
function(alpha0,alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out){
.e1 <- exp(tau0)
.e2 <- x/2
.e4 <- exp(-(x/.e1))
.e5 <- 1 + .e4
.e7 <- exp(alpha0 + z_alpha %*% t(alpha_k))
.e8 <- 1/.e5
.e9 <- exp(beta)
.e12 <- matrix(.e2) %*% out$nodes + .e2
.e14 <- exp(-(.e12/.e1))
.e15 <- .e7 - 1
.e16 <- 1 + .e14
.e17 <- 1/.e16
.e18 <- 1 - .e8
.e19 <- .e18^.e9
.e20 <- .e8^.e15
.e21 <- 1 - .e17
.e22 <- .e21^.e9
.e24 <- .e19 * .e20 + pophaz * exp(alpha_pophaz)
.e25 <- exp(t(log(.e17))%*%.e15)
.e26 <- .e22 %*% .e25
.e28 <- d * .e19 * .e20
.e31 <- (log1p(.e14) * out$weights)*c(.e26)*c(.e7)
.e33 <- .e28 * .e7 * log1p(.e4)
.e34 <- .e7 - 2
.e35 <- .e9 - 1
.e36 <- z_alpha %*% rep(1, ncol(z_alpha))
cbind(alpha0 = .e33/.e24 + x * rowSums(x = -.e31)/2,
alpha_k = .e33 * .e36/.e24 + x * rowSums(x = -(.e31 * c(.e36)))/2,
beta = -(.e28 * .e9 * log(.e18)/.e24 - x * rowSums(x = c(.e26 * .e9) * log(.e21) * out$weights)/2),
tau0 = -(x * (d * (.e18^.e35 * .e20 * .e9 - .e19 * .e8^.e34 * .e15) * .e4/(.e24 * .e5^2 * .e1) -
rowSums(x = (.e21^.e35 * c(.e25) * .e9 -
.e22 * .e17^c(.e34) * c(.e15)) * .e14 * .e12 * out$weights/(.e16^2 * .e1))/2)))
}
adtneh_ll_fn_deriv_i2 <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,out)
{
x2 <- ifelse((x > exp(tau0)),exp(tau0), x)
.e1 <- x2/2
.e3 <- matrix(.e1) %*% out$nodes
.e4 <- exp(tau0)
.e5 <- .e3 + .e1
.e7 <- exp(-(.e5/.e4))
.e8 <- .e3 + x/2
.e11 <- exp(-(.e8/.e4))
.e12 <- exp(alpha0 + z_alpha %*% t(alpha_k))
.e13 <- 1 + .e11
.e14 <- exp(beta)
.e15 <- 1 - 1/.e13
.e16 <- 1 + .e7
.e17 <- 1/.e16
.e18 <- .e12 - 1
.e19 <- .e15^.e14
.e20 <- .e17^c(.e18)
.e21 <- .e19 * .e20
.e24 <- .e21 * c(.e12) * log1p(.e7) * out$weights
cbind(alpha0 = x2 * rowSums(x = -.e24)/2,
alpha_k = x2 * rowSums(x = -(.e24 * c(z_alpha %*% rep(1, ncol(z_alpha)))))/2,
beta = x2 * rowSums(x = .e21 * .e14 * log(.e15) * out$weights)/2,
tau0 = x2 * rowSums(x = (.e15^(.e14 - 1) * .e20 * .e11 * .e14 * .e8/.e13^2 - .e19 * .e17^(c(.e12 - 2)) * .e7 * c(.e18) * .e5/.e16^2) * out$weights/.e4)/2)
}
adtneh_ll_fn_deriv_i <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out) {
gradi <- rbind(
adtneh_ll_fn_deriv_i1(
alpha0,
alpha_k,
beta,
tau0,
x[which((x <= tau0))],
d[which((x <= tau0))],
pophaz[which((x <= tau0))],
cumpophaz[which((x <= tau0))],
alpha_pophaz[which((x <= tau0))],
matrix(z_alpha[which((x <= tau0)),], ncol = n_z_alpha),
n_z_alpha,
n_z_tau,
out
),
adtneh_ll_fn_deriv_i2(
alpha0, alpha_k,
beta,
tau0,
x[which((x > tau0))],
d[which((x > tau0))],
pophaz[which((x > tau0))],
cumpophaz[which((x > tau0))],
alpha_pophaz[which((x > tau0))],
matrix(z_alpha[which((x > tau0)),], ncol = n_z_alpha) ,
n_z_alpha,
n_z_tau,out)
)
return(gradi)
}
adtneh_ll_fn_deriv_all <- function(theta) {
alpha0 <- theta[1]
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
alpha <- (alpha0 + z_alpha %*% alpha_k)
tau <- (tau0)
u <- (1/(1 + exp(-x / (tau))))
alpha_pophaz <- c(matrix(1, ncol = 1, nrow = length(pophaz)) %*% 0)
pophaz <- c(pophaz)
cumpophaz <- c(cumpophaz)
gradval1 <- adtneh_ll_fn_deriv_i(alpha0, alpha_k, beta, tau0, x, d,
pophaz, cumpophaz,
alpha_pophaz, z_alpha,
n_z_alpha, n_z_tau, out)
gradval <- colSums(gradval1, na.rm = TRUE)
return(gradval)
}
} else if (n_z_alpha == 0 & n_z_tau == 0) {
out <- statmod::gauss.quad(20, "legendre")
adtneh_ll_fn_deriv_i1 <-
function(alpha0,alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out) {
.e1 <- exp(tau0)
.e2 <- x/2
.e3 <- exp(alpha0)
.e4 <- exp(beta)
.e5 <- x/.e1
.e6 <- .e3 - 1
.e9 <- matrix(.e2) %*% out$nodes + .e2
.e10 <- .e9/.e1
.e11 <- 1 + .e5
.e12 <- .e11^.e4
.e13 <- .e5^.e6
.e14 <- 1 - .e10
.e15 <- .e10^.e6
.e16 <- .e14^.e4
.e18 <- .e12 * .e13 + pophaz * exp(alpha_pophaz)
.e19 <- .e15 * .e16
.e20 <- d * .e12
.e21 <- .e3 - 2
.e22 <- .e4 - 1
cbind(alpha0 = -(.e20 * .e3 * (log(x) - tau0) * .e13/.e18 -
x * rowSums(x = .e19 * .e3 * (log(.e9) - tau0) * out$weights)/2),
beta = -(.e20 * .e4 * log1p(.e5) * .e13/.e18 - x *
rowSums(x = .e19 * .e4 * log(.e14) * out$weights)/2),
tau0 = x * (d * (.e11^.e22 * .e4 * .e13 + .e12 * .e6 * .e5^.e21)/(.e18 * .e1) +
rowSums(x = (.e15 * .e14^.e22 * .e4 - .e10^.e21 * .e16 * .e6) * .e9 * out$weights/.e1)/2))
}
adtneh_ll_fn_deriv_i2 <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,out) {
.e1 <- exp(tau0)*rep(1,length(x))
.e2 <- .e1/2
.e4 <- matrix(.e2) %*% out$nodes
.e5 <- .e2 + .e4
.e6 <- .e5/.e1
.e7 <- exp(alpha0)
.e8 <- exp(beta)
.e9 <- 1 - .e6
.e10 <- .e7 - 1
.e11 <- .e6^.e10
.e12 <- .e9^.e8
.e13 <- .e11 * .e12
cbind(alpha0 = .e1 * rowSums(x = .e13 * .e7 * (log(.e5) -
tau0) * out$weights)/2,
beta = .e1 * rowSums(x = .e13 * .e8 * log(.e9) * out$weights)/2,
tau0 = (0.5 * rowSums(.e13 * out$weights) +
rowSums(x = (.e6^(.e7 - 2) * .e12 * .e10 -
.e11 * .e9^(.e8 - 1) * .e8) *
(matrix(0.5 * .e1, nrow = , ncol = ,
byrow = , dimnames = ) %*% out$nodes - .e4) *
out$weights/.e1)/2) * .e1)
}
adtneh_ll_fn_deriv_i <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out) {
gradi <- rbind(
adtneh_ll_fn_deriv_i1(
alpha0,
beta,
tau0,
x[which((x <= exp(tau0)))],
d[which((x <= exp(tau0)))],
pophaz[which((x <= exp(tau0)))],
cumpophaz[which((x <= exp(tau0)))],
alpha_pophaz[which((x <= exp(tau0)))],
n_z_alpha,
n_z_tau,
out
),
adtneh_ll_fn_deriv_i2(
alpha0,
beta,
tau0,
x[which((x > exp(tau0)))],
d[which((x > exp(tau0)))],
pophaz[which((x > exp(tau0)))],
cumpophaz[which((x > exp(tau0)))],
alpha_pophaz[which((x > exp(tau0)))],
n_z_alpha,
n_z_tau,
out
)
)
return(gradi)
}
adtneh_ll_fn_deriv_all <- function(theta
) {
alpha0 <- theta[1]
beta <- (theta[2])
tau0 <- theta[3]
alpha <- (alpha0)
tau <- (tau0)
alpha_pophaz <- c(matrix(1, ncol = 1, nrow = length(pophaz)) %*% 0)
pophaz <- c(pophaz)
cumpophaz <- c(cumpophaz)
gradval1 <- adtneh_ll_fn_deriv_i(alpha0, beta, tau0, x, d,
pophaz, cumpophaz,
alpha_pophaz,
n_z_alpha, n_z_tau, out)
gradval <- colSums(gradval1, na.rm = TRUE)
return(gradval)
}
}
}
adtneh_ll_fn <- function(alpha0 = NULL,
alpha_k = NULL,
beta = NULL,
tau0 = NULL,
tau_z = NULL,
alpha = NULL) {
if (pophaz.alpha) {
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0, alpha)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
theta <- c(alpha0, alpha_k, beta, tau0, alpha)
}else if (n_z_alpha == 0 & n_z_tau > 0) {
theta <- c(alpha0, beta, tau0, tau_z, alpha)
}else if (n_z_alpha > 0 & n_z_tau > 0) {
theta <- c(alpha0, alpha_k, beta, tau0, tau_z, alpha)
}
}else{
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
theta <- c(alpha0, alpha_k, beta, tau0)
}else if (n_z_alpha == 0 & n_z_tau > 0) {
theta <- c(alpha0, beta, tau0, tau_z)
}else if (n_z_alpha > 0 & n_z_tau > 0) {
theta <- c(alpha0, alpha_k, beta, tau0, tau_z)
}
}
if (anyNA(theta))
return(1e10)
if (pophaz.alpha) {
alpha_pophaz <- matrix(1,
ncol = 1,
nrow = length(pophaz)) %*% theta[n_z_alpha + 2 + n_z_tau + 1 + 1]
}else{
alpha_pophaz <- matrix(1, ncol = 1,
nrow = length(pophaz)) %*% 0
}
n_z_tau <- ncol(z_tau)
n_z_alpha <- ncol(z_alpha)
n_z_tau_ad <- n_z_tau - 1
n_z_alpha_ad <- n_z_alpha - 1
alpha0 <- theta[1]
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
tau_z <- theta[(n_z_alpha + 2 + 1 + 1):(n_z_alpha + 2 + n_z_tau + 1)]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau > 0 & n_z_alpha == 0) {
beta <- (theta[2])
tau0 <- theta[2 + 1]
tau_z <- theta[(2 + 1 + 1):(2 + n_z_tau + 1)]
if (any(alpha0 <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha == 0) {
tau0 <- theta[n_z_alpha + 2 + 1]
if (alpha0 <= 0 | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}
if (anyNA(sum_minusloklik)) {
return(1e10)
} else {
return(sum_minusloklik)
}
}
adtneh_ll_fn3 <- function(alpha0, alpha_k, beta,
tau0, tau_z, alpha = NULL) {
if (pophaz.alpha) {
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0, alpha)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0, alpha)
}else if (n_z_alpha == 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, beta, tau0, tau_z, alpha)
}else if (n_z_alpha > 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0, tau_z, alpha)
}
}else{
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(parse(text = paste("c(",
paste(paste("alpha", 1:n_z_alpha,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0)
} else if (n_z_alpha == 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, beta, tau0, tau_z)
} else if (n_z_alpha > 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0, tau_z)
}
}
n_z_tau <- ncol(z_tau)
n_z_alpha <- ncol(z_alpha)
n_z_tau_ad <- n_z_tau - 1
n_z_alpha_ad <- n_z_alpha - 1
if (anyNA(theta))
return(1e10)
if (pophaz.alpha) {
alpha_pophaz <- matrix(1,
ncol = 1,
nrow = length(pophaz)) %*% theta[n_z_alpha + 2 + n_z_tau + 1 + 1]
}else{
alpha_pophaz <- matrix(1, ncol = 1,
nrow = length(pophaz)) %*% 0
}
alpha0 <- theta[1]
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
tau_z <- theta[(n_z_alpha + 2 + 1 + 1):(n_z_alpha + 2 + n_z_tau + 1)]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau > 0 & n_z_alpha == 0) {
beta <- (theta[2])
tau0 <- theta[2 + 1]
tau_z <- theta[(2 + 1 + 1):(2 + n_z_tau + 1)]
if (any(alpha0 <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha == 0) {
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
if (alpha0 <= 0 | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}
if (anyNA(sum_minusloklik)) {
return(1e10)
} else {
return(sum_minusloklik)
}
}
if (!is.null(theta_lower) && !is.null(theta_upper)) {
start_val <- gen_start_values(nvalues = 10000,
min = theta_lower,
max = theta_upper,
npar = length(theta_init))
} else if (NCD == "iter") {
stop("add bounds for initial values")
}
if (optim_func == "optim") {
if (method_opt == "L-BFGS-B") {
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
tempgr <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn_deriv_all(val)[i])
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
gr = tempgr,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_init[i] + rnorm(1))/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}else if ( NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
tempgr <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn_deriv_all(val)[i])
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
gr = tempgr,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
} else if (NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par,
TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}
} else {
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
} else if (NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}
else{
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
} else if ( NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}}
}else if (optim_func == "bbmle") {
if (method_opt == "L-BFGS-B" |
method_opt == "Rvmmin" |
method_opt == "Rcgmin" |
method_opt == "nlminb" |
method_opt == "bobyqa" ) {
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else {
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
tempgr <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn_deriv_all(val)[i])
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
gr = tempgr,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_init[i] + rnorm(1))/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else {
NCD <- ncoor_des
if (pophaz.alpha) {
alpha <- theta[length(theta)]
}else {
alpha <- 1
}
if (is.null(NCD)) {
par_adtneh2 <- genpar_adtneh2(theta, z_tau, z_alpha)
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else{
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau > 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}
add.arguments <- function(f, params) {
t <- paste("arg <- alist(",
paste(sapply(1:length(params),
function(i)
paste(names(params)[i], "=", sep = "")),
collapse = ","),
")", sep = "")
formals(f) <- eval(parse(text = t))
f
}
adtneh_ll_fn2noalpha2 <- add.arguments(adtneh_ll_fn3, start)
bbmle::parnames(adtneh_ll_fn2noalpha2) <- names(theta2)
start2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",(start)[i], sep = "")),
collapse = ","),
")", sep = "")
start3 <- eval(parse(text = start2))
theta_lower2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_lower)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_upper2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_upper)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_lower3 <- eval(parse(text = theta_lower2))
theta_upper3 <- eval(parse(text = theta_upper2))
optimized <- bbmle::mle2(
start = start3,
minuslogl = adtneh_ll_fn2noalpha2,
fixed = fixed,
optimizer = optimizer,
method = method_opt,
lower = unlist(theta_lower3),
upper = unlist(theta_upper3),
control = list(maxit = maxit_opt,
trace = trace))
} else {
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}
}
else
{ if (pophaz.alpha) {
alpha <- theta[length(theta)]
}else {
alpha <- 1
}
if (method_opt == "all.methods") {
NCD <- ncoor_des
if (is.null(NCD)) {
par_adtneh2 <- genpar_adtneh2(theta, z_tau, z_alpha)
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else{
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau > 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}
add.arguments <- function(f, params) {
t <- paste("arg <- alist(",
paste(sapply(1:length(params),
function(i)
paste(names(params)[i], "=", sep = "")),
collapse = ","),
")", sep = "")
formals(f) <- eval(parse(text = t))
f
}
adtneh_ll_fn2noalpha2 <- add.arguments(adtneh_ll_fn3, start)
bbmle::parnames(adtneh_ll_fn2noalpha2) <- names(theta2)
start2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",(start)[i], sep = "")),
collapse = ","),
")", sep = "")
start3 <- eval(parse(text = start2))
theta_lower2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_lower)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_upper2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_upper)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_lower3 <- eval(parse(text = theta_lower2))
theta_upper3 <- eval(parse(text = theta_upper2))
optimized2 <- optimx::optimx(
par = unlist(start3),
fn = adtneh_ll_fn,
lower = unlist(theta_lower3),
upper = unlist(theta_upper3),
control = list(all.methods = TRUE,
maximize = FALSE,
trace = trace,
save.failures = TRUE))
optimized <- bbmle::mle2(
start = start3,
minuslogl = adtneh_ll_fn2noalpha2,
fixed = fixed,
lower = unlist(theta_lower3),
upper = unlist(theta_upper3),
control = list(all.methods = TRUE,
maximize = FALSE,
trace = trace,
save.failures = TRUE))
} else {
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else{
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else{
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else{
par_adtneh2 <- genpar_adtneh2(theta, z_tau, z_alpha)
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else{
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau > 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}
add.arguments <- function(f, params) {
t <- paste("arg <- alist(",
paste(sapply(1:length(params),
function(i)
paste(names(params)[i], "=", sep = "")),
collapse = ","),
")", sep = "")
formals(f) <- eval(parse(text = t))
f
}
adtneh_ll_fn2noalpha2 <- add.arguments(adtneh_ll_fn3, start)
bbmle::parnames(adtneh_ll_fn2noalpha2) <- names(theta2)
if (pophaz.alpha) {
alpha <- theta[length(theta)]
}else {
alpha <- 1
}
start2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",(start)[i], sep = "")),
collapse = ","),
")", sep = "")
start3 <- eval(parse(text = start2))
theta_lower2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_lower)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_upper2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_upper)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_lower3 <- eval(parse(text = theta_lower2))
theta_upper3 <- eval(parse(text = theta_upper2))
optimized <- bbmle::mle2(
start = start3,
minuslogl = adtneh_ll_fn2noalpha2,
fixed = fixed,
optimizer = optimizer,
method = method_opt,
control = list(maxit = maxit_opt,
trace = trace))
}
}
}
}
if (inherits(optimized, "try-error")) {
stop("The function cannot be evaluated using these initial parameters!")
}
if (optim_func == "bbmle") {
par <- try(optimized@details$par, TRUE)
AIC <- try(2* optimized@details$value + 2*length(par), TRUE)
loglik <- try(-as.numeric(bbmle::logLik(optimized)) , TRUE)
evaluations <- try(optimized@details$counts, TRUE)
iterations <- try(evaluations, TRUE)
convergence <- try(optimized@details$convergence, TRUE)
message <- try(optimized@details$message, TRUE)
par_adtneh2 <- genpar_adtneh2(par, z_tau, z_alpha)
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[c(which(stringr::str_detect(names(par_adtneh2),
"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[c(which(stringr::str_detect(names(par_adtneh2),
"tau_z")))]
varcov_star <- optimized@vcov
if (any(inherits(optimized@vcov, "try-error"))) {
varcov_star <- try(solve(optimized@details$hessian), TRUE)
}
}else{
par <- try(optimized$par, TRUE)
AIC <- try(2*optimized$value + 2*length(par), TRUE)
loglik <- try(-optimized$value, TRUE)
evaluations <- try(optimized$counts, TRUE)
iterations <- try(evaluations, TRUE)
convergence <- try(optimized$convergence, TRUE)
message <- try(optimized$message, TRUE)
SD <- try(numDeriv::hessian(adtneh_ll_fn, par), TRUE)
if (hessian_varcov) {
varcov_star <- try(solve(SD), TRUE)
} else{
varcov_star <- try(adtneh_vartheta(z_tau,
z_alpha,
time = x,
theta = par,
d,
pophaz), TRUE)
}
if (any(inherits(varcov_star, "try-error"))) {
varcov_star <- try(solve(optimized$hessian), TRUE)
}
}
std_err_star <- try(sqrt(diag(varcov_star)), TRUE)
thetaest <- matrix(par, nrow = 1)
n_z_tau <- ncol(z_tau)
n_z_tau_ad <- n_z_tau - 1
if (pophaz.alpha) {
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
'pophaz.alpha')}else
if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0',
'pophaz.alpha')
}
}else
{
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0')
}else if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0')
}
}
thetaest2 <- thetaest
id_beta <- which(colnames(thetaest) == 'beta')
id_tau0 <- which(stringr::str_detect(colnames(thetaest),"tau"))
if (n_z_tau > 0) {
thetaest2[c(id_tau0)] <-
exp(c(thetaest[c(id_tau0)][1], thetaest[c(id_tau0)][1] + thetaest[c(id_tau0)[-c(1)]]))
} else{
thetaest2[id_tau0] <- exp(thetaest[id_tau0])
}
thetaest2[id_beta] <- exp(thetaest[id_beta]) + 1
if (pophaz.alpha) {
id_pophaz.alpha <- which(colnames(thetaest) == 'pophaz.alpha')
thetaest2[id_pophaz.alpha] <- exp(thetaest[id_pophaz.alpha])
idalpha_beta <- which(stringr::str_detect(colnames(thetaest), "alpha"))[-c(id_pophaz.alpha)]
if (n_z_alpha > 0) {
thetaest2[idalpha_beta] <- exp(c(thetaest[c(idalpha_beta)][1], thetaest[c(idalpha_beta)][1] + thetaest[c(idalpha_beta)[-c(1)]]))
}else{
thetaest2[idalpha_beta] <- exp(thetaest[c(idalpha_beta)])
}
}else{
idalpha_beta <- which(stringr::str_detect(colnames(thetaest), "alpha"))
if (n_z_alpha > 0) {
thetaest2[idalpha_beta] <- exp(c(thetaest[c(idalpha_beta)][1], thetaest[c(idalpha_beta)][1] + thetaest[c(idalpha_beta)[-c(1)]]))
}else {
thetaest2[idalpha_beta] <- exp(c(thetaest[c(idalpha_beta)]))
}
}
Dthetaest2 <- thetaest2
Dthetaest2[-id_beta] <- 1
Dthetaest2[id_beta] <- exp(thetaest[id_beta])
Dthetaest2[id_tau0] <- thetaest2[id_tau0]
if (pophaz.alpha) {
Dthetaest2[id_pophaz.alpha] <- exp(thetaest[id_pophaz.alpha])
}
varcov <- try(c(Dthetaest2^2)*varcov_star,TRUE)
std_err <- try((sqrt(diag(varcov))), TRUE)
if (inherits(std_err, "try-error")) {
}else{
std_err <- try(matrix(std_err, ncol = length(std_err)), TRUE)
if (pophaz.alpha) {
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
'pophaz.alpha')
}else
if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0',
'pophaz.alpha')
}
}else{
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0')
}else if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0')
}
}
colnames(thetaest) <- try(paste0(colnames(thetaest),"*"),TRUE)
std_err_star <- matrix(std_err_star, ncol = length(std_err_star))
colnames(std_err_star) <- try(paste0(colnames(thetaest),"*"), TRUE)
}
if (is.null(ncoor_des)) {
iter_coords <- 0
}else if (is.numeric(ncoor_des)) {
iter_coords <- ncoor_des
}else if (ncoor_des == "iter") {
iter_coords <- j
}
if (method_opt == "all.methods") {
return(
list(iter_coords = iter_coords,
coefficients = thetaest,
estimates = thetaest2,
loglik = loglik,
iterations = iterations,
evaluations = evaluations,
convergence = convergence,
message = message,
varcov = varcov,
varcov_star = varcov_star,
std_err = std_err,
std_err_star = std_err_star,
AIC = AIC,
optimized2 = optimized2
)
)
}else{
return(
list(iter_coords = iter_coords,
coefficients = thetaest,
estimates = thetaest2,
loglik = loglik,
iterations = iterations,
evaluations = evaluations,
convergence = convergence,
message = message,
varcov = varcov,
varcov_star = varcov_star,
std_err = std_err,
std_err_star = std_err_star,
AIC = AIC
)
)
}
}
Try the curesurv package in your browser
Any scripts or data that you put into this service are public.
curesurv documentation built on April 12, 2025, 2:21 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fit.opt.maxim.adtneh2bis
fit.opt.maxim.adtneh2bis <- function(x = x,
d = d,
z_tau = z_tau,
z_alpha = z_alpha,
theta_init = theta_init,
theta_lower = theta_lower,
theta_upper = theta_upper,
pophaz = pophaz,
cumpophaz = cumpophaz,
method_opt = method_opt,
pophaz.alpha = pophaz.alpha,
maxit_opt = maxit_opt,
gradient = gradient,
hessian_varcov = hessian_varcov,
ncoor_des = ncoor_des,
optim_func = optim_func,
iter_eps = iter_eps,
trace = trace,
optim_fixed = optim_fixed,
optimizer = optimizer) {
fixed <- optim_fixed
theta <- theta2 <- theta_init
n_z_tau <- ncol(z_tau)
n_z_alpha <- ncol(z_alpha)
n_z_tau_ad <- n_z_tau - 1
n_z_alpha_ad <- n_z_alpha - 1
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if (pophaz.alpha) {
if (n_z_tau >= 1 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"),
'alpha')
} else if (n_z_tau >= 1 & n_z_alpha == 0) {
names(theta) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
names(theta2) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"),
'alpha')
} else if (n_z_tau == 0 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
'pophaz.alpha')
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha ,sep = "_"),
'beta',
'tau0',
'alpha')
}
else if (n_z_tau == 0 & n_z_alpha == 0) {
names(theta) <- c('alpha0','beta','tau0', 'pophaz.alpha')
names(theta2) <- c('alpha0','beta','tau0', 'alpha')
}
}else{
if (n_z_tau >= 1 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"))
}else if (n_z_tau == 0 & n_z_alpha >= 1) {
names(theta) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0')
names(theta2) <- c('alpha0',
paste('alpha', 1:n_z_alpha,sep = "_"),
'beta',
'tau0')
} else if (n_z_tau >= 1 & n_z_alpha == 0) {
names(theta) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
names(theta2) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, sep = "_"))
} else if (n_z_tau == 0 & n_z_alpha == 0) {
names(theta) <- c('alpha0','beta','tau0')
names(theta2) <- c('alpha0','beta','tau0')
}
}
if (gradient == TRUE) {
if (n_z_alpha > 0 & n_z_tau > 0) {
adtneh_ll_fn_deriv_i <- Deriv::Deriv( ~ -sum(
d * log(c(exp(alpha_pophaz)) * pophaz +
ifelse((x <= (tau)),
(((1/(1+exp(-x / (tau)))))^ ((alpha) - 1) * (1 - ((1/(1+exp(-x / (tau)))))) ^ ((beta - 1))),
0)) -
ifelse((x <= (tau)),
(tau) * beta(alpha, beta) * stats::pbeta(u, alpha, beta),
(tau) * beta(alpha, beta) * stats::pbeta(1, alpha, beta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)),
c("alpha0","alpha_k","beta","tau0","tau_z"),
combine = "cbind",
function(alpha0, alpha_k, beta, tau0, tau_z, x, d, pophaz, cumpophaz,
alpha_pophaz, z_alpha, z_tau, n_z_alpha, n_z_tau)NULL)
adtneh_ll_fn_deriv_all <- function(theta
) {
tau0 <- theta[n_z_alpha + 2 + 1]
alpha0 <- theta[1]
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- theta[n_z_alpha + 2]
tau0 <- theta[n_z_alpha + 2 + 1]
tau_z <- theta[(n_z_alpha + 2 + 1 + 1):(n_z_alpha + 2 + n_z_tau + 1)]
alpha <- (alpha0 + z_alpha %*% alpha_k)
tau <- (tau0 + z_tau %*% tau_z)
u <- (1 / (1 + exp(-x / (tau))))
alpha_pophaz <- c(matrix(1, ncol = 1, nrow = length(pophaz)) %*% 0)
pophaz <- c(pophaz)
cumpophaz <- c(cumpophaz)
gradval <- colSums(adtneh_ll_fn_deriv_i(alpha0,
alpha_k,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz,
z_alpha,
n_z_alpha,
n_z_tau))
names(gradval) <- NULL
return(gradval)
}
}else if (n_z_alpha == 0 & n_z_tau > 0) {
} else if (n_z_alpha > 0 & n_z_tau == 0) {
out <- statmod::gauss.quad(20, "legendre")
adtneh_ll_fn_deriv_i1 <-
function(alpha0,alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out){
.e1 <- exp(tau0)
.e2 <- x/2
.e4 <- exp(-(x/.e1))
.e5 <- 1 + .e4
.e7 <- exp(alpha0 + z_alpha %*% t(alpha_k))
.e8 <- 1/.e5
.e9 <- exp(beta)
.e12 <- matrix(.e2) %*% out$nodes + .e2
.e14 <- exp(-(.e12/.e1))
.e15 <- .e7 - 1
.e16 <- 1 + .e14
.e17 <- 1/.e16
.e18 <- 1 - .e8
.e19 <- .e18^.e9
.e20 <- .e8^.e15
.e21 <- 1 - .e17
.e22 <- .e21^.e9
.e24 <- .e19 * .e20 + pophaz * exp(alpha_pophaz)
.e25 <- exp(t(log(.e17))%*%.e15)
.e26 <- .e22 %*% .e25
.e28 <- d * .e19 * .e20
.e31 <- (log1p(.e14) * out$weights)*c(.e26)*c(.e7)
.e33 <- .e28 * .e7 * log1p(.e4)
.e34 <- .e7 - 2
.e35 <- .e9 - 1
.e36 <- z_alpha %*% rep(1, ncol(z_alpha))
cbind(alpha0 = .e33/.e24 + x * rowSums(x = -.e31)/2,
alpha_k = .e33 * .e36/.e24 + x * rowSums(x = -(.e31 * c(.e36)))/2,
beta = -(.e28 * .e9 * log(.e18)/.e24 - x * rowSums(x = c(.e26 * .e9) * log(.e21) * out$weights)/2),
tau0 = -(x * (d * (.e18^.e35 * .e20 * .e9 - .e19 * .e8^.e34 * .e15) * .e4/(.e24 * .e5^2 * .e1) -
rowSums(x = (.e21^.e35 * c(.e25) * .e9 -
.e22 * .e17^c(.e34) * c(.e15)) * .e14 * .e12 * out$weights/(.e16^2 * .e1))/2)))
}
adtneh_ll_fn_deriv_i2 <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,out)
{
x2 <- ifelse((x > exp(tau0)),exp(tau0), x)
.e1 <- x2/2
.e3 <- matrix(.e1) %*% out$nodes
.e4 <- exp(tau0)
.e5 <- .e3 + .e1
.e7 <- exp(-(.e5/.e4))
.e8 <- .e3 + x/2
.e11 <- exp(-(.e8/.e4))
.e12 <- exp(alpha0 + z_alpha %*% t(alpha_k))
.e13 <- 1 + .e11
.e14 <- exp(beta)
.e15 <- 1 - 1/.e13
.e16 <- 1 + .e7
.e17 <- 1/.e16
.e18 <- .e12 - 1
.e19 <- .e15^.e14
.e20 <- .e17^c(.e18)
.e21 <- .e19 * .e20
.e24 <- .e21 * c(.e12) * log1p(.e7) * out$weights
cbind(alpha0 = x2 * rowSums(x = -.e24)/2,
alpha_k = x2 * rowSums(x = -(.e24 * c(z_alpha %*% rep(1, ncol(z_alpha)))))/2,
beta = x2 * rowSums(x = .e21 * .e14 * log(.e15) * out$weights)/2,
tau0 = x2 * rowSums(x = (.e15^(.e14 - 1) * .e20 * .e11 * .e14 * .e8/.e13^2 - .e19 * .e17^(c(.e12 - 2)) * .e7 * c(.e18) * .e5/.e16^2) * out$weights/.e4)/2)
}
adtneh_ll_fn_deriv_i <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out) {
gradi <- rbind(
adtneh_ll_fn_deriv_i1(
alpha0,
alpha_k,
beta,
tau0,
x[which((x <= tau0))],
d[which((x <= tau0))],
pophaz[which((x <= tau0))],
cumpophaz[which((x <= tau0))],
alpha_pophaz[which((x <= tau0))],
matrix(z_alpha[which((x <= tau0)),], ncol = n_z_alpha),
n_z_alpha,
n_z_tau,
out
),
adtneh_ll_fn_deriv_i2(
alpha0, alpha_k,
beta,
tau0,
x[which((x > tau0))],
d[which((x > tau0))],
pophaz[which((x > tau0))],
cumpophaz[which((x > tau0))],
alpha_pophaz[which((x > tau0))],
matrix(z_alpha[which((x > tau0)),], ncol = n_z_alpha) ,
n_z_alpha,
n_z_tau,out)
)
return(gradi)
}
adtneh_ll_fn_deriv_all <- function(theta) {
alpha0 <- theta[1]
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
alpha <- (alpha0 + z_alpha %*% alpha_k)
tau <- (tau0)
u <- (1/(1 + exp(-x / (tau))))
alpha_pophaz <- c(matrix(1, ncol = 1, nrow = length(pophaz)) %*% 0)
pophaz <- c(pophaz)
cumpophaz <- c(cumpophaz)
gradval1 <- adtneh_ll_fn_deriv_i(alpha0, alpha_k, beta, tau0, x, d,
pophaz, cumpophaz,
alpha_pophaz, z_alpha,
n_z_alpha, n_z_tau, out)
gradval <- colSums(gradval1, na.rm = TRUE)
return(gradval)
}
} else if (n_z_alpha == 0 & n_z_tau == 0) {
out <- statmod::gauss.quad(20, "legendre")
adtneh_ll_fn_deriv_i1 <-
function(alpha0,alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out) {
.e1 <- exp(tau0)
.e2 <- x/2
.e3 <- exp(alpha0)
.e4 <- exp(beta)
.e5 <- x/.e1
.e6 <- .e3 - 1
.e9 <- matrix(.e2) %*% out$nodes + .e2
.e10 <- .e9/.e1
.e11 <- 1 + .e5
.e12 <- .e11^.e4
.e13 <- .e5^.e6
.e14 <- 1 - .e10
.e15 <- .e10^.e6
.e16 <- .e14^.e4
.e18 <- .e12 * .e13 + pophaz * exp(alpha_pophaz)
.e19 <- .e15 * .e16
.e20 <- d * .e12
.e21 <- .e3 - 2
.e22 <- .e4 - 1
cbind(alpha0 = -(.e20 * .e3 * (log(x) - tau0) * .e13/.e18 -
x * rowSums(x = .e19 * .e3 * (log(.e9) - tau0) * out$weights)/2),
beta = -(.e20 * .e4 * log1p(.e5) * .e13/.e18 - x *
rowSums(x = .e19 * .e4 * log(.e14) * out$weights)/2),
tau0 = x * (d * (.e11^.e22 * .e4 * .e13 + .e12 * .e6 * .e5^.e21)/(.e18 * .e1) +
rowSums(x = (.e15 * .e14^.e22 * .e4 - .e10^.e21 * .e16 * .e6) * .e9 * out$weights/.e1)/2))
}
adtneh_ll_fn_deriv_i2 <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,out) {
.e1 <- exp(tau0)*rep(1,length(x))
.e2 <- .e1/2
.e4 <- matrix(.e2) %*% out$nodes
.e5 <- .e2 + .e4
.e6 <- .e5/.e1
.e7 <- exp(alpha0)
.e8 <- exp(beta)
.e9 <- 1 - .e6
.e10 <- .e7 - 1
.e11 <- .e6^.e10
.e12 <- .e9^.e8
.e13 <- .e11 * .e12
cbind(alpha0 = .e1 * rowSums(x = .e13 * .e7 * (log(.e5) -
tau0) * out$weights)/2,
beta = .e1 * rowSums(x = .e13 * .e8 * log(.e9) * out$weights)/2,
tau0 = (0.5 * rowSums(.e13 * out$weights) +
rowSums(x = (.e6^(.e7 - 2) * .e12 * .e10 -
.e11 * .e9^(.e8 - 1) * .e8) *
(matrix(0.5 * .e1, nrow = , ncol = ,
byrow = , dimnames = ) %*% out$nodes - .e4) *
out$weights/.e1)/2) * .e1)
}
adtneh_ll_fn_deriv_i <-
function(alpha0, alpha_k=NULL,
beta,
tau0,
x,
d,
pophaz,
cumpophaz,
alpha_pophaz, z_alpha=NULL,
n_z_alpha,
n_z_tau,
out) {
gradi <- rbind(
adtneh_ll_fn_deriv_i1(
alpha0,
beta,
tau0,
x[which((x <= exp(tau0)))],
d[which((x <= exp(tau0)))],
pophaz[which((x <= exp(tau0)))],
cumpophaz[which((x <= exp(tau0)))],
alpha_pophaz[which((x <= exp(tau0)))],
n_z_alpha,
n_z_tau,
out
),
adtneh_ll_fn_deriv_i2(
alpha0,
beta,
tau0,
x[which((x > exp(tau0)))],
d[which((x > exp(tau0)))],
pophaz[which((x > exp(tau0)))],
cumpophaz[which((x > exp(tau0)))],
alpha_pophaz[which((x > exp(tau0)))],
n_z_alpha,
n_z_tau,
out
)
)
return(gradi)
}
adtneh_ll_fn_deriv_all <- function(theta
) {
alpha0 <- theta[1]
beta <- (theta[2])
tau0 <- theta[3]
alpha <- (alpha0)
tau <- (tau0)
alpha_pophaz <- c(matrix(1, ncol = 1, nrow = length(pophaz)) %*% 0)
pophaz <- c(pophaz)
cumpophaz <- c(cumpophaz)
gradval1 <- adtneh_ll_fn_deriv_i(alpha0, beta, tau0, x, d,
pophaz, cumpophaz,
alpha_pophaz,
n_z_alpha, n_z_tau, out)
gradval <- colSums(gradval1, na.rm = TRUE)
return(gradval)
}
}
}
adtneh_ll_fn <- function(alpha0 = NULL,
alpha_k = NULL,
beta = NULL,
tau0 = NULL,
tau_z = NULL,
alpha = NULL) {
if (pophaz.alpha) {
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0, alpha)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
theta <- c(alpha0, alpha_k, beta, tau0, alpha)
}else if (n_z_alpha == 0 & n_z_tau > 0) {
theta <- c(alpha0, beta, tau0, tau_z, alpha)
}else if (n_z_alpha > 0 & n_z_tau > 0) {
theta <- c(alpha0, alpha_k, beta, tau0, tau_z, alpha)
}
}else{
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
theta <- c(alpha0, alpha_k, beta, tau0)
}else if (n_z_alpha == 0 & n_z_tau > 0) {
theta <- c(alpha0, beta, tau0, tau_z)
}else if (n_z_alpha > 0 & n_z_tau > 0) {
theta <- c(alpha0, alpha_k, beta, tau0, tau_z)
}
}
if (anyNA(theta))
return(1e10)
if (pophaz.alpha) {
alpha_pophaz <- matrix(1,
ncol = 1,
nrow = length(pophaz)) %*% theta[n_z_alpha + 2 + n_z_tau + 1 + 1]
}else{
alpha_pophaz <- matrix(1, ncol = 1,
nrow = length(pophaz)) %*% 0
}
n_z_tau <- ncol(z_tau)
n_z_alpha <- ncol(z_alpha)
n_z_tau_ad <- n_z_tau - 1
n_z_alpha_ad <- n_z_alpha - 1
alpha0 <- theta[1]
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
tau_z <- theta[(n_z_alpha + 2 + 1 + 1):(n_z_alpha + 2 + n_z_tau + 1)]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau > 0 & n_z_alpha == 0) {
beta <- (theta[2])
tau0 <- theta[2 + 1]
tau_z <- theta[(2 + 1 + 1):(2 + n_z_tau + 1)]
if (any(alpha0 <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha == 0) {
tau0 <- theta[n_z_alpha + 2 + 1]
if (alpha0 <= 0 | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}
if (anyNA(sum_minusloklik)) {
return(1e10)
} else {
return(sum_minusloklik)
}
}
adtneh_ll_fn3 <- function(alpha0, alpha_k, beta,
tau0, tau_z, alpha = NULL) {
if (pophaz.alpha) {
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0, alpha)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0, alpha)
}else if (n_z_alpha == 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, beta, tau0, tau_z, alpha)
}else if (n_z_alpha > 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0, tau_z, alpha)
}
}else{
if (n_z_alpha == 0 & n_z_tau == 0) {
theta <- c(alpha0, beta, tau0)
} else if (n_z_alpha > 0 & n_z_tau == 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(parse(text = paste("c(",
paste(paste("alpha", 1:n_z_alpha,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0)
} else if (n_z_alpha == 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, beta, tau0, tau_z)
} else if (n_z_alpha > 0 & n_z_tau > 0) {
if (optim_func == "bbmle") {
alpha_k <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
tau_z <- eval(
parse(text = paste("c(",
paste(paste("tau", 1:n_z_tau,
sep = "_"),
collapse = ","),
")", sep = "")))
}
theta <- c(alpha0, alpha_k, beta, tau0, tau_z)
}
}
n_z_tau <- ncol(z_tau)
n_z_alpha <- ncol(z_alpha)
n_z_tau_ad <- n_z_tau - 1
n_z_alpha_ad <- n_z_alpha - 1
if (anyNA(theta))
return(1e10)
if (pophaz.alpha) {
alpha_pophaz <- matrix(1,
ncol = 1,
nrow = length(pophaz)) %*% theta[n_z_alpha + 2 + n_z_tau + 1 + 1]
}else{
alpha_pophaz <- matrix(1, ncol = 1,
nrow = length(pophaz)) %*% 0
}
alpha0 <- theta[1]
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
tau_z <- theta[(n_z_alpha + 2 + 1 + 1):(n_z_alpha + 2 + n_z_tau + 1)]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau > 0 & n_z_alpha == 0) {
beta <- (theta[2])
tau0 <- theta[2 + 1]
tau_z <- theta[(2 + 1 + 1):(2 + n_z_tau + 1)]
if (any(alpha0 <= 0) | any(tau0 + z_tau %*% tau_z <= 0)) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha_k <- theta[2:(n_z_alpha + 1)]
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
if (any(alpha0 + z_alpha %*% alpha_k <= 0) | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}else if (n_z_tau == 0 & n_z_alpha == 0) {
beta <- (theta[n_z_alpha + 2])
tau0 <- theta[n_z_alpha + 2 + 1]
if (alpha0 <= 0 | tau0 <= 0) {
sum_minusloklik <- 1e10
}else {
sum_minusloklik <- -sum(d * log(c(exp(alpha_pophaz)) * pophaz +
lexc_ad2bis(z_tau, z_alpha, x, theta)) -
cumLexc_ad2bis(z_tau, z_alpha, x, theta) -
cumpophaz * c(exp(alpha_pophaz)), na.rm = TRUE)
}
}
if (anyNA(sum_minusloklik)) {
return(1e10)
} else {
return(sum_minusloklik)
}
}
if (!is.null(theta_lower) && !is.null(theta_upper)) {
start_val <- gen_start_values(nvalues = 10000,
min = theta_lower,
max = theta_upper,
npar = length(theta_init))
} else if (NCD == "iter") {
stop("add bounds for initial values")
}
if (optim_func == "optim") {
if (method_opt == "L-BFGS-B") {
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
tempgr <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn_deriv_all(val)[i])
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
gr = tempgr,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_init[i] + rnorm(1))/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}else if ( NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
tempgr <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn_deriv_all(val)[i])
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
gr = tempgr,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
} else if (NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par,
TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}
} else {
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
} else if (NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}
else{
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else if (is.numeric(NCD)) {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
} else if ( NCD == "iter") {
optimized_theta <- try(stats::optim(par = theta,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
if (any(is.na(optimized_theta) || (inherits(optimized_theta, "try-error")))) {
init.cd <- theta
} else {
init.cd <- optimized_theta
}
j = 0
mydiff <- matrix(data = 1,
nrow = 10000,
ncol = 1)
loglik_val <- matrix(data = 0,
nrow = 10000,
ncol = 1)
theta_val <- matrix(data = 0,
nrow = 10000,
ncol = length(theta))
while (abs(mydiff[j + 1, ]) > iter_eps && (j = j + 1) < 10000) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.valpar <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500)), TRUE)
new.val <- try(new.valpar$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- init.cd[i] + 1
}
init.cd <- replace(init.cd, i, new.val )
}
theta_val[j, ] <- theta
theta_val[j + 1, ] <- try(init.cd, TRUE)
loglik_val[j + 1,] <- try(c(new.valpar$value), TRUE)
mydiff[j + 1, ] <- try(c(loglik_val[j + 1, ] -
loglik_val[j, ]), TRUE)
message("difference between loglikelihood at iter = ",
j + 1,
"equal",
mydiff[j + 1, ],
". The j+1 th logliklihood equal",
loglik_val[j + 1, ],
".\n")
if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] > 0) {
loglik_val[(j + 1), ] = 10000
mydiff[(j + 1), ] = 10000
} else if (abs(loglik_val[(j + 1), ]) == Inf |
abs(loglik_val[(j + 1), ]) == .Machine$double.xmax &
loglik_val[(j + 1), ] < 0) {
loglik_val[(j + 1), ] = -10000
mydiff[(j + 1), ] = -10000
}
if (j >= 2) {
par(mfrow = c(1, 2))
plot(2:(j + 1),
loglik_val[2:(j + 1),],
type = "l",
col = 2,
ylab = "loglik",
xlab = "iter")
graphics::grid()
plot(2:(j + 1),
abs(mydiff[2:(j + 1),]),
type = "l",
col = "green",
ylab = "loglik_abs_diff",
xlab = "iter")
graphics::grid()
par(mfrow = c(1, 1))
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}}
}else if (optim_func == "bbmle") {
if (method_opt == "L-BFGS-B" |
method_opt == "Rvmmin" |
method_opt == "Rcgmin" |
method_opt == "nlminb" |
method_opt == "bobyqa" ) {
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else {
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
tempgr <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn_deriv_all(val)[i])
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
gr = tempgr,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_init[i] + rnorm(1))/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else {
NCD <- ncoor_des
if (pophaz.alpha) {
alpha <- theta[length(theta)]
}else {
alpha <- 1
}
if (is.null(NCD)) {
par_adtneh2 <- genpar_adtneh2(theta, z_tau, z_alpha)
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else{
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau > 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}
add.arguments <- function(f, params) {
t <- paste("arg <- alist(",
paste(sapply(1:length(params),
function(i)
paste(names(params)[i], "=", sep = "")),
collapse = ","),
")", sep = "")
formals(f) <- eval(parse(text = t))
f
}
adtneh_ll_fn2noalpha2 <- add.arguments(adtneh_ll_fn3, start)
bbmle::parnames(adtneh_ll_fn2noalpha2) <- names(theta2)
start2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",(start)[i], sep = "")),
collapse = ","),
")", sep = "")
start3 <- eval(parse(text = start2))
theta_lower2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_lower)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_upper2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_upper)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_lower3 <- eval(parse(text = theta_lower2))
theta_upper3 <- eval(parse(text = theta_upper2))
optimized <- bbmle::mle2(
start = start3,
minuslogl = adtneh_ll_fn2noalpha2,
fixed = fixed,
optimizer = optimizer,
method = method_opt,
lower = unlist(theta_lower3),
upper = unlist(theta_upper3),
control = list(maxit = maxit_opt,
trace = trace))
} else {
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}
}
else
{ if (pophaz.alpha) {
alpha <- theta[length(theta)]
}else {
alpha <- 1
}
if (method_opt == "all.methods") {
NCD <- ncoor_des
if (is.null(NCD)) {
par_adtneh2 <- genpar_adtneh2(theta, z_tau, z_alpha)
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else{
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau > 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}
add.arguments <- function(f, params) {
t <- paste("arg <- alist(",
paste(sapply(1:length(params),
function(i)
paste(names(params)[i], "=", sep = "")),
collapse = ","),
")", sep = "")
formals(f) <- eval(parse(text = t))
f
}
adtneh_ll_fn2noalpha2 <- add.arguments(adtneh_ll_fn3, start)
bbmle::parnames(adtneh_ll_fn2noalpha2) <- names(theta2)
start2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",(start)[i], sep = "")),
collapse = ","),
")", sep = "")
start3 <- eval(parse(text = start2))
theta_lower2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_lower)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_upper2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_upper)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_lower3 <- eval(parse(text = theta_lower2))
theta_upper3 <- eval(parse(text = theta_upper2))
optimized2 <- optimx::optimx(
par = unlist(start3),
fn = adtneh_ll_fn,
lower = unlist(theta_lower3),
upper = unlist(theta_upper3),
control = list(all.methods = TRUE,
maximize = FALSE,
trace = trace,
save.failures = TRUE))
optimized <- bbmle::mle2(
start = start3,
minuslogl = adtneh_ll_fn2noalpha2,
fixed = fixed,
lower = unlist(theta_lower3),
upper = unlist(theta_upper3),
control = list(all.methods = TRUE,
maximize = FALSE,
trace = trace,
save.failures = TRUE))
} else {
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
hessian = TRUE,
lower = theta_lower[i],
upper = theta_upper[i],
control = list(
maxit = maxit_opt, trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
method = method_opt,
hessian = TRUE,
lower = theta_lower,
upper = theta_upper,
control = list(maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else{
if (gradient) {
NCD <- ncoor_des
if (is.null(NCD)) {
optimized <- stats::optim(par = theta,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))
} else{
init.cd <- theta
for (j in 1:NCD) {
message("cycle = ", j, "\n")
for (i in 1:length(init.cd)) {
tempf <- Vectorize(function(par){
val <- replace(init.cd, i, par)
return(adtneh_ll_fn(val))
})
new.val <- try(stats::optim(par = init.cd[i],
fn = tempf,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500))$par, TRUE)
message("new initial value of parameter number ",i," is equal to: ",new.val,"\n")
if (anyNA(new.val) | inherits(new.val, "try-error")) {
new.val <- (theta_lower[i] + theta_upper[i])/2
}
init.cd <- replace(init.cd, i, new.val )
}
}
if (any(is.na(init.cd))) {
initG <- theta
} else {
initG <- init.cd
}
optimized <- try(stats::optim(par = initG ,
fn = adtneh_ll_fn,
gr = adtneh_ll_fn_deriv_all,
method = method_opt,
control = list(pgtol = 1e-15,
maxit = maxit_opt,
trace = trace,
REPORT = 500)), TRUE)
}
}else{
par_adtneh2 <- genpar_adtneh2(theta, z_tau, z_alpha)
if (n_z_tau > 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else{
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau > 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"tau_z")))]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
tau_z = as.numeric(tau_z))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha > 0) {
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[
c(which(stringr::str_detect(names(par_adtneh2),"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
alpha_k = as.numeric(alpha_k),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}else if (n_z_tau == 0 & n_z_alpha == 0) {
alpha0 <- par_adtneh2["alpha0"]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
if (pophaz.alpha) {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0),
alpha = as.numeric(alpha))
}else {
start <- c(alpha0 = as.numeric(alpha0),
beta = as.numeric(beta),
tau0 = as.numeric(tau0))
}
names(start) <- names(theta2)
}
add.arguments <- function(f, params) {
t <- paste("arg <- alist(",
paste(sapply(1:length(params),
function(i)
paste(names(params)[i], "=", sep = "")),
collapse = ","),
")", sep = "")
formals(f) <- eval(parse(text = t))
f
}
adtneh_ll_fn2noalpha2 <- add.arguments(adtneh_ll_fn3, start)
bbmle::parnames(adtneh_ll_fn2noalpha2) <- names(theta2)
if (pophaz.alpha) {
alpha <- theta[length(theta)]
}else {
alpha <- 1
}
start2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",(start)[i], sep = "")),
collapse = ","),
")", sep = "")
start3 <- eval(parse(text = start2))
theta_lower2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_lower)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_upper2 <- paste("list(",
paste(sapply(1:length(start),
function(i)
paste(names(start)[i], "=",
(theta_upper)[i], sep = "")),
collapse = ","),
")", sep = "")
theta_lower3 <- eval(parse(text = theta_lower2))
theta_upper3 <- eval(parse(text = theta_upper2))
optimized <- bbmle::mle2(
start = start3,
minuslogl = adtneh_ll_fn2noalpha2,
fixed = fixed,
optimizer = optimizer,
method = method_opt,
control = list(maxit = maxit_opt,
trace = trace))
}
}
}
}
if (inherits(optimized, "try-error")) {
stop("The function cannot be evaluated using these initial parameters!")
}
if (optim_func == "bbmle") {
par <- try(optimized@details$par, TRUE)
AIC <- try(2* optimized@details$value + 2*length(par), TRUE)
loglik <- try(-as.numeric(bbmle::logLik(optimized)) , TRUE)
evaluations <- try(optimized@details$counts, TRUE)
iterations <- try(evaluations, TRUE)
convergence <- try(optimized@details$convergence, TRUE)
message <- try(optimized@details$message, TRUE)
par_adtneh2 <- genpar_adtneh2(par, z_tau, z_alpha)
alpha0 <- par_adtneh2["alpha0"]
alpha_k <- par_adtneh2[c(which(stringr::str_detect(names(par_adtneh2),
"alpha_k")))]
beta <- par_adtneh2["beta"]
tau0 <- par_adtneh2["tau0"]
tau_z <- par_adtneh2[c(which(stringr::str_detect(names(par_adtneh2),
"tau_z")))]
varcov_star <- optimized@vcov
if (any(inherits(optimized@vcov, "try-error"))) {
varcov_star <- try(solve(optimized@details$hessian), TRUE)
}
}else{
par <- try(optimized$par, TRUE)
AIC <- try(2*optimized$value + 2*length(par), TRUE)
loglik <- try(-optimized$value, TRUE)
evaluations <- try(optimized$counts, TRUE)
iterations <- try(evaluations, TRUE)
convergence <- try(optimized$convergence, TRUE)
message <- try(optimized$message, TRUE)
SD <- try(numDeriv::hessian(adtneh_ll_fn, par), TRUE)
if (hessian_varcov) {
varcov_star <- try(solve(SD), TRUE)
} else{
varcov_star <- try(adtneh_vartheta(z_tau,
z_alpha,
time = x,
theta = par,
d,
pophaz), TRUE)
}
if (any(inherits(varcov_star, "try-error"))) {
varcov_star <- try(solve(optimized$hessian), TRUE)
}
}
std_err_star <- try(sqrt(diag(varcov_star)), TRUE)
thetaest <- matrix(par, nrow = 1)
n_z_tau <- ncol(z_tau)
n_z_tau_ad <- n_z_tau - 1
if (pophaz.alpha) {
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
'pophaz.alpha')}else
if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0',
'pophaz.alpha')
}
}else
{
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(thetaest) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0')
}else if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(thetaest) <- c('alpha0',
'beta',
'tau0')
}
}
thetaest2 <- thetaest
id_beta <- which(colnames(thetaest) == 'beta')
id_tau0 <- which(stringr::str_detect(colnames(thetaest),"tau"))
if (n_z_tau > 0) {
thetaest2[c(id_tau0)] <-
exp(c(thetaest[c(id_tau0)][1], thetaest[c(id_tau0)][1] + thetaest[c(id_tau0)[-c(1)]]))
} else{
thetaest2[id_tau0] <- exp(thetaest[id_tau0])
}
thetaest2[id_beta] <- exp(thetaest[id_beta]) + 1
if (pophaz.alpha) {
id_pophaz.alpha <- which(colnames(thetaest) == 'pophaz.alpha')
thetaest2[id_pophaz.alpha] <- exp(thetaest[id_pophaz.alpha])
idalpha_beta <- which(stringr::str_detect(colnames(thetaest), "alpha"))[-c(id_pophaz.alpha)]
if (n_z_alpha > 0) {
thetaest2[idalpha_beta] <- exp(c(thetaest[c(idalpha_beta)][1], thetaest[c(idalpha_beta)][1] + thetaest[c(idalpha_beta)[-c(1)]]))
}else{
thetaest2[idalpha_beta] <- exp(thetaest[c(idalpha_beta)])
}
}else{
idalpha_beta <- which(stringr::str_detect(colnames(thetaest), "alpha"))
if (n_z_alpha > 0) {
thetaest2[idalpha_beta] <- exp(c(thetaest[c(idalpha_beta)][1], thetaest[c(idalpha_beta)][1] + thetaest[c(idalpha_beta)[-c(1)]]))
}else {
thetaest2[idalpha_beta] <- exp(c(thetaest[c(idalpha_beta)]))
}
}
Dthetaest2 <- thetaest2
Dthetaest2[-id_beta] <- 1
Dthetaest2[id_beta] <- exp(thetaest[id_beta])
Dthetaest2[id_tau0] <- thetaest2[id_tau0]
if (pophaz.alpha) {
Dthetaest2[id_pophaz.alpha] <- exp(thetaest[id_pophaz.alpha])
}
varcov <- try(c(Dthetaest2^2)*varcov_star,TRUE)
std_err <- try((sqrt(diag(varcov))), TRUE)
if (inherits(std_err, "try-error")) {
}else{
std_err <- try(matrix(std_err, ncol = length(std_err)), TRUE)
if (pophaz.alpha) {
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"),
'pophaz.alpha')
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
'pophaz.alpha')
}else
if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0',
'pophaz.alpha')
}
}else{
if (n_z_tau > 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau > 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0',
paste('tau', 1:n_z_tau, colnames(z_tau), sep = "_"))
}else if (n_z_tau == 0 & n_z_alpha > 0) {
colnames(std_err) <- c('alpha0',
paste('alpha', 1:n_z_alpha, colnames(z_alpha) ,sep = "_"),
'beta',
'tau0')
}else if (n_z_tau == 0 & n_z_alpha == 0) {
colnames(std_err) <- c('alpha0',
'beta',
'tau0')
}
}
colnames(thetaest) <- try(paste0(colnames(thetaest),"*"),TRUE)
std_err_star <- matrix(std_err_star, ncol = length(std_err_star))
colnames(std_err_star) <- try(paste0(colnames(thetaest),"*"), TRUE)
}
if (is.null(ncoor_des)) {
iter_coords <- 0
}else if (is.numeric(ncoor_des)) {
iter_coords <- ncoor_des
}else if (ncoor_des == "iter") {
iter_coords <- j
}
if (method_opt == "all.methods") {
return(
list(iter_coords = iter_coords,
coefficients = thetaest,
estimates = thetaest2,
loglik = loglik,
iterations = iterations,
evaluations = evaluations,
convergence = convergence,
message = message,
varcov = varcov,
varcov_star = varcov_star,
std_err = std_err,
std_err_star = std_err_star,
AIC = AIC,
optimized2 = optimized2
)
)
}else{
return(
list(iter_coords = iter_coords,
coefficients = thetaest,
estimates = thetaest2,
loglik = loglik,
iterations = iterations,
evaluations = evaluations,
convergence = convergence,
message = message,
varcov = varcov,
varcov_star = varcov_star,
std_err = std_err,
std_err_star = std_err_star,
AIC = AIC
)
)
}
}
Try the curesurv package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.