Plot update

R/model-timeseries.R

@@ -25,13 +25,13 @@ ztimeseries$methods(
                               fd  = "First Differences: E(Y|X1) - E(Y|X)",
                               acf = "Autocorrelation Function",
                               ev.shortrun = "Expected Values Immediately Resulting from Shock",
-                              ev.longrun = "Long Run Expected Values after Innovation",
+                              ev.longrun = "Long Run Expected Values after Persistant Shift",
                               pv.shortrun = "Predicted Values Immediately Resulting from Shock",
                               pv.longrun = "Long Run Predicted Values after Innovation",
                               evseries.shock = "Expected Values Over Time from Shock",
-                              evseries.innovation ="Expected Values Over Time from Innovation",
+                              evseries.innovation ="Expected Values Over Time from Persistant Shift",
                               pvseries.shock = "Predicted Values Over Time from Shock",
-                              pvseries.innovation ="Predicted Values Over Time from Innovation")
+                              pvseries.innovation ="Predicted Values Over Time from Persistant Shift")
   }
 )
 

R/plots.R

@@ -21,7 +21,7 @@
 #' @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) {
+simulations.plot <-function(y, y1=NULL, xlab=NULL, ylab=NULL, 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)))
@@ -49,6 +49,14 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
                 # Create a sequence of names
                 nameseq <- paste("Y=", min(y):max(y), sep="")
 
+                # Asign default axis labels
+                if (is.null(ylab)) {
+                    ylab <- "Frequency"
+                }
+                if (is.null(xlab)) {
+                    xlab <- "Y Value"
+                }
+
                 # 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)
@@ -74,12 +82,32 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
                 all.names <- names(y)
             }
 
+            # Assign Default Axis Labels
+            if (is.null(ylab)) {
+                ylab <- "Frequency"
+            }
+            if (is.null(xlab)) {
+                xlab <- "Value"
+            }
+
             # 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)){
+            # Assign Default Graph Labels
+            if (is.null(ylab)) {
+                ylab <- "Density"
+            }
+            if (is.null(xlab)) {
+                if (main == "First Differences: E(Y|X1) - E(Y|X)") {
+                    xlab <- "Difference in Outcome"
+                } else {
+                    xlab <- "Outcome"
+                }
+            }
+
             if(ncol(as.matrix(y))>1){
                 ncoly <- ncol(y)
                 hold.dens <- list()
@@ -124,6 +152,7 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
                 }
 
             }else{
+
                 den.y <- density(y)
                 output <- plot(den.y, xlab=xlab, ylab=ylab, main=main, col=line.col[1])
                 if(!identical(col[1],"n")){
@@ -145,7 +174,8 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
             }
 
             yseq<-min(c(y,y1)):max(c(y,y1))
-            nameseq<- paste("Y=",yseq,sep="")
+            nameseqX <- paste("(Y|X)=",yseq,sep="")
+            nameseqX1 <- paste("(Y|X1)=",yseq,sep="")
             n.y<-length(yseq)
 
             colors<-rev(heat.colors(n.y^2))
@@ -176,24 +206,29 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
                 }
             }
 
-            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)
+            axis(side=1,labels=nameseqX, at=seq(0,1,length=n.y), cex.axis=1, las=1)
+            axis(side=2,labels=nameseqX1, 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
+                    n.y1 <- ncol(y1)
+                    nameseqX <- 1:n.y
+                    nameseqX1 <- 1:n.y1
                 }else{
-                    nameseq <- colnames(y)
+                    nameseqX <- colnames(y)
+                    nameseqX1 <- colnames(y1)
                 }
             }else{
-                nameseq <- names(y)
+                nameseqX <- names(y)
+                nameseqX1 <- names(y1)
             }
 
             n.y <- ncol(y)
@@ -220,7 +255,7 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
             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 )
+            image(z=comp, axes=FALSE, col=colors, zlim=c(min(comp),max(comp)),main=main, xlab = xlab, ylab = ylab)
 
             locations.x<-seq(from=0,to=1,length=nrow(comp))
             locations.y<-locations.x
@@ -231,13 +266,26 @@ simulations.plot <-function(y, y1=NULL, xlab="", ylab="", main="", col=NULL, lin
                 }
             }
 
