Nothing
##' Likelihood profiles for fitted strategic models
##'
##' Calculate profile likelihood to assess convergence of a model.
##'
##' Likelihood profiling can help determine if a model fit failed to reach a
##' global maximum, which can be an issue (especially for the
##' \code{\link{ultimatum}} model). The process of profiling is as follows: a
##' parameter selected to be profiled is fixed at certain values spaced around
##' its originally estimated value, while the log-likelihood is maximized with
##' respect to the other parameters in the model. For models with large numbers
##' of observations or parameters, profiling may take a long time, as \eqn{p
##' \times (2s + 1)}{p * (2s + 1)} models will be fit (p: number of parameters;
##' s: number of steps).
##'
##' The function will issue a warning if a model fit in profiling has a
##' log-likelihood exceeding that of the original model. This means the
##' original fit failed to reach a global maximum, and any inferences based on
##' the fitted model are invalid. If this occurs, refit the model, passing the
##' \code{profile.game} output to the fitting function's \code{profile} argument
##' (as in the example below). The new fit will use the coefficients from the
##' profile fit with the highest log-likelihood as starting values.
##'
##' The function is based loosely on \code{\link{profile.glm}} in the \pkg{MASS}
##' package. However, that function focuses on the calculation of exact
##' confidence intervals for regression coefficients, whereas this one is for
##' diagnosing non-convergence. Future versions of the \pkg{games} package may
##' incorporate the confidence interval functionality as well.
##' @param fitted a fitted model of class \code{game}.
##' @param which integer vector giving the indices of the parameters to be
##' profiled. The default is to use all parameters. Parameters that were held
##' fixed in the original fitting are ignored if selected.
##' @param steps number of steps to take (in each direction) from the original
##' value for each parameter to be profiled.
##' @param dist distance the last step should be from the original parameter
##' value (in terms of standard errors if \code{use.se} is \code{TRUE};
##' absolute value otherwise). Should be a numeric vector of length equal to 1
##' or \code{length(coef(fitted))}.
##' @param use.se logical: whether \code{dist} refers to standard errors
##' @param report logical: whether to print status bar (for complex models or
##' those with many parameters, profiling can be lengthy)
##' @param ... other arguments to be passed to the fitting function (see
##'ode{\link{maxLik}}).
##' @return A list of data frames, each containing the estimated coefficients
##' across the profiled values for a particular parameter. The first column of
##' each data frame is the log-likelihood for the given fits. The returned
##' object is of class \code{c("profile.game", "profile")}.
##' @seealso \code{\link{plot.profile.game}} for plotting profiled likelihoods
##' @export
##' @author Brenton Kenkel (\email{brenton.kenkel@@gmail.com})
##' @example inst/examples/profile.game.r
profile.game <- function(fitted, which = 1:p, steps = 5, dist = 3, use.se =
TRUE, report = TRUE, ...)
{
## get the regressors from the original model
mf <- match(c("subset", "na.action"), names(fitted$call), 0L)
mf <- fitted$call[c(1L, mf)]
mf$formula <- fitted$formulas
mf$data <- fitted$model
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
regr <- list()
for (i in seq_len(length(fitted$formulas)[2]))
regr[[i]] <- model.matrix(fitted$formulas, data = mf, rhs = i)
## other info needed to refit the original model: coefficients (for starting
## values), fixed parameters, outcome, link, and type. also will retrieve
## the standard errors to determine values to profile at
cf <- fitted$coefficients
p <- length(cf)
if (length(dim(fitted$y))) { ## for ultimatum objects (if offer and accept
## both included)
y <- fitted$y[, 1]
} else {
y <- as.numeric(fitted$y) ## use as.numeric() because y is stored as a
## factor in a game object
}
fixed <- fitted$fixed
link <- fitted$link
type <- fitted$type
method <- fitted$convergence$method
if (length(dist) == 1L)
dist <- rep(dist, length(cf))
if (use.se) {
se <- sqrt(diag(fitted$vcov))
} else {
se <- rep(1, length(cf))
}
## special stuff for ultimatum models
if (inherits(fitted, "ultimatum")) {
minOffer <- fitted$minOffer
maxOffer <- fitted$maxOffer
acc <- if (length(dim(fitted$y))) fitted$y[, 2] else NULL
offerOnly <- fitted$outcome == "offer"
offertol <- fitted$offertol
names(regr) <- c("X", "Z")
} else {
minOffer <- maxOffer <- acc <- offerOnly <- offertol <- NULL
}
## ## Special stuff for egame12cor models
## if (inherits(fitted, "egame12cor")) {
## gridsize <- fitted$gridsize
## } else {
## gridsize <- NULL
## }
## retrieve the appropriate log-likelihood and gradient, based on the type
## of model. this relies on the class attribute of the object being ordered
## correctly, which will generally be the case if users don't mess with the
## objects created, but in the future it might be better to use something
## like the "family" attribute of glm objects
logLik <- switch(class(fitted)[2],
egame12 = logLik12,
## egame12cor = logLik12cor,
egame122 = logLik122,
egame123 = logLik123,
ultimatum = logLikUlt)
logLikGrad <- switch(class(fitted)[2],
egame12 = logLikGrad12,
## egame12cor = NULL,
egame122 = logLikGrad122,
egame123 = logLikGrad123,
ultimatum = logLikGradUlt)
if (!fitted$convergence$gradient)
logLikGrad <- NULL
## looping over each parameter value
##
## the code in the loop is loosely based on that of MASS:::profile.glm,
## written by D.M. Bates and W.N. Venables, licensed under the GPL. the
## difference is that theirs uses a nice adaptive algorithm to go
## approximately the right distance in the specified number of steps,
## whereas mine is dumb and depends entirely on the supplied values of
## 'steps' and 'dist'
didNotConverge <- FALSE
ans <- list()
k <- 0
if (report) {
cat("\nEstimating likelihood profiles...\n")
pb <- txtProgressBar(min = 1, max = length(which) * (2*steps + 1))
}
for (i in which) {
if (fixed[i]) # skip fixed parameters
next
## calculate the new parameter values to maximize at
thisAns <- matrix(nrow = 2*steps, ncol = length(cf) + 1)
thisAns <- data.frame(thisAns)
names(thisAns) <- c("logLik", names(cf))
cfvals <- seq(cf[i] - dist[i]*se[i], cf[i] + dist[i]*se[i], length.out =
2*steps + 1)
cfvals <- cfvals[-(steps + 1)] # prevent refitting from original
# coefficients (causes the
# non-convergence warning to be issued
# inappropriately)
fvec <- fixed
fvec[i] <- TRUE
## inner loop: refitting while fixing parameter i at the j'th element of
## 'cfvals'
for (j in seq_along(cfvals)) {
sval <- cf
sval[i] <- cfvals[j]
results <-
tryCatch(maxLik(logLik = logLik, grad = logLikGrad, start =
sval, fixed = fvec, method = method, y = y, regr
= regr, link = link, type = type, acc = acc,
minOffer=minOffer, maxOffer = maxOffer, offerOnly = offerOnly,
offertol = offertol, ...),
error = identity)
if (inherits(results, "error")) {
thisAns[j, ] <- NA
} else {
thisAns[j, ] <- c(results$max, results$estimate)
}
k <- k + 1
setTxtProgressBar(pb, k)
}
## check if any of the profiled fits yield a log-likelihood higher than
## that of the the original fitted model
if (isTRUE(any(thisAns[, 1] > sum(fitted$log.likelihood))))
didNotConverge <- TRUE
ans[[names(cf)[i]]] <- thisAns
}
if (report)
cat("\n")
if (didNotConverge)
warning("some profiled fits have higher log-likelihood than original fit; refit the model using \"profile\" option")
attr(ans, "original.fit") <- fitted
class(ans) <- c("profile.game", "profile")
return(ans)
}
##' Plot profiles of strategic model log-likelihoods
##'
##' Plot output of \code{\link{profile.game}}.
##'
##' This method provides plots for a quick assessment of whether \code{game}
##' models have failed to converge to a global maximum. For each parameter
##' profiled (see \code{\link{profile.game}} for details of the profiling
##' process), a spline interpolation of the log-likelihood profile is provided,
##' with an "x" marking the value at the original parameter estimate.
##'
##' Sometimes the plot will seem to indicate that the original fit did not reach
##' the global maximum, even though \code{profile.game} did not issue the
##' non-covergence warning. This is an artifact of the interpolation, which can
##' be confirmed by re-running \code{plot.profile.game} with the argument
##' \code{show.pts = TRUE}.
##' @param x an object of class \code{profile.game}, typically created by
##' running \code{\link{profile.game}} on a fitted \code{game} model
##' @param show.pts logical: plot a point for the log-likelihood of each
##' profiled model?
##' @param ... other arguments, currently ignored
##' @return \code{x}, invisibly
##' @method plot profile.game
##' @seealso \code{\link{profile.game}}
##' @export
##' @author Brenton Kenkel
plot.profile.game <- function(x, show.pts = FALSE, ...)
{
## this function is based largely on MASS:::plot.profile, written by
## D.M. Bates and W.N. Venables, licensed under GPL
nm <- names(x)
origcf <- coef(attr(x, "original.fit"))
origll <- sum(attr(x, "original.fit")$log.likelihood)
## saves the user's original graphical parameters (which are restored after
## the new plot is made) and sets new ones to display all plots at once
nr <- ceiling(sqrt(length(nm)))
oldpar <- par(mfrow = c(nr, nr))
on.exit(par(oldpar))
for (nam in nm) {
xvals <- x[[nam]][, nam]
yvals <- x[[nam]][, "logLik"]
## make sure interpolation passes through original fit
xvals <- c(xvals, origcf[nam])
yvals <- c(yvals, origll)
splineVals <- spline(xvals, yvals)
ylim <- c(min(yvals, splineVals$y), max(yvals, splineVals$y))
plot(xvals, yvals, xlab = nam, ylab = "log-likelihood", type = "n",
ylim = ylim)
## plot an "x" at the main model estimate
points(xvals[length(xvals)], yvals[length(yvals)], pch = 4)
## plot other points (if desired)
if (show.pts) {
points(xvals[-length(xvals)], yvals[-length(xvals)])
}
## spline fit for the likelihood curve
lines(splineVals$x, splineVals$y)
}
invisible(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.