#' Plots the observed, predict and an independent variable from a
#' regression model.
#'
#' @param Regression.object The model.
#' @param predictor.number The index of the predictor to be used in the plot,
#' where 1 means the first predictor.
#' @details With numeric predicts, other variables are assumed to have an
#' average effect. With factors, the first level is from the intercept, and
#' the other levels are the intercept plus the parameter (thus, if you modify
#' the contrats, this plot will not display the predicted values correctly).
#' @importFrom flipFormat Equation
#' @importFrom graphics abline segments
#' @export
PredictionPlot <- function(Regression.object, predictor.number = NULL)
{
if (!inherits(Regression.object, "Regression"))
stop("'PredictionPlot' requires a regression object.")
# Working out which variable to use.
coefs <- coef(Regression.object)
coefs.n <- length(coefs)
if (is.null(predictor.number))
{
if (coefs.n > 2)
warning("Only the first predictor variable is shown on this plot")
predictor.number <- 1
}
else if (!predictor.number %in% coefs.n)
warningPredictorVariableDoesNotExist()
# Identifying the variable to use
variable.number <- predictor.number + 1
y.name <- names(Regression.object$model)[1]
x.name <- names(Regression.object$model)[variable.number]
x <- Regression.object$model[[variable.number]]#model.matrix(Regression.object)[, 2]
y <- Regression.object$model[[1]]
# Creating the plot.
equation <- Equation(Regression.object)
plt <- plot(x, y, xlab = x.name, ylab = y.name, main = equation)
# Plotting the line of best fit.
if (is.factor(x))
{
# Identifying the parameter to use
n.parameters <- unlist(lapply(x, nlevels))
factors <- n.parameters != 0
n.parameters[!factors] <- 1
n.parameters[factors] <- n.parameters[factors] - 1
cum.parameters <- cumsum(n.parameters)
par.indices <- c(1, (cum.parameters[predictor.number] + 1):cum.parameters[predictor.number + 1])
pars <- coefs[par.indices]
pars[-1] <- pars[-1] + pars[1]
n.pars <- n.parameters[variable.number] + 1
for (i in 1:n.pars)
segments(i - 0.5, pars[i], i + 0.5, pars[i], col = "red", lwd = "2", lty = 2)
return(plt)
}
b <- coefs[[x.name]]
x.mean <- mean(x)
y.mean <- mean(y)
a <- y.mean - b * x.mean
plt <- abline(a, b, col = "red", lwd = 2, lty = 2)
plt
}
#' Produce an effects plot of a regression model
#'
#' @param model A model of class \code{Regression}
#' @param max.factor.label.length The maximum length of any factor label, as plotted on the
#' x-axis for categorical variables. Longer labels are truncated, possibly with an integer
#' suffix to ensure uniquueness.
#' @param y.axis.title The title of the y-axis for each plot. Defaults to the outcome label.
#' @export
EffectsPlot <- function(model,
max.factor.label.length = NULL,
y.axis.title = NULL)
{
if (any(aliased.var <- model$summary$aliased))
model <- updateAliasedModel(model)
## model$formula <- formula(model$original)
## model$summary$aliased <- model$summary$aliased[!model$summary$aliased]
effects <- allEffects(model)
relevant.coefs <- !isDummyVariable(names(effects))
effects <- effects[relevant.coefs]
limitEffectsLabels <- function(ef, max.len){
if (ef$variables[[1]]$is.factor)
{
old.levs <- ef$variables[[1]]$levels
new.levs <- sapply(old.levs, function(x) substr(x, 1, min(max.len, nchar(x))))
new.levs <- make.unique(new.levs, sep = "")
indices <- match(levels(ef$x[[1]]), old.levs)
levels(ef$x[[1]]) <- new.levs[indices]
ef$variables[[1]]$levels <- new.levs
return(ef)
}
else
return(ef)
}
ylab <- y.axis.title
if (is.null(ylab) || nchar(ylab) == 0)
{
ylab <- model$outcome.label
if (model$type == "Binary Logit")
ylab <- paste0("probability(", ylab, ")")
}
y.axis <- list(lab = ylab)
if (model$type == "Binary Logit")
y.axis$lim <- c(0, 1)
if (!is.null(max.factor.label.length))
{
effects <- lapply(effects, limitEffectsLabels, max.factor.label.length)
class(effects) <- "efflist"
}
type <- if (model$type %in% c("Ordered Logit", "Multinomial Logit"))
"probability"
else
"response"
# Top 4 most important, ordered by increasing p-value from anova
anova.obj <- model$anova[relevant.coefs, ]
p.values <- anova.obj[, ncol(anova.obj)]
p.values <- p.values[!rownames(anova.obj) == "Residuals"]
effects <- effects[order(p.values)]
## effects <- effects[1:min(4, length(effects))]
plot(effects,
type = type,
axes = list(x = list(rug = FALSE),
y = y.axis))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.