-            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)
+            axis(side=1,labels=nameseqX, at=seq(0,1,length=n.y), cex.axis=1, las=1)
+            axis(side=2,labels=nameseqX1, 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)){
+            ## Assign Default Axis Titles
+            if (is.null(ylab)) {
+                ylab <- "Density"
+            }
+            if (is.null(xlab)) {
+                if (main == "First Differences: E(Y|X1) - E(Y|X)") {
+                    xlab <- "Difference in Outcome"
+                } else {
+                    xlab <- "Outcome"
+                }
+
+            }
+
 
             if(is.null(col)){
                 semi.col.x <-rgb(142,229,238,150,maxColorValue=255)
@@ -593,10 +641,10 @@ ci.plot <- function(obj, qi = "ev", var = NULL, ..., main = NULL, sub = NULL,
         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]], 
+            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))
@@ -610,7 +658,7 @@ ci.plot <- function(obj, qi = "ev", var = NULL, ..., main = NULL, sub = NULL,
                 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))

vignettes/zelig_graphing.Rmd

@@ -0,0 +1,89 @@
+---
+title: "Graphing With Zelig"
+date: "`r Sys.Date()`"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{zelig-graphing}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+---
+
+```{r setup, include=FALSE}
+knitr::opts_knit$set(
+    stop_on_error = 2L
+)
+knitr::opts_chunk$set(
+    fig.height = 11,
+    fig.width = 7
+)
+
+options(cite = FALSE)
+```
+
+## Graphing Syntax:
+
+
+```{r, eval=FALSE}
+z.out <- zelig(Y ~ X, data=your.data, model="your choice", ...) # Run your model
+x.out <- setx(z.out, X = value) # counterfactual for simulation
+x.out <- setx1(x.out, X = other.value) # Optional Counterfactual: For time series and first differences
+s.out <- sim(x.out) # Simulation required for plotting
+plot(s.out) # Graph your simulation
+```
+
+
+## General Overview
+
+Graphing in Zelig works in concert with the simulation and counterfactual analysis functionality included in Zelig. Zelig provides a number graphs to represent the results of the simulations. However, as simulation is not univerally part of the standard statistical workflow, interpretation of some the graphs produced by Zelig is not always straightforward.  
+
+## Types of Graphs
+
+#### Value Distribution
+
+The majority of graphs produced by the Zelig `plot` function are estimates of the probability distributions of simulated model statistics. The probabilty distributions are estimated by the `density` function in the `R` `stats` package. This implements a kernel density estimator. Due to this, the Zelig graphs aren't true probability distributions, but they are reasonable approximations of them.
+
+<img src="Density.png" height="300" />
+
+An example of a Zelig generated distribution plot from a Poisson model.
+
+#### Histograms:
+
+Zelig's histograms represent the relative frequencies of outputs based on user input counterfactuals and the simulated model parameters. The are output when a model's dependent variable has discrete, rather than continuous, output.
+
+<img src="Histogram.png" height="300"/>
+
+An example of a Zelig generated histogram from a Poisson model.
+
+#### Discrete Value Comparisons:
+
+When comparing the results of counterfactuals for discrete valued models, Zelig outputs Contingency Tables of the following type. They visually represent a heatmap of the model output for the two generated counterfactuals.
+
+<img src="Logit Comparison.png" height="300"/>
+
+Contingency Table generated from a Logistic Regression Model
+
+<img src="Poisson Comparison.png" height="300"/>
+
+Contingency Table generated from a Poisson Regression Model
+
+#### Time Series:
+
+Time Series models in Zelig produce a different series of graphics then the other forms of models. In order to be able to utilize the plotting features of Zelig Time series models, the  `setx1` method must have been called prior to running `sim` on the model.
+
+<img src="Correlelogram.png" height="300"/>
+
+The autocorrelation graph, showing the correlation between current values and past values of the model outcome.
+
+<img src="Shock trend.png" height="300"/>
+
+A graph showing the trend in the outcome variable if the underlying variables change to the values set in `setx1` from their means or the values specified in `setx` for a single period.
+
+<img src="Shift trend.png" height="300"/>
+
+A graph showing the trend in the outcome variable if the underlying variables persistently change to the values set in `setx1` from their means or the values specified in `setx`.
+
+
+<!--- ### Customizing Zelig Graphs --->
+
+