R/mhurdle.lnl.R In mhurdle: Multiple Hurdle Tobit Models

Defines functions mhurdle.lnl

```mhurdle.lnl <- function(param, X1, X2, X3, X4, y, gradient = FALSE,
fitted = FALSE, dist = NULL, corr = FALSE,
robust = TRUE){

#----------------------------------------------------------------
# 1. Declare the transformation functions for the parameters and
#----------------------------------------------------------------

# if robust, compute the transformation of the parameter and its
# gradient in order to get the structural parameter, ie, 0 is
# entered for sigma and is actually log(sigma), so that sigma =
# exp(0) = 1

if (robust){
frho <- function(x) atan(x) * 2 / pi
grho <- function(x) 2 / pi / (1 +  x ^ 2)
fmu <- function(x) exp(x)
gmu <- function(x) exp(x)
fsd <- function(x) exp(x)
gsd <- function(x) exp(x)
}
else{
frho <- function(x) x
grho <- function(x) rep(1, length(x))
fmu <- function(x) x
gmu <- function(x) 1
fsd <- function(x) x
gsd <- function(x) 1
}

#-----------------------------------------------------
# 2. Extract the coefficients and compute the relevant
# crossproducts
#-----------------------------------------------------

##.........................................................
## 2.a  Dummies for the existing equations hi and number of
##  coefficients Ki
##.........................................................

h1 <- ! is.null(X1) ;  K1 <- ifelse(h1, ncol(X1), 0)
h3 <- ! is.null(X3) ;  K3 <- ifelse(h3, ncol(X3), 0)
h4 <- ! is.null(X4) ;  K4 <- ifelse(h4, ncol(X4), 0)
K2 <- ncol(X2)

##.............................
##   2.b Equation 1 (selection)
##.............................

if (h1){
beta1 <- param[1:K1]
bX1 <- as.numeric(crossprod(t(X1), beta1))
Phi1 <- pnorm(bX1) ; phi1 <- dnorm(bX1)
Phi1 <- sanitize(Phi1, string = c("Phi1", "mhurdle.lnl"),
verbal = FALSE, replace = TRUE)
}
else{
bX1 <- Inf ; beta1 <- NULL
Phi1 <- 1 ; phi1 <- 0;
}

##.........................
##   2.c Equation 2 (level)
##.........................

beta2 <- param[(K1 + 1):(K1 + K2)]
bX2 <- as.numeric(crossprod(t(X2), beta2))

##............................
##   2.d Equation 3 (purchase)
##............................

if (h3){
beta3 <- param[(K1 + K2 + 1):(K1 + K2 + K3)]
bX3 <- as.numeric(crossprod(t(X3), beta3))
Phi3 <- pnorm(bX3) ; phi3 <- dnorm(bX3)
Phi3 <- sanitize(Phi3, string = c("Phi3", "mhurdle.lnl"),
verbal = FALSE, replace = TRUE)
}
else{
bX3 <- Inf ; beta3 <- NULL
Phi3 <- 1 ; phi3 <- 0
}

##.........................
##   2.e Standard deviation
##.........................

sd <- param[K1 + K2 + K3 + 1]
sd <- fsd(sd)

##.......................................
##   2.f Equation 4 (heterosckedasticity)
##.......................................

if (h4){
beta4 <- param[(K1 + K2 + K3 + 2):(K1 + K2 + K3 + K4 + 1)]
bX4 <- as.numeric(crossprod(t(X4), beta4))
pbX4 <- pnorm(bX4)
#PKG sigma <- sd * pbX4
}
else{
sigma <- sd ; pbX4 <- 1
}

##...............................
##   2.g Correlation coefficients
##...............................

# KR is the length of rho, either 1 (h1 or h3) or 3 (h1 and h3)
KR <- ifelse(corr, h1 + h3 + h1 * h3, 0)
if (corr){
posrho <- (K1 + K2 + K3 + K4 + 2):(K1 + K2 + K3 + K4 + 1 + KR)
rho <- param[posrho]
# In case of only one correlation coefficient, use the whole
# vector with only one non-null component
if (h1 & ! h3) rho <- c(rho, 0, 0)
if (! h1 & h3) rho <- c(0, 0, rho)
rho <- frho(rho)
if (h1 & h3){
# In case of a trivariate normal distribution, check the
# joint relation of the three coefficients
rho3 <- function(rho) rho[1] * rho[2] + c(-1, 1) *
sqrt(1 + rho[1] ^ 2 * rho[2] ^ 2 - rho[1] ^ 2 - rho[2] ^ 2)
if (rho[3] < rho3(rho)[1]) rho[3] <- rho3(rho)[1] + 1E-04
if (rho[3] > rho3(rho)[2]) rho[3] <- rho3(rho)[2] - 1E-04
}
}
else rho <- rep(0, 3)

##........................................................
##   2.h Transformation coefficient (Box-Cox and ish only)
##........................................................

if (dist %in% c("bc", "bc2", "ihs")) lambda <- param[K1 + K2 + K3 + 1 + K4 + KR + 1]
if (dist %in% c("bc", "bc2")) sgn <- sign(lambda) else sgn <- + 1

##...............................................
##   2.i Position parameter (ln and Box-Cox only)
##...............................................

if (dist %in% c("ln2", "bc2")){
posmu <- K1 + K2 + K3 + 1 + K4 + KR + ifelse(dist == "ln2", 1, 2)
mu <- fmu(param[posmu])
}

#----------------------------------------------
# 3. Compute the elements of the log likelihood
#----------------------------------------------

##................................................
##   3.a. Transformation of the dependent variable
##................................................

Ty <- switch(dist,
"ln"  = log2(y) + log(Phi3),
"ln2" = log2(y * Phi3 + mu),
"bc"  = (exp(lambda * log(y * Phi3)) - 1) / lambda,
"bc2" = (exp(lambda * log(y * Phi3 + mu)) - 1) / lambda,
"ihs" = log(lambda * y * Phi3 + sqrt(1 + (lambda  * y * Phi3) ^ 2)) / lambda,
y * Phi3
)

##................................
##   3.b Logarithm of the jacobian
##................................

lnJ <- switch(dist,
"ln"  = - log2(y),
"ln2" = log(Phi3) - log(mu + Phi3 * y),
"bc"  = (lambda - 1) * log2(y) + lambda * log(Phi3),
"bc2" = (lambda - 1) * log(Phi3 * y + mu) + log(Phi3),
"ihs" = - 0.5 * log(1 + (lambda * Phi3 * y) ^ 2) + log(Phi3),
log(Phi3)
)

##............................................
##   3.c  Residual of the consumption equation
##............................................

resid <- (Ty - bX2)
# problem with bc and lambda < 0, for y = 0, resid = -inf and
# log(dnorm(resid)) = -inf
resid[y == 0] <- 0

##..............................................................
##  3.d Computation of the (opposite of) the range of z and of z
##  for y = 0
##..............................................................

# The minimum value of z is - Inf, except for bc/bc2 with lambda >
# 0 and for tn (left truncation)
mzn <- + Inf
if (dist %in% c("bc", "bc2") && lambda > 0) mzn <- (1 / lambda + bX2) / sigma
if (dist == "tn") mzn <- bX2 / sigma

# The maximum value of z is + Inf, except for bc/bc2 with lambda <
# 0 (right truncation)
mzx <- - Inf
if (dist %in% c("bc", "bc2") && lambda < 0) mzx <- (1 / lambda + bX2) / sigma

# (T(0) - b2x2) / sigma
mz0 <- + Inf
if (dist == "bc2") mz0 <- ( - (mu ^ lambda - 1) / lambda + bX2) / sigma
if (dist == "bc" && lambda > 0)  mz0 <- (1 / lambda + bX2) / sigma
if (dist == "ln2") mz0 <- (bX2 - log(mu)) / sigma
if (dist %in% c("n", "tn")) mz0 <- bX2 / sigma

##.................................................
## 3.e Trivariate and bivariate cummulative normals
##.................................................

Pr123A <- PHI3(bX1, mz0, bX3, rho)
Pr123B <- PHI3(bX1, mzx, bX3, rho)
if (h1) arg1 <- (bX1 + rho[1] * resid / sigma) / sqrt(1 - rho[1] ^ 2) else arg1 <- Inf
if (h3) arg3 <- (bX3 + rho[3] * resid / sigma) / sqrt(1 - rho[3] ^ 2) else arg3 <- Inf
Pr13 <- PHI2(arg1, arg3,
(rho[2] - rho[1] * rho[3]) / sqrt(1 - rho[1] ^ 2) / sqrt(1 - rho[3] ^ 2)
)
Pr13\$f <- sanitize(Pr13\$f, string = c("Pr13", "mhurdle.ln"),
verbal = FALSE, replace = TRUE)
# PI is the correction of the truncature
PI <- punorm(mzn) - punorm(mzx)
##......................................
## 3.e Computation of the log-likelihhod
##......................................

Pplus <- (Pr123A\$f - Pr123B\$f) / PI
Numerator <- PI - Pr123A\$f + Pr123B\$f
lnL_null <- log(Numerator) - log(PI)
lnL_one <- log(PI - Numerator) - log(PI)
lnL_null[y != 0] <- 0
lnL_pos <-
- log(sigma) +
dnorm(resid / sigma, log = TRUE) +
log(Pr13\$f) +
lnJ - log(PI)
lnL_pos[y == 0] <- 0
lnL <- lnL_null * (y == 0) + lnL_pos * (y != 0)
if (any(is.na(lnL) | any(is.infinite(lnL)))){
lnL[is.na(lnL)] <- 0
warnings("infinite or missing values of lnLi")
}
attr(lnL, "parts") <- c(lnLNull = sum(lnL_null[y == 0]),
lnLOne  = sum(lnL_one[y != 0]),
lnLPos  = sum(lnL_pos[y != 0]))

#-------------------------------
# 4. Computation of the gradient
#-------------------------------

##......................................
##   4.a Derivatives respective to beta1
##......................................

if (h1){
lnL_beta1 <- (y == 0) * ( (- Pr123A\$d1 + Pr123B\$d1) / Numerator) +
(y != 0) * ( Pr13\$d1 / sqrt(1 - rho[1] ^ 2) / Pr13\$f)
}

##......................................
##   4.b Derivatives respective to beta2
##......................................

PI2 <-  (dnorm(mzn) - dnorm(mzx) ) / sigma
lnL_beta2 <- (y == 0) * ( (PI2 - Pr123A\$d2 / sigma + Pr123B\$d2 / sigma) / Numerator -
PI2 / PI) +
(y != 0) * ( resid / sigma ^ 2 -
( Pr13\$d1 * rho[1] / sigma / sqrt(1 - rho[1] ^ 2) +
Pr13\$d2 * rho[3] / sigma / sqrt(1 - rho[3] ^ 2) ) / Pr13\$f -
PI2 / PI)

##......................................
##   4.c Derivatives respective to beta3
##......................................

if (h3){
# derivative of T(Phi3 y) with bX3
Ty_bX3 <- switch(dist,
"ln" = mills(bX3),
"ln2" = y * phi3 / (mu + y * Phi3),
"bc" = exp(lambda * log(y * Phi3)) * mills(bX3),
"bc2" = exp((lambda - 1) * log(y * Phi3 + mu)) * phi3 * y,
"ihs" = y * phi3 / sqrt( 1 + (lambda * y * Phi3) ^ 2),
y * phi3
)
# derivative of lnJ with bX3
lnJ_bX3 <- switch(dist,
"ln" = 0,
"ln2" = mills(bX3) - y * phi3 / (mu + y * Phi3),
"bc" = lambda * mills(bX3),
"bc2" = (lambda - 1) * phi3 * y / (Phi3 * y + mu) + mills(bX3),
"ihs" = - phi3 * Phi3 * lambda ^ 2 * y ^ 2 /
(1 + (lambda * y * Phi3) ^ 2) + mills(bX3),
mills(bX3)
)

lnL_beta3 <- (y == 0) * (- Pr123A\$d3 + Pr123B\$d3) / Numerator +
(y != 0) * (- resid / sigma ^ 2 * Ty_bX3 +
( Pr13\$d1 * Ty_bX3 * rho[1] / sigma / sqrt(1 - rho[1] ^ 2) +
Pr13\$d2 * (1 + Ty_bX3 * rho[3] / sigma) /
sqrt(1 - rho[3] ^ 2) ) / Pr13\$f +
lnJ_bX3
)
}

##   4.d  Derivative respective to sigma
PI_sigma <- 0
Pr123B_sigma <- 0
Pr123A_sigma <- ifelse(is.infinite(mz0), 0, Pr123A\$d2 * (- mz0 / sigma))
if (dist == "tn" | (dist %in% c("bc", "bc2") && lambda > 0)) PI_sigma <- PI_sigma + dnorm(mzn) * (- mzn / sigma)
if (dist %in% c("bc", "bc2") && lambda < 0){
PI_sigma <- PI_sigma - dnorm(mzx) * (- mzx / sigma)
Pr123B_sigma <- Pr123B\$d2 * (- mzx / sigma)
}
lnL_sigma <- (y == 0) * ( (PI_sigma - Pr123A_sigma + Pr123B_sigma) / Numerator - PI_sigma / PI) +
(y != 0) * (- 1 / sigma + resid ^ 2 / sigma ^ 3 +
( - Pr13\$d1 * rho[1] * resid / sigma ^ 2 / sqrt(1 - rho[1] ^ 2) -
Pr13\$d2 * rho[3] * resid / sigma ^ 2 /
sqrt(1 - rho[3] ^ 2))  / Pr13\$f -
PI_sigma / PI)

##......................................
##   4.e Derivatives respective to beta4
##......................................

if (! is.null(X4)){
}

##....................................
##   4.f Derivatives respective to rho
##....................................

if (corr){
if (h1 & ! h3){
Pr123A\$dr <- cbind(Pr123A\$dr, 0, 0)
Pr123B\$dr <- cbind(Pr123B\$dr, 0, 0)
}
if (! h1 & h3){
Pr123A\$dr <- cbind(0, 0, Pr123A\$dr)
Pr123B\$dr <- cbind(0, 0, Pr123B\$dr)
}
if ( (h1 & h3) & !is.matrix(Pr123A\$dr) ) Pr123A\$dr <- cbind(0, Pr123A\$dr, 0)

if (length(Pr123B\$dr) == 1 && Pr123B\$dr == 0) Pr123B\$dr <- cbind(0, 0, 0)
lnL_rho <- c()
if (h1){
Drho12 <- (rho[1] * rho[2] - rho[3]) /
(1 - rho[1] ^ 2) ^ 1.5  / sqrt(1 - rho[3] ^ 2)
lnL_rho12 <- (y == 0) * (- Pr123A\$dr[, 1] + Pr123B\$dr[, 1]) / Numerator  +
(y != 0) * ( Pr13\$d1 * (resid / sigma + rho[1] / (1 - rho[1] ^ 2) *
(bX1 + rho[1] * resid / sigma) ) /
(1 - rho[1] ^ 2) ^ 0.5 +
Pr13\$dr * Drho12) / Pr13\$f
lnL_rho <- cbind(lnL_rho, rho12 = lnL_rho12 * gradientrho[1])
}
if (h1 & h3){
Drho13 <- 1 / sqrt(1 - rho[1] ^ 2) / sqrt(1 - rho[3] ^ 2)
lnL_rho13 <- (y == 0) * (- Pr123A\$dr[, 2] + Pr123B\$dr[, 2]) / Numerator  +
(y != 0) * ( Pr13\$dr * Drho13) / Pr13\$f
lnL_rho <- cbind(lnL_rho, rho13 = lnL_rho13 * gradientrho[2])
}
if (h3){
Drho23 <- (rho[3] * rho[2] - rho[1]) /
(1 - rho[3] ^ 2) ^ 1.5  / sqrt(1 - rho[1] ^ 2)
lnL_rho23 <- (y == 0) * (- Pr123A\$dr[, 3] + Pr123B\$dr[, 3]) / Numerator  +
(y != 0) * ( Pr13\$d2 * (resid / sigma + rho[3] / (1 - rho[3] ^ 2) *
(bX3 + rho[3] * resid / sigma) ) /
(1 - rho[3] ^ 2) ^ 0.5 +
Pr13\$dr * Drho23) / Pr13\$f
lnL_rho <- cbind(lnL_rho, rho23 = lnL_rho23 * gradientrho[3])
}
}

##......................................
##   4.g Derivative respective to lambda
##......................................

if (dist %in% c("bc", "bc2")){
lnJ_lambda <- switch(dist,
"bc"  = log2(y) + log(Phi3),
"bc2" = log(Phi3 * y + mu),
"ihs" = - lambda * y ^ 2 * Phi3 ^ 2 / (1 + (lambda * y * Phi3) ^ 2)
)
Ty_lambda <- (log(Phi3 * y + mu) * (Phi3 * y + mu) ^ lambda - Ty) / lambda
if (mu == 0) T0_lambda <- ( 1 / lambda ^ 2 * (lambda > 0) + 0 * (lambda < 0))
else T0_lambda <- (log(mu) * mu ^ lambda - (mu ^ lambda - 1) / lambda) / lambda
Tymax_lambda <- 0 * (lambda > 0) + (1 / lambda ^ 2) * (lambda < 0)
PI_lambda <- - sign(lambda) * dnorm( (bX2 + 1 / lambda) / sigma ) / (sigma * lambda ^ 2)
lnL_lambda <- vector(mode = "numeric", length = length(y))
lnL_lambda[y == 0] <- ( (PI_lambda - Pr123A\$d2 * (- T0_lambda /  sigma) + Pr123B\$d2 *
(- Tymax_lambda / sigma)) / Numerator - PI_lambda / PI)[y == 0]
lnL_lambda[y > 0] <- ( (- resid / sigma ^ 2 +
(Pr13\$d1 * rho[1] / sqrt(1 - rho[1] ^ 2) / sigma +
Pr13\$d2 * rho[3] / sqrt(1 - rho[3] ^ 2) / sigma) / Pr13\$f ) *
Ty_lambda + lnJ_lambda - PI_lambda / PI)[y > 0]
}
if (dist == "ihs"){
Ty_lambda <- (y * Phi3) / lambda / sqrt(1 + (lambda * y * Phi3) ^ 2) - Ty / lambda
lnL_lambda <- vector(mode = "numeric", length = length(y))
}

##..................................
##   4.h Derivative respective to mu
##..................................
if (dist == "bc2"){
Ty_mu <- (Phi3 * y + mu) ^ (lambda - 1)
T0_mu <- mu ^ (lambda - 1)
lnJ_mu <- (lambda - 1) / (Phi3 * y + mu)
lnL_mu <- vector(mode = "numeric", length = length(y))
lnL_mu[y == 0] <- (- Pr123A\$d2 * (- T0_mu / sigma) / Numerator)[y == 0]
lnL_mu[y > 0] <- (( - resid / sigma ^ 2 +
(Pr13\$d1 * rho[1] / sqrt(1 - rho[1] ^ 2) / sigma +
Pr13\$d2 * rho[3] / sqrt(1 - rho[3] ^ 2) / sigma) / Pr13\$f
) * Ty_mu + lnJ_mu)[y != 0]
}
if (dist == "ln2"){
Ty_mu <- 1 / (mu + Phi3 * y)
T0_mu <- 1 / mu
lnJ_mu <- - 1 / (mu + Phi3 * y)
lnL_mu <- vector(mode = "numeric", length = length(y))
lnL_mu[y == 0] <- ( - Pr123A\$d2 * (- T0_mu / sigma ) / Numerator)[y == 0]
lnL_mu[y != 0] <- (( - resid / sigma ^ 2 +
(Pr13\$d1 * rho[1] / sqrt(1 - rho[1] ^ 2) / sigma +
Pr13\$d2 * rho[3] / sqrt(1 - rho[3] ^ 2) / sigma) / Pr13\$f
) * Ty_mu + lnJ_mu)[y != 0]
}
# for ln2, when mu = 0, lnL_mu is 0 / 0 = 0 for y = 0, the NaN
# is correctly replaced by 0
}

#------------------------------------
# 5. Computation of the fitted values
#------------------------------------

if (fitted){
if (dist %in% c("bc", "bc2")){
if (length(mzn) == 1) mzn <- rep(mzn, length(y))
if (length(mzx) == 1) mzx <- rep(mzx, length(y))
if (length(Phi3) == 1) Phi3 <- rep(Phi3, length(y))
if (length(Phi1) == 1) Phi1 <- rep(Phi1, length(y))

resid <- function(z, index){
switch(dist,
"ln"  = log2(z) + log(Phi3[index]),
"ln2" = log2(z * Phi3[index] + mu),
"bc"  = (exp(lambda * log(z * Phi3[index])) - 1) / lambda,
"bc2" = (exp(lambda * log(z * Phi3[index] + mu)) - 1) / lambda,
"ihs" = log(lambda * z * Phi3[index] + sqrt(1 + (lambda  * z * Phi3[index]) ^ 2)) / lambda,
z * Phi3[index]
) - bX2[index]
}
lnJ <- function(z, index){
switch(dist,
"ln"  = - log2(z),
"ln2" = log(Phi3[index]) - log(mu + Phi3[index] * z),
"bc"  = (lambda - 1) * log2(z) + lambda * log(Phi3[index]),
"bc2" = (lambda - 1) * log(Phi3[index] * z + mu) + log(Phi3[index]),
"ihs" = - 0.5 * log(1 + (lambda * Phi3[index] * z) ^ 2) + log(Phi3[index]),
log(Phi3[index])
)
}

arg1 <- function(z, index) if (h1) (bX1[index] + rho[1] * resid(z, index) / sigma) / sqrt(1 - rho[1] ^ 2) else Inf
arg3 <- function(z, index) if (h3) (bX3[index] + rho[3] * resid(z, index) / sigma) / sqrt(1 - rho[3] ^ 2) else Inf
Pr13 <- function(z, index) PHI2(arg1(z, index), arg3(z, index),
(rho[2] - rho[1] * rho[3]) / sqrt(1 - rho[1] ^ 2) / sqrt(1 - rho[3] ^ 2))\$f
# Pbnorm plante et pas PHI2
PI <- function(index) punorm(mzn[index]) - punorm(mzx[index])
fplus <- function(z, index){
x <- z * exp(- log(sigma) + dnorm(resid(z, index) / sigma, log = TRUE) + log(Pr13(z, index)) + lnJ(z, index) - log(PI(index)))
x
}
if (dist %in% c("bc", "bc2") && lambda < 0) maxint <- max(y) * 3
else maxint <- +Inf
E <- sapply(seq_len(length(y)), function(i) integrate(function(x) fplus(x, i), 0, maxint)\$value)
}
if (dist %in% c("ln", "ln2")){
if (h1) arg1 <- bX1 + rho[1] * sigma else arg1 <- + Inf
if (h3) arg3 <- bX3 + rho[3] * sigma else arg3 <- + Inf
E <- exp(bX2 + 0.5 * sigma ^ 2) / Phi3 * ptnorm(arg1, mz0 + sigma, arg3, rho)
if (dist == "ln2") E <- E - mu * ptnorm(bX1, mz0, bX3, rho)/ Phi3
}
if (dist %in% c("n", "tn")){
phi2 <- dnorm(bX2 / sigma)
phi1 <- dnorm(bX1)
phi3 <- dnorm(bX3)
Pr13 <- pbnorm((bX1 - rho[1] * bX2 / sigma) / sqrt(1 - rho[1] ^ 2),
(bX3 - rho[3] * bX2 / sigma) / sqrt(1 - rho[3] ^ 2),
(rho[2] - rho[1] * rho[3]) / sqrt( (1 - rho[1] ^ 2) * (1 - rho[3] ^ 2)))
Pr23 <- pbnorm((bX2 / sigma - rho[1] * bX1) / sqrt(1 - rho[1] ^ 2),
(bX3         - rho[2] * bX1) / sqrt(1 - rho[2] ^ 2),
(rho[3] - rho[1] * rho[2]) / sqrt( (1 - rho[1] ^ 2) * (1 - rho[2] ^ 2)))

Pr12 <- pbnorm((bX1         - rho[2] * bX3) / sqrt(1 - rho[2] ^ 2),
(bX2 / sigma - rho[3] * bX3) / sqrt(1 - rho[3] ^ 2),
(rho[1] - rho[2] * rho[3]) / sqrt(1 - rho[2] ^ 2) / sqrt(1 - rho[3] ^ 2)
)
E <- bX2 / Phi3 + sigma / (Phi3 * Pplus) * (phi2 * Pr13 +
rho[1] * phi1 * Pr23 / sqrt( (1 - rho[1] ^ 2) * (1 - rho[3] ^ 2)) +
rho[3] * phi3 * Pr12 / sqrt( (1 - rho[2] ^ 2) * (1 - rho[3] ^ 2)))

E <- bX2 / Phi3 + sigma / (Phi3 * Pplus) * (phi2 * Pr13 + rho[1] * phi1 * Pr23 + rho[3] * phi3 * Pr12)

}
if (! is.null(attr(y, "geomean"))){
E <- E * attr(y, "geomean")
}
attr(lnL, "fitted") <- cbind(E = E,
Ep = E / Pplus,
p = Pplus)
}
lnL
}
```

Try the mhurdle package in your browser

Any scripts or data that you put into this service are public.

mhurdle documentation built on Dec. 11, 2021, 9:21 a.m.