Version 1.37

- zero-augmented gamma improved
- generic extract.samples added, as well as method for map2stan
- converted map and map2stan to use alist constructions, allowing <- in
formulas
- flist_untag() supports alist to formula list conversion. used
internally.
- DIC() works on both map and map2stan models
- DIC.map2stan() does what it sounds like it does
- lots of improvements to map2stan, mainly supporting planned features
NYI
- renaming PCI to PI
- added check_index() and coerce_index() for handling index variables
in multilevel models
- renamed homeworkch2 to homeworkch3
- bug fix to man page for col.alpha
- added Rinder data

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.36
-Date: 2014-03-04
+Version: 1.37
+Date: 2014-03-14
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>
 Depends:

R/colors.r

@@ -17,4 +17,6 @@ col.dist <- function( x , mu=0 , sd=1 , col="slateblue" ) {
     cols <- sapply( x , function(z) exp(-(z-mu)^2/sd) )
     cols <- sapply( cols , function(z) col.alpha(col,z) )
     cols
-}
+}
+
+rangi2 <- col.desat( "blue" , 0.5 )

R/distributions.r

@@ -95,15 +95,6 @@ rgampois <- function( n , mu , scale ) {
     rnbinom( n , size=shape , prob=prob )
 }
 
-###
-# zero-inflated gamma functions
-dzerogamma <- function( x , prob , mu , scale , log=TRUE ) {
-    p1 <- dgamma2( x , mu=mu , s=scale , log=FALSE )
-    p <- ifelse( x==0 , prob , (1-prob)*p1 )
-    if ( log==TRUE ) p <- log(p)
-    p
-}
-
 # laplace (double exponential)
 dlaplace <- function(x,location=0,lambda=1,log=FALSE) {
     # f(y) = (1/(2b)) exp( -|y-a|/b )
@@ -183,6 +174,15 @@ rzipois <- function(n,p,lambda) {
     return(y)
 }
 
+# zero-augmented gamma2 distribution
+dzagamma2 <- function(x,prob,mu,scale,log=FALSE) {
+    K <- as.data.frame(cbind(x=x,prob=prob,mu=mu,scale=scale))
+    llg <- dgamma( x , shape=mu/scale , scale=scale , log=TRUE )
+    ll <- ifelse( K$x==0 , log(K$prob) , log(1-K$prob) + llg )
+    if ( log==FALSE ) ll <- exp(ll)
+    ll
+}
+
 # zero-inflated binomial distribution
 # this is slow, but accurate
 dzibinom <- function(x,p_zero,size,prob,log=FALSE) {

R/map-class.r

@@ -47,7 +47,7 @@ setMethod("show", "map", function(object){
         print( object@formula[[i]] )
     }
 
-    cat("\nEstimates:\n")
+    cat("\nMAP values:\n")
     print(coef(object))
 
     cat("\nLog-likelihood: ")
@@ -62,4 +62,26 @@ setMethod("summary", "map", function(object){
 
     precis(object)
 
-})
+})
+
+setGeneric("extract.samples",
+function( object , n=10000 , clean.names=TRUE , ... ) {
+    require(MASS)
+    mu <- 0
+    if ( class(object)[1] %in% c("mer","bmer","glmerMod","lmerMod") ) {
+        mu <- fixef(object)
+    } else {
+        mu <- xcoef(object)
+    }
+    result <- as.data.frame( mvrnorm( n=n , mu=mu , Sigma=vcov(object) ) )
+    if ( clean.names==TRUE ) {
+        # convert (Intercept) to Intercept
+        for ( i in 1:ncol(result) ) {
+            if ( colnames(result)[i] == "(Intercept)" ) {
+                colnames(result)[i] <- "Intercept"
+            }
+        }
+    }
+    result
+}
+)

R/map.r

@@ -1,5 +1,18 @@
 # MAP - maximum a posteriori estimation
 
+# utility function to convert alist() construction with <- tagged linear model to regular ~ variety
+flist_untag <- function(flist) {
+    for ( i in 1:length(flist) ) {
+        if ( class(flist[[i]])=="<-" ) {
+            # tagged formula, so convert to ~ formula expression
+            flist[[i]][[1]] <- as.name("~")
+        }
+        flist[[i]] <- eval(flist[[i]])
+    }
+    as.list(flist)
+}
+
+# main event
 map <- function( flist , start , data , method="BFGS" , hessian=TRUE , debug=FALSE , ... ) {
 
     ########################################
@@ -19,6 +32,7 @@ map <- function( flist , start , data , method="BFGS" , hessian=TRUE , debug=FAL
     }
 
     flist.orig <- flist
+    flist <- flist_untag(flist)
 
     ########################################
     # check for common issues

R/map2stan-class.r

@@ -12,8 +12,10 @@ setMethod("coef", "map2stan", function(object) {
     object@coef
 })
 
-extract.samples <- function(object) {
-    p <- extract(object@stanfit)
+setMethod("extract.samples","map2stan",
+function(object) {
+    require(rstan)
+    p <- rstan::extract(object@stanfit)
     # get rid of dev and lp__
     p[['dev']] <- NULL
     p[['lp__']] <- NULL
@@ -23,6 +25,7 @@ extract.samples <- function(object) {
     }
     return(p)
 }
+)
 
 plotchains <- function(object , pars=names(object@start) , ...) {
     if ( class(object)=="map2stan" )
@@ -63,8 +66,15 @@ function (object, ...)
 
 DIC <- function( object , n=1000 ) {
     if ( class(object)=="map2stan") {
-        val <- attr(object,"DIC")
-        attr(val,"pD") <- attr(object,"pD")
+        if ( !is.null(attr(object,"DIC")) ) {
+            val <- attr(object,"DIC")
+            attr(val,"pD") <- attr(object,"pD")
+        } else {
+            # must compute
+            v <- DIC.map2stan(object)
+            val <- v[1]
+            attr(val,"pD") <- v[2]
+        }
     }
     if ( class(object)=="map" ) {
         post <- sample.qa.posterior( object , n=n )
@@ -82,6 +92,35 @@ DIC <- function( object , n=1000 ) {
     return(val)
 }
 
+DIC.map2stan <- function( object ) {
+    fit <- object@stanfit
+    # compute DIC
+    dev.post <- extract(fit, "dev", permuted = TRUE, inc_warmup = FALSE)
+    dbar <- mean( dev.post$dev )
+    # to compute dhat, need to feed parameter averages back into compiled stan model
+    post <- extract( fit )
+    Epost <- list()
+    for ( i in 1:length(post) ) {
+        dims <- length( dim( post[[i]] ) )
+        name <- names(post)[i]
+        if ( name!="lp__" & name!="dev" ) {
+            if ( dims==1 ) {
+                Epost[[ name ]] <- mean( post[[i]] )
+            } else {
+                Epost[[ name ]] <- apply( post[[i]] , 2:dims , mean )
+            }
+        }
+    }#i
+
+    # push expected values back through model and fetch deviance
+    #message("Taking one more sample now, at expected values of parameters, in order to compute DIC")
+    fit2 <- stan( fit=fit , init=list(Epost) , data=object@data , pars="dev" , chains=1 , iter=1 , refresh=-1 )
+    dhat <- as.numeric( extract(fit2,"dev") )
+    pD <- dbar - dhat
+    dic <- dbar + pD
+    return( c( dic , pD ) )
+}
+
 setMethod("show", "map2stan", function(object){
 
     cat("map2stan model fit\n")
@@ -154,6 +193,7 @@ resample <- function( object , iter=1e4 , warmup=1000 , chains=1 , cores=1 , ...
 
     result <- object
     result@stanfit <- fit
+    attr(result,"DIC") <- NULL # clear out any old DIC calculation
     return(result)
 }