Nothing
#' Estimate Instrumental Variable Probit model by Maximum Likelihood.
#'
#' Estimation of Probit model with one endogenous and continuous variable by Maximum Likelihood.
#'
#' @name ivpml
#' @param formula a symbolic description of the model of the form \code{y ~ x | z} where \code{y} is the binary dependent variable, \code{x} includes the exogenous and the endogenous continuous variable, and \code{z} is the complete set of instruments.
#' @param data the data of class \code{data.frame}.
#' @param messages if \code{TRUE}, then additional messages for the estimation procedure are displayed.
#' @param ... arguments passed to \code{maxLik}.
#' @param x,object an object of class \code{ivpml}.
#' @param digits the number of digits.
#' @param eigentol the standard errors are only calculated if the ratio of the smallest and largest eigenvalue of the Hessian matrix is less than \code{eigentol}. Otherwise the Hessian is treated as singular.
#' @param newdata optionally, a data frame in which to look for variables with which to predict.
#' @param type the type of prediction required. The default, \code{type = xb}, is on the linear prediction. If \code{type = pr}, the predicted probabilities of a positive outcome is returned. Finally, if \code{type = stdp} the standard errors of the linear predictions for each individual is returned.
#' @param asf if \code{TRUE}, the average structural function is used. This option is not allowed with \code{xb} or \code{stdp}.
#' @author Mauricio Sarrias.
#' @details
#'
#' The IV probit for cross-sectional data has the following structure:
#'
#' \deqn{
#' y_{1i}^* = x_i^\top\beta + \gamma y_{2i}+ \epsilon_i,
#' }
#' with
#' \deqn{
#' y_{2i} = z_i^\top\delta + \upsilon_i,
#' }
#' where \eqn{y_{1i}^*} is the latent (unobserved) dependent variable for individual \eqn{i = 1,...,N};
#' \eqn{y_{2i}} is the endogenous continuous variable; \eqn{z_i} is the vector of exogenous variables
#' which also includes the instruments for \eqn{y_{2i}};
#' and \eqn{(\epsilon, \upsilon)} are normal jointly distributed.
#'
#' The model is estimated using the \code{maxLik} function from \code{\link[maxLik]{maxLik}} package using
#' analytic gradient.
#'
#' @references
#' Greene, W. H. (2012). Econometric Analysis. 7 edition. Prentice Hall.
#' @examples
#' \donttest{
#' # Data
#' library("AER")
#' data("PSID1976")
#' PSID1976$lfp <- as.numeric(PSID1976$participation == "yes")
#' PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0,
#' levels = c(FALSE, TRUE),
#' labels = c("no", "yes")))
#'
#' # IV probit model by MLE
#' # (nwincome is endogenous and heducation is the additional instrument)
#' PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000)
#' fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age +
#' youngkids + oldkids + nwincome |
#' education + experience + I(experience^2) + age +
#' youngkids + oldkids + heducation,
#' data = PSID1976)
#' summary(fiml.probit)
#' }
#' @import Formula maxLik stats
#' @keywords models
#' @export
ivpml <- function(formula, data, messages = TRUE, ...){
callT <- match.call(expand.dots = TRUE)
callF <- match.call(expand.dots = FALSE)
nframe <- length(sys.calls())
# ============================-
# 1. Model frame ----
# ============================-
mf <- callT
m <- match(c("formula", "data", "subset", "na.action"), names(mf), 0)
mf <- mf[c(1, m)]
f <- Formula(formula)
if (length(f)[2L] < 2L) stop("No instruments in the second part of formula")
mf$formula <- f
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
# ============================-
# 2. Variables----
# ============================-
y1 <- model.response(mf)
y.var <- f[[2]]
X <- model.matrix(f, data = mf, rhs = 1)
Z <- model.matrix(f, data = mf, rhs = 2)
y2 <- X[, !(colnames(X) %in% colnames(Z)), drop = FALSE]
if (ncol(Z) < ncol(X)) stop("Model underidentified: check your variables in formula")
if (messages && (ncol(Z) == ncol(X))) cat("\nEstimating a just identified model....\n")
if (messages && (ncol(Z) > ncol(X))) cat("\nEstimating an overidentified model....\n")
end.var <- colnames(y2)
instruments <- colnames(Z)
if (length(end.var) > 1L) stop("ivpml only works with one endogenous variable")
# Check dependent variable
if (!all(y1 %in% c( 0, 1, TRUE, FALSE))){
stop( "all dependent variables must be either 0, 1, TRUE, or FALSE")
}
if (!is.numeric(y1)) y1 <- as.numeric(y1)
# ============================-
# 3. Initial values----
# ============================-
# These are different as those in STATA for the first equation
if (messages) cat("\nObtaining starting values from probit and linear model...\n")
probit <- glm.fit(X, y1, family = binomial("probit"))
linear <- lm.fit(Z, y2)
sigma2 <- sum(resid(linear) ^ 2) / linear$df.residual
lnsigma <- log(sqrt(sigma2))
rho <- cor(resid(linear), resid(probit))
athrho <- 0.5 * log((1 + rho)/(1 - rho))
start <- c(coef(probit), coef(linear), lnsigma, athrho)
nam_prob <- paste(y.var, colnames(X), sep = ":")
nam_lin <- paste(end.var, colnames(Z), sep = ":")
names(start) <- c(nam_prob, nam_lin, "lnsigma", "atanhrho")
# ============================-
# 4. Optimization----
# ============================-
if (is.null(callT$method)) callT$method <- 'nr'
opt <- callT
m <- match(c("print.level", "ftol", "tol", "reltol",
"gradtol", "steptol", "lambdatol", "qrtol",
"iterlim", "fixed", "activePar", "method", "control", "constraints", "gradient", "hessian"),
names(opt), 0L)
opt <- opt[c(1L, m)]
opt$start <- start
opt[[1]] <- as.name('maxLik')
opt$logLik <- as.name('lnbinary_iv')
opt[c('y1', 'y2', 'X', 'Z')] <- list(as.name('y1'), as.name('y2'), as.name('X'), as.name('Z'))
out <- eval(opt, sys.frame(which = nframe))
# ============================-
# 5. Save results----
# ============================-
out$end.var <- end.var
out$y.var <- y.var
out$instruments <- instruments
out$formula <- f
out$mf <- mf
out$call <- callT
#class(out) <- c("ivpml", "maxLik", class(out))
class(out) <- c("ivpml", class(out)) # JSS reviewer
return(out)
}
############################---
# S3 method for ivpml class ----
#############################---
#' @rdname ivpml
#' @method terms ivpml
#' @export
terms.ivpml <- function(x, ...){
formula(x$formula)
}
#' @rdname ivpml
#' @method model.matrix ivpml
#' @import stats
#' @export
model.matrix.ivpml <- function(object, ...){
f <- object$formula
mf <- object$mf
X <- model.matrix(f, data = mf, rhs = 1)
Z <- model.matrix(f, data = mf, rhs = 2)
y2 <- X[, !(colnames(X) %in% colnames(Z)), drop = FALSE]
out <- list(X = X, Z = Z, y2 = y2)
return(out)
}
#' @rdname ivpml
#' @method estfun ivpml
#' @importFrom sandwich estfun
#' @export estfun.ivpml
estfun.ivpml <- function(x, ...){
class(x) <- c("maxLik", "maxim")
estfun(x, ...)
}
#' @rdname ivpml
#' @method bread ivpml
#' @importFrom sandwich bread
#' @export bread.ivpml
bread.ivpml <- function(x, ...){
class(x) <- c("maxLik", "maxim")
bread(x, ...)
}
# #' @rdname ivpml
# #' @import stats
# #' @method AIC ivpml
# #' @export
# AIC.ivpml <- function(object, k = 2, ...){
# -2*logLik(object) + k * length(coef(object))
# }
#
# #' @rdname ivpml
# #' @import stats
# #' @method BIC ivpml
# #' @export
# BIC.ivpml <- function(object, ...){
# AIC(object, k = log(nrow(object$gradientObs)), ...)
# }
#' @rdname ivpml
#' @method vcov ivpml
#' @import stats
#' @export
vcov.ivpml <- function(object, ...){
class(object) <- c("maxLik", "maxim")
vcov(object, ...)
}
#' @rdname ivpml
#' @import stats
df.residual.ivpml <- function(object, ...){
return(row(object$gradientObs) - length(coef(object)))
}
#' @rdname ivpml
#' @export
coef.ivpml <- function(object, ...){
class(object) <- c("maxLik", "maxim")
coef(object, ...)
}
#' @rdname ivpml
#' @export
logLik.ivpml <- function(object, ...){
structure(object$maximum, df = length(coef(object)), nobs = nrow(object$gradientObs), class = "logLik")
}
#' @rdname ivpml
#' @method print ivpml
#' @import stats
#' @export
print.ivpml <- function(x, ...){
cat("Maximum Likelihood estimation\n")
cat(maximType(x), ", ", nIter(x), " iterations\n", sep = "")
cat("Return code ", returnCode(x), ": ", returnMessage(x),
"\n", sep = "")
if (!is.null(x$estimate)) {
cat("Log-Likelihood:", x$maximum)
cat(" (", sum(activePar(x)), " free parameter(s))\n",
sep = "")
cat("Estimate(s):", x$estimate, "\n")
}
}
#' @rdname ivpml
#' @method summary ivpml
#' @import stats
#' @importFrom miscTools stdEr
#' @export
summary.ivpml <- function(object, eigentol = 1e-12, ...){
result <- object$maxim
nParam <- length(coef(object))
activePar <- activePar(object)
if ((object$code < 100) & !is.null(coef(object))) {
t <- coef(object)/stdEr(object, eigentol = eigentol)
p <- 2 * pnorm(-abs(t))
t[!activePar(object)] <- NA
p[!activePar(object)] <- NA
results <- cbind(Estimate = coef(object), `Std. error` = stdEr(object,
eigentol = eigentol), `z value` = t, `Pr(> z)` = p)
}
else {
results <- NULL
}
summary <- list(maximType = object$type, iterations = object$iterations,
returnCode = object$code, returnMessage = object$message,
loglik = object$maximum, estimate = results, fixed = !activePar,
NActivePar = sum(activePar), constraints = object$constraints,
end.var = object$end.var, instruments = object$instruments)
class(summary) <- "summary.ivpml"
summary
}
#' @rdname ivpml
#' @method print summary.ivpml
#' @import stats
#' @export
print.summary.ivpml <- function(x, digits = max(3, getOption("digits") - 2),
...){
cat("--------------------------------------------\n")
cat("Maximum Likelihood estimation of IV Probit model \n")
cat(maximType(x), ", ", nIter(x), " iterations\n", sep = "")
cat("Return code ", returnCode(x), ": ", returnMessage(x),
"\n", sep = "")
if (!is.null(x$estimate)) {
cat("Log-Likelihood:", x$loglik, "\n")
cat(x$NActivePar, " free parameters\n")
cat("Estimates:\n")
printCoefmat(x$estimate, digits = digits)
cat("\nInstrumented:", x$end.var)
cat("\nInstruments:", x$instruments, "\n")
cat("Wald test of exogeneity (corr = 0): chi2", round(x$estimate["atanhrho",3]^2, 2),
"with 1 df. Prob > chi2 = ", round(pchisq(x$estimate["atanhrho",3]^2,1, lower.tail = FALSE), 4), "\n")
}
if (!is.null(x$constraints)) {
cat("\nWarning: constrained likelihood estimation.",
"Inference is probably wrong\n")
cat("Constrained optimization based on", x$constraints$type,
"\n")
if (!is.null(x$constraints$code))
cat("Return code:", x$constraints$code, "\n")
if (!is.null(x$constraints$message))
cat(x$constraints$message, "\n")
cat(x$constraints$outer.iterations, " outer iterations, barrier value",
x$constraints$barrier.value, "\n")
}
cat("--------------------------------------------\n")
}
############################-
# Effects and other functions ----
#############################-
#' @rdname ivpml
#' @method predict ivpml
#' @export
predict.ivpml <- function(object, newdata = NULL,
type = c("xb", "pr", "stdp"),
asf = TRUE,
...){
# xb: linear prediction excluding endogeneity; the default
# pr: probability of a positive outcome
# stdp: standard error of the linear prediction
# asf: average structural function; the default. This option
# is not allowed with xb or stdp
type <- match.arg(type)
mf <- if (is.null(newdata)) object$mf else newdata
f <- object$formula
X <- model.matrix(f, data = mf, rhs = 1)
Z <- model.matrix(f, data = mf, rhs = 2)
y2 <- X[, !(colnames(X) %in% colnames(Z)), drop = FALSE]
K <- ncol(X)
P <- ncol(Z)
param <- coef(object)
beta <- param[1L:K]
delta <- param[(K + 1L):(K + P)]
lnsigma <- param[(K + P + 1L)]
atanhrho <- tail(param, n = 1L)
sigma <- exp(lnsigma)
rho <- tanh(atanhrho)
index1 <- crossprod(t(X), beta)
index2 <- crossprod(t(Z), delta)
mi <- (index1 + (rho / sigma) * (y2 - index2)) / sqrt(1 - rho ^ 2)
if (type == "pr"){
if (asf){
out <- as.vector(pnorm(mi))
} else {
out <- as.vector(pnorm(index1))
}
}
if (type == "xb"){
out <- as.vector(index1)
}
if (type == "stdp"){
out <- c()
for (i in 1:nrow(X)){
out <- c(out, sqrt(X[i, , drop = F]%*% vcov(object)[1:K, 1:K] %*% t(X[i, , drop = F])))
}
}
return(out)
}
#' Get Model Summaries for use with "mtable" for objects of class ivpml
#'
#' A generic function to collect coefficients and summary statistics from a \code{ivpml} object. It is used in \code{mtable}
#'
#' @param obj a \code{ivpml} object,
#' @param alpha level of the confidence intervals,
#' @param ... further arguments,
#'
#' @details For more details see package \pkg{memisc}.
#' @import stats
#' @importFrom memisc getSummary
#' @export
getSummary.ivpml <- function(obj, alpha = 0.05, ...){
s <- summary(obj)$estimate
f <- obj$formula
mf <- obj$mf
X <- model.matrix(f, data = mf, rhs = 1)
Z <- model.matrix(f, data = mf, rhs = 2)
end.var <- obj$end.var
y.var <- obj$y.var
cf.eq1 <- s[rownames(s) %in% c(paste(y.var, colnames(X), sep = ":"), "lnsigma", "atanhrho"), ]
cf.eq2 <- s[rownames(s) %in% paste(end.var, colnames(Z), sep = ":"), ]
cval <- qnorm(1 - alpha/2)
cf.eq1 <- cbind(cf.eq1, cf.eq1[, 1] - cval * cf.eq1[, 2], cf.eq1[, 1] + cval * cf.eq1[, 2])
cf.eq2 <- cbind(cf.eq2, cf.eq2[, 1] - cval * cf.eq2[, 2], cf.eq2[, 1] + cval * cf.eq2[, 2])
rownames(cf.eq1) <- c(colnames(X), "lnsigma", "atanhrho")
rownames(cf.eq2) <- colnames(Z)
all.vars <- unique(c(colnames(X), colnames(Z), "lnsigma", "atanhrho"))
# Make Table
coef <- array(dim = c(length(all.vars), 6, 2),
dimnames = list(all.vars, c("est", "se", "stat", "p", "lwr", "upr"), c(y.var, end.var)))
coef[rownames(cf.eq1),,1] <- cf.eq1
coef[rownames(cf.eq2),,2] <- cf.eq2
# Statistics
sumstat <- c(logLik = logLik(obj), deviance = NA, AIC = AIC(obj), BIC = BIC(obj), N = nrow(obj$gradientObs),
LR = NA, df = NA, p = NA, Aldrich.Nelson = NA, McFadden = NA, Cox.Snell = NA)
list(coef = coef, sumstat = sumstat, contrasts = obj$contrasts, xlevels = obj$xlevels, call = obj$call)
}
#### Effects
#' Get Model Summaries for use with "mtable" for objects of class effect.ivpml
#'
#' A generic function to collect coefficients and summary statistics from a \code{effect.ivpml} object. It is used in \code{mtable}
#'
#' @param obj an \code{effect.ivpml} object,
#' @param alpha level of the confidence intervals,
#' @param ... further arguments,
#'
#' @details For more details see package \pkg{memisc}.
#' @import stats
#' @importFrom memisc getSummary
#' @export
getSummary.effect.ivpml <- function(obj, alpha = 0.05, ...){
cf <- summary(obj)$CoefTable
cval <- qnorm(1 - alpha/2)
coef <- cbind(cf, cf[, 1] - cval * cf[, 2], cf[, 1] + cval * cf[, 2])
dim(coef) <- c(dim(coef)[1], dim(coef)[2], 1)
dimnames(coef) <- list(rownames(cf), c("est", "se", "stat", "p", "lwr", "upr"), all.vars(obj$formula)[1])
#colnames(coef) <- c("est", "se", "stat", "p", "lwr", "upr")
# Statistics
sumstat <- c(logLik = obj$maximum, deviance = NA, AIC = NA, BIC = NA, N = nrow(obj$gradientObs),
LR = NA, df = NA, p = NA, Aldrich.Nelson = NA, McFadden = NA, Cox.Snell = NA)
list(coef = coef, sumstat = sumstat, contrasts = NULL, xlevels = NULL, call = obj$call)
}
mdydx.ivpml <- function(coeff, object, asf){
f <- object$formula
mf <- object$mf
X <- model.matrix(f, data = mf, rhs = 1)
Z <- model.matrix(f, data = mf, rhs = 2)
y2 <- X[, !(colnames(X) %in% colnames(Z)), drop = FALSE]
bhat <- coeff
sigma <- exp(bhat["lnsigma"])
rho <- tanh(bhat["atanhrho"])
r <- (rho / sigma)
sr <- sqrt(1 - rho^2)
K <- ncol(X)
P <- ncol(Z)
beta <- bhat[1L:K]
delta <- bhat[(K + 1L):(K + P)]
# Obtain classes
all.var <- all.vars(formula(f, rhs = 1))[-1L]
classes <- rep("numeric", length(all.var))
class.mf <- attributes(terms(mf))[["dataClasses"]][-1L]
classes[paste0("factor(", all.var, ")") %in% names(class.mf)] <- class.mf[names(class.mf) %in% paste0("factor(", all.var, ")")]
names(beta) <- colnames(X)
## Compute marginal effects
mes <- c()
mes.name <- c()
for (k in 1:length(all.var)){
if (classes[k] == "numeric") {
xb <- crossprod(t(X), beta)
upsilon <- y2 - crossprod(t(Z), delta)
mi <- (xb + r * upsilon) / sr
bk <- make.inter.num(all.var[k], names(beta), beta, X)
if (asf) {
me <- dnorm(mi) * bk / sr
} else {
me <- dnorm(xb) * bk
}
mes <- cbind(mes, me)
mes.name <- c(mes.name, all.var[k])
}
if (classes[k] == "factor"){
levs <- attributes(mf[, paste0("factor(", all.var[k], ")")])$levels
levs <- levs[-1L]
## Make P0
beta.temp <- beta
vb <- make.inter.factor(all.var[k], names(beta), levs)
if (any(vb$names %in% names(beta))) beta.temp[names(beta) %in% vb$names] <- 0
xb <- crossprod(t(X), beta.temp)
upsilon <- y2 - crossprod(t(Z), delta)
mi <- (xb + r * upsilon) / sr
p0 <- pnorm(mi)
for (j in 1:length(levs)){
## Make P1
Xtemp <- X
if (any(vb$names %in% names(beta))) Xtemp[, names(beta) %in% vb$names] <- 0
vbj <- make.inter.factor(all.var[k], names(beta), levs[j])
if (any(vbj$names %in% names(beta))) Xtemp[, names(beta) %in% vbj$names] <- X[, vbj$names.inte]
if (vbj$names[1] %in% names(beta)) Xtemp[, names(beta) %in% vbj$names[1]] <- 1
xbt <- crossprod(t(Xtemp), beta)
mit <- (xbt + r * upsilon) / sr
p1 <- pnorm(mit)
me <- p1 - p0
mes <- cbind(mes, me)
mes.name <- c(mes.name, paste0("factor(",all.var[k],")",levs[j], sep = ""))
}
}
}
colnames(mes) <- mes.name
mes <- colMeans(mes)
return(mes)
}
#### Likelihood function ----
#' @importFrom utils tail
lnbinary_iv <- function(param, y1, y2, X, Z, gradient = TRUE, hessian = TRUE){
# Likelihood function for binary Probit IV model with cross-sectional
# data and one continuous endogenous variable
F <- pnorm
f <- dnorm
#ff <- function(x) -x * dnorm(x)
K <- ncol(X)
P <- ncol(Z)
beta <- param[1L:K]
delta <- param[(K + 1L):(K + P)]
lnsigma <- param[(K + P + 1L)]
atanhrho <- tail(param, n = 1L)
sigma <- exp(lnsigma)
rho <- tanh(atanhrho)
index1 <- crossprod(t(X), beta)
index2 <- crossprod(t(Z), delta)
q <- 2 * y1 - 1
ai <- as.numeric(q * (index1 + (rho / sigma) * (y2 - index2)) / sqrt(1 - rho ^ 2))
P1 <- F(ai)
bi <- as.numeric((y2 - index2) / sigma)
P2 <- (1 / sigma) * f(bi)
Pi <- pmax(P1 * P2, .Machine$double.eps)
Li <- log(Pi)
if (gradient){
m <- function(x) f(x) / F(x)
gb <- X * (m(ai) * q / sqrt(1 - rho^2))
gd <- Z * (- m(ai) * q * (rho / sigma) / sqrt(1 - rho^2) + (bi / sigma))
glnsigma <- -(m(ai) * q * rho / sqrt(1 - rho ^ 2)) * bi + bi^2 - 1
gathrho <- m(ai) * q * (index1 * rho + bi) / sqrt(sech(atanhrho)^2)
attr(Li,'gradient') <- cbind(gb, gd, glnsigma, gathrho)
}
if (hessian){
h <- function(x) - x * m(x) - m(x)^2
H <- matrix(0, K + P + 2, K + P + 2)
cqr <- q / sqrt(1 - rho^2)
bb <- h(ai) * cqr^2
H[1:K, 1:K] <- crossprod(bb * X, X)
bd <- - h(ai) * cqr^2 * (rho / sigma)
H[1:K, (K + 1):(K + P)] <- crossprod(bd * X, Z)
bls <- -h(ai) * cqr^2 * rho * bi
H[1:K, (K + P + 1)] <- colSums(bls * X)
bt <- h(ai) * cqr * (q * (index1 * rho + bi) / sqrt(sech(atanhrho)^2))
H[1:K, (K + P + 2)] <- colSums(as.numeric(bt) * X)
H[(K + 1):(K + P), 1:K] <- t(H[1:K, (K + 1):(K + P)])
cqr2 <- q * (rho / sigma) / sqrt(1 - rho^2)
dd <- h(ai) * cqr2^2 - (1 / sigma^2)
H[(K + 1):(K + P), (K + 1):(K + P)] <- crossprod(dd * Z, Z)
dls <- (bi / sigma)* (h(ai)* (q * rho / sqrt(1 - rho^2))^2 -2)
H[(K + 1):(K + P), (K + P + 1)] <- colSums(dls * Z)
dt <- -h(ai) * cqr2 * (q * (index1 * rho + bi) / sqrt(sech(atanhrho)^2))
H[(K + 1):(K + P), (K + P + 2)] <- colSums(as.numeric(dt) * Z)
H[(K + P + 1), 1:K]<- t(H[1:K, (K + P + 1)])
H[(K + P + 1), (K + 1):(K + P)] <- t(H[(K + 1):(K + P), (K + P + 1)])
H[K + P + 1, K + P + 1] <- sum(h(ai) * (q * rho * bi / sqrt(1 - rho^2))^2 + m(ai) * (q * rho * bi / sqrt(1 - rho^2)) - 2 * bi^2)
H[K + P + 1, K + P + 2] <- sum(- bi * (h(ai) * (q * rho / sqrt(1 - rho^2)) * (q * (index1 * rho + bi) / sqrt(sech(atanhrho)^2)) + m(ai) * (q / sqrt(sech(atanhrho)^2))))
H[K + P + 2, 1:K] <- t(H[1:K, (K + P + 2)])
H[K + P + 2, (K + 1):(K + P)] <- t(H[(K + 1):(K + P), (K + P + 2)])
H[K + P + 2, K + P + 1] <- H[K + P + 1, K + P + 2]
H[K + P + 2, K + P + 2] <- sum(h(ai) * (q * (index1 * rho + bi) / sqrt(sech(atanhrho)^2))^2 + q * m(ai) * (index1 + bi * rho)/sqrt(sech(atanhrho)^2))
attr(Li, 'hessian') <- H
}
return(Li)
}
sech <- function(x) 1 / cosh(x)
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.