Nothing
mnl.loglikelihood <- function(x, y, coefficients) {
U <- mnl.utilities(x, coefficients)
P <- mnl.probabilities(U)
sum(y*log(P))
}
mnl.utilities <- function(x, coefficients, npars=dim(x)[2], nobs=dim(x)[1], nalts=dim(x)[3]) {
pars3d <- aperm(array(coefficients, c(npars, nobs, nalts)), c(2,1,3))
U <- x*pars3d
apply(U, c(1,3), sum)
}
mnl.probabilities <- function(U) {
eU <- exp(U)
sum.eU <- apply(eU, 1, sum)
eU/sum.eU
}
do.estimate.mlogit <- function(x, y, specification=NULL, method='BHHH', ...) {
# x is in 3d format
# y is either a vector of indices of the chosen alternatives, or a 0-1 matrix (nobs x nalts)
P <- 0
mnl.logLik <- function(param) {
U <- mnl.utilities(x, param, npars, nobs, nalts)
P <<- mnl.probabilities(U)
sum(y*log(P))
}
mnl.logLik.sum <- function(param) {
sum(mnl.logLik(param))
}
mnl.grad <- function(param) {
d <- y - P
apply(aperm(array(d, c(nobs,nalts,npars)), c(1,3,2))*x, c(1,2), sum)
}
mnl.hess <- function(param) {
g <- mnl.grad(param)
-t(g)%*%g
}
npars <- dim(x)[2]
nobs <- dim(x)[1]
nalts <- dim(x)[3]
params <- rep(0, npars)
ll0 <- mnl.logLik(params)
s <- Sys.time()
maxlik <- maxLik(mnl.logLik, grad=mnl.grad, hess=mnl.hess,
start=params,
method=method, ...
)
time <- Sys.time() - s
coefs <- maxlik$estimate
hessian <- maxlik$hessian
names(coefs) <- rownames(hessian) <- colnames(hessian) <- dimnames(x)[[2]]
U <- mnl.utilities(x, coefs, npars, nobs, nalts)
P <- mnl.probabilities(U)
fitted.values <- P
residuals <- y - fitted.values
g <- maxlik$gradient
conv <- (t(g)%*%solve(-hessian))%*%g
result <- list(coefficients=coefs, logLik=maxlik$maximum, logLik0=ll0,
aic=-2 * maxlik$maximum + 2 * npars,
bic=-2 * maxlik$maximum + npars * log(nobs),
iter=maxlik$iterations, hessian=hessian,
gradient=g,
fitted.values=fitted.values,
residuals=residuals,
specification = specification,
convergence = conv,
method=method,
time=time,
code = maxlik$code,
message=maxlik$message,
last.step=maxlik$last.spep
)
attr(result, 'class') <- 'mnl'
invisible(result)
}
"estimate.mlogit" <- function (object, ...) UseMethod("estimate.mlogit")
"estimate.mlogit.formula" <- function(f, data, method='BHHH', choices=NULL, base.choice=1,
varying = NULL, sep='.', ...) {
mnlspec <- mnl.spec(f, data=data, choices=choices, base.choice=base.choice,
varying=varying, sep=sep)
estimate.mlogit(mnlspec, data, method=method, ...)
}
"estimate.mlogit.mnl.spec" <- function(object, data, method='BHHH', ...) {
mnl.data <- get.mnl.data(data, spec=object)
invisible(do.estimate.mlogit(mnl.data$x, mnl.data$y, specification=object, method=method, ...))
}
"estimate.mlogit.bic.mlogit" <- function(object, ...) {
estimate.mlogit(object$bma.specifications, ...)
}
"estimate.mlogit.list" <- function(object, data, verbose=TRUE, ...) {
result <- list()
for (ispec in 1:length(object)){
if(verbose) cat('\n\nEstimating specification #', ispec,
'\n===============================\n')
result[[ispec]] <- estimate.mlogit(object[[ispec]], data, ...)
if(verbose) print(summary(result[[ispec]]))
}
invisible(result)
}
get.mnl.data <- function(x, y=NULL, spec, intercept.prefix='asc.') {
# Returns a list with components:
# x - 3d array (nobs x nvars x nchoices)
# y - 0-1 response matrix (nobs x nchoices)
data <- list()
data$y <- y
nchoices <- length(spec$choices)
nx <- nrow(x)
if(is.null(y)) data$y <- x[,spec$response]
unique.y <- unique(data$y)
if(!all(is.element(unique.y, 1:nchoices))) { # replace values of y by its index within choices
ynew <- rep(NA, length(data$y))
for(ichoice in 1:nchoices)
ynew[data$y==spec$choices[ichoice]] <- ichoice
keep.idx <- !is.na(ynew)
data$y <- ynew[keep.idx]
x <- x[keep.idx,]
nx <- nrow(x)
}
# create a response matrix from a vector
ym <- matrix(0, ncol=nchoices, nrow=nx)
index <- matrix(c(1:nx,data$y), ncol=2)
ym[index] <- 1
data$y <- ym
other.choices <- spec$choices[-spec$base.choice]
other.choices.idx <- (1:nchoices)[-spec$base.choice]
lvars <- sum(spec$same.coefs) + sum(spec$variable.used[,!spec$same.coefs])
col.names <- colnames(x)
data$x <- array(NA, c(nx, lvars, nchoices))
var.count <- 0
x.col.names <- c()
for (var.sh in colnames(spec$variable.used)) {
if (spec$same.coefs[var.sh]) { #variables with the same coefficient
var.count <- var.count + 1
if (is.element(var.sh, col.names)){ # is the variable non-alternative specific name in the data?
data$x[,var.count,] <- x[,var.sh]
} else { # get alternative specific name from the data
which.choice <- (1:nchoices)[is.element(spec$full.var.names[,var.sh], col.names)]
alt.spec.vars <- spec$full.var.names[which.choice, var.sh]
d <- x[,alt.spec.vars]
if(is.element(spec$base.choice, which.choice)) {
# generate differences
d <- d - x[,spec$full.var.names[spec$base.choice, var.sh]]
}
data$x[,var.count,which.choice] <- as.matrix(d)
}
x.col.names <- c(x.col.names, var.sh)
} else { # different coefficients for different alternatives
var.full <- spec$full.var.names[,var.sh]
var.base <- spec$full.var.names[spec$base.choice, var.sh]
get.difference <- FALSE
if (is.element(var.base, col.names)) get.difference <- TRUE
for (ichoice in other.choices.idx) {
if(spec$variable.used[ichoice,var.sh]) {
var.count <- var.count + 1
if (is.element(var.full[ichoice], col.names)) {#is the variable alt.-specific name in the data?
# generate differences
d <- x[,var.full[ichoice]]
if (get.difference) d <- d - x[,var.base]
data$x[,var.count,ichoice] <- d
} else data$x[,var.count,ichoice] <- x[,var.sh]
data$x[,var.count,-ichoice] <- 0
x.col.names <- c(x.col.names, var.full[ichoice])
}
}
}
}
if(any(spec$intercepts)) {
int.array <- array(0, c(nx, sum(spec$intercepts), nchoices))
dimnames(int.array)[[3]] <- spec$choices
j<-1
for (ichoice in (1:nchoices)[spec$intercepts]) {
int.array[,j,(1:nchoices)[ichoice]] <- 1
j <- j+1
}
data$x <- abind(int.array, data$x, along=2)
x.col.names <- c(paste(intercept.prefix, spec$choices[spec$intercepts], sep=''), x.col.names)
}
data$x[,,spec$base.choice] <- 0
dimnames(data$x) <- list(1:nx, x.col.names, dimnames(data$x)[[3]])
return(data)
}
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.