Nothing
#' Plot effect stars for vglm objects.
#'
#' @aliases effectstars.vglm
#' @export
#' @export effectstars.vglm
#' @importFrom VGAM coef
#' @importFrom miscTools insertCol
#' @description Plots effect stars for \code{\link[VGAM]{vglm}}-objects. In particular, the method
#' works for multinomial logit models created by family \code{\link[VGAM]{multinomial}} and for
#' models with ordinal response like \code{\link[VGAM]{sratio}}, \code{\link[VGAM]{cratio}},
#' \code{\link[VGAM]{cumulative}} or \code{\link[VGAM]{acat}}. \cr
#'
#' For more details on plotting effect stars see \code{\link{effectstars}}.
#'
#' @param x A \code{\link[VGAM]{vglm}}-object.
#' @param p.values Should the p-values of the single coefficients be included in the labels? Default
#' is \code{FALSE}.
#' @param symmetric Should the parameters be transformed to parameters with symmetric (sum-to-zero)
#' side constraints instead of using reference levels. Default is \code{TRUE} for
#' \code{\link[VGAM]{multinomial}}-models. If the \code{\link[VGAM]{multinomial}}-model
#' contains object-specific covariates (\code{xij} argument from \code{\link[VGAM]{vglm.control}}) symmetric side
#' constraints are not possible. In ordinal response models,
#' no side constraints are needed and the option is obsolete.
#' @param plot.parallel Should parallel parameters (equal over all response categories) be
#' represented by effect stars. Default is \code{FALSE}.
#' @param ... further arguments for generic function \code{\link{effectstars}}.
#' @seealso \code{\link{effectstars}} \code{\link{effectstars.DIFlasso}}
#' @author Gunther Schauberger \cr
#' \email{gunther.schauberger@tum.de} \cr \cr
#' \url{https://www.sg.tum.de/epidemiologie/team/schauberger/}
#'
#' @references Tutz, G. and Schauberger, G. (2013): \emph{Visualization of Categorical Response Models -
#' from Data Glyphs to Parameter Glyphs}, Journal of Computational and Graphical Statistics 22(1), 156--177.
#' \cr \cr Gerhard Tutz (2012): \emph{Regression for Categorical Data}, Cambridge University Press
#' @examples
#' \dontrun{
#' ############################################
#' ### Examples for multinomial logit model
#' ############################################
#'
#' ### German election data
#' data(election)
#' library(VGAM)
#' m_elect <- vglm(Partychoice ~ Gender + West + Age + Union + Highschool + Unemployment
#' + Pol.Interest + Democracy + Religion, family = multinomial(), data = election)
#'
#' effectstars(m_elect)
#'
#' # include p.values
#' effectstars(m_elect, p.values = TRUE)
#'
#' ### German election data with category-specific covariates
#'
#' data(election)
#' election[,13:16] <- election[,13:16] - election[,12]
#' election[,18:21] <- election[,18:21] - election[,17]
#' election[,23:26] <- election[,23:26] - election[,22]
#' election[,28:31] <- election[,28:31] - election[,27]
#'
#' election$Social <- election$Social_SPD
#' election$Immigration <- election$Immigration_SPD
#' election$Nuclear <- election$Nuclear_SPD
#' election$Left_Right <- election$Left_Right_SPD
#'
#' m.all <- vglm(Partychoice ~ Social + Immigration + Nuclear + Left_Right + Age +
#' Religion + Democracy + Pol.Interest + Unemployment + Highschool + Union + West +
#' Gender, data = election,
#' family = multinomial(parallel = TRUE~-1 + Social + Immigration +
#' Nuclear + Left_Right, refLevel = 1),
#' xij = list(Social ~ Social_SPD + Social_FDP + Social_Greens + Social_Left,
#' Immigration ~ Immigration_SPD + Immigration_FDP +
#' Immigration_Greens + Immigration_Left,
#' Nuclear ~ Nuclear_SPD + Nuclear_FDP +
#' Nuclear_Greens + Nuclear_Left,
#' Left_Right ~ Left_Right_SPD + Left_Right_FDP +
#' Left_Right_Greens + Left_Right_Left),
#' form2 = ~Social + Immigration + Nuclear + Left_Right + Age +
#' Religion + Democracy + Pol.Interest + Unemployment + Highschool + Union + West +
#' Gender + Social_SPD + Social_FDP + Social_Greens + Social_Left +
#' Immigration_SPD + Immigration_FDP + Immigration_Greens + Immigration_Left +
#' Nuclear_SPD + Nuclear_FDP + Nuclear_Greens + Nuclear_Left +
#' Left_Right_SPD + Left_Right_FDP + Left_Right_Greens + Left_Right_Left
#' )
#'
#' effectstars(m.all, symmetric = FALSE, p.values = TRUE)
#' summary(m.all)
#'
#'
#' ### Chilean plebiscite data
#' data(plebiscite)
#' m_chile <- vglm(Vote ~ ., family = multinomial(), data = plebiscite)
#' effectstars(m_chile)
#'
#' # choose fixed circle sizes and use reference category instead of symmetric side constraints
#' effectstars(m_chile, symmetric = FALSE, fixed = TRUE)
#'
#' ############################################
#' ### Examples for ordinal data
#' ############################################
#'
#' ### Munich insolvency data
#' data(insolvency)
#' insolvency$Age <- scale(insolvency$Age)
#'
#' my_formula <- Insolvency ~ Age + Gender
#'
#' m_acat <- vglm(my_formula, data = insolvency,family = acat())
#' m_cratio <- vglm(my_formula, data = insolvency,family = cratio())
#' m_sratio <- vglm(my_formula, data = insolvency,family = sratio())
#' m_cumulative <- vglm(my_formula, data = insolvency,family = cumulative())
#'
#' summary(m_acat)
#' effectstars(m_acat, p.values = TRUE)
#'
#' summary(m_cratio)
#' effectstars(m_cratio, p.values = TRUE)
#'
#' summary(m_sratio)
#' effectstars(m_sratio, p.values = TRUE)
#'
#' summary(m_cumulative)
#' effectstars(m_cumulative, p.values = TRUE)
#' }
effectstars.vglm <- function(x, p.values = FALSE, symmetric = TRUE, plot.parallel = FALSE, ...){
if(!(attributes(x@family)$vfamily[1] %in% c("multinomial","cumulative","sratio","acat","cratio"))){
stop("Currently effect stars are are only implemented for family multinomial
and for the ordinal families cumulative, sratio, cratio and acat")
}
if(attributes(x@family)$vfamily[1] == "multinomial"){
# check if xij is specified
if(!all(dim(x@Xm2)==c(0,0))){
# check if symmetric=TRUE is specified
if(symmetric){
warning("If argument xij was specified in vglm, symmetric side constraints can not be used.
Option 'symmetric' is set 'symmetric = FALSE'!")
symmetric <- FALSE
}
}
coef.se <- coef.se.multinomial(x = x, symmetric = symmetric)
}else{
coef.se <- coef.se.ordinal(x = x)
}
coefs <- coef.se$coefficients
se <- coef.se$stderr
which.parallel <- coef.se$which.parallel
is_parallel <- coef.se$is_parallel
p <- nrow(coefs)
k <- ncol(coefs)
# calculate pvalues
pvalues <- 2*( 1 - pnorm(abs(coefs/se)))
# create printable pvalues
pval.print <- pvalues
is.NAN <- is.nan(pvalues)
pval.print[pvalues < 5e-04] <- 0
pval.print[pvalues >= 5e-04 & pvalues < 0.001] <- 0.001
pval.print <- format(pval.print, digits = 1, nsmall = 3)
# are labels in command?
plot_labels <- NULL
if(!is.null(match.call()$labels)){
labels <- plot_labels <- eval(match.call()$labels)
}else{
# category names as labels for coefs, se, and pvalues
labels <- plot_labels <- (x@misc$ynames)[1:k]
}
# rownames of coefficients, including categories and intercations
complete.names <- x@misc$colnames.x
complete.names[complete.names == "(Intercept)"] <- "Intercept"
# allnames contains all names of all covariates
if(all.vars(x@misc$formula)[2]=="."){
allnames <- attributes(terms(x@misc$formula,data=get(x@misc$dataname)))$term.labels
}else{
allnames <- attributes(terms(x@misc$formula))$term.labels
}
if(all.vars(x@misc$formula)[2]=="."){
allnames <- c("Intercept",allnames)}
else{
if(attributes(terms(x@misc$formula))$intercept==1){
allnames <- c("Intercept",allnames)
}
}
allnames <- unique(c(allnames,unlist(strsplit(allnames,":"))))
##############
# initialize factor_cat, contains all factor levels of covariates if covariates are factors
factor_cat <- complete.names
# after deleting covariate names from complete.names only the categories are left
for(i in 1:length(allnames)){
factor_cat <- sub(allnames[i],"",factor_cat,fixed = TRUE)
}
for(i in 1:length(factor_cat)){
if(substr(factor_cat[i], 1, 1)==":"){
factor_cat[i] <- substr(factor_cat[i], 2,nchar(factor_cat[i]))
}
if(substr(factor_cat[i], nchar(factor_cat[i]), nchar(factor_cat[i]))==":"){
factor_cat[i] <- substr(factor_cat[i], 1,nchar(factor_cat[i])-1)
}}
# create plot_subs, either by subs from call or with factor levels from factor_cat
if(is.null(match.call()$subs)){
if(!all(factor_cat==rep("",p))){
plot_subs <- factor_cat
plot_subs[factor_cat!=""] <- paste0("(",factor_cat[factor_cat!=""],")")
}
}else{
plot_subs <- match.call()$subs
}
# create plot_names, either by names from call or by removing factor_cat from complete.names
if(!is.null(match.call()$names)){
plot_names <- match.call()$names
}else{
plot_names <- complete.names
# get all factor labels, also from interactions
individ.facs <- unlist(strsplit(factor_cat,":"))
unique_facs <- c(factor_cat, individ.facs)
for(i in 1:length(unique_facs)){
if(substr(unique_facs[i], 1, 1)==":"){
unique_facs[i] <- substr(unique_facs[i], 2,nchar(unique_facs[i]))
}
if(substr(unique_facs[i], nchar(unique_facs[i]), nchar(unique_facs[i]))==":"){
unique_facs[i] <- substr(unique_facs[i], 1,nchar(unique_facs[i])-1)
}
}
unique_facs <- unique(unique_facs[unique_facs!=""])
# remove everything from unique.facs from the covariate names
if(length(unique_facs)>0){
unique_facs <- unique_facs[order(nchar(unique_facs),decreasing = TRUE)]
for(i in 1:length(unique_facs)){
plot_names <- sub(unique_facs[i],"",plot_names,fixed = TRUE)
}
}
}
# name cols and rows of coefs, se and pvalues
colnames(coefs) <- colnames(se) <- colnames(pvalues) <- labels
rownames(coefs) <- rownames(se) <- rownames(pvalues) <- complete.names
# if p.values, create a label matrix containing pvalues
if(p.values){
plot_labels <- rep(plot_labels,each = p)
plot_labels[c(!is.NAN)] <- paste0(plot_labels[!is.NAN], "\n (", pval.print[!is.NAN], ")")
plot_labels <- matrix(plot_labels,ncol=k)
}
## remove everything corresponding to parallel covariates if necessary
plot_coefs <- coefs
if(!plot.parallel){
if(!all(!is_parallel)){
if(is.null(match.call()$names)){
plot_names <- plot_names[-which.parallel]
}
plot_coefs <- plot_coefs[-which.parallel,,drop=FALSE]
if(is.matrix(plot_labels)){
plot_labels <- plot_labels[-which.parallel,,drop=FALSE]
}
plot_subs <- plot_subs[-which.parallel]
}
}
# create a list of arguments for function effect stars
arg.list <- list(x = exp(plot_coefs), names = plot_names, subs = plot_subs, labels = plot_labels,
control = star.ctrl(), cols = NULL, fixed = FALSE, scale = 1)
if(!is.null(match.call()$control)){
arg.list$control <- match.call()$control
}
if(!is.null(match.call()$cols)){
arg.list$cols <- match.call()$cols
}
if(!is.null(match.call()$fixed)){
arg.list$fixed <- match.call()$fixed
}
if(!is.null(match.call()$scale)){
arg.list$scale <- match.call()$scale
}
# plot effect stars
do.call(effectstars, arg.list)
# return invisible stuff
invisible(list(model = x, coefficients = coefs, stderr = se, pvalues = pvalues))
}
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.