1#' Plot Quantities of Interest in a Zelig-fashion 2#' 3#' Various graph generation for different common types of simulated results from 4#' Zelig 5#' @usage simulations.plot(y, y1=NULL, xlab="", ylab="", main="", col=NULL, line.col=NULL, 6#' axisnames=TRUE) 7#' @param y A matrix or vector of simulated results generated by Zelig, to be 8#' graphed. 9#' @param y1 For comparison of two sets of simulated results at different 10#' choices of covariates, this should be an object of the same type and 11#' dimension as y. If no comparison is to be made, this should be NULL. 12#' @param xlab Label for the x-axis. 13#' @param ylab Label for the y-axis. 14#' @param main Main plot title. 15#' @param col A vector of colors. Colors will be used in turn as the graph is 16#' built for main plot objects. For nominal/categorical data, this colors 17#' renders as the bar color, while for numeric data it renders as the background 18#' color. 19#' @param line.col A vector of colors. Colors will be used in turn as the graph is 20#' built for line color shading of plot objects. 21#' @param axisnames a character-vector, specifying the names of the axes 22#' @return nothing 23#' @author James Honaker 24simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, line.col=NULL, axisnames=TRUE) { 25 26 binarytest <- function(j){ 27 if(!is.null(attr(j,"levels"))) return(identical( sort(levels(j)),c(0,1))) 28 return(FALSE) 29 } 30 31 32 33 ## Univariate Plots ## 34 if(is.null(y1)){ 35 36 if (is.null(col)) 37 col <- rgb(100,149,237,maxColorValue=255) 38 39 if (is.null(line.col)) 40 line.col <- "black" 41 42 # Integer Values 43 if ((length(unique(y))<11 & all(as.integer(y) == y)) | is.factor(y) | is.character(y)) { 44 45 if(is.factor(y) | is.character(y)){ 46 y <- as.numeric(y) 47 } 48 49 # Create a sequence of names 50 nameseq <- paste("Y=", min(y):max(y), sep="") 51 52 # Set the heights of the barplots. 53 # Note that tablar requires that all out values are greater than zero. 54 # So, we subtract the min value (ensuring everything is at least zero) 55 # then add 1 56 bar.heights <- tabulate(y - min(y) + 1) / length(y) 57 58 # Barplot with (potentially) some zero columns 59 output <- barplot(bar.heights, xlab=xlab, ylab=ylab, main=main, col=col[1], 60 axisnames=axisnames, names.arg=nameseq) 61 62 # Vector of 1's and 0's 63 } else if(ncol(as.matrix(y))>1 & binarytest(y) ){ 64 65 n.y <- nrow(y) 66 # Precedence is names > colnames > 1:n 67 if(is.null(names(y))){ 68 if(is.null(colnames(y) )){ 69 all.names <- 1:n.y 70 }else{ 71 all.names <- colnames(y) 72 } 73 }else{ 74 all.names <- names(y) 75 } 76 77 # Barplot with (potentially) some zero columns 78 output <- barplot( apply(y,2,sum)/n.y, xlab=xlab, ylab=ylab, main=main, col=col[1], 79 axisnames=axisnames, names.arg=all.names) 80 81 # Continuous Values 82 } else if(is.numeric(y)){ 83 if(ncol(as.matrix(y))>1){ 84 ncoly <- ncol(y) 85 hold.dens <- list() 86 ymax <- xmax <- xmin <- rep(0,ncol(y)) 87 for(i in 1:ncoly){ 88 hold.dens[[i]] <- density(y[,i]) 89 ymax[i] <- max(hold.dens[[i]]$y) 90 xmax[i] <- max(hold.dens[[i]]$x) 91 xmin[i] <- min(hold.dens[[i]]$x) 92 } 93 shift <- 0:ncoly 94 all.xlim <- c(min(xmin), max(xmax)) 95 all.ylim <- c(0,ncoly) 96 97 # Precedence is names > colnames > 1:n 98 if(is.null(names(y))){ 99 if(is.null(colnames(y) )){ 100 all.names <- 1:ncoly 101 }else{ 102 all.names <- colnames(y) 103 } 104 }else{ 105 all.names <- names(y) 106 } 107 shrink <- 0.9 108 for(i in 1:ncoly ){ 109 if(i<ncoly){ 110 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") 111 if(!identical(col[1],"n")){ 112 polygon(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], col=col[1]) 113 } 114 abline(h=shift[i+1]) 115 text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4) 116 par(new=TRUE) 117 }else{ 118 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") 119 if(!identical(col[1],"n")){ 120 polygon(hold.dens[[i]]$x, shrink*hold.dens[[i]]$y/ymax[i] + shift[i], col=col[1]) 121 } 122 text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4) 123 } 124 } 125 126 }else{ 127 den.y <- density(y) 128 output <- plot(den.y, xlab=xlab, ylab=ylab, main=main, col=line.col[1]) 129 if(!identical(col[1],"n")){ 130 polygon(den.y$x, den.y$y, col=col[1]) 131 } 132 } 133 } 134 135 ## Comparison Plots ## 136 137 }else{ 138 139 # Integer - Plot and shade a matrix 140 if(( length(unique(y))<11 & all(as.integer(y) == y) ) | is.factor(y) | is.character(y)){ 141 142 if(is.factor(y) | is.character(y)){ 143 y <- as.numeric(y) 144 y1 <- as.numeric(y1) 145 } 146 147 yseq<-min(c(y,y1)):max(c(y,y1)) 148 nameseq<- paste("Y=",yseq,sep="") 149 n.y<-length(yseq) 150 151 colors<-rev(heat.colors(n.y^2)) 152 lab.colors<-c("black","white") 153 comp<-matrix(NA,nrow=n.y,ncol=n.y) 154 155 for(i in 1:n.y){ 156 for(j in 1:n.y){ 157 flag<- y==yseq[i] & y1==yseq[j] 158 comp[i,j]<-mean(flag) 159 } 160 } 161 162 old.pty<-par()$pty 163 old.mai<-par()$mai 164 165 par(pty="s") 166 par(mai=c(0.3,0.3,0.3,0.1)) 167 168 image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main ) 169 170 locations.x<-seq(from=0,to=1,length=nrow(comp)) 171 locations.y<-locations.x 172 173 for(m in 1:n.y){ 174 for(n in 1:n.y){ 175 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]) 176 } 177 } 178 179 axis(side=1,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=1) 180 axis(side=2,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=3) 181 box() 182 par(pty=old.pty,mai=old.mai) 183 ## Two Vectors of 1's and 0's 184 }else if( ncol(as.matrix(y))>1 & binarytest(y) & ncol(as.matrix(y1))>1 & binarytest(y1) ) { 185 186 # Everything in this section assumes ncol(y)==ncol(y1) 187 188 # Precedence is names > colnames > 1:n 189 if(is.null(names(y))){ 190 if(is.null(colnames(y) )){ 191 nameseq <- 1:n.y 192 }else{ 193 nameseq <- colnames(y) 194 } 195 }else{ 196 nameseq <- names(y) 197 } 198 199 n.y <- ncol(y) 200 yseq <- 1:n.y 201 202 y <- y %*% yseq 203 y1 <- y1 %*% yseq 204 205 ## FROM HERE ON -- Replicates above. Should address more generically 206 colors<-rev(heat.colors(n.y^2)) 207 lab.colors<-c("black","white") 208 comp<-matrix(NA,nrow=n.y,ncol=n.y) 209 210 for(i in 1:n.y){ 211 for(j in 1:n.y){ 212 flag<- y==yseq[i] & y1==yseq[j] 213 comp[i,j]<-mean(flag) 214 } 215 } 216 217 old.pty<-par()$pty 218 old.mai<-par()$mai 219 220 par(pty="s") 221 par(mai=c(0.3,0.3,0.3,0.1)) 222 223 image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main ) 224 225 locations.x<-seq(from=0,to=1,length=nrow(comp)) 226 locations.y<-locations.x 227 228 for(m in 1:n.y){ 229 for(n in 1:n.y){ 230 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]) 231 } 232 } 233 234 axis(side=1,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=1) 235 axis(side=2,labels=nameseq, at=seq(0,1,length=n.y), cex.axis=1, las=3) 236 box() 237 par(pty=old.pty,mai=old.mai) 238 239 ## Numeric - Plot two densities on top of each other 240 }else if(is.numeric(y) & is.numeric(y1)){ 241 242 if(is.null(col)){ 243 semi.col.x <-rgb(142,229,238,150,maxColorValue=255) 244 semi.col.x1<-rgb(255,114,86,150,maxColorValue=255) 245 col<-c(semi.col.x,semi.col.x1) 246 }else if(length(col)<2){ 247 col<-c(col,col) 248 } 249 250 if(ncol(as.matrix(y))>1){ 251 shrink <- 0.9 252 ncoly <- ncol(y) # Assumes columns of y match cols y1. Should check or enforce. 253 # Precedence is names > colnames > 1:n 254 if(is.null(names(y))){ 255 if(is.null(colnames(y) )){ 256 all.names <- 1:ncoly 257 }else{ 258 all.names <- colnames(y) 259 } 260 }else{ 261 all.names <- names(y) 262 } 263 264 hold.dens.y <- hold.dens.y1 <- list() 265 ymax <- xmax <- xmin <- rep(0,ncoly) 266 for(i in 1:ncoly){ 267 hold.dens.y[[i]] <- density(y[,i]) 268 hold.dens.y1[[i]] <- density(y1[,i], bw=hold.dens.y[[i]]$bw) 269 ymax[i] <- max(hold.dens.y[[i]]$y, hold.dens.y1[[i]]$y) 270 xmax[i] <- max(hold.dens.y[[i]]$x, hold.dens.y1[[i]]$x) 271 xmin[i] <- min(hold.dens.y[[i]]$x, hold.dens.y1[[i]]$x) 272 } 273 all.xlim <- c(min(xmin), max(xmax)) 274 all.ylim <- c(0,ncoly) 275 shift <- 0:ncoly 276 for(i in 1:ncoly ){ 277 if(i<ncoly){ 278 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") 279 par(new=TRUE) 280 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") 281 282 if(!identical(col[1],"n")){ 283 polygon(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], col=col[1]) 284 } 285 if(!identical(col[2],"n")){ 286 polygon(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], col=col[2]) 287 } 288 abline(h=shift[i+1]) 289 text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4) 290 par(new=TRUE) 291 }else{ 292 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") 293 par(new=TRUE) 294 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") 295 296 if(!identical(col[1],"n")){ 297 polygon(hold.dens.y[[i]]$x, shrink*hold.dens.y[[i]]$y/ymax[i] + shift[i], col=col[1]) 298 } 299 if(!identical(col[2],"n")){ 300 polygon(hold.dens.y1[[i]]$x, shrink*hold.dens.y1[[i]]$y/ymax[i] + shift[i], col=col[2]) 301 } 302 text(x=all.xlim[1], y=(shift[i] + shift[i+1])/2, labels=all.names[i], pos=4) 303 } 304 } 305 }else{ 306 den.y<-density(y) 307 den.y1<-density(y1,bw=den.y$bw) 308 309 all.xlim<-c(min(c(den.y$x,den.y1$x)),max(c(den.y$x,den.y1$x))) 310 all.ylim<-c(min(c(den.y$y,den.y1$y)),max(c(den.y$y,den.y1$y))) 311 312 output<-plot(den.y,xlab=xlab,ylab=ylab,main=main,col=col[1],xlim=all.xlim,ylim=all.ylim) 313 par(new=TRUE) 314 output<-plot(den.y1,xlab=xlab,ylab=ylab,main="",col=col[2],xlim=all.xlim,ylim=all.ylim) 315 316 if(!identical(col[1],"n")){ 317 polygon(den.y$x,den.y$y,col=col[1]) 318 } 319 if(!identical(col[2],"n")){ 320 polygon(den.y1$x,den.y1$y,col=col[2]) 321 } 322 } 323 } 324 } 325} 326 327 328 329 330 331 332#' Default Plot Design For Zelig Model QI's 333#' 334#' @usage qi.plot(obj, ...) 335#' @param obj A reference class zelig5 object 336#' @param ... Parameters to be passed to the `truehist' function which is 337#' implicitly called for numeric simulations 338#' @author James Honaker with panel layouts from Matt Owen 339 340qi.plot <- function (obj, ...) { 341 # Save old state 342 old.par <- par(no.readonly=T) 343 344 if(is_timeseries(obj)){ 345 par(mfcol=c(3,1)) 346 if(obj$bsetx & !obj$bsetx1) { 347 ## If only setx and not setx1 were called on the object 348 zeligACFplot(obj$get_qi("acf", xvalue="x")) 349 } 350 else{ 351 zeligACFplot(obj$get_qi("acf", xvalue="x1")) 352 } 353 ci.plot(obj, qi="pvseries.shock") 354 ci.plot(obj, qi="pvseries.innovation") 355 return() 356 } 357 358 # Determine whether two "Expected Values" qi's exist 359 both.ev.exist <- (length(obj$sim.out$x$ev)>0) & (length(obj$sim.out$x1$ev)>0) 360 # Determine whether two "Predicted Values" qi's exist 361 both.pv.exist <- (length(obj$sim.out$x$pv)>0) & (length(obj$sim.out$x1$pv)>0) 362 363 color.x <- rgb(242, 122, 94, maxColorValue=255) 364 color.x1 <- rgb(100, 149, 237, maxColorValue=255) 365 # Interpolation of the above colors in rgb color space: 366 color.mixed <- rgb(t(round((col2rgb(color.x) + col2rgb(color.x1))/2)), maxColorValue=255) 367 368 if (! ("x" %in% names(obj$sim.out))) { 369 return(par(old.par)) 370 } else if (! ("x1" %in% names(obj$sim.out))) { 371 372 373 panels <- matrix(1:2, 2, 1) 374 375 # The plotting device: 376 # 377 # +-----------+ 378 # | 1 | 379 # +-----------+ 380 # | 2 | 381 # +-----------+ 382 } else { 383 panels <- matrix(c(1:5, 5), ncol=2, nrow=3, byrow = TRUE) 384 385 # the plotting device: 386 # 387 # +-----+-----+ 388 # | 1 | 2 | 389 # +-----+-----+ 390 # | 3 | 4 | 391 # +-----+-----+ 392 # | 5 | 393 # +-----------+ 394 395 panels <- if (xor(both.ev.exist, both.pv.exist)) 396 rbind(panels, c(6, 6)) 397 398 # the plotting device: 399 # 400 # +-----+-----+ 401 # | 1 | 2 | 402 # +-----+-----+ 403 # | 3 | 4 | 404 # +-----+-----+ 405 # | 5 | 406 # +-----------+ 407 # | 6 | 408 # +-----------+ 409 410 else if (both.ev.exist && both.pv.exist) 411 rbind(panels, c(6, 7)) 412 else 413 panels 414 415 # the plotting device: 416 # 417 # +-----+-----+ 418 # | 1 | 2 | 419 # +-----+-----+ 420 # | 3 | 4 | 421 # +-----+-----+ 422 # | 5 | 423 # +-----+-----+ 424 # | 6 | 7 | 425 # +-----+-----+ 426 } 427 428 layout(panels) 429 430 titles <- obj$setx.labels 431 432 # Plot each simulation 433 if(length(obj$sim.out$x$pv)>0) 434 simulations.plot(obj$get_qi(qi="pv", xvalue="x"), main = titles$pv, col = color.x, line.col = "black") 435 436 if(length(obj$sim.out$x1$pv)>0) 437 simulations.plot(obj$get_qi(qi="pv", xvalue="x1"), main = titles$pv1, col = color.x1, line.col = "black") 438 439 if(length(obj$sim.out$x$ev)>0) 440 simulations.plot(obj$get_qi(qi="ev", xvalue="x"), main = titles$ev, col = color.x, line.col = "black") 441 442 if(length(obj$sim.out$x1$ev)>0) 443 simulations.plot(obj$get_qi(qi="ev", xvalue="x1"), main = titles$ev1, col = color.x1, line.col = "black") 444 445 if(length(obj$sim.out$x1$fd)>0) 446 simulations.plot(obj$get_qi(qi="fd", xvalue="x1"), main = titles$fd, col = color.mixed, line.col = "black") 447 448 if(both.pv.exist) 449 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") 450 451 if(both.ev.exist) 452 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") 453 454 455 # Restore old state 456 par(old.par) 457 458 # Return old parameter invisibly 459 invisible(old.par) 460} 461 462 463 464#' Method for plotting qi simulations across a range within a variable, with confidence intervals 465#' 466#' @param obj A reference class zelig5 object 467#' @param qi a character-string specifying the quantity of interest to plot 468#' @param var The variable to be used on the x-axis. Default is the variable 469#' across all the chosen values with smallest nonzero variance 470#' @param ... Parameters to be passed to the `truehist' function which is 471#' implicitly called for numeric simulations 472#' @param main a character-string specifying the main heading of the plot 473#' @param sub a character-string specifying the sub heading of the plot 474#' @param xlab a character-string specifying the label for the x-axis 475#' @param ylab a character-string specifying the label for the y-axis 476#' @param xlim Limits to the x-axis 477#' @param ylim Limits to the y-axis 478#' @param legcol ``legend color'', an valid color used for plotting the line 479#' colors in the legend 480#' @param col a valid vector of colors of at least length 3 to use to color the 481#' confidence intervals 482#' @param leg ``legend position'', an integer from 1 to 4, specifying the 483#' position of the legend. 1 to 4 correspond to ``SE'', ``SW'', ``NW'', and 484#' ``NE'' respectively. Setting to 0 or ``n'' turns off the legend. 485#' @param legpos ``legend type'', exact coordinates and sizes for legend. 486#' Overrides argment ``leg.type'' 487#' @param ci vector of length three of confidence interval levels to draw. 488#' @param discont optional point of discontinuity along the x-axis at which 489#' to interupt the graph 490#' @return the current graphical parameters. This is subject to change in future 491#' implementations of Zelig 492#' @author James Honaker 493#' @usage ci.plot(obj, qi="ev", var=NULL, ..., main = NULL, sub = 494#' NULL, xlab = NULL, ylab = NULL, xlim = NULL, ylim = 495#' NULL, legcol="gray20", col=NULL, leg=1, legpos= 496#' NULL, ci = c(80, 95, 99.9), discont=NULL) 497#' @export 498 499ci.plot <- function(obj, qi = "ev", var = NULL, ..., main = NULL, sub = NULL, 500 xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, 501 legcol = "gray20", col = NULL, leg = 1, legpos = NULL, 502 ci = c(80, 95, 99.9), discont = NULL) { 503 504 is_zelig(obj) 505 if(!is_timeseries(obj)) is_simsrange(obj$sim.out) 506 msg <- 'Simulations for more than one fitted observation are required.' 507 is_length_not_1(obj$sim.out$range, msg = msg) 508 if (!is.null(obj$sim.out$range1)) { 509 is_length_not_1(obj$sim.out$range1, msg) 510 if (length(obj$sim.out$range) != length(obj$sim.out$range1)) 511 stop('The two fitted data ranges are not the same length.', 512 call. = FALSE) 513 } 514 515 ########################### 516 #### Utility Functions #### 517 # Define function to cycle over range list and extract correct qi's 518 ## CAN THESE NOW BE REPLACED WITH THE GETTER METHODS? 519 520 extract.sims <- function(obj, qi) { 521 d <- length(obj$sim.out$range) 522 k <- length(obj$sim.out$range[[1]][qi][[1]][[1]]) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED 523 hold <- matrix(NA, nrow = k, ncol = d) 524 for (i in 1:d) { 525 hold[, i] <- obj$sim.out$range[[i]][qi][[1]][[1]] # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED 526 } 527 return(hold) 528 } 529 530 extract.sims1 <- function(obj, qi) { 531 # Should find better architecture for alternate range sims 532 d <- length(obj$sim.out$range1) 533 k <- length(obj$sim.out$range1[[1]][qi][[1]][[1]]) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED 534 hold <- matrix(NA, nrow = k, ncol = d) 535 for (i in 1:d) { 536 hold[, i] <- obj$sim.out$range1[[i]][qi][[1]][[1]] # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED 537 } 538 return(hold) 539 } 540 541 # Define functions to compute confidence intervals CAN WE MERGE THESE TOGETHER SO AS NOT TO 542 # HAVE TO SORT TWICE? 543 ci.upper <- function(x, alpha) { 544 pos <- max(round((1 - (alpha/100)) * length(x)), 1) 545 return(sort(x)[pos]) 546 } 547 548 ci.lower <- function(x, alpha) { 549 pos <- max(round((alpha/100) * length(x)), 1) 550 return(sort(x)[pos]) 551 } 552 553 ########################### 554 555 if(length(ci)<3){ 556 ci<-rep(ci,3) 557 } 558 if(length(ci)>3){ 559 ci<-ci[1:3] 560 } 561 ci<-sort(ci) 562 563 ## Timeseries: 564 if(is_timeseries(obj)){ 565 #xmatrix<- ## Do we need to know the x in which the shock/innovation occcured? For secondary graphs, titles, legends? 566 xname <- "Time" 567 qiseries <- c("pvseries.shock","pvseries.innovation","evseries.shock","evseries.innovation") 568 if (!qi %in% qiseries){ 569 cat(paste("Error: For Timeseries models, argument qi must be one of ", paste(qiseries, collapse=" or ") ,".\n", sep="") ) 570 return() 571 } 572 if (obj$bsetx & !obj$bsetx1) { 573 ## If setx has been called and setx1 has not been called 574 ev<-t( obj$get_qi(qi=qi, xvalue="x") ) # NOTE THE NECESSARY TRANSPOSE. Should we more clearly standardize this? 575 } else { 576 ev<-t( obj$get_qi(qi=qi, xvalue="x1") ) # NOTE THE NECESSARY TRANSPOSE. Should we more clearly standardize this? 577 } 578 d<-ncol(ev) 579 xseq<-1:d 580 ev1 <- NULL # Maybe want to add ability to overlay another graph? 581 582 # Define xlabel 583 if (is.null(xlab)) 584 xlab <- xname 585 if (is.null(ylab)){ 586 if(qi %in% c("pvseries.shock", "pvseries.innovation")) 587 ylab<- as.character(obj$setx.labels["pv"]) 588 if(qi %in% c("evseries.shock", "evseries.innovation")) 589 ylab<- as.character(obj$setx.labels["ev"]) 590 } 591 592 if (is.null(main)) 593 main <- as.character(obj$setx.labels[qi]) 594 if (is.null(discont)) 595 discont <- 22.5 # NEED TO SET AUTOMATICALLY 596 597 ## Everything Else: 598 }else{ 599 d <- length(obj$sim.out$range) 600 601 if (d < 1) { 602 return() # Should add warning 603 } 604 num_cols <- length(obj$setx.out$range[[1]]$mm[[1]] ) 605 xmatrix <- matrix(NA,nrow = d, ncol = num_cols) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED 606 for (i in 1:d){ 607 xmatrix[i,] <- matrix(obj$setx.out$range[[i]]$mm[[1]], 608 ncol = num_cols) # THAT IS A LONG PATH THAT MAYBE SHOULD BE CHANGED 609 } 610 611 if (d == 1 && is.null(var)) { 612 warning("Must specify the `var` parameter when plotting the confidence interval of an unvarying model. Plotting nothing.") 613 return(invisible(FALSE)) 614 } 615 616 xvarnames <- names(as.data.frame( obj$setx.out$range[[1]]$mm[[1]])) # MUST BE A BETTER WAY/PATH TO GET NAMES 617 618 if(is.character(var)){ 619 if( !(var %in% xvarnames ) ){ 620 warning("Specified variable for confidence interval plot is not in estimated model. Plotting nothing.") 621 return(invisible(FALSE)) 622 } 623 } 624 625 if (is.null(var)) { 626 # Determine x-axis variable based on variable with unique fitted values equal to the number of scenarios 627 length_unique <- function(x) length(unique(x)) 628 var.unique <- apply(xmatrix, 2, length_unique) 629 var.seq <- 1:ncol(xmatrix) 630 position <- var.seq[var.unique == d] 631 if (length(position) > 1) { 632 position <- position[1] # arbitrarily pick the first variable if more than one 633 message(sprintf('%s chosen as the x-axis variable. Use the var argument to specify a different variable.', xvarnames[position])) 634 } 635 } else { 636 if(is.numeric(var)){ 637 position <- var 638 }else if(is.character(var)){ 639 position <- grep(var,xvarnames) 640 } 641 } 642 position <- min(position) 643 xseq <- xmatrix[,position] 644 xname <- xvarnames[position] 645 # Define xlabel 646 if (is.null(xlab)) 647 xlab <- paste("Range of",xname) 648 649 # Use "qi" argument to select quantities of interest and set labels 650 ev1<-NULL 651 if(!is.null(obj$sim.out$range1)){ 652 ev1<-extract.sims1(obj,qi=qi) 653 } 654 ev<-extract.sims(obj,qi=qi) 655 if (is.null(ylab)){ 656 ylab <- as.character(obj$setx.labels[qi]) 657 } 658 659 } 660 661 662 663 664 # 665 k<-ncol(ev) 666 n<-nrow(ev) 667 668 # 669 if(is.null(col)){ 670 myblue1<-rgb( 100, 149, 237, alpha=50, maxColorValue=255) 671 myblue2<-rgb( 152, 245, 255, alpha=50, maxColorValue=255) 672 myblue3<-rgb( 191, 239, 255, alpha=70, maxColorValue=255) 673 myred1 <-rgb( 237, 149, 100, alpha=50, maxColorValue=255) 674 myred2 <-rgb( 255, 245, 152, alpha=50, maxColorValue=255) 675 myred3 <-rgb( 255, 239, 191, alpha=70, maxColorValue=255) 676 677 col<-c(myblue1,myblue2,myblue3,myred1,myred2,myred3) 678 }else{ 679 if(length(col)<6){ 680 col<-rep(col,6)[1:6] 681 } 682 } 683 684 # Define function to numerically extract summaries of distributions from set of all simulated qi's 685 form.history <- function (k,xseq,results,ci=c(80,95,99.9)){ 686 687 history<-matrix(NA, nrow=k,ncol=8) 688 for (i in 1:k) { 689 v <- c( 690 xseq[i], 691 median(results[,i]), 692 693 ci.upper(results[,i],ci[1]), 694 ci.lower(results[,i],ci[1]), 695 696 ci.upper(results[,i],ci[2]), 697 ci.lower(results[,i],ci[2]), 698 699 ci.upper(results[,i],ci[3]), 700 ci.lower(results[,i],ci[3]) 701 ) 702 703 history[i, ] <- v 704 } 705 if (k == 1) { 706 left <- c( 707 xseq[1]-.5, 708 median(results[,1]), 709 710 ci.upper(results[,1],ci[1]), 711 ci.lower(results[,1],ci[1]), 712 713 ci.upper(results[,1],ci[2]), 714 ci.lower(results[,1],ci[2]), 715 716 ci.upper(results[,1],ci[3]), 717 ci.lower(results[,1],ci[3]) 718 ) 719 right <- c( 720 xseq[1]+.5, 721 median(results[,1]), 722 723 ci.upper(results[,1],ci[1]), 724 ci.lower(results[,1],ci[1]), 725 726 ci.upper(results[,1],ci[2]), 727 ci.lower(results[,1],ci[2]), 728 729 ci.upper(results[,1],ci[3]), 730 ci.lower(results[,1],ci[3]) 731 ) 732 v <- c( 733 xseq[1], 734 median(results[,1]), 735 736 ci.upper(results[,1],ci[1]), 737 ci.lower(results[,1],ci[1]), 738 739 ci.upper(results[,1],ci[2]), 740 ci.lower(results[,1],ci[2]), 741 742 ci.upper(results[,1],ci[3]), 743 ci.lower(results[,1],ci[3]) 744 ) 745 history <- rbind(left, v, right) 746 } 747 748 return(history) 749 } 750 751 history<- form.history(k,xseq,ev,ci) 752 if(!is.null(ev1)){ 753 history1<- form.history(k,xseq,ev1,ci) 754 }else{ 755 history1<-NULL 756 } 757 758 # This is for small sets that have been duplicated so as to have observable volume 759 if(k==1){ 760 k<-3 761 } 762 763 # Specify x-axis length 764 all.xlim <- if (is.null(xlim)) 765 c(min(c(history[, 1],history1[, 1])),max(c(history[, 1],history1[, 1]))) 766 else 767 xlim 768 769 770 # Specify y-axis length 771 all.ylim <-if (is.null(ylim)) 772 c(min(c(history[, -1],history1[, -1])),max(c(history[, -1],history1[, -1]))) 773 else 774 ylim 775 776 777 # Define y label 778 if (is.null(ylab)) 779 ylab <- "Expected Values: E(Y|X)" 780 781 782 ## This is the plot 783 784 par(bty="n") 785 centralx<-history[,1] 786 centraly<-history[,2] 787 788 789 if(is.null(discont)){ 790 gotok <- k 791 }else{ 792 gotok <- sum(xseq < discont) 793 if((gotok<2) | (gotok>(k-2))){ 794 cat("Warning: Discontinuity is located at edge or outside the range of x-axis.\n") 795 gotok<-k 796 discont<-NULL 797 } 798 if(gotok<k){ 799 gotokp1<- gotok+1 800 centralx<-c(centralx[1:gotok], NA, centralx[gotok+1:length(centralx)]) 801 centraly<-c(centraly[1:gotok], NA, centraly[gotok+1:length(centraly)]) 802 } 803 } 804 805 plot(x=centralx, y=centraly, type="l", xlim=all.xlim, ylim=all.ylim, main = main, sub = sub, xlab=xlab, ylab=ylab) 806 807 polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,7],history[gotok:1,8]),col=col[3],border="white") 808 polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,5],history[gotok:1,6]),col=col[2],border="gray90") 809 polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,3],history[gotok:1,4]),col=col[1],border="gray60") 810 polygon(c(history[1:gotok,1],history[gotok:1,1]),c(history[1:gotok,7],history[gotok:1,8]),col=NA,border="white") 811 812 if(!is.null(discont)){ 813 polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,7],history[k:gotokp1,8]),col=col[3],border="white") 814 polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,5],history[k:gotokp1,6]),col=col[2],border="gray90") 815 polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,3],history[k:gotokp1,4]),col=col[1],border="gray60") 816 polygon(c(history[gotokp1:k,1],history[k:gotokp1,1]),c(history[gotokp1:k,7],history[k:gotokp1,8]),col=NA,border="white") 817 abline(v=discont, lty=5, col="grey85") 818 } 819 820 if(!is.null(ev1)){ 821 822 lines(x=history1[1:gotok, 1], y=history1[1:gotok, 2], type="l") 823 if(!is.null(discont)){ 824 lines(x=history1[gotokp1:k, 1], y=history1[gotokp1:k, 2], type="l") 825 } 826 827 polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,7],history1[gotok:1,8]),col=col[6],border="white") 828 polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,5],history1[gotok:1,6]),col=col[5],border="gray90") 829 polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,3],history1[gotok:1,4]),col=col[4],border="gray60") 830 polygon(c(history1[1:gotok,1],history1[gotok:1,1]),c(history1[1:gotok,7],history1[gotok:1,8]),col=NA,border="white") 831 832 if(!is.null(discont)){ 833 polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,7],history1[k:gotokp1,8]),col=col[6],border="white") 834 polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,5],history1[k:gotokp1,6]),col=col[5],border="gray90") 835 polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,3],history1[k:gotokp1,4]),col=col[4],border="gray60") 836 polygon(c(history1[gotokp1:k,1],history1[k:gotokp1,1]),c(history1[gotokp1:k,7],history1[k:gotokp1,8]),col=NA,border="white") 837 } 838 } 839 840 ## This is the legend 841 if((leg != "n") & (leg != 0)){ 842 843 if(is.null(legpos)){ 844 if(leg==1){ 845 legpos<-c(.91,.04,.2,.05) 846 }else if(leg==2){ 847 legpos<-c(.09,.04,.2,.05) 848 }else if(leg==3){ 849 legpos<-c(.09,.04,.8,.05) 850 }else{ 851 legpos<-c(.91,.04,.8,.05) 852 } 853 } 854 855 lx<-min(all.xlim)+ legpos[1]*(max(all.xlim)- min(all.xlim)) 856 hx<-min(all.xlim)+ (legpos[1]+legpos[2])*(max(all.xlim)- min(all.xlim)) 857 858 deltax<-(hx-lx)*.1 859 860 my<-min(all.ylim) +legpos[3]*min(max(all.ylim) - min(all.ylim)) 861 dy<-legpos[4]*(max(all.ylim) - min(all.ylim)) 862 863 864 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) 865 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) 866 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) 867 lines(c(lx-5*deltax,lx),c(my,my),col="white",lwd=3) 868 lines(c(lx-5*deltax,lx),c(my,my),col=legcol) 869 lines(c(lx,hx),c(my,my)) 870 871 polygon(c(lx,lx,hx,hx),c(my-3*dy,my+3*dy,my+3*dy,my-3*dy),col=col[3],border="white") 872 polygon(c(lx,lx,hx,hx),c(my-2*dy,my+2*dy,my+2*dy,my-2*dy),col=col[2],border="gray90") 873 polygon(c(lx,lx,hx,hx),c(my-1*dy,my+1*dy,my+1*dy,my-1*dy),col=col[1],border="gray60") 874 polygon(c(lx,lx,hx,hx),c(my-3*dy,my+3*dy,my+3*dy,my-3*dy),col=NA,border="white") 875 876 text(lx,my,labels="median",pos=2,cex=0.5,col=legcol) 877 text(lx,my+2*dy,labels=paste("ci",ci[2],sep=""),pos=2,cex=0.5,col=legcol) 878 text(hx,my+1*dy,labels=paste("ci",ci[1],sep=""),pos=4,cex=0.5,col=legcol) 879 text(hx,my+3*dy,labels=paste("ci",ci[3],sep=""),pos=4,cex=0.5,col=legcol) 880 } 881 882} 883 884#' Receiver Operator Characteristic Plots 885#' 886#' The 'rocplot' command generates a receiver operator characteristic plot to 887#' compare the in-sample (default) or out-of-sample fit for two logit or probit 888#' regressions. 889#' 890#' @usage 891#' rocplot(z1, z2, 892#' cutoff = seq(from=0, to=1, length=100), lty1="solid", 893#' lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"), 894#' col1=par("col"), col2=par("col"), 895#' main="ROC Curve", 896#' xlab = "Proportion of 1's Correctly Predicted", 897#' ylab="Proportion of 0's Correctly Predicted", 898#' plot = TRUE, 899#' ... 900#' ) 901#' 902#' @param z1 first model 903#' @param z2 second model 904#' @param cutoff A vector of cut-off values between 0 and 1, at which to 905#' evaluate the proportion of 0s and 1s correctly predicted by the first and 906#' second model. By default, this is 100 increments between 0 and 1 907#' inclusive 908#' @param lty1 the line type of the first model (defaults to 'line') 909#' @param lty2 the line type of the second model (defaults to 'dashed') 910#' @param lwd1 the line width of the first model (defaults to 1) 911#' @param lwd2 the line width of the second model (defaults to 1) 912#' @param col1 the color of the first model (defaults to 'black') 913#' @param col2 the color of the second model (defaults to 'black') 914#' @param main a title for the plot (defaults to "ROC Curve") 915#' @param xlab a label for the X-axis 916#' @param ylab a lavel for the Y-axis 917#' @param plot whether to generate a plot to the selected device 918#' @param \dots additional parameters to be passed to the plot 919#' @return if plot is TRUE, rocplot simply generates a plot. Otherwise, a list 920#' with the following is produced: 921#' \item{roc1}{a matrix containing a vector of x-coordinates and 922#' y-coordinates corresponding to the number of ones and zeros correctly 923#' predicted for the first model.} 924#' \item{roc2}{a matrix containing a vector of x-coordinates and 925#' y-coordinates corresponding to the number of ones and zeros correctly 926#' predicted for the second model.} 927#' \item{area1}{the area under the first ROC curve, calculated using 928#' Reimann sums.} 929#' \item{area2}{the area under the second ROC curve, calculated using 930#' Reimann sums.} 931#' @export 932#" @author Kosuke Imai and Olivia Lau 933rocplot <- function(z1, z2, 934 cutoff = seq(from=0, to=1, length=100), lty1="solid", 935 lty2="dashed", lwd1=par("lwd"), lwd2=par("lwd"), 936 col1=par("col"), col2=par("col"), 937 main="ROC Curve", 938 xlab = "Proportion of 1's Correctly Predicted", 939 ylab="Proportion of 0's Correctly Predicted", 940 plot = TRUE, 941 ...) { 942 y1 <- z1$data[as.character(z1$formula[[2]])] 943 y2 <- z2$data[as.character(z2$formula[[2]])] 944 fitted1 <- fitted(z1)[[1]] 945 fitted2 <- fitted(z2)[[1]] 946 roc1 <- roc2 <- matrix(NA, nrow = length(cutoff), ncol = 2) 947 colnames(roc1) <- colnames(roc2) <- c("ones", "zeros") 948 for (i in 1:length(cutoff)) { 949 roc1[i,1] <- mean(fitted1[y1==1] >= cutoff[i]) 950 roc2[i,1] <- mean(fitted2[y2==1] >= cutoff[i]) 951 roc1[i,2] <- mean(fitted1[y1==0] < cutoff[i]) 952 roc2[i,2] <- mean(fitted2[y2==0] < cutoff[i]) 953 } 954 if (plot) { 955 plot(0:1, 0:1, type = "n", xaxs = "i", yaxs = "i", 956 main=main, xlab=xlab, ylab=ylab, ...) 957 lines(roc1, lty = lty1, lwd = lwd1, col=col1) 958 lines(roc2, lty = lty2, lwd = lwd2, col=col2) 959 abline(1, -1, lty = "dotted") 960 } 961 else { 962 area1 <- area2 <- array() 963 for (i in 2:length(cutoff)) { 964 area1[i-1] <- (roc1[i,2] - roc1[(i-1),2]) * roc1[i,1] 965 area2[i-1] <- (roc2[i,2] - roc2[(i-1),2]) * roc2[i,1] 966 } 967 return(list(roc1 = roc1, 968 roc2 = roc2, 969 area1 = sum(na.omit(area1)), 970 area2 = sum(na.omit(area2)))) 971 } 972} 973 974 975#' Plot Autocorrelation Function from Zelig QI object 976#' @keywords internal 977 978 979zeligACFplot <- 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, ...){ 980 981 x <- z$expected.acf 982 ci.x <- z$ci.acf 983 984 if(omitzero){ 985 x<-x[2:length(x)] 986 ci.x$ci.upper <- ci.x$ci.upper[2:length(ci.x$ci.upper)] 987 ci.x$ci.lower <- ci.x$ci.lower[2:length(ci.x$ci.lower)] 988 } 989 990 if(is.null(ylim)){ 991 ylim<-c(min( c(ci.x$ci.lower, 0, x) ), max( c(ci.x$ci.upper, 0 , x) )) 992 993 } 994 if(is.null(col)){ 995 col <- rgb(100,149,237,maxColorValue=255) 996 } 997 998 bout <- barplot(x, col=col, main=main, xlab=xlab, ylab=ylab, ylim=ylim, ...) 999 1000 n <- length(x) 1001 xseq <- as.vector(bout) 1002 NAseq <- rep(NA, n) 1003 1004 xtemp <- cbind( xseq-epsilon, xseq+epsilon, NAseq) 1005 xtemp <- as.vector(t(xtemp)) 1006 ytemp <- cbind(ci.x$ci.upper, ci.x$ci.upper, NAseq) 1007 ytemp <- as.vector(t(ytemp)) 1008 lines(x=xtemp ,y=ytemp, col=barcol) 1009 1010 ytemp <- cbind(ci.x$ci.lower, ci.x$ci.lower, NAseq) 1011 ytemp <- as.vector(t(ytemp)) 1012 lines(x=xtemp ,y=ytemp, col=barcol) 1013 1014 xtemp <- cbind( xseq, xseq, NAseq) 1015 xtemp <- as.vector(t(xtemp)) 1016 ytemp <- cbind(ci.x$ci.upper, ci.x$ci.lower, NAseq) 1017 ytemp <- as.vector(t(ytemp)) 1018 lines(x=xtemp ,y=ytemp, col=barcol) 1019} 1020