Nothing
R/plots.R
In Zelig: Everyone's Statistical Software
Defines functions
zeligACFplot
rocplot
ci.plot
simulations.plot
Documented in
ci.plot
rocplot
simulations.plot
zeligACFplot
#' Plot Quantities of Interest in a Zelig-fashion
#'
#' Various graph generation for different common types of simulated results from
#' Zelig
#' @usage simulations.plot(y, y1=NULL, xlab="", ylab="", main="", col=NULL, line.col=NULL,
#' axisnames=TRUE)
#' @param y A matrix or vector of simulated results generated by Zelig, to be
#' graphed.
#' @param y1 For comparison of two sets of simulated results at different
#' choices of covariates, this should be an object of the same type and
#' dimension as y. If no comparison is to be made, this should be NULL.
#' @param xlab Label for the x-axis.
#' @param ylab Label for the y-axis.
#' @param main Main plot title.
#' @param col A vector of colors. Colors will be used in turn as the graph is
#' built for main plot objects. For nominal/categorical data, this colors
#' renders as the bar color, while for numeric data it renders as the background
#' color.
#' @param line.col A vector of colors. Colors will be used in turn as the graph is
#' built for line color shading of plot objects.
#' @param axisnames a character-vector, specifying the names of the axes
#' @return nothing
#' @author James Honaker
simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, line.col=NULL, axisnames=TRUE) {
binarytest <- function(j){
if(!is.null(attr(j,"levels"))) return(identical( sort(levels(j)),c(0,1)))
return(FALSE)
}
## Univariate Plots ##
if(is.null(y1)){
if (is.null(col))
col <- rgb(100,149,237,maxColorValue=255)
if (is.null(line.col))
line.col <- "black"
# Integer Values
if ((length(unique(y))<11 & all(as.integer(y) == y)) | is.factor(y) | is.character(y)) {
if(is.factor(y) | is.character(y)){
y <- as.numeric(y)
}
# Create a sequence of names
nameseq <- paste("Y=", min(y):max(y), sep="")
# Set the heights of the barplots.
# Note that tablar requires that all out values are greater than zero.
# So, we subtract the min value (ensuring everything is at least zero)
# then add 1
bar.heights <- tabulate(y - min(y) + 1) / length(y)
# Barplot with (potentially) some zero columns
output <- barplot(bar.heights, xlab=xlab, ylab=ylab, main=main, col=col[1],
axisnames=axisnames, names.arg=nameseq)
# Vector of 1's and 0's
} else if(ncol(as.matrix(y))>1 & binarytest(y) ){
n.y <- nrow(y)
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
all.names <- 1:n.y
}else{
all.names <- colnames(y)
}
}else{
all.names <- names(y)
}
# Barplot with (potentially) some zero columns
output <- barplot( apply(y,2,sum)/n.y, xlab=xlab, ylab=ylab, main=main, col=col[1],
axisnames=axisnames, names.arg=all.names)
# Continuous Values
} else if(is.numeric(y)){
if(ncol(as.matrix(y))>1){
ncoly <- ncol(y)
hold.dens <- list()
ymax <- xmax <- xmin <- rep(0,ncol(y))
for(i in 1:ncoly){
hold.dens[[i]] <- density(y[,i])
ymax[i] <- max(hold.dens[[i]]$y)
xmax[i] <- max(hold.dens[[i]]$x)
xmin[i] <- min(hold.dens[[i]]$x)
}
shift <- 0:ncoly
all.xlim <- c(min(xmin), max(xmax))
all.ylim <- c(0,ncoly)
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
all.names <- 1:ncoly
}else{
all.names <- colnames(y)
}
}else{
all.names <- names(y)
}
shrink <- 0.9
for(i in 1:ncoly ){
if(i<ncoly){
output <- plot(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], xaxt="n", yaxt="n", xlab="", ylab="", main="", col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], col=col[1])
}
abline(h=shift[i+1])
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
par(new=TRUE)
}else{
output <- plot(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], yaxt="n", xlab=xlab, ylab=ylab, main=main, col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], col=col[1])
}
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
}
}
}else{
den.y <- density(y)
output <- plot(den.y, xlab=xlab, ylab=ylab, main=main, col=line.col[1])
if(!identical(col[1],"n")){
polygon(den.y$x, den.y$y, col=col[1])
}
}
}
## Comparison Plots ##
}else{
# Integer - Plot and shade a matrix
if(( length(unique(y))<11 & all(as.integer(y) == y) ) | is.factor(y) | is.character(y)){
if(is.factor(y) | is.character(y)){
y <- as.numeric(y)
y1 <- as.numeric(y1)
}
yseq<-min(c(y,y1)):max(c(y,y1))
nameseq<- paste("Y=",yseq,sep="")
n.y<-length(yseq)
colors<-rev(heat.colors(n.y^2))
lab.colors<-c("black","white")
comp<-matrix(NA,nrow=n.y,ncol=n.y)
for(i in 1:n.y){
for(j in 1:n.y){
flag<- y==yseq[i] & y1==yseq[j]
comp[i,j]<-mean(flag)
}
}
old.pty<-par()$pty
old.mai<-par()$mai
par(pty="s")
par(mai=c(0.3,0.3,0.3,0.1))
image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main )
locations.x<-seq(from=0,to=1,length=nrow(comp))
locations.y<-locations.x
for(m in 1:n.y){
for(n in 1:n.y){
text(x=locations.x[m],y=locations.y[n],labels=paste(round(100*comp[m,n])), col=lab.colors[(comp[m,n]> ((max(comp)-min(comp))/2) )+1])
}
}
axis(side=1,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=1)
axis(side=2,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=3)
box()
par(pty=old.pty,mai=old.mai)
## Two Vectors of 1's and 0's
}else if( ncol(as.matrix(y))>1 & binarytest(y) & ncol(as.matrix(y1))>1 & binarytest(y1) ) {
# Everything in this section assumes ncol(y)==ncol(y1)
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
nameseq <- 1:n.y
}else{
nameseq <- colnames(y)
}
}else{
nameseq <- names(y)
}
n.y <- ncol(y)
yseq <- 1:n.y
y <- y %*% yseq
y1 <- y1 %*% yseq
## FROM HERE ON -- Replicates above. Should address more generically
colors<-rev(heat.colors(n.y^2))
lab.colors<-c("black","white")
comp<-matrix(NA,nrow=n.y,ncol=n.y)
for(i in 1:n.y){
for(j in 1:n.y){
flag<- y==yseq[i] & y1==yseq[j]
comp[i,j]<-mean(flag)
}
}
old.pty<-par()$pty
old.mai<-par()$mai
par(pty="s")
par(mai=c(0.3,0.3,0.3,0.1))
image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main )
locations.x<-seq(from=0,to=1,length=nrow(comp))
locations.y<-locations.x
for(m in 1:n.y){
for(n in 1:n.y){
text(x=locations.x[m],y=locations.y[n],labels=paste(round(100*comp[m,n])), col=lab.colors[(comp[m,n]> ((max(comp)-min(comp))/2) )+1])
}
}
axis(side=1,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=1)
axis(side=2,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=3)
box()
par(pty=old.pty,mai=old.mai)
## Numeric - Plot two densities on top of each other
}else if(is.numeric(y) & is.numeric(y1)){
if(is.null(col)){
semi.col.x <-rgb(142,229,238,150,maxColorValue=255)
semi.col.x1<-rgb(255,114,86,150,maxColorValue=255)
col<-c(semi.col.x,semi.col.x1)
}else if(length(col)<2){
col<-c(col,col)
}
if(ncol(as.matrix(y))>1){
shrink <- 0.9
ncoly <- ncol(y) # Assumes columns of y match cols y1. Should check or enforce.
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
all.names <- 1:ncoly
}else{
all.names <- colnames(y)
}
}else{
all.names <- names(y)
}
hold.dens.y <- hold.dens.y1 <- list()
ymax <- xmax <- xmin <- rep(0,ncoly)
for(i in 1:ncoly){
hold.dens.y[[i]] <- density(y[,i])
hold.dens.y1[[i]] <- density(y1[,i], bw=hold.dens.y[[i]]$bw)
ymax[i] <- max(hold.dens.y[[i]]$y, hold.dens.y1[[i]]$y)
xmax[i] <- max(hold.dens.y[[i]]$x, hold.dens.y1[[i]]$x)
xmin[i] <- min(hold.dens.y[[i]]$x, hold.dens.y1[[i]]$x)
}
all.xlim <- c(min(xmin), max(xmax))
all.ylim <- c(0,ncoly)
shift <- 0:ncoly
for(i in 1:ncoly ){
if(i<ncoly){
output <- plot(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], xaxt="n", yaxt="n", xlab="", ylab="", main="", col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
par(new=TRUE)
output <- plot(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], xaxt="n", yaxt="n", xlab="", ylab="", main="", col=line.col[2], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], col=col[1])
}
if(!identical(col[2],"n")){
polygon(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], col=col[2])
}
abline(h=shift[i+1])
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
par(new=TRUE)
}else{
output <- plot(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], yaxt="n", xlab=xlab, ylab=ylab, main=main, col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
par(new=TRUE)
output <- plot(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], yaxt="n", xlab=xlab, ylab=ylab, main=main, col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], col=col[1])
}
if(!identical(col[2],"n")){
polygon(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], col=col[2])
}
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
}
}
}else{
den.y<-density(y)
den.y1<-density(y1,bw=den.y$bw)
all.xlim<-c(min(c(den.y$x,den.y1$x)),max(c(den.y$x,den.y1$x)))
all.ylim<-c(min(c(den.y$y,den.y1$y)),max(c(den.y$y,den.y1$y)))
output<-plot(den.y,xlab=xlab,ylab=ylab,main=main,col=col[1],xlim=all.xlim,ylim=all.ylim)
par(new=TRUE)
output<-plot(den.y1,xlab=xlab,ylab=ylab,main="",col=col[2],xlim=all.xlim,ylim=all.ylim)
if(!identical(col[1],"n")){
polygon(den.y$x,den.y$y,col=col[1])
}
if(!identical(col[2],"n")){
polygon(den.y1$x,den.y1$y,col=col[2])
}
}
}
}
}
#' Default Plot Design For Zelig Model QI's
#'
#' @usage qi.plot(obj, ...)
#' @param obj A reference class zelig5 object
#' @param ... Parameters to be passed to the `truehist' function which is
#' implicitly called for numeric simulations
#' @author James Honaker with panel layouts from Matt Owen
qi.plot <- function (obj, ...) {
# Save old state
old.par <- par(no.readonly=T)
if(is_timeseries(obj)){
par(mfcol=c(3,1))
if(obj$bsetx & !obj$bsetx1) {
## If only setx and not setx1 were called on the object
zeligACFplot(obj$get_qi("acf", xvalue="x"))
}
else{
zeligACFplot(obj$get_qi("acf", xvalue="x1"))
}
ci.plot(obj, qi="pvseries.shock")
ci.plot(obj, qi="pvseries.innovation")
return()
}
# Determine whether two "Expected Values" qi's exist
both.ev.exist <- (length(obj$sim.out$x$ev)>0) & (length(obj$sim.out$x1$ev)>0)
# Determine whether two "Predicted Values" qi's exist
both.pv.exist <- (length(obj$sim.out$x$pv)>0) & (length(obj$sim.out$x1$pv)>0)
color.x <- rgb(242, 122, 94, maxColorValue=255)
color.x1 <- rgb(100, 149, 237, maxColorValue=255)
# Interpolation of the above colors in rgb color space:
color.mixed <- rgb(t(round((col2rgb(color.x) + col2rgb(color.x1))/2)), maxColorValue=255)
if (! ("x" %in% names(obj$sim.out))) {
return(par(old.par))
} else if (! ("x1" %in% names(obj$sim.out))) {
panels <- matrix(1:2, 2, 1)
# The plotting device:
#
# +-----------+
# | 1 |
# +-----------+
# | 2 |
# +-----------+
} else {
panels <- matrix(c(1:5, 5), ncol=2, nrow=3, byrow = TRUE)
# the plotting device:
#
# +-----+-----+
# | 1 | 2 |
# +-----+-----+
# | 3 | 4 |
# +-----+-----+
# | 5 |
# +-----------+
panels <- if (xor(both.ev.exist, both.pv.exist))
rbind(panels, c(6, 6))
# the plotting device:
#
# +-----+-----+
# | 1 | 2 |
# +-----+-----+
# | 3 | 4 |
# +-----+-----+
# | 5 |
# +-----------+
# | 6 |
# +-----------+
else if (both.ev.exist && both.pv.exist)
rbind(panels, c(6, 7))
else
panels
# the plotting device:
#
# +-----+-----+
# | 1 | 2 |
# +-----+-----+
# | 3 | 4 |
# +-----+-----+
# | 5 |
# +-----+-----+
# | 6 | 7 |
# +-----+-----+
}
layout(panels)
titles <- obj$setx.labels
# Plot each simulation
if(length(obj$sim.out$x$pv)>0)
simulations.plot(obj$get_qi(qi="pv", xvalue="x"), main = titles$pv, col = color.x, line.col = "black")
if(length(obj$sim.out$x1$pv)>0)
simulations.plot(obj$get_qi(qi="pv", xvalue="x1"), main = titles$pv1, col = color.x1, line.col = "black")
if(length(obj$sim.out$x$ev)>0)
simulations.plot(obj$get_qi(qi="ev", xvalue="x"), main = titles$ev, col = color.x, line.col = "black")
if(length(obj$sim.out$x1$ev)>0)
simulations.plot(obj$get_qi(qi="ev", xvalue="x1"), main = titles$ev1, col = color.x1, line.col = "black")
if(length(obj$sim.out$x1$fd)>0)
simulations.plot(obj$get_qi(qi="fd", xvalue="x1"), main = titles$fd, col = color.mixed, line.col = "black")
if(both.pv.exist)
simulations.plot(y=obj$get_qi(qi="pv", xvalue="x"), y1=obj$get_qi(qi="pv", xvalue="x1"), main = "Comparison of Y|X and Y|X1", col = paste(c(color.x, color.x1), "80", sep=""), line.col = "black")
if(both.ev.exist)
simulations.plot(y=obj$get_qi(qi="ev", xvalue="x"), y1=obj$get_qi(qi="ev", xvalue="x1"), main = "Comparison of E(Y|X) and E(Y|X1)", col = paste(c(color.x, color.x1), "80", sep=""), line.col = "black")
# Restore old state
par(old.par)
# Return old parameter invisibly
invisible(old.par)
}
#' Method for plotting qi simulations across a range within a variable, with confidence intervals
#'
#' @param obj A reference class zelig5 object
#' @param qi a character-string specifying the quantity of interest to plot
#' @param var The variable to be used on the x-axis. Default is the variable
#' across all the chosen values with smallest nonzero variance
#' @param ... Parameters to be passed to the `truehist' function which is
#' implicitly called for numeric simulations
#' @param main a character-string specifying the main heading of the plot
#' @param sub a character-string specifying the sub heading of the plot
#' @param xlab a character-string specifying the label for the x-axis
#' @param ylab a character-string specifying the label for the y-axis
#' @param xlim Limits to the x-axis
#' @param ylim Limits to the y-axis
#' @param legcol ``legend color'', an valid color used for plotting the line
#' colors in the legend
#' @param col a valid vector of colors of at least length 3 to use to color the
#' confidence intervals
#' @param leg ``legend position'', an integer from 1 to 4, specifying the
#' position of the legend. 1 to 4 correspond to ``SE'', ``SW'', ``NW'', and
#' ``NE'' respectively. Setting to 0 or ``n'' turns off the legend.
#' @param legpos ``legend type'', exact coordinates and sizes for legend.
#' Overrides argment ``leg.type''
#' @param ci vector of length three of confidence interval levels to draw.
#' @param discont optional point of discontinuity along the x-axis at which
#' to interupt the graph
#' @return the current graphical parameters. This is subject to change in future
#' implementations of Zelig
#' @author James Honaker
#' @usage ci.plot(obj, qi="ev", var=NULL, ..., main = NULL, sub =
#' NULL, xlab = NULL, ylab = NULL, xlim = NULL, ylim =
#' NULL, legcol="gray20", col=NULL, leg=1, legpos=
#' NULL, ci = c(80, 95, 99.9), discont=NULL)
#' @export
ci.plot <- function(obj, qi = "ev", var = NULL, ..., main = NULL, sub = NULL,
xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL,
legcol = "gray20", col = NULL, leg = 1, legpos = NULL,
ci = c(80, 95, 99.9), discont = NULL) {
is_zelig(obj)
if(!is_timeseries(obj)) is_simsrange(obj$sim.out)
msg <- 'Simulations for more than one fitted observation are required.'
is_length_not_1(obj$sim.out$range, msg = msg)
if (!is.null(obj$sim.out$range1)) {
is_length_not_1(obj$sim.out$range1, msg)
if (length(obj$sim.out$range) != length(obj$sim.out$range1))
stop('The two fitted data ranges are not the same length.',
call. = FALSE)
}
###########################
#### Utility Functions ####
# Define function to cycle over range list and extract correct qi's
## CAN THESE NOW BE REPLACED WITH THE GETTER METHODS?
extract.sims <- function(obj, qi) {
d <- length(obj$sim.out$range)
k <- length(obj$sim.out$range[[1]][qi][[1]][[1]]) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
hold <- matrix(NA, nrow = k, ncol = d)
for (i in 1:d) {
hold[, i] <- obj$sim.out$range[[i]][qi][[1]][[1]] # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
}
return(hold)
}
extract.sims1 <- function(obj, qi) {
# Should find better architecture for alternate range sims
d <- length(obj$sim.out$range1)
k <- length(obj$sim.out$range1[[1]][qi][[1]][[1]]) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
hold <- matrix(NA, nrow = k, ncol = d)
for (i in 1:d) {
hold[, i] <- obj$sim.out$range1[[i]][qi][[1]][[1]] # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
}
return(hold)
}
# Define functions to compute confidence intervals CAN WE MERGE THESE TOGETHER SO AS NOT TO
# HAVE TO SORT TWICE?
ci.upper <- function(x, alpha) {
pos <- max(round((1 - (alpha/100)) * length(x)), 1)
return(sort(x)[pos])
}
ci.lower <- function(x, alpha) {
pos <- max(round((alpha/100) * length(x)), 1)
return(sort(x)[pos])
}
###########################
if(length(ci)<3){
ci<-rep(ci,3)
}
if(length(ci)>3){
ci<-ci[1:3]
}
ci<-sort(ci)
## Timeseries:
if(is_timeseries(obj)){
#xmatrix<- ## Do we need to know the x in which the shock/innovation occcured? For secondary graphs, titles, legends?
xname <- "Time"
qiseries <- c("pvseries.shock","pvseries.innovation","evseries.shock","evseries.innovation")
if (!qi %in% qiseries){
cat(paste("Error: For Timeseries models, argument qi must be one of ", paste(qiseries, collapse=" or ") ,".\n", sep="") )
return()
}
if (obj$bsetx & !obj$bsetx1) {
## If setx has been called and setx1 has not been called
ev<-t( obj$get_qi(qi=qi, xvalue="x") ) # NOTE THE NECESSARY TRANSPOSE. Should we more clearly standardize this?
} else {
ev<-t( obj$get_qi(qi=qi, xvalue="x1") ) # NOTE THE NECESSARY TRANSPOSE. Should we more clearly standardize this?
}
d<-ncol(ev)
xseq<-1:d
ev1 <- NULL # Maybe want to add ability to overlay another graph?
# Define xlabel
if (is.null(xlab))
xlab <- xname
if (is.null(ylab)){
if(qi %in% c("pvseries.shock", "pvseries.innovation"))
ylab<- as.character(obj$setx.labels["pv"])
if(qi %in% c("evseries.shock", "evseries.innovation"))
ylab<- as.character(obj$setx.labels["ev"])
}
if (is.null(main))
main <- as.character(obj$setx.labels[qi])
if (is.null(discont))
discont <- 22.5 # NEED TO SET AUTOMATICALLY
## Everything Else:
}else{
d <- length(obj$sim.out$range)
if (d < 1) {
return() # Should add warning
}
num_cols <- length(obj$setx.out$range[[1]]$mm[[1]] )
xmatrix <- matrix(NA,nrow = d, ncol = num_cols) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
for (i in 1:d){
xmatrix[i,] <- matrix(obj$setx.out$range[[i]]$mm[[1]],
ncol = num_cols) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
}
if (d == 1 && is.null(var)) {
warning("Must specify the `var` parameter when plotting the confidence interval of an unvarying model. Plotting nothing.")
return(invisible(FALSE))
}
xvarnames <- names(as.data.frame( obj$setx.out$range[[1]]$mm[[1]])) # MUST BE A BETTER WAY/PATH TO GET NAMES
if(is.character(var)){
if( !(var %in% xvarnames ) ){
warning("Specified variable for confidence interval plot is not in estimated model. Plotting nothing.")
return(invisible(FALSE))
}
}
if (is.null(var)) {
# Determine x-axis variable based on variable with unique fitted values equal to the number of scenarios
length_unique <- function(x) length(unique(x))
var.unique <- apply(xmatrix, 2, length_unique)
var.seq <- 1:ncol(xmatrix)
position <- var.seq[var.unique == d]
if (length(position) > 1) {
position <- position[1] # arbitrarily pick the first variable if more than one
message(sprintf('%s chosen as the x-axis variable. Use the var argument to specify a different variable.', xvarnames[position]))
}
} else {
if(is.numeric(var)){
position <- var
}else if(is.character(var)){
position <- grep(var,xvarnames)
}
}
position <- min(position)
xseq <- xmatrix[,position]
xname <- xvarnames[position]
# Define xlabel
if (is.null(xlab))
xlab <- paste("Range of",xname)
# Use "qi" argument to select quantities of interest and set labels
ev1<-NULL
if(!is.null(obj$sim.out$range1)){
ev1<-extract.sims1(obj,qi=qi)
}
ev<-extract.sims(obj,qi=qi)
if (is.null(ylab)){
ylab <- as.character(obj$setx.labels[qi])
}
}
#
k<-ncol(ev)
n<-nrow(ev)
#
if(is.null(col)){
myblue1<-rgb( 100, 149, 237, alpha=50, maxColorValue=255)
myblue2<-rgb( 152, 245, 255, alpha=50, maxColorValue=255)
myblue3<-rgb( 191, 239, 255, alpha=70, maxColorValue=255)
myred1 <-rgb( 237, 149, 100, alpha=50, maxColorValue=255)
myred2 <-rgb( 255, 245, 152, alpha=50, maxColorValue=255)
myred3 <-rgb( 255, 239, 191, alpha=70, maxColorValue=255)
col<-c(myblue1,myblue2,myblue3,myred1,myred2,myred3)
}else{
if(length(col)<6){
col<-rep(col,6)[1:6]
}
}
# Define function to numerically extract summaries of distributions from set of all simulated qi's
form.history <- function (k,xseq,results,ci=c(80,95,99.9)){
history<-matrix(NA, nrow=k,ncol=8)
for (i in 1:k) {
v <- c(
xseq[i],
median(results[,i]),
ci.upper(results[,i],ci[1]),
ci.lower(results[,i],ci[1]),
ci.upper(results[,i],ci[2]),
ci.lower(results[,i],ci[2]),
ci.upper(results[,i],ci[3]),
ci.lower(results[,i],ci[3])
)
history[i, ] <- v
}
if (k == 1) {
left <- c(
xseq[1]-.5,
median(results[,1]),
ci.upper(results[,1],ci[1]),
ci.lower(results[,1],ci[1]),
ci.upper(results[,1],ci[2]),
ci.lower(results[,1],ci[2]),
ci.upper(results[,1],ci[3]),
ci.lower(results[,1],ci[3])
)
right <- c(
xseq[1]+.5,
median(results[,1]),
ci.upper(results[,1],ci[1]),
ci.lower(results[,1],ci[1]),
ci.upper(results[,1],ci[2]),
ci.lower(results[,1],ci[2]),
ci.upper(results[,1],ci[3]),
ci.lower(results[,1],ci[3])
)
v <- c(
xseq[1],
median(results[,1]),
ci.upper(results[,1],ci[1]),
ci.lower(results[,1],ci[1]),
ci.upper(results[,1],ci[2]),
ci.lower(results[,1],ci[2]),
ci.upper(results[,1],ci[3]),
ci.lower(results[,1],ci[3])
)
history <- rbind(left, v, right)
}
return(history)
}
history<- form.history(k,xseq,ev,ci)
if(!is.null(ev1)){
history1<- form.history(k,xseq,ev1,ci)
}else{
history1<-NULL
}
# This is for small sets that have been duplicated so as to have observable volume
if(k==1){
k<-3
}
# Specify x-axis length
all.xlim <- if (is.null(xlim))
c(min(c(history[, 1],history1[, 1])),max(c(history[, 1],history1[, 1])))
else
xlim
# Specify y-axis length
all.ylim <-if (is.null(ylim))
c(min(c(history[, -1],history1[, -1])),max(c(history[, -1],history1[, -1])))
else
ylim
# Define y label
if (is.null(ylab))
ylab <- "Expected Values: E(Y|X)"
## This is the plot
par(bty="n")
centralx<-history[,1]
centraly<-history[,2]
if(is.null(discont)){
gotok <- k
}else{
gotok <- sum(xseq < discont)
if((gotok<2) | (gotok>(k-2))){
cat("Warning: Discontinuity is located at edge or outside the range of x-axis.\n")
gotok<-k
discont<-NULL
}
if(gotok<k){
gotokp1<- gotok+1
centralx<-c(centralx[1:gotok], NA, centralx[gotok+1:length(centralx)])
centraly<-c(centraly[1:gotok], NA, centraly[gotok+1:length(centraly)])
}
}
plot(x=centralx, y=centraly, type="l", xlim=all.xlim, ylim=all.ylim, main = main, sub = sub, xlab=xlab, ylab=ylab)
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,7],history[gotok:1,8]),col=col[3],border="white")
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,5],history[gotok:1,6]),col=col[2],border="gray90")
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,3],history[gotok:1,4]),col=col[1],border="gray60")
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,7],history[gotok:1,8]),col=NA,border="white")
if(!is.null(discont)){
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,7],history[k:gotokp1,8]),col=col[3],border="white")
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,5],history[k:gotokp1,6]),col=col[2],border="gray90")
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,3],history[k:gotokp1,4]),col=col[1],border="gray60")
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,7],history[k:gotokp1,8]),col=NA,border="white")
abline(v=discont, lty=5, col="grey85")
}
if(!is.null(ev1)){
lines(x=history1[1:gotok, 1], y=history1[1:gotok, 2], type="l")
if(!is.null(discont)){
lines(x=history1[gotokp1:k, 1], y=history1[gotokp1:k, 2], type="l")
}
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,7],history1[gotok:1,8]),col=col[6],border="white")
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,5],history1[gotok:1,6]),col=col[5],border="gray90")
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,3],history1[gotok:1,4]),col=col[4],border="gray60")
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,7],history1[gotok:1,8]),col=NA,border="white")
if(!is.null(discont)){
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,7],history1[k:gotokp1,8]),col=col[6],border="white")
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,5],history1[k:gotokp1,6]),col=col[5],border="gray90")
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,3],history1[k:gotokp1,4]),col=col[4],border="gray60")
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,7],history1[k:gotokp1,8]),col=NA,border="white")
}
}
## This is the legend
if((leg != "n") & (leg != 0)){
if(is.null(legpos)){
if(leg==1){
legpos<-c(.91,.04,.2,.05)
}else if(leg==2){
legpos<-c(.09,.04,.2,.05)
}else if(leg==3){
legpos<-c(.09,.04,.8,.05)
}else{
legpos<-c(.91,.04,.8,.05)
}
}
lx<-min(all.xlim)+ legpos[1]*(max(all.xlim)- min(all.xlim))
hx<-min(all.xlim)+ (legpos[1]+legpos[2])*(max(all.xlim)- min(all.xlim))
deltax<-(hx-lx)*.1
my<-min(all.ylim) +legpos[3]*min(max(all.ylim) - min(all.ylim))
dy<-legpos[4]*(max(all.ylim) - min(all.ylim))
lines(c(hx+deltax,hx+2*deltax,hx+2*deltax,hx+deltax),c(my+3*dy,my+3*dy,my-3*dy,my-3*dy),col=legcol)
lines(c(hx+3*deltax,hx+4*deltax,hx+4*deltax,hx+3*deltax),c(my+1*dy,my+1*dy,my-1*dy,my-1*dy),col=legcol)
lines(c(lx-deltax,lx-2*deltax,lx-2*deltax,lx-deltax),c(my+2*dy,my+2*dy,my-2*dy,my-2*dy),col=legcol)
lines(c(lx-5*deltax,lx),c(my,my),col="white",lwd=3)
lines(c(lx-5*deltax,lx),c(my,my),col=legcol)
lines(c(lx,hx),c(my,my))
polygon(c(lx,lx,hx,hx),c(my-3*dy,my+3*dy,my+3*dy,my-3*dy),col=col[3],border="white")
polygon(c(lx,lx,hx,hx),c(my-2*dy,my+2*dy,my+2*dy,my-2*dy),col=col[2],border="gray90")
polygon(c(lx,lx,hx,hx),c(my-1*dy,my+1*dy,my+1*dy,my-1*dy),col=col[1],border="gray60")
polygon(c(lx,lx,hx,hx),c(my-3*dy,my+3*dy,my+3*dy,my-3*dy),col=NA,border="white")
text(lx,my,labels="median",pos=2,cex=0.5,col=legcol)
text(lx,my+2*dy,labels=paste("ci",ci[2],sep=""),pos=2,cex=0.5,col=legcol)
text(hx,my+1*dy,labels=paste("ci",ci[1],sep=""),pos=4,cex=0.5,col=legcol)
text(hx,my+3*dy,labels=paste("ci",ci[3],sep=""),pos=4,cex=0.5,col=legcol)
}
}
#' Receiver Operator Characteristic Plots
#'
#' The 'rocplot' command generates a receiver operator characteristic plot to
#' compare the in-sample (default) or out-of-sample fit for two logit or probit
#' regressions.
#'
#' @usage
#' rocplot(z1, z2,
#' cutoff = seq(from=0, to=1, length=100), lty1="solid",
#' lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"),
#' col1=par("col"), col2=par("col"),
#' main="ROC Curve",
#' xlab = "Proportion of 1's Correctly Predicted",
#' ylab="Proportion of 0's Correctly Predicted",
#' plot = TRUE,
#' ...
#' )
#'
#' @param z1 first model
#' @param z2 second model
#' @param cutoff A vector of cut-off values between 0 and 1, at which to
#' evaluate the proportion of 0s and 1s correctly predicted by the first and
#' second model. By default, this is 100 increments between 0 and 1
#' inclusive
#' @param lty1 the line type of the first model (defaults to 'line')
#' @param lty2 the line type of the second model (defaults to 'dashed')
#' @param lwd1 the line width of the first model (defaults to 1)
#' @param lwd2 the line width of the second model (defaults to 1)
#' @param col1 the color of the first model (defaults to 'black')
#' @param col2 the color of the second model (defaults to 'black')
#' @param main a title for the plot (defaults to "ROC Curve")
#' @param xlab a label for the X-axis
#' @param ylab a lavel for the Y-axis
#' @param plot whether to generate a plot to the selected device
#' @param \dots additional parameters to be passed to the plot
#' @return if plot is TRUE, rocplot simply generates a plot. Otherwise, a list
#' with the following is produced:
#' \item{roc1}{a matrix containing a vector of x-coordinates and
#' y-coordinates corresponding to the number of ones and zeros correctly
#' predicted for the first model.}
#' \item{roc2}{a matrix containing a vector of x-coordinates and
#' y-coordinates corresponding to the number of ones and zeros correctly
#' predicted for the second model.}
#' \item{area1}{the area under the first ROC curve, calculated using
#' Reimann sums.}
#' \item{area2}{the area under the second ROC curve, calculated using
#' Reimann sums.}
#' @export
#" @author Kosuke Imai and Olivia Lau
rocplot <- function(z1, z2,
cutoff = seq(from=0, to=1, length=100), lty1="solid",
lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"),
col1=par("col"), col2=par("col"),
main="ROC Curve",
xlab = "Proportion of 1's Correctly Predicted",
ylab="Proportion of 0's Correctly Predicted",
plot = TRUE,
...) {
y1 <- z1$data[as.character(z1$formula[[2]])]
y2 <- z2$data[as.character(z2$formula[[2]])]
fitted1 <- fitted(z1)[[1]]
fitted2 <- fitted(z2)[[1]]
roc1 <- roc2 <- matrix(NA, nrow = length(cutoff), ncol = 2)
colnames(roc1) <- colnames(roc2) <- c("ones", "zeros")
for (i in 1:length(cutoff)) {
roc1[i,1] <- mean(fitted1[y1==1] >= cutoff[i])
roc2[i,1] <- mean(fitted2[y2==1] >= cutoff[i])
roc1[i,2] <- mean(fitted1[y1==0] < cutoff[i])
roc2[i,2] <- mean(fitted2[y2==0] < cutoff[i])
}
if (plot) {
plot(0:1, 0:1, type = "n", xaxs = "i", yaxs = "i",
main=main, xlab=xlab, ylab=ylab, ...)
lines(roc1, lty = lty1, lwd = lwd1, col=col1)
lines(roc2, lty = lty2, lwd = lwd2, col=col2)
abline(1, -1, lty = "dotted")
}
else {
area1 <- area2 <- array()
for (i in 2:length(cutoff)) {
area1[i-1] <- (roc1[i,2] - roc1[(i-1),2]) * roc1[i,1]
area2[i-1] <- (roc2[i,2] - roc2[(i-1),2]) * roc2[i,1]
}
return(list(roc1 = roc1,
roc2 = roc2,
area1 = sum(na.omit(area1)),
area2 = sum(na.omit(area2))))
}
}
#' Plot Autocorrelation Function from Zelig QI object
#' @keywords internal
zeligACFplot <- function(z, omitzero=FALSE, barcol="black", epsilon=0.1, col=NULL, main="Autocorrelation Function", xlab="Period", ylab="Correlation of Present Shock with Future Outcomes", ylim=NULL, ...){
x <- z$expected.acf
ci.x <- z$ci.acf
if(omitzero){
x<-x[2:length(x)]
ci.x$ci.upper <- ci.x$ci.upper[2:length(ci.x$ci.upper)]
ci.x$ci.lower <- ci.x$ci.lower[2:length(ci.x$ci.lower)]
}
if(is.null(ylim)){
ylim<-c(min( c(ci.x$ci.lower, 0, x) ), max( c(ci.x$ci.upper, 0 , x) ))
}
if(is.null(col)){
col <- rgb(100,149,237,maxColorValue=255)
}
bout <- barplot(x, col=col, main=main, xlab=xlab, ylab=ylab, ylim=ylim, ...)
n <- length(x)
xseq <- as.vector(bout)
NAseq <- rep(NA, n)
xtemp <- cbind( xseq-epsilon, xseq+epsilon, NAseq)
xtemp <- as.vector(t(xtemp))
ytemp <- cbind(ci.x$ci.upper, ci.x$ci.upper, NAseq)
ytemp <- as.vector(t(ytemp))
lines(x=xtemp ,y=ytemp, col=barcol)
ytemp <- cbind(ci.x$ci.lower, ci.x$ci.lower, NAseq)
ytemp <- as.vector(t(ytemp))
lines(x=xtemp ,y=ytemp, col=barcol)
xtemp <- cbind( xseq, xseq, NAseq)
xtemp <- as.vector(t(xtemp))
ytemp <- cbind(ci.x$ci.upper, ci.x$ci.lower, NAseq)
ytemp <- as.vector(t(ytemp))
lines(x=xtemp ,y=ytemp, col=barcol)
}
Try the Zelig package in your browser
Any scripts or data that you put into this service are public.
Zelig documentation built on Jan. 8, 2021, 2:26 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions zeligACFplot rocplot ci.plot simulations.plot
Documented in ci.plot rocplot simulations.plot zeligACFplot
#' Plot Quantities of Interest in a Zelig-fashion
#'
#' Various graph generation for different common types of simulated results from
#' Zelig
#' @usage simulations.plot(y, y1=NULL, xlab="", ylab="", main="", col=NULL, line.col=NULL,
#' axisnames=TRUE)
#' @param y A matrix or vector of simulated results generated by Zelig, to be
#' graphed.
#' @param y1 For comparison of two sets of simulated results at different
#' choices of covariates, this should be an object of the same type and
#' dimension as y. If no comparison is to be made, this should be NULL.
#' @param xlab Label for the x-axis.
#' @param ylab Label for the y-axis.
#' @param main Main plot title.
#' @param col A vector of colors. Colors will be used in turn as the graph is
#' built for main plot objects. For nominal/categorical data, this colors
#' renders as the bar color, while for numeric data it renders as the background
#' color.
#' @param line.col A vector of colors. Colors will be used in turn as the graph is
#' built for line color shading of plot objects.
#' @param axisnames a character-vector, specifying the names of the axes
#' @return nothing
#' @author James Honaker
simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, line.col=NULL, axisnames=TRUE) {
binarytest <- function(j){
if(!is.null(attr(j,"levels"))) return(identical( sort(levels(j)),c(0,1)))
return(FALSE)
}
## Univariate Plots ##
if(is.null(y1)){
if (is.null(col))
col <- rgb(100,149,237,maxColorValue=255)
if (is.null(line.col))
line.col <- "black"
# Integer Values
if ((length(unique(y))<11 & all(as.integer(y) == y)) | is.factor(y) | is.character(y)) {
if(is.factor(y) | is.character(y)){
y <- as.numeric(y)
}
# Create a sequence of names
nameseq <- paste("Y=", min(y):max(y), sep="")
# Set the heights of the barplots.
# Note that tablar requires that all out values are greater than zero.
# So, we subtract the min value (ensuring everything is at least zero)
# then add 1
bar.heights <- tabulate(y - min(y) + 1) / length(y)
# Barplot with (potentially) some zero columns
output <- barplot(bar.heights, xlab=xlab, ylab=ylab, main=main, col=col[1],
axisnames=axisnames, names.arg=nameseq)
# Vector of 1's and 0's
} else if(ncol(as.matrix(y))>1 & binarytest(y) ){
n.y <- nrow(y)
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
all.names <- 1:n.y
}else{
all.names <- colnames(y)
}
}else{
all.names <- names(y)
}
# Barplot with (potentially) some zero columns
output <- barplot( apply(y,2,sum)/n.y, xlab=xlab, ylab=ylab, main=main, col=col[1],
axisnames=axisnames, names.arg=all.names)
# Continuous Values
} else if(is.numeric(y)){
if(ncol(as.matrix(y))>1){
ncoly <- ncol(y)
hold.dens <- list()
ymax <- xmax <- xmin <- rep(0,ncol(y))
for(i in 1:ncoly){
hold.dens[[i]] <- density(y[,i])
ymax[i] <- max(hold.dens[[i]]$y)
xmax[i] <- max(hold.dens[[i]]$x)
xmin[i] <- min(hold.dens[[i]]$x)
}
shift <- 0:ncoly
all.xlim <- c(min(xmin), max(xmax))
all.ylim <- c(0,ncoly)
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
all.names <- 1:ncoly
}else{
all.names <- colnames(y)
}
}else{
all.names <- names(y)
}
shrink <- 0.9
for(i in 1:ncoly ){
if(i<ncoly){
output <- plot(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], xaxt="n", yaxt="n", xlab="", ylab="", main="", col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], col=col[1])
}
abline(h=shift[i+1])
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
par(new=TRUE)
}else{
output <- plot(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], yaxt="n", xlab=xlab, ylab=ylab, main=main, col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], col=col[1])
}
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
}
}
}else{
den.y <- density(y)
output <- plot(den.y, xlab=xlab, ylab=ylab, main=main, col=line.col[1])
if(!identical(col[1],"n")){
polygon(den.y$x, den.y$y, col=col[1])
}
}
}
## Comparison Plots ##
}else{
# Integer - Plot and shade a matrix
if(( length(unique(y))<11 & all(as.integer(y) == y) ) | is.factor(y) | is.character(y)){
if(is.factor(y) | is.character(y)){
y <- as.numeric(y)
y1 <- as.numeric(y1)
}
yseq<-min(c(y,y1)):max(c(y,y1))
nameseq<- paste("Y=",yseq,sep="")
n.y<-length(yseq)
colors<-rev(heat.colors(n.y^2))
lab.colors<-c("black","white")
comp<-matrix(NA,nrow=n.y,ncol=n.y)
for(i in 1:n.y){
for(j in 1:n.y){
flag<- y==yseq[i] & y1==yseq[j]
comp[i,j]<-mean(flag)
}
}
old.pty<-par()$pty
old.mai<-par()$mai
par(pty="s")
par(mai=c(0.3,0.3,0.3,0.1))
image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main )
locations.x<-seq(from=0,to=1,length=nrow(comp))
locations.y<-locations.x
for(m in 1:n.y){
for(n in 1:n.y){
text(x=locations.x[m],y=locations.y[n],labels=paste(round(100*comp[m,n])), col=lab.colors[(comp[m,n]> ((max(comp)-min(comp))/2) )+1])
}
}
axis(side=1,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=1)
axis(side=2,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=3)
box()
par(pty=old.pty,mai=old.mai)
## Two Vectors of 1's and 0's
}else if( ncol(as.matrix(y))>1 & binarytest(y) & ncol(as.matrix(y1))>1 & binarytest(y1) ) {
# Everything in this section assumes ncol(y)==ncol(y1)
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
nameseq <- 1:n.y
}else{
nameseq <- colnames(y)
}
}else{
nameseq <- names(y)
}
n.y <- ncol(y)
yseq <- 1:n.y
y <- y %*% yseq
y1 <- y1 %*% yseq
## FROM HERE ON -- Replicates above. Should address more generically
colors<-rev(heat.colors(n.y^2))
lab.colors<-c("black","white")
comp<-matrix(NA,nrow=n.y,ncol=n.y)
for(i in 1:n.y){
for(j in 1:n.y){
flag<- y==yseq[i] & y1==yseq[j]
comp[i,j]<-mean(flag)
}
}
old.pty<-par()$pty
old.mai<-par()$mai
par(pty="s")
par(mai=c(0.3,0.3,0.3,0.1))
image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main )
locations.x<-seq(from=0,to=1,length=nrow(comp))
locations.y<-locations.x
for(m in 1:n.y){
for(n in 1:n.y){
text(x=locations.x[m],y=locations.y[n],labels=paste(round(100*comp[m,n])), col=lab.colors[(comp[m,n]> ((max(comp)-min(comp))/2) )+1])
}
}
axis(side=1,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=1)
axis(side=2,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=3)
box()
par(pty=old.pty,mai=old.mai)
## Numeric - Plot two densities on top of each other
}else if(is.numeric(y) & is.numeric(y1)){
if(is.null(col)){
semi.col.x <-rgb(142,229,238,150,maxColorValue=255)
semi.col.x1<-rgb(255,114,86,150,maxColorValue=255)
col<-c(semi.col.x,semi.col.x1)
}else if(length(col)<2){
col<-c(col,col)
}
if(ncol(as.matrix(y))>1){
shrink <- 0.9
ncoly <- ncol(y) # Assumes columns of y match cols y1. Should check or enforce.
# Precedence is names > colnames > 1:n
if(is.null(names(y))){
if(is.null(colnames(y) )){
all.names <- 1:ncoly
}else{
all.names <- colnames(y)
}
}else{
all.names <- names(y)
}
hold.dens.y <- hold.dens.y1 <- list()
ymax <- xmax <- xmin <- rep(0,ncoly)
for(i in 1:ncoly){
hold.dens.y[[i]] <- density(y[,i])
hold.dens.y1[[i]] <- density(y1[,i], bw=hold.dens.y[[i]]$bw)
ymax[i] <- max(hold.dens.y[[i]]$y, hold.dens.y1[[i]]$y)
xmax[i] <- max(hold.dens.y[[i]]$x, hold.dens.y1[[i]]$x)
xmin[i] <- min(hold.dens.y[[i]]$x, hold.dens.y1[[i]]$x)
}
all.xlim <- c(min(xmin), max(xmax))
all.ylim <- c(0,ncoly)
shift <- 0:ncoly
for(i in 1:ncoly ){
if(i<ncoly){
output <- plot(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], xaxt="n", yaxt="n", xlab="", ylab="", main="", col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
par(new=TRUE)
output <- plot(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], xaxt="n", yaxt="n", xlab="", ylab="", main="", col=line.col[2], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], col=col[1])
}
if(!identical(col[2],"n")){
polygon(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], col=col[2])
}
abline(h=shift[i+1])
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
par(new=TRUE)
}else{
output <- plot(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], yaxt="n", xlab=xlab, ylab=ylab, main=main, col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
par(new=TRUE)
output <- plot(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], yaxt="n", xlab=xlab, ylab=ylab, main=main, col=line.col[1], xlim=all.xlim, ylim=all.ylim, type="l")
if(!identical(col[1],"n")){
polygon(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], col=col[1])
}
if(!identical(col[2],"n")){
polygon(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], col=col[2])
}
text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4)
}
}
}else{
den.y<-density(y)
den.y1<-density(y1,bw=den.y$bw)
all.xlim<-c(min(c(den.y$x,den.y1$x)),max(c(den.y$x,den.y1$x)))
all.ylim<-c(min(c(den.y$y,den.y1$y)),max(c(den.y$y,den.y1$y)))
output<-plot(den.y,xlab=xlab,ylab=ylab,main=main,col=col[1],xlim=all.xlim,ylim=all.ylim)
par(new=TRUE)
output<-plot(den.y1,xlab=xlab,ylab=ylab,main="",col=col[2],xlim=all.xlim,ylim=all.ylim)
if(!identical(col[1],"n")){
polygon(den.y$x,den.y$y,col=col[1])
}
if(!identical(col[2],"n")){
polygon(den.y1$x,den.y1$y,col=col[2])
}
}
}
}
}
#' Default Plot Design For Zelig Model QI's
#'
#' @usage qi.plot(obj, ...)
#' @param obj A reference class zelig5 object
#' @param ... Parameters to be passed to the `truehist' function which is
#' implicitly called for numeric simulations
#' @author James Honaker with panel layouts from Matt Owen
qi.plot <- function (obj, ...) {
# Save old state
old.par <- par(no.readonly=T)
if(is_timeseries(obj)){
par(mfcol=c(3,1))
if(obj$bsetx & !obj$bsetx1) {
## If only setx and not setx1 were called on the object
zeligACFplot(obj$get_qi("acf", xvalue="x"))
}
else{
zeligACFplot(obj$get_qi("acf", xvalue="x1"))
}
ci.plot(obj, qi="pvseries.shock")
ci.plot(obj, qi="pvseries.innovation")
return()
}
# Determine whether two "Expected Values" qi's exist
both.ev.exist <- (length(obj$sim.out$x$ev)>0) & (length(obj$sim.out$x1$ev)>0)
# Determine whether two "Predicted Values" qi's exist
both.pv.exist <- (length(obj$sim.out$x$pv)>0) & (length(obj$sim.out$x1$pv)>0)
color.x <- rgb(242, 122, 94, maxColorValue=255)
color.x1 <- rgb(100, 149, 237, maxColorValue=255)
# Interpolation of the above colors in rgb color space:
color.mixed <- rgb(t(round((col2rgb(color.x) + col2rgb(color.x1))/2)), maxColorValue=255)
if (! ("x" %in% names(obj$sim.out))) {
return(par(old.par))
} else if (! ("x1" %in% names(obj$sim.out))) {
panels <- matrix(1:2, 2, 1)
# The plotting device:
#
# +-----------+
# | 1 |
# +-----------+
# | 2 |
# +-----------+
} else {
panels <- matrix(c(1:5, 5), ncol=2, nrow=3, byrow = TRUE)
# the plotting device:
#
# +-----+-----+
# | 1 | 2 |
# +-----+-----+
# | 3 | 4 |
# +-----+-----+
# | 5 |
# +-----------+
panels <- if (xor(both.ev.exist, both.pv.exist))
rbind(panels, c(6, 6))
# the plotting device:
#
# +-----+-----+
# | 1 | 2 |
# +-----+-----+
# | 3 | 4 |
# +-----+-----+
# | 5 |
# +-----------+
# | 6 |
# +-----------+
else if (both.ev.exist && both.pv.exist)
rbind(panels, c(6, 7))
else
panels
# the plotting device:
#
# +-----+-----+
# | 1 | 2 |
# +-----+-----+
# | 3 | 4 |
# +-----+-----+
# | 5 |
# +-----+-----+
# | 6 | 7 |
# +-----+-----+
}
layout(panels)
titles <- obj$setx.labels
# Plot each simulation
if(length(obj$sim.out$x$pv)>0)
simulations.plot(obj$get_qi(qi="pv", xvalue="x"), main = titles$pv, col = color.x, line.col = "black")
if(length(obj$sim.out$x1$pv)>0)
simulations.plot(obj$get_qi(qi="pv", xvalue="x1"), main = titles$pv1, col = color.x1, line.col = "black")
if(length(obj$sim.out$x$ev)>0)
simulations.plot(obj$get_qi(qi="ev", xvalue="x"), main = titles$ev, col = color.x, line.col = "black")
if(length(obj$sim.out$x1$ev)>0)
simulations.plot(obj$get_qi(qi="ev", xvalue="x1"), main = titles$ev1, col = color.x1, line.col = "black")
if(length(obj$sim.out$x1$fd)>0)
simulations.plot(obj$get_qi(qi="fd", xvalue="x1"), main = titles$fd, col = color.mixed, line.col = "black")
if(both.pv.exist)
simulations.plot(y=obj$get_qi(qi="pv", xvalue="x"), y1=obj$get_qi(qi="pv", xvalue="x1"), main = "Comparison of Y|X and Y|X1", col = paste(c(color.x, color.x1), "80", sep=""), line.col = "black")
if(both.ev.exist)
simulations.plot(y=obj$get_qi(qi="ev", xvalue="x"), y1=obj$get_qi(qi="ev", xvalue="x1"), main = "Comparison of E(Y|X) and E(Y|X1)", col = paste(c(color.x, color.x1), "80", sep=""), line.col = "black")
# Restore old state
par(old.par)
# Return old parameter invisibly
invisible(old.par)
}
#' Method for plotting qi simulations across a range within a variable, with confidence intervals
#'
#' @param obj A reference class zelig5 object
#' @param qi a character-string specifying the quantity of interest to plot
#' @param var The variable to be used on the x-axis. Default is the variable
#' across all the chosen values with smallest nonzero variance
#' @param ... Parameters to be passed to the `truehist' function which is
#' implicitly called for numeric simulations
#' @param main a character-string specifying the main heading of the plot
#' @param sub a character-string specifying the sub heading of the plot
#' @param xlab a character-string specifying the label for the x-axis
#' @param ylab a character-string specifying the label for the y-axis
#' @param xlim Limits to the x-axis
#' @param ylim Limits to the y-axis
#' @param legcol ``legend color'', an valid color used for plotting the line
#' colors in the legend
#' @param col a valid vector of colors of at least length 3 to use to color the
#' confidence intervals
#' @param leg ``legend position'', an integer from 1 to 4, specifying the
#' position of the legend. 1 to 4 correspond to ``SE'', ``SW'', ``NW'', and
#' ``NE'' respectively. Setting to 0 or ``n'' turns off the legend.
#' @param legpos ``legend type'', exact coordinates and sizes for legend.
#' Overrides argment ``leg.type''
#' @param ci vector of length three of confidence interval levels to draw.
#' @param discont optional point of discontinuity along the x-axis at which
#' to interupt the graph
#' @return the current graphical parameters. This is subject to change in future
#' implementations of Zelig
#' @author James Honaker
#' @usage ci.plot(obj, qi="ev", var=NULL, ..., main = NULL, sub =
#' NULL, xlab = NULL, ylab = NULL, xlim = NULL, ylim =
#' NULL, legcol="gray20", col=NULL, leg=1, legpos=
#' NULL, ci = c(80, 95, 99.9), discont=NULL)
#' @export
ci.plot <- function(obj, qi = "ev", var = NULL, ..., main = NULL, sub = NULL,
xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL,
legcol = "gray20", col = NULL, leg = 1, legpos = NULL,
ci = c(80, 95, 99.9), discont = NULL) {
is_zelig(obj)
if(!is_timeseries(obj)) is_simsrange(obj$sim.out)
msg <- 'Simulations for more than one fitted observation are required.'
is_length_not_1(obj$sim.out$range, msg = msg)
if (!is.null(obj$sim.out$range1)) {
is_length_not_1(obj$sim.out$range1, msg)
if (length(obj$sim.out$range) != length(obj$sim.out$range1))
stop('The two fitted data ranges are not the same length.',
call. = FALSE)
}
###########################
#### Utility Functions ####
# Define function to cycle over range list and extract correct qi's
## CAN THESE NOW BE REPLACED WITH THE GETTER METHODS?
extract.sims <- function(obj, qi) {
d <- length(obj$sim.out$range)
k <- length(obj$sim.out$range[[1]][qi][[1]][[1]]) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
hold <- matrix(NA, nrow = k, ncol = d)
for (i in 1:d) {
hold[, i] <- obj$sim.out$range[[i]][qi][[1]][[1]] # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
}
return(hold)
}
extract.sims1 <- function(obj, qi) {
# Should find better architecture for alternate range sims
d <- length(obj$sim.out$range1)
k <- length(obj$sim.out$range1[[1]][qi][[1]][[1]]) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
hold <- matrix(NA, nrow = k, ncol = d)
for (i in 1:d) {
hold[, i] <- obj$sim.out$range1[[i]][qi][[1]][[1]] # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
}
return(hold)
}
# Define functions to compute confidence intervals CAN WE MERGE THESE TOGETHER SO AS NOT TO
# HAVE TO SORT TWICE?
ci.upper <- function(x, alpha) {
pos <- max(round((1 - (alpha/100)) * length(x)), 1)
return(sort(x)[pos])
}
ci.lower <- function(x, alpha) {
pos <- max(round((alpha/100) * length(x)), 1)
return(sort(x)[pos])
}
###########################
if(length(ci)<3){
ci<-rep(ci,3)
}
if(length(ci)>3){
ci<-ci[1:3]
}
ci<-sort(ci)
## Timeseries:
if(is_timeseries(obj)){
#xmatrix<- ## Do we need to know the x in which the shock/innovation occcured? For secondary graphs, titles, legends?
xname <- "Time"
qiseries <- c("pvseries.shock","pvseries.innovation","evseries.shock","evseries.innovation")
if (!qi %in% qiseries){
cat(paste("Error: For Timeseries models, argument qi must be one of ", paste(qiseries, collapse=" or ") ,".\n", sep="") )
return()
}
if (obj$bsetx & !obj$bsetx1) {
## If setx has been called and setx1 has not been called
ev<-t( obj$get_qi(qi=qi, xvalue="x") ) # NOTE THE NECESSARY TRANSPOSE. Should we more clearly standardize this?
} else {
ev<-t( obj$get_qi(qi=qi, xvalue="x1") ) # NOTE THE NECESSARY TRANSPOSE. Should we more clearly standardize this?
}
d<-ncol(ev)
xseq<-1:d
ev1 <- NULL # Maybe want to add ability to overlay another graph?
# Define xlabel
if (is.null(xlab))
xlab <- xname
if (is.null(ylab)){
if(qi %in% c("pvseries.shock", "pvseries.innovation"))
ylab<- as.character(obj$setx.labels["pv"])
if(qi %in% c("evseries.shock", "evseries.innovation"))
ylab<- as.character(obj$setx.labels["ev"])
}
if (is.null(main))
main <- as.character(obj$setx.labels[qi])
if (is.null(discont))
discont <- 22.5 # NEED TO SET AUTOMATICALLY
## Everything Else:
}else{
d <- length(obj$sim.out$range)
if (d < 1) {
return() # Should add warning
}
num_cols <- length(obj$setx.out$range[[1]]$mm[[1]] )
xmatrix <- matrix(NA,nrow = d, ncol = num_cols) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
for (i in 1:d){
xmatrix[i,] <- matrix(obj$setx.out$range[[i]]$mm[[1]],
ncol = num_cols) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED
}
if (d == 1 && is.null(var)) {
warning("Must specify the `var` parameter when plotting the confidence interval of an unvarying model. Plotting nothing.")
return(invisible(FALSE))
}
xvarnames <- names(as.data.frame( obj$setx.out$range[[1]]$mm[[1]])) # MUST BE A BETTER WAY/PATH TO GET NAMES
if(is.character(var)){
if( !(var %in% xvarnames ) ){
warning("Specified variable for confidence interval plot is not in estimated model. Plotting nothing.")
return(invisible(FALSE))
}
}
if (is.null(var)) {
# Determine x-axis variable based on variable with unique fitted values equal to the number of scenarios
length_unique <- function(x) length(unique(x))
var.unique <- apply(xmatrix, 2, length_unique)
var.seq <- 1:ncol(xmatrix)
position <- var.seq[var.unique == d]
if (length(position) > 1) {
position <- position[1] # arbitrarily pick the first variable if more than one
message(sprintf('%s chosen as the x-axis variable. Use the var argument to specify a different variable.', xvarnames[position]))
}
} else {
if(is.numeric(var)){
position <- var
}else if(is.character(var)){
position <- grep(var,xvarnames)
}
}
position <- min(position)
xseq <- xmatrix[,position]
xname <- xvarnames[position]
# Define xlabel
if (is.null(xlab))
xlab <- paste("Range of",xname)
# Use "qi" argument to select quantities of interest and set labels
ev1<-NULL
if(!is.null(obj$sim.out$range1)){
ev1<-extract.sims1(obj,qi=qi)
}
ev<-extract.sims(obj,qi=qi)
if (is.null(ylab)){
ylab <- as.character(obj$setx.labels[qi])
}
}
#
k<-ncol(ev)
n<-nrow(ev)
#
if(is.null(col)){
myblue1<-rgb( 100, 149, 237, alpha=50, maxColorValue=255)
myblue2<-rgb( 152, 245, 255, alpha=50, maxColorValue=255)
myblue3<-rgb( 191, 239, 255, alpha=70, maxColorValue=255)
myred1 <-rgb( 237, 149, 100, alpha=50, maxColorValue=255)
myred2 <-rgb( 255, 245, 152, alpha=50, maxColorValue=255)
myred3 <-rgb( 255, 239, 191, alpha=70, maxColorValue=255)
col<-c(myblue1,myblue2,myblue3,myred1,myred2,myred3)
}else{
if(length(col)<6){
col<-rep(col,6)[1:6]
}
}
# Define function to numerically extract summaries of distributions from set of all simulated qi's
form.history <- function (k,xseq,results,ci=c(80,95,99.9)){
history<-matrix(NA, nrow=k,ncol=8)
for (i in 1:k) {
v <- c(
xseq[i],
median(results[,i]),
ci.upper(results[,i],ci[1]),
ci.lower(results[,i],ci[1]),
ci.upper(results[,i],ci[2]),
ci.lower(results[,i],ci[2]),
ci.upper(results[,i],ci[3]),
ci.lower(results[,i],ci[3])
)
history[i, ] <- v
}
if (k == 1) {
left <- c(
xseq[1]-.5,
median(results[,1]),
ci.upper(results[,1],ci[1]),
ci.lower(results[,1],ci[1]),
ci.upper(results[,1],ci[2]),
ci.lower(results[,1],ci[2]),
ci.upper(results[,1],ci[3]),
ci.lower(results[,1],ci[3])
)
right <- c(
xseq[1]+.5,
median(results[,1]),
ci.upper(results[,1],ci[1]),
ci.lower(results[,1],ci[1]),
ci.upper(results[,1],ci[2]),
ci.lower(results[,1],ci[2]),
ci.upper(results[,1],ci[3]),
ci.lower(results[,1],ci[3])
)
v <- c(
xseq[1],
median(results[,1]),
ci.upper(results[,1],ci[1]),
ci.lower(results[,1],ci[1]),
ci.upper(results[,1],ci[2]),
ci.lower(results[,1],ci[2]),
ci.upper(results[,1],ci[3]),
ci.lower(results[,1],ci[3])
)
history <- rbind(left, v, right)
}
return(history)
}
history<- form.history(k,xseq,ev,ci)
if(!is.null(ev1)){
history1<- form.history(k,xseq,ev1,ci)
}else{
history1<-NULL
}
# This is for small sets that have been duplicated so as to have observable volume
if(k==1){
k<-3
}
# Specify x-axis length
all.xlim <- if (is.null(xlim))
c(min(c(history[, 1],history1[, 1])),max(c(history[, 1],history1[, 1])))
else
xlim
# Specify y-axis length
all.ylim <-if (is.null(ylim))
c(min(c(history[, -1],history1[, -1])),max(c(history[, -1],history1[, -1])))
else
ylim
# Define y label
if (is.null(ylab))
ylab <- "Expected Values: E(Y|X)"
## This is the plot
par(bty="n")
centralx<-history[,1]
centraly<-history[,2]
if(is.null(discont)){
gotok <- k
}else{
gotok <- sum(xseq < discont)
if((gotok<2) | (gotok>(k-2))){
cat("Warning: Discontinuity is located at edge or outside the range of x-axis.\n")
gotok<-k
discont<-NULL
}
if(gotok<k){
gotokp1<- gotok+1
centralx<-c(centralx[1:gotok], NA, centralx[gotok+1:length(centralx)])
centraly<-c(centraly[1:gotok], NA, centraly[gotok+1:length(centraly)])
}
}
plot(x=centralx, y=centraly, type="l", xlim=all.xlim, ylim=all.ylim, main = main, sub = sub, xlab=xlab, ylab=ylab)
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,7],history[gotok:1,8]),col=col[3],border="white")
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,5],history[gotok:1,6]),col=col[2],border="gray90")
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,3],history[gotok:1,4]),col=col[1],border="gray60")
polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,7],history[gotok:1,8]),col=NA,border="white")
if(!is.null(discont)){
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,7],history[k:gotokp1,8]),col=col[3],border="white")
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,5],history[k:gotokp1,6]),col=col[2],border="gray90")
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,3],history[k:gotokp1,4]),col=col[1],border="gray60")
polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,7],history[k:gotokp1,8]),col=NA,border="white")
abline(v=discont, lty=5, col="grey85")
}
if(!is.null(ev1)){
lines(x=history1[1:gotok, 1], y=history1[1:gotok, 2], type="l")
if(!is.null(discont)){
lines(x=history1[gotokp1:k, 1], y=history1[gotokp1:k, 2], type="l")
}
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,7],history1[gotok:1,8]),col=col[6],border="white")
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,5],history1[gotok:1,6]),col=col[5],border="gray90")
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,3],history1[gotok:1,4]),col=col[4],border="gray60")
polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,7],history1[gotok:1,8]),col=NA,border="white")
if(!is.null(discont)){
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,7],history1[k:gotokp1,8]),col=col[6],border="white")
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,5],history1[k:gotokp1,6]),col=col[5],border="gray90")
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,3],history1[k:gotokp1,4]),col=col[4],border="gray60")
polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,7],history1[k:gotokp1,8]),col=NA,border="white")
}
}
## This is the legend
if((leg != "n") & (leg != 0)){
if(is.null(legpos)){
if(leg==1){
legpos<-c(.91,.04,.2,.05)
}else if(leg==2){
legpos<-c(.09,.04,.2,.05)
}else if(leg==3){
legpos<-c(.09,.04,.8,.05)
}else{
legpos<-c(.91,.04,.8,.05)
}
}
lx<-min(all.xlim)+ legpos[1]*(max(all.xlim)- min(all.xlim))
hx<-min(all.xlim)+ (legpos[1]+legpos[2])*(max(all.xlim)- min(all.xlim))
deltax<-(hx-lx)*.1
my<-min(all.ylim) +legpos[3]*min(max(all.ylim) - min(all.ylim))
dy<-legpos[4]*(max(all.ylim) - min(all.ylim))
lines(c(hx+deltax,hx+2*deltax,hx+2*deltax,hx+deltax),c(my+3*dy,my+3*dy,my-3*dy,my-3*dy),col=legcol)
lines(c(hx+3*deltax,hx+4*deltax,hx+4*deltax,hx+3*deltax),c(my+1*dy,my+1*dy,my-1*dy,my-1*dy),col=legcol)
lines(c(lx-deltax,lx-2*deltax,lx-2*deltax,lx-deltax),c(my+2*dy,my+2*dy,my-2*dy,my-2*dy),col=legcol)
lines(c(lx-5*deltax,lx),c(my,my),col="white",lwd=3)
lines(c(lx-5*deltax,lx),c(my,my),col=legcol)
lines(c(lx,hx),c(my,my))
polygon(c(lx,lx,hx,hx),c(my-3*dy,my+3*dy,my+3*dy,my-3*dy),col=col[3],border="white")
polygon(c(lx,lx,hx,hx),c(my-2*dy,my+2*dy,my+2*dy,my-2*dy),col=col[2],border="gray90")
polygon(c(lx,lx,hx,hx),c(my-1*dy,my+1*dy,my+1*dy,my-1*dy),col=col[1],border="gray60")
polygon(c(lx,lx,hx,hx),c(my-3*dy,my+3*dy,my+3*dy,my-3*dy),col=NA,border="white")
text(lx,my,labels="median",pos=2,cex=0.5,col=legcol)
text(lx,my+2*dy,labels=paste("ci",ci[2],sep=""),pos=2,cex=0.5,col=legcol)
text(hx,my+1*dy,labels=paste("ci",ci[1],sep=""),pos=4,cex=0.5,col=legcol)
text(hx,my+3*dy,labels=paste("ci",ci[3],sep=""),pos=4,cex=0.5,col=legcol)
}
}
#' Receiver Operator Characteristic Plots
#'
#' The 'rocplot' command generates a receiver operator characteristic plot to
#' compare the in-sample (default) or out-of-sample fit for two logit or probit
#' regressions.
#'
#' @usage
#' rocplot(z1, z2,
#' cutoff = seq(from=0, to=1, length=100), lty1="solid",
#' lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"),
#' col1=par("col"), col2=par("col"),
#' main="ROC Curve",
#' xlab = "Proportion of 1's Correctly Predicted",
#' ylab="Proportion of 0's Correctly Predicted",
#' plot = TRUE,
#' ...
#' )
#'
#' @param z1 first model
#' @param z2 second model
#' @param cutoff A vector of cut-off values between 0 and 1, at which to
#' evaluate the proportion of 0s and 1s correctly predicted by the first and
#' second model. By default, this is 100 increments between 0 and 1
#' inclusive
#' @param lty1 the line type of the first model (defaults to 'line')
#' @param lty2 the line type of the second model (defaults to 'dashed')
#' @param lwd1 the line width of the first model (defaults to 1)
#' @param lwd2 the line width of the second model (defaults to 1)
#' @param col1 the color of the first model (defaults to 'black')
#' @param col2 the color of the second model (defaults to 'black')
#' @param main a title for the plot (defaults to "ROC Curve")
#' @param xlab a label for the X-axis
#' @param ylab a lavel for the Y-axis
#' @param plot whether to generate a plot to the selected device
#' @param \dots additional parameters to be passed to the plot
#' @return if plot is TRUE, rocplot simply generates a plot. Otherwise, a list
#' with the following is produced:
#' \item{roc1}{a matrix containing a vector of x-coordinates and
#' y-coordinates corresponding to the number of ones and zeros correctly
#' predicted for the first model.}
#' \item{roc2}{a matrix containing a vector of x-coordinates and
#' y-coordinates corresponding to the number of ones and zeros correctly
#' predicted for the second model.}
#' \item{area1}{the area under the first ROC curve, calculated using
#' Reimann sums.}
#' \item{area2}{the area under the second ROC curve, calculated using
#' Reimann sums.}
#' @export
#" @author Kosuke Imai and Olivia Lau
rocplot <- function(z1, z2,
cutoff = seq(from=0, to=1, length=100), lty1="solid",
lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"),
col1=par("col"), col2=par("col"),
main="ROC Curve",
xlab = "Proportion of 1's Correctly Predicted",
ylab="Proportion of 0's Correctly Predicted",
plot = TRUE,
...) {
y1 <- z1$data[as.character(z1$formula[[2]])]
y2 <- z2$data[as.character(z2$formula[[2]])]
fitted1 <- fitted(z1)[[1]]
fitted2 <- fitted(z2)[[1]]
roc1 <- roc2 <- matrix(NA, nrow = length(cutoff), ncol = 2)
colnames(roc1) <- colnames(roc2) <- c("ones", "zeros")
for (i in 1:length(cutoff)) {
roc1[i,1] <- mean(fitted1[y1==1] >= cutoff[i])
roc2[i,1] <- mean(fitted2[y2==1] >= cutoff[i])
roc1[i,2] <- mean(fitted1[y1==0] < cutoff[i])
roc2[i,2] <- mean(fitted2[y2==0] < cutoff[i])
}
if (plot) {
plot(0:1, 0:1, type = "n", xaxs = "i", yaxs = "i",
main=main, xlab=xlab, ylab=ylab, ...)
lines(roc1, lty = lty1, lwd = lwd1, col=col1)
lines(roc2, lty = lty2, lwd = lwd2, col=col2)
abline(1, -1, lty = "dotted")
}
else {
area1 <- area2 <- array()
for (i in 2:length(cutoff)) {
area1[i-1] <- (roc1[i,2] - roc1[(i-1),2]) * roc1[i,1]
area2[i-1] <- (roc2[i,2] - roc2[(i-1),2]) * roc2[i,1]
}
return(list(roc1 = roc1,
roc2 = roc2,
area1 = sum(na.omit(area1)),
area2 = sum(na.omit(area2))))
}
}
#' Plot Autocorrelation Function from Zelig QI object
#' @keywords internal
zeligACFplot <- function(z, omitzero=FALSE, barcol="black", epsilon=0.1, col=NULL, main="Autocorrelation Function", xlab="Period", ylab="Correlation of Present Shock with Future Outcomes", ylim=NULL, ...){
x <- z$expected.acf
ci.x <- z$ci.acf
if(omitzero){
x<-x[2:length(x)]
ci.x$ci.upper <- ci.x$ci.upper[2:length(ci.x$ci.upper)]
ci.x$ci.lower <- ci.x$ci.lower[2:length(ci.x$ci.lower)]
}
if(is.null(ylim)){
ylim<-c(min( c(ci.x$ci.lower, 0, x) ), max( c(ci.x$ci.upper, 0 , x) ))
}
if(is.null(col)){
col <- rgb(100,149,237,maxColorValue=255)
}
bout <- barplot(x, col=col, main=main, xlab=xlab, ylab=ylab, ylim=ylim, ...)
n <- length(x)
xseq <- as.vector(bout)
NAseq <- rep(NA, n)
xtemp <- cbind( xseq-epsilon, xseq+epsilon, NAseq)
xtemp <- as.vector(t(xtemp))
ytemp <- cbind(ci.x$ci.upper, ci.x$ci.upper, NAseq)
ytemp <- as.vector(t(ytemp))
lines(x=xtemp ,y=ytemp, col=barcol)
ytemp <- cbind(ci.x$ci.lower, ci.x$ci.lower, NAseq)
ytemp <- as.vector(t(ytemp))
lines(x=xtemp ,y=ytemp, col=barcol)
xtemp <- cbind( xseq, xseq, NAseq)
xtemp <- as.vector(t(xtemp))
ytemp <- cbind(ci.x$ci.upper, ci.x$ci.lower, NAseq)
ytemp <- as.vector(t(ytemp))
lines(x=xtemp ,y=ytemp, col=barcol)
}
Try the Zelig package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.