1.393 - on path to 1.4

- Added Hurricanes data - map and map2stan no longer require start values --- they will sample from priors if no explicit inits are given - Added some basic error trapping to map and map2stan - compare supports WAIC now - map now allows vector parameters, via an abstraction layer that translates to scalar for optim and then back again - precis suppresses vector/matrix parameters by default; optional depth parameter displays them - added postcheck, which computes and plots posterior validation check for each case in data (uses link and sim internally)
rmcelreath · Jul 11, 2014 · 77f4d13 · 77f4d13
1 parent 23e15c3
commit 77f4d13
Show file tree

Hide file tree

Showing 16 changed files with 1,448 additions and 635 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.392
+Version: 1.393
 Date: 2014-05-01
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>

diff --git a/R/WAIC-methods.R b/R/WAIC-methods.R
@@ -4,7 +4,7 @@
 
 WAIC <- function( object , n=1000 , refresh=0.1 , ... ) {
 
-    if ( class(object)!="map2stan" ) stop("Requires map2stan fit")
+    if ( !(class(object)%in%c("map2stan")) ) stop("Requires map2stan fit")
 
     if ( !is.null(attr(object,"WAIC")) ) {
         # already have it stored in object, so just return it

diff --git a/R/compare.r b/R/compare.r
@@ -8,7 +8,7 @@ compare.show <- function( object ) {
 setMethod( "show" , "compareIC" , function(object) compare.show(object) )
 
 # new compare function, defaulting to DIC
-compare <- function( ... , n=1e3 , sort="DIC" ) {
+compare <- function( ... , n=1e3 , sort="DIC" , WAIC=FALSE ) {
     # retrieve list of models
     L <- list(...)
     if ( is.list(L[[1]]) && length(L)==1 )
@@ -18,17 +18,30 @@ compare <- function( ... , n=1e3 , sort="DIC" ) {
     mnames <- match.call()
     mnames <- as.character(mnames)[2:(length(L)+1)]
 
-    DIC.list <- lapply( L , function(z) DIC( z , n=n ) )
-    pD.list <- sapply( DIC.list , function(x) attr(x,"pD") )
+    if ( WAIC==FALSE ) {
+        DIC.list <- lapply( L , function(z) DIC( z , n=n ) )
+        pD.list <- sapply( DIC.list , function(x) attr(x,"pD") )
+    } else {
+        # use WAIC instead of DIC
+        WAIC.list <- lapply( L , function(z) WAIC( z , n=n ) )
+        pD.list <- sapply( WAIC.list , function(x) attr(x,"pWAIC") )
+        DIC.list <- WAIC.list
+    }
+
     DIC.list <- unlist(DIC.list)
 
     dDIC <- DIC.list - min( DIC.list )
     w.DIC <- ICweights( DIC.list )
 
-    result <- data.frame( DIC=DIC.list , pD=pD.list , dDIC=dDIC , weight=w.DIC )
+    if ( WAIC==FALSE )
+        result <- data.frame( DIC=DIC.list , pD=pD.list , dDIC=dDIC , weight=w.DIC )
+    else
+        result <- data.frame( WAIC=DIC.list , pWAIC=pD.list , dWAIC=dDIC , weight=w.DIC )
+
     rownames(result) <- mnames
 
     if ( !is.null(sort) ) {
+        if ( WAIC==TRUE & sort=="DIC" ) sort <- "WAIC"
         result <- result[ order( result[[sort]] ) , ]
     }
 
@@ -37,8 +50,8 @@ compare <- function( ... , n=1e3 , sort="DIC" ) {
 
 # plot method for compareIC results shows deviance in and expected deviance out of sample, for each model, ordered top-to-bottom by rank
 setMethod("plot" , "compareIC" , function(x,y,...) {
-    dev_in <- x@output$DIC - x@output$pD
-    dev_out <- x@output$DIC
+    dev_in <- x@output[[1]] - x@output[[2]]
+    dev_out <- x@output[[1]]
     n <- length(dev_in)
     dotchart( dev_in[n:1] , labels=rownames(x@output)[n:1] , xlab="deviance" , pch=16 , xlim=c(min(dev_in),max(dev_out)) , ... )
     points( dev_out[n:1] , 1:n )
@@ -208,4 +221,76 @@ m3 <- map(
 
 plot(x)
 
+# now map2stan
+
+m0 <- map2stan(
+    alist(
+        pulled_left ~ dbinom(1,theta),
+        logit(theta) <- a,
+        a ~ dnorm(0,1)
+    ) ,
+    data=d,
+    start=list(a=0)
+)
+
+m1 <- map2stan(
+    alist(
+        pulled_left ~ dbinom(1,theta),
+        logit(theta) <- a + bp*prosoc_left,
+        a ~ dnorm(0,1),
+        bp ~ dnorm(0,1)
+    ) ,
+    data=d,
+    start=list(a=0,bp=0)
+)
+
+m2 <- map2stan(
+    alist(
+        pulled_left ~ dbinom(1,theta),
+        logit(theta) <- a + bp*prosoc_left + bpc*condition*prosoc_left,
+        a ~ dnorm(0,1),
+        bp ~ dnorm(0,1),
+        bpc ~ dnorm(0,1)
+    ) ,
+    data=d,
+    start=list(a=0,bp=0,bpc=0)
+)
+
+m3 <- map2stan(
+    alist(
+        pulled_left ~ dbinom(1,theta),
+        logit(theta) <- a + bp*prosoc_left + bc*condition + bpc*condition*prosoc_left,
+        a ~ dnorm(0,1),
+        bp ~ dnorm(0,1),
+        bc ~ dnorm(0,1),
+        bpc ~ dnorm(0,1)
+    ) ,
+    data=d,
+    start=list(a=0,bp=0,bc=0,bpc=0)
+)
+
+m4 <- map2stan(
+    alist(
+        pulled_left ~ dbinom(1,theta),
+        logit(theta) <- a + aj + bp*prosoc_left + bc*condition + bpc*condition*prosoc_left,
+        a ~ dnorm(0,1),
+        aj[actor] ~ dnorm(0,sigma_actor),
+        bp ~ dnorm(0,1),
+        bc ~ dnorm(0,1),
+        bpc ~ dnorm(0,1),
+        sigma_actor ~ dcauchy(0,1)
+    ) ,
+    data=d,
+    start=list(a=0,bp=0,bc=0,bpc=0,sigma_actor=1,aj=rep(0,7))
+)
+
+( x1 <- compare(m0,m1,m2,m3,m4,WAIC=FALSE) )
+
+( x2 <- compare(m0,m1,m2,m3,m4,WAIC=TRUE) )
+
+plot(x1)
+
+plot(x2)
+
+
 }
diff --git a/R/map-class.r b/R/map-class.r
@@ -102,4 +102,40 @@ function( object , n=10000 , clean.names=TRUE , ... ) {
     }
     result
 }
-)
+)
+
+setMethod("extract.samples", "map",
+function(object,n=1e4,...){
+    require(MASS)
+    mu <- object@coef
+    result <- as.data.frame( mvrnorm( n=n , mu=mu , Sigma=vcov(object) ) )
+    # convert vector parameters to vectors in list
+    veclist <- attr(object,"veclist")
+    name_head <- function(aname) strsplit( aname , "[" , fixed=TRUE )[[1]][1]
+    name_index <- function(aname) as.numeric(regmatches( aname , regexec( "\\[(.+)\\]" , aname ) )[[1]][2])
+    if ( length(veclist) > 0 ) {
+        new_result <- list()
+        # copy non-vector samples into new list
+        for ( i in 1:length(result) ) {
+            if ( !( name_head(names(result)[i]) %in% names(veclist) ) ) {
+                new_result[[ names(result)[i] ]] <- result[[i]]
+            }
+        }#i
+        # now build vectors out of paramters with [n] in name
+        for ( i in 1:length(veclist) ) {
+            # empty matrix with parameters on cols and samples on rows
+            # so n-by-m, where m in number of pars in vector and n is number of samples
+            new_matrix <- matrix( 0 , ncol=veclist[[i]]$n , nrow=n )
+            for ( j in 1:length(result) ) {
+                if ( name_head(names(result)[j]) == names(veclist)[i] ) {
+                    the_index <- name_index( names(result)[j] )
+                    new_matrix[,the_index] <- result[[j]]
+                }
+            }
+            new_result[[ names(veclist)[i] ]] <- new_matrix
+        }#i
+        result <- new_result
+    }
+    # return result
+    result
+})
diff --git a/R/map-predict.r b/R/map-predict.r
@@ -49,10 +49,21 @@ function( fit , data , n=1000 , post , probs=NULL , refresh=0.1 , flatten=TRUE ,
             }
 
             # make environment
-            e <- list( as.list(data) , as.list(post[s,]) )
+            init <- list() # holds one row of samples across all params
+            for ( j in 1:length(post) ) {
+                par_name <- names(post)[ j ]
+                dims <- dim( post[[par_name]] )
+                # scalar
+                if ( length(dims)<2 ) init[[par_name]] <- post[[par_name]][s]
+                # vector
+                if ( length(dims)==2 ) init[[par_name]] <- post[[par_name]][s,]
+                # matrix
+                if ( length(dims)==3 ) init[[par_name]] <- post[[par_name]][s,,]
+            }#j
+            e <- list( as.list(data) , as.list(init) )
             e <- unlist( e , recursive=FALSE )
             value[s,] <- eval(parse(text=lm),envir=e)
-        }
+        }#s
         link_out[[i]] <- value
         names(link_out)[i] <- parout
     }