1.54 - Titanic Iceberg

- bugfix for map2stan imputation combined with certain kinds of varying effects. Added a test for this. - man pages for pairs.map2stan and plot.map2stan - plot.map2stan (trace plots) now recognize n_cols and max_rows arguments for arranging plot grid - WAIC method for map2stan defaults to use all non-warmup samples returned by Stan
rmcelreath · Apr 14, 2015 · ff3a7e8 · ff3a7e8
1 parent 9704c7d
commit ff3a7e8
Show file tree

Hide file tree

Showing 11 changed files with 153 additions and 39 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.53
-Date: 2015-04-03
+Version: 1.54
+Date: 2015-04-13
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>
 Imports: coda, MASS, mvtnorm

diff --git a/NAMESPACE b/NAMESPACE
@@ -17,5 +17,7 @@ importFrom(parallel, makeCluster)
 importFrom(parallel, clusterExport)
 importFrom(mvtnorm, dmvnorm)
 importFrom(mvtnorm, rmvnorm)
+export(dmvnorm)
+export(rmvnorm)
 exportClasses(map2stan)
 exportClasses(map)
diff --git a/R/map2stan-class.r b/R/map2stan-class.r
@@ -14,24 +14,32 @@ setMethod("coef", "map2stan", function(object) {
 })
 
 setMethod("extract.samples","map2stan",
-function(object) {
+function(object,n,...) {
     #require(rstan)
-    p <- rstan::extract(object@stanfit)
+    p <- rstan::extract(object@stanfit,...)
     # get rid of dev and lp__
     p[['dev']] <- NULL
     p[['lp__']] <- NULL
     # get rid of those ugly dimnames
     for ( i in 1:length(p) ) {
         attr(p[[i]],"dimnames") <- NULL
     }
+    if ( !missing(n) ) {
+        for ( i in 1:length(p) ) {
+            n_dims <- length( dim(p[[i]]) )
+            if ( n_dims==1 ) p[[i]] <- p[[i]][1:n]
+            if ( n_dims==2 ) p[[i]] <- p[[i]][1:n,]
+            if ( n_dims==3 ) p[[i]] <- p[[i]][1:n,,]
+        }
+    }
     return(p)
 }
 )
 
 setMethod("extract.samples","stanfit",
-function(object) {
+function(object,...) {
     #require(rstan)
-    p <- rstan::extract(object)
+    p <- rstan::extract(object,...)
     # get rid of dev and lp__
     #p[['dev']] <- NULL
     #p[['lp__']] <- NULL
@@ -236,13 +244,16 @@ setMethod("plot" , "map2stan" , function(x,y,...) {
 
 setMethod("pairs" , "map2stan" , function(x, n=500 , alpha=0.7 , cex=0.7 , pch=16 , adj=1 , pars , ...) {
     #require(rstan)
-    posterior <- extract.samples(x)
-    if ( !missing(pars) ) {
-        # select out named parameters
-        p <- list()
-        for ( k in pars ) p[[k]] <- posterior[[k]]
-        posterior <- p
-    }
+    if ( missing(pars) )
+        posterior <- extract.samples(x)
+    else
+        posterior <- extract.samples(x,pars=pars)
+    #if ( !missing(pars) ) {
+    #    # select out named parameters
+    #    p <- list()
+    #    for ( k in pars ) p[[k]] <- posterior[[k]]
+    #    posterior <- p
+    #}
     panel.dens <- function(x, ...) {
         usr <- par("usr"); on.exit(par(usr))
         par(usr = c(usr[1:2], 0, 1.5) )
@@ -272,13 +283,16 @@ setMethod("pairs" , "map2stan" , function(x, n=500 , alpha=0.7 , cex=0.7 , pch=1
 # my trace plot function
 #rethink_palette <- c("#5BBCD6","#F98400","#F2AD00","#00A08A","#FF0000")
 rethink_palette <- c("#8080FF","#F98400","#F2AD00","#00A08A","#FF0000")
-tracerplot <- function( object , col=rethink_palette , alpha=1 , bg=gray(0.7,0.5) , ask=TRUE , window , n_cols=2 , ... ) {
+tracerplot <- function( object , pars , col=rethink_palette , alpha=1 , bg=gray(0.7,0.5) , ask=TRUE , window , n_cols=2 , max_rows=5 , ... ) {
     chain.cols <- col
 
     if ( class(object)!="map2stan" ) stop( "requires map2stan fit" )
 
     # get all chains, not mixed, from stanfit
-    post <- extract(object@stanfit,permuted=FALSE,inc_warmup=TRUE)
+    if ( missing(pars) )
+        post <- extract(object@stanfit,permuted=FALSE,inc_warmup=TRUE)
+    else
+        post <- extract(object@stanfit,pars=pars,permuted=FALSE,inc_warmup=TRUE)
 
     # names
     dimnames <- attr(post,"dimnames")
@@ -296,8 +310,8 @@ tracerplot <- function( object , col=rethink_palette , alpha=1 , bg=gray(0.7,0.5
     n_rows_per_page <- n_rows
     paging <- FALSE
     n_pages <- 1
-    if ( n_rows_per_page > 5 ) {
-        n_rows_per_page <- 5
+    if ( n_rows_per_page > max_rows ) {
+        n_rows_per_page <- max_rows
         n_pages <- ceiling(n_pars/(n_cols*n_rows_per_page))
         paging <- TRUE
     }

diff --git a/R/map2stan-templates.r b/R/map2stan-templates.r
@@ -71,9 +71,14 @@ map2stan.templates <- list(
                     # vector, so constuct vector data type in Stan code
                     # and add vector to transformed parameters
                     pars <- list()
+                    n <- length(k[[1]])-1
                     for ( i in 2:(n+1) ) pars[[i-1]] <- as.character(k[[1]][[i]])
                     vname <- concat( "Mu_" , paste(pars,collapse="") )
 
+                    # get tpars from parent
+                    m_tpars1 <- get( "m_tpars1" , envir=parenvir )
+                    m_tpars2 <- get( "m_tpars2" , envir=parenvir )
+
                     # add declaration to transformed parameters
                     m_tpars1 <- concat( m_tpars1 , indent , "vector[" , n , "] " , vname , ";\n" )
                     # add transformation

diff --git a/R/map2stan.r b/R/map2stan.r
@@ -900,21 +900,22 @@ map2stan <- function( flist , data , start , pars , constraints=list() , types=l
             d[[ N_txt ]] <- N
 
             # mark grouping variable used
-            #fp[['used_predictors']] <- listappend( fp[['used_predictors']] , list(var=vprior$group,N=length(d[[vprior$group]]) ) )
-            # check likelihoods for matching length and use that N_name
-            if ( length(fp[['likelihood']])>0 ) {
-                groupN <- length(d[[vprior$group]])
-                for ( j in 1:length(fp[['likelihood']]) ) {
-                    if ( fp[['likelihood']][[j]]$N_cases == groupN ) {
-                        #fp[['used_predictors']] <- listappend( fp[['used_predictors']] , list(var=vprior$group,N=fp[['likelihood']][[j]]$N_name ) )
-                        fp[['used_predictors']][[vprior$group]] <- list( var=vprior$group , N=fp[['likelihood']][[j]]$N_name )
-                    }
-                }#j
-            } else {
-                # just add raw integer length
-                #fp[['used_predictors']] <- listappend( fp[['used_predictors']] , list(var=vprior$group,N=length(d[[vprior$group]]) ) )
-                fp[['used_predictors']][[vprior$group]] <- list( var=vprior$group , N=length(d[[vprior$group]]) )
-            }
+            # check if already there
+            if ( is.null(fp[['used_predictors']][[vprior$group]]) ) {
+                # check likelihoods for matching length and use that N_name
+                if ( length(fp[['likelihood']])>0 ) {
+                    groupN <- length(d[[vprior$group]])
+                    for ( j in 1:length(fp[['likelihood']]) ) {
+                        if ( fp[['likelihood']][[j]]$N_cases == groupN ) {
+                            fp[['used_predictors']][[vprior$group]] <- list( var=vprior$group , N=fp[['likelihood']][[j]]$N_name )
+                            break
+                        }
+                    }#j
+                } else {
+                    # just add raw integer length
+                    fp[['used_predictors']][[vprior$group]] <- list( var=vprior$group , N=length(d[[vprior$group]]) )
+                }
+            }# is.null
 
             # check for explicit start value
             for ( k in vprior$pars_out )

diff --git a/R/postcheck.r b/R/postcheck.r
@@ -1,6 +1,6 @@
 # graphical posterior validation checks for map and map2stan
 
-postcheck <- function( fit , prob=0.89 , window=20 , n=1000 , col=rangi2 , ... ) {
+postcheck <- function( fit , x , prob=0.89 , window=20 , n=1000 , col=rangi2 , ... ) {
 
     undot <- function( astring ) {
         astring <- gsub( "." , "_" , astring , fixed=TRUE )
@@ -20,6 +20,15 @@ postcheck <- function( fit , prob=0.89 , window=20 , n=1000 , col=rangi2 , ... )
     if ( class(fit)=="map2stan" ) outcome <- undot(outcome)
     y <- fit@data[[ outcome ]]
 
+    # check for x-axis variable definition
+    if ( !missing(x) ) {
+        if ( x %in% names(fit@data) ) {
+            x_var <- fit@data[[x]]
+            r <- range(x_var)*c(0.9,1.1)
+            x_seq <- seq(from=r[1],to=r[2],length.out=30)
+        }
+    }
+
     # compute posterior predictions for each case
 
     mu <- apply( pred , 2 , mean )

diff --git a/R/sim.r b/R/sim.r
@@ -159,7 +159,7 @@ function( fit , data , n=1000 , post , ll=FALSE , refresh=0.1 , replace=list() ,
 )
 
 setMethod("sim", "map2stan",
-function( fit , data , n , post , refresh=0.1 , replace=list() , ... ) {
+function( fit , data , n=1000 , post , refresh=0.1 , replace=list() , ... ) {
 
     ########################################
     # check arguments
@@ -171,14 +171,16 @@ function( fit , data , n , post , refresh=0.1 , replace=list() , ... ) {
         #data <- as.data.frame(data)
     }
 
+    if ( n==0 ) {
+        ptemp <- extract.samples(fit)
+        n <- dim(ptemp[[1]])[1]
+    }
+
     if ( missing(post) ) 
         post <- extract.samples(fit,n=n)
 
-    if ( missing(n) )
-        n <- dim(post[[1]])[1] # in case n=0 was passed
-
     # get linear model values from link
-    # use our posterior samples, so later parameters have right correlation structure
+    # use our posterior samples, so later parameters have right correlation structure with link values
     # don't flatten result, so we end up with a named list, even if only one element
     pred <- link( fit , data=data , n=n , post=post , flatten=FALSE , refresh=refresh , replace=replace , ... )
 

diff --git a/R/zWAIC-methods.R b/R/zWAIC-methods.R
@@ -9,8 +9,9 @@ function( object , n=1000 , refresh=0.1 , pointwise=FALSE , ... ) {
 }
 )
 
+# by default uses all samples returned by Stan; indicated by n=0
 setMethod("WAIC", "map2stan",
-function( object , n=1000 , refresh=0.1 , pointwise=FALSE , ... ) {
+function( object , n=0 , refresh=0.1 , pointwise=FALSE , ... ) {
 
     if ( !(class(object)%in%c("map2stan")) ) stop("Requires map2stan fit")
 
@@ -24,6 +25,9 @@ function( object , n=1000 , refresh=0.1 , pointwise=FALSE , ... ) {
                 attr(new_waic,"pWAIC") <- sum(unlist(attr(old_waic,"pWAIC")))
                 attr(new_waic,"se") <- attr(old_waic,"se")
                 return( new_waic )
+            } else {
+                # return pointwise
+                return( old_waic )
             }
         } else {
             # old waic not pointwise

diff --git a/man/pairs.map2stan.Rd b/man/pairs.map2stan.Rd
@@ -0,0 +1,36 @@
+\name{pairs.map2stan}
+\alias{pairs.map2stan}
+\alias{pairs.map}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Trace plot of map2stan chains}
+\description{
+  Shows \code{pairs} plots for Stan samples produced by a \code{map} or \code{map2stan} fit. 
+}
+\usage{
+pairs(x, n=500 , alpha=0.7 , cex=0.7 , pch=16 , adj=1 , pars , ...)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{x}{map or map2stan model fit}
+  \item{n}{Number of samples to show in scatter plots}
+  \item{alpha}{alpha transparency of scatter plots}
+  \item{cex}{Character expansion factor for scatter plots}
+  \item{pch}{Point type for scatter plots}
+  \item{adj}{\code{adj} argument for \code{\link{density}} estimate}
+  \item{pars}{Character vector of parameters to display}
+  \item{...}{Additional plot arguments}
+}
+\details{
+  This is the default pairs plot method for both \code{\link{map}} and \code{\link{map2stan}} model fits. 
+}
+\value{
+}
+\references{}
+\author{Richard McElreath}
+\seealso{}
+\examples{
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ }
+
diff --git a/man/plot.map2stan.Rd b/man/plot.map2stan.Rd
@@ -0,0 +1,37 @@
+\name{plot.map2stan}
+\alias{plot.map2stan}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Trace plot of map2stan chains}
+\description{
+  Shows trace plots for Stan samples produced by a \code{map2stan} fit, annotated by number of effective samples. Automatic paging and adjustable number of columns and rows per page.
+}
+\usage{
+plot( object , pars , col=rethink_palette , alpha=1 , bg=gray(0.7,0.5) , ask=TRUE , window , n_cols=2 , max_rows=5 , ... )
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{object}{map2stan model fit}
+  \item{pars}{Character vector of parameters to display}
+  \item{col}{Vector of colors to use for chains. Recycled.}
+  \item{alpha}{alpha transparency of chains}
+  \item{bg}{Background fill for warmup samples}
+  \item{ask}{Whether to pause for paging. Default is \code{TRUE}.}
+  \item{window}{Range of samples to display. Default is all samples.}
+  \item{n_cols}{Number of columns per page}
+  \item{max_rows}{Maximum number of rows per page}
+  \item{...}{Additional arguments to pass to plot commands}
+}
+\details{
+  This is the default trace plot method for \code{\link{map2stan}} model fits. 
+}
+\value{
+}
+\references{}
+\author{Richard McElreath}
+\seealso{}
+\examples{
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ }
+
diff --git a/man/sim.Rd b/man/sim.Rd
@@ -1,5 +1,9 @@
 \name{sim}
 \alias{sim}
+\alias{sim-methods}
+\alias{sim,map2stan-method}
+\alias{sim,map-method}
+\alias{sim,lm-method}
 \title{Simulate posterior observations}
 \description{
   Simulates posterior observations for \code{map} and \code{map2stan} model fits.