Nothing
#' 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)
}
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.