Nothing
one_equation_model <- function(X2, y, dist = NULL, start = NULL, check_gradient = FALSE, robust = TRUE){
if (dist %in% c("ln", "tn", "bc")) stop("the specified model doesn't allow zero observations")
if (dist == "n"){
start_length <- ncol(X2) + 1
if (is.null(start)) start <- c(rep(0.1, ncol(X2)), 1)
names(start) <- c(colnames(X2), "sigma")
# the tobit model, efficiently estimated using Survreg from
# the survival package
tobit.lnl <- function(param, fitted = FALSE){
beta <- param[1 : (length(param) - 1)]
sigma <- param[length(param)]
bX <- as.numeric(crossprod(t(X2), beta))
lnL <- (y == 0) * pnorm(- bX / sigma, log.p = TRUE) +
(y > 0) * (- log(sigma) + dnorm( (y - bX) / sigma, log = TRUE))
gbeta <- - (y == 0) * mills(- bX / sigma) / sigma + (y > 0) * (y - bX) / sigma ^ 2
gsigma <- (y == 0) * mills(- bX / sigma) * bX / sigma ^ 2 +
(y > 0) * (- 1 / sigma + (y - bX) ^ 2 / sigma ^ 3)
gradi <- cbind(gbeta * X2, sigma = gsigma)
attr(lnL, "gradient") <- gradi
lnL
}
if (check_gradient) return(compare_gradient(tobit.lnl, start))
tobit <- survival::survreg(survival::Surv(y, y > 0, type = "left") ~ X2 - 1, dist = "gaussian")
beta <- coef(tobit)
sigma <- tobit$scale
bX <- as.numeric(crossprod(t(X2), beta))
resid <- y - bX
z <- tobit.lnl(c(beta, sigma))
gradi <- attr(z, "gradient")
result <- list(coefficients = c(beta, sigma = sigma),
vcov = solve(crossprod(gradi)),
logLik = logLik(tobit),
gradient = gradi,
fitted = cbind(1 - pnorm(bX / sigma), bX + sigma * mills(bX / sigma)))
}
if (dist == "ln2"){
start_length <- ncol(X2) + 2
if (is.null(start)){
thelm <- lm(log(y + 1) ~ X2 - 1)
sigma <- summary(thelm)$sigma
beta <- coef(thelm)
start <- c(beta, sigma = sigma, alpha = 1)
}
else if (length(start) != start_length) stop("the starting values vector has a wrong length")
ltobit.lnl <- function(param, fitted = FALSE){
beta <- param[1 : (length(param) - 2)]
sigma <- param[length(param) - 1]
alpha <- param[length(param)]
bX <- as.numeric(crossprod(t(X2), beta))
lnL <- (y == 0) * pnorm( (log( alpha) - bX) / sigma, log.p = TRUE) +
(y > 0) * ( dnorm( (log(y + alpha) - bX) / sigma, log = TRUE) -
log(y + alpha) - log(sigma)
)
mu <- mills( (log(alpha) - bX) / sigma)
gbeta <- - (y == 0) * mu / sigma +
(y > 0) * (log(y + alpha) - bX) / sigma ^ 2
gsigma <- - (y == 0) * mu * (log(alpha) - bX) / sigma ^ 2 +
(y > 0) * ( - 1 / sigma + (log(alpha + y) - bX) ^ 2 / sigma ^ 3)
galpha <- (y == 0) * mu / (alpha * sigma) -
(y > 0) * ( 1 / (y + alpha) + ( log(y + alpha) - bX) / (sigma ^ 2 * (y + alpha)))
attr(lnL, "gradient") <- cbind(gbeta * X2, gsigma, galpha)
if (fitted){
P0 <- pnorm( (log(alpha) - bX) / sigma)
Ep <- exp(bX + 0.5 * sigma ^ 2) * pnorm( (bX - log(alpha)) / sigma + sigma) / (1 - P0)
attr(lnL, "fitted") <- cbind(P0, Ep)
}
lnL
}
if (check_gradient) return(compare_gradient(ltobit.lnl, start))
result <- maxLik::maxLik(ltobit.lnl, start = start, method = "bhhh", print.level = 0)
fitted <- attr(ltobit.lnl(coef(result), fitted = TRUE), "fitted")
result <- list(coefficients = coef(result),
vcov = vcov(result),
logLik = logLik(result),
gradient = result$gradientObs,
fitted = fitted)
}
if (dist == "bc2"){
if (is.null(start)) start <- c(rep(0.1, ncol(X2)), 0.1, 1)
names(start) <- c(colnames(X2), "alpha", "sigma")
f2 <- function(x) boxcox.lnl(x, X2, y, gradient = TRUE, truncated = FALSE, robust = FALSE)
if (check_gradient) return(compare_gradient(f2, start))
result <- boxcox.fit(X2, y, robust = robust, start = start,
check_gradient = check_gradient, truncated = FALSE)
}
coef.names <- list(h1 = NULL,
h2 = colnames(X2),
h3 = NULL,
sd = "sd",
h4 = NULL,
corr = NULL,
tr = NULL,
pos = NULL)
if (dist == "bc2") coef.names$tr <- "tr"
if (dist %in% c("bc2", "ln2")) coef.names$pos <- "pos"
result <- list(coefficients = result$coefficients,
vcov = result$vcov,
fitted.values = result$fitted,
logLik = result$logLik,
gradient = result$gradient,
model = NULL,
formula = NULL,
coef.names = coef.names,
call = NULL
)
class(result) <- c("mhurdle","maxLik")
result
return(result)
}
# Compute the estimation of hurdle models in the cases where it can be
# done using two independent estimations (a binomial logit model and a
# normal/log-normal/truncated linear model). This is relevant for
# uncorrelated models with selection
two_parts_model <- function(X1, X2, y, dist = NULL, start = NULL, check_gradient = FALSE){
if (is.null(start)) start <- c(rep(0.1, ncol(X1) + ncol(X2)), 1)
names(start) <- c(paste("h1", colnames(X1), sep = "."),
paste("h2", colnames(X2), sep = "."),
"sd")
yp <- y[y > 0]
q <- 2 * (as.numeric(y > 0)) - 1
lnl.probit <- function(x){
bX1 <- as.numeric(crossprod(x, t(X1)))
lnL <- pnorm(q * bX1, log.p = TRUE)
attr(lnL, "gradient") <- q * mills(q * bX1) * X1
attr(lnL, "hessian") <- crossprod(q ^ 2 * dmills(q * bX1) * X1, X1)
attr(lnL, "parts") <- c(lnLNull = sum(pnorm(- bX1[y == 0], log.p = TRUE)),
lnLOne = sum(pnorm( bX1[y > 0], log.p = TRUE)))
lnL
}
lnl.ln <- function(x){
sigma <- x[length(x)]
beta <- x[- length(x)]
bX2 <- as.numeric(crossprod(beta, t(X2)))
resid <- log2(y) - bX2
lnL <- (- log(sigma) - log2(y) + dnorm(resid / sigma, log = TRUE)) * (y > 0)
lnL_beta <- resid / sigma ^ 2
lnL_sigma <- resid ^ 2 / sigma ^ 3 - 1 / sigma
attr(lnL, "gradient") <- cbind(lnL_beta * X2, sigma = lnL_sigma) * (y > 0)
attr(lnL, "hessian") <- bdiag(- 1 / sigma ^ 2 * crossprod(X2[y > 0, ]),
- 3 * sum( (y - bX2)[y > 0] ^ 2) / sigma ^ 4 + 1 / sigma ^ 2)
attr(lnL, "hessian") <- rbind(cbind( - 1 / sigma ^ 2 * crossprod(X2[y > 0, ]), - 2 * apply((resid * X2)[y > 0, ], 2, sum)),
c(- 2 * apply((resid * X2)[y > 0, ], 2, sum), 1 / sigma ^ 2 - 3 * sum(resid[y > 0] ^ 2)))
# attr(lnL, "hessian") <- bdiag(- 1 / sigma ^ 2 * crossprod(X2[y > 0, ]), 1 / sigma ^ 2 - 3 * sum(resid[y > 0] ^ 2))
lnL
}
lnl.boxcox <- function(x){
boxcox.lnl(x, X2, yp, truncated = TRUE, robust = FALSE, gradient = TRUE)
}
lnl.tn <- function(x){
sigma <- x[length(x)]
beta <- x[- length(x)]
bX2 <- as.numeric(crossprod(beta, t(X2)))
resid <- y - bX2
lnL <- - log(sigma) + dnorm(resid / sigma, log = TRUE) - pnorm(bX2 / sigma, log.p = TRUE)
lnL <- lnL * (y > 0)
mu <- mills( bX2 / sigma)
dmu <- dmills( bX2 / sigma)
lnL_beta <- resid / sigma ^ 2 - mu / sigma
lnL_sigma <- - 1 / sigma + resid ^ 2 / sigma ^ 3 + mu * bX2 / sigma ^2
lnL_beta_beta <- crossprod(- (bX2 + dmu) / sigma ^ 2 * X2, X2)
lnL_beta_sigma <- - 2 * apply(resid / sigma + mu + dmu * bX2 / sigma * X2, 2, sum) / sigma ^ 2
lnL_sigma_sigma <- sum(1 - 3 * resid ^ 2 / sigma ^ 2 - dmu * bX2 / sigma ^ 2 - 2 * mu * bX2 / sigma) / sigma ^ 2
attr(lnL, "gradient") <- cbind(lnL_beta * X2, sigma = lnL_sigma) * (y > 0)
attr(lnL, "hessian") <- rbind(cbind(lnL_beta_beta * crossprod(X2), lnL_beta_sigma),
c(lnL_beta_sigma, lnL_sigma_sigma))
lnL
}
lnl.tn.olsen <- function(x, gradient = TRUE, hessian = TRUE){
theta <- x[ length(x)]
gamma <- x[- length(x)]
gX2 <- as.numeric(crossprod(gamma, t(X2)))
lnL <- dnorm(theta * y - gX2, log = TRUE) + log(theta) - pnorm(gX2, log.p = TRUE)
lnL <- - sum(lnL * (y > 0))
lnL_gamma <- - mills(gX2) + (theta * y - gX2)
lnL_theta <- 1 / theta - y * (theta * y - gX2)
lnL_gamma_gamma <- - crossprod( (1 + dmills(gX2)) * X2 * (y > 0), X2 * (y > 0))
lnL_gamma_theta <- apply(y * X2 * (y > 0), 2, sum)
lnL_theta_theta <- - sum( (1 / theta ^ 2 + y ^ 2) * (y > 0))
attr(lnL, "gradient") <- - apply(cbind(lnL_gamma * X2, theta = lnL_theta) * (y > 0), 2, sum)
attr(lnL, "hessian") <- - rbind(cbind(lnL_gamma_gamma, lnL_gamma_theta),
c(lnL_gamma_theta, lnL_theta_theta))
lnL
}
if (dist == "tn") f_second <- function(x) lnl.tn(x)
if (dist == "ln") f_second <- function(x) lnl.ln(x)
if (dist == "bc") f_second <- function(x) lnl.boxcox(x)
fs_second <- function(x) sum(as.numeric(f_second(x)))
fs_first <- function(x) sum(as.numeric(lnl.probit(x)))
lnl_two_parts <- function(x){
coefs_one <- x[1:ncol(X1)]
coefs_two <- x[(ncol(X1) + 1):(ncol(X1) + ncol(X2) + 1)]
one <- lnl.probit(coefs_one)
two <- f_second(coefs_two)
lnL <- as.numeric(one) + as.numeric(two)
grad_one <- attr(one, "gradient")
grad_two <- attr(two, "gradient")
grad <- cbind(grad_one, grad_two)
attr(lnL, "gradient") <- grad
attr(lnL, "parts") <- c(lnLNull = unname(attr(one, "parts")["lnLNull"]),
lnLOne = unname(attr(one, "parts")["lnLOne"]),
lnLPos = sum(as.numeric(two)) + unname(attr(one, "parts")["lnLOne"]))
lnL
}
if (check_gradient) return(compare_gradient(lnl_two_parts, start))
first <- glm(y != 0 ~ X1 - 1, family = binomial(link = "probit"))
coefs_first <- coef(first)
names(coefs_first) <- paste("h1.", colnames(X1), sep = "")
first <- lnl.probit(coefs_first)
L.null <- as.numeric(first)
lp <- as.numeric(X1 %*% coefs_first)
gradient_first <- attr(first, "gradient")
hessian_first <- attr(first, "hessian")
# Computation of the likelihood for positive observations for log-normal distribution
if (dist == "ln"){
second <- lm(log(yp) ~ X2[y > 0, ] - 1)
beta2 <- coef(second)
sigma <- sqrt(deviance(second) / nobs(second))
coefs_second <- c(beta2, sigma)
names(coefs_second) <- c(paste("h2.", colnames(X2), sep = ""), "sigma")
}
else{
if (dist == "tn"){
start <- rep(0.1, ncol(X2) + 1)
# la version olsen ne marche pas pour une raison inconnue
# za <- nlm(lnl.tn.olsen, start, print.level = 2L)
second <- truncreg::truncreg(yp ~ X2[y > 0, ] - 1)
}
if (dist == "bc") second <- boxcoxreg(yp ~ X2[y > 0, ] - 1, method = "bhhh", print.level = 0)
coefs_second <- coef(second)
names(coefs_second) <- c(paste("h2.", colnames(X2), sep = ""), sigma = sigma)
}
gradient_second <- attr(f_second(coefs_second), "gradient")
hessian_second <- attr(f_second(coefs_second), "hessian")
gradient <- cbind(gradient_first, gradient_second)
colnames(gradient) <- names(c(coefs_first, coefs_second))
L.pos <- logLik(second)
coef.names <- list(h1 = colnames(X1),
h2 = colnames(X2),
h3 = NULL,
sd = "sd",
h4 = NULL,
corr = NULL,
tr = NULL,
pos = NULL)
if (dist == "bc") coef.names$tr <- "tr"
logLik <- lnl_two_parts(c(coefs_first, coefs_second))
logLik <- structure(as.numeric(logLik), df = length(c(coefs_first, coefs_second)),
parts = attr(logLik, "parts"),
nobs = length(y), class = "logLik")
result <- list(coefficients = c(coefs_first, coefs_second),
fitted.values = NULL,
logLik = logLik,
gradient = gradient,
hessian = bdiag(hessian_first, hessian_second),
model = NULL,
formula = NULL,
coef.names = coef.names,
call = NULL
)
result$dpar <- - sum(diag(solve(result$hessian, crossprod(result$gradient))))
class(result) <- c("mhurdle","maxLik")
result
}
newton <- function(fun, coefs, trace, direction = c("min", "max"), ...){
if (trace){
cat("Initial values of the coefficients:\n")
}
direction <- match.arg(direction)
i <- 1
eps <- 10
while (abs(eps) > 1E-07){
f <- fun(coefs, gradient = TRUE, hessian = TRUE, ...)
g <- attr(f, "gradient")
if (is.matrix(g)) g <- apply(g, 2, sum)
h <- attr(f, "hessian")
if (direction == "max"){
f <- - f
g <- - g
h <- - h
}
lambda <- 1
newcoefs <- coefs - as.numeric(solve(h, g))
as_scalar <- function(x) sum(as.numeric(x))
while (as_scalar(- fun(newcoefs, ...)) > as_scalar(f)){
lambda <- lambda / 2
if(trace) cat(paste("function is increasing, lambda set to:", lambda, "\n"))
newcoefs <- coefs - lambda * as.numeric(solve(h, g))
}
eps <- as.numeric(crossprod(solve(h, g), g))
if (trace) cat(paste("iteration:", i, "criteria:", round(eps, 5), "\n"))
i <- i + 1
if (i > 500) stop("max iter reached")
coefs <- newcoefs
}
coefs
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.