bug patch: plot(compare(.))

- patch to fix a plotting problem with compare's plot method.
- Added Rhat version number 4 to precis output
- fix log inversion bug with dlkjcorr

R/compare.r

@@ -122,7 +122,7 @@ compare <- function( ... , n=1e3 , sort="WAIC" , func=WAIC , WAIC=TRUE , refresh
 
 # 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,xlim,SE=TRUE,dSE=TRUE,weights=FALSE,...) {
-    dev_in <- x[[1]] - x[[2]]*2
+    dev_in <- x[[1]] - x[[5]]*2 # criterion - penalty*2
     dev_out <- x[[1]]
     if ( !is.null(x[['SE']]) ) devSE <- x[['SE']]
     dev_out_lower <- dev_out - devSE
@@ -137,7 +137,8 @@ setMethod("plot" , "compareIC" , function(x,y,xlim,SE=TRUE,dSE=TRUE,weights=FALS
     if ( missing(xlim) ) {
         xlim <- c(min(dev_in),max(dev_out))
         if ( SE==TRUE & !is.null(x[['SE']]) ) {
-            xlim <- c(min(dev_in),max(dev_out_upper))
+            xlim[1] <- min(dev_in,dev_out_lower)
+            xlim[2] <- max(dev_out_upper)
         }
     }
     main <- colnames(x)[1]

R/distributions.r

@@ -196,10 +196,10 @@ rlaplace <- function(n,location=0,lambda=1) {
 }
 
 
-# onion method correlation matrix
+# onion method correlation matrix - omits normalization constant
 dlkjcorr <- function( x , eta=1 , log=TRUE ) {
     ll <- det(x)^(eta-1)
-    if ( log==FALSE ) ll <- exp(ll)
+    if ( log==TRUE ) ll <- log(ll)
     return(ll)
 }
 

R/precis.r

@@ -150,6 +150,12 @@ precis_format <- function( result , depth , sort , decreasing ) {
         result <- result[o,]
     }
 
+    # label Rhat with version
+    rhat_col <- which( colnames(result)=="Rhat" )
+    if ( !is.null(rhat_col) ) {
+        colnames(result)[rhat_col] <- "Rhat4"
+    }
+
     return(result)
 }
 

man/precis.Rd

@@ -26,7 +26,7 @@ precis( model , depth=1 , pars , ci=TRUE , prob=0.89 ,
   Can also provide expected value, standard deviation, and HPDI columns for a data frame.
 }
 \value{
-  A data frame with a row for each parameter.
+  A data frame with a row for each parameter. The n_eff and Rhat columns are calculated by rstan. Rhat4 indicates Rhat as defined in Gelman et al 2013 "Bayesian Data Analysis". This is different from the classic 1992 Gelman and Rubin definition, as Rhat4 uses more information from multiple chains.
 }
 \references{}
 \author{Richard McElreath}

tests/rethinking_tests/test_ulam1.R

@@ -1,4 +1,4 @@
-# test_dir('rethinking/tests/ulam_tests',reporter="summary")
+# test_dir('rethinking/tests/rethinking_tests',filter="ulam1")
 
 context('ulam')
 library(rethinking)
@@ -58,4 +58,214 @@ test_that("Pareto regression post match",
     expect_equiv_eps( round(coef(mp),2) , 1.01 )
 )
 
+# beta
+
+test_that("ulam dbeta2 code",{
+    y <- rbeta2( 100 , 0.3 , 3 )
+    x <- rnorm( 100 )
+    z <- ulam(
+        alist(
+            y ~ dbeta2( p , theta ),
+            logit(p) <- a + b*x,
+            a ~ normal(0,1.5),
+            b ~ normal(0,0.5),
+            theta ~ exponential(1)
+        ), data=list(y=y,x=x) , sample=FALSE )
+    expect_known_hash( z$model , "f2210bdc47" )
+})
+
+test_that("ulam dbeta2 no prior error",{
+    y <- rbeta2( 100 , 0.3 , 3 )
+    expect_error(
+        z <- ulam(
+        alist(
+            y ~ dbeta2( p , theta ),
+            logit(p) <- a,
+            #a ~ normal(0,1.5),
+            theta ~ exponential(1)
+        ), data=list(y=y) , sample=TRUE )
+    , "Do all parameters have priors?" )
+})
+
+# binomial tests
+
+data( UCBadmit )
+UCBadmit$male <- as.integer( ifelse( UCBadmit$applicant.gender=="male" , 1 , 0 ) )
+UCBadmit$dept <- rep( 1:6 , each=2 )
+UCBadmit$ag <- UCBadmit$applicant.gender
+UCBadmit$applicant.gender <- NULL
+
+test_that("quap to ulam binomial",{
+    z <- quap(
+        alist(
+            admit ~ dbinom(applications,p),
+            logit(p) <- a,
+            a ~ dnorm(0,4)
+        ), data=UCBadmit )
+    zz <- ulam( z , sample=FALSE )
+    expect_known_hash( zz$model , "93f10cbc50" )
+})
+
+test_that("ulam constraints",{
+    z <- ulam(
+        alist(
+            admit ~ binomial(applications,p),
+            logit(p) <- a,
+            a ~ normal(0,4)
+        ),
+        constraints=list(a="lower=0"),
+        data=UCBadmit , sample=FALSE )
+    expect_known_hash( z$model , "ea06a8ca36" )
+})
+
+test_that("ulam binomial log_lik",{
+    z <- ulam(
+        alist(
+            admit ~ binomial(applications,p),
+            logit(p) <- a + b*male,
+            a ~ normal(0,1),
+            b ~ normal(0,0.5)
+        ), data=UCBadmit , sample=TRUE , log_lik=TRUE )
+    expect_known_hash( z@model , "097208dacb" )
+    expect_equivalent( length(WAIC(z)) , 4 )
+})
+
+# test to make sure constant '1' in binomial here isn't [i]'d in log_lik code
+data( chimpanzees )
+test_that("ulam binomial size 1 log_lik test",{
+    z <- ulam(
+        alist(
+            y ~ binomial(1,p),
+            logit(p) <- a,
+            a ~ normal(0,4)
+        ), data=list(y = chimpanzees$pulled_left) , sample=FALSE , log_lik=TRUE )
+    expect_known_hash( z$model , "1eb0e4e07f" )
+})
+
+# continuous missing data
+
+UCBadmit$x <- rnorm(12)
+UCBadmit$x[1:2] <- NA
+UCBadmit$x[2] <- NA
+
+# explicit code with merge_missing
+test_that("ulam continuous missing data 1",{
+    z2b2 <- ulam(
+        alist(
+            admit ~ binomial(applications,p),
+            logit(p) <- a + b*male + bx*x_merge,
+            x_merge ~ normal( 0 , 1 ),
+            x_merge <- merge_missing( x , x_impute ),
+            a ~ normal(0,4),
+            b ~ normal(0,1),
+            bx ~ normal(0,1)
+        ), data=UCBadmit , sample=TRUE , log_lik=FALSE )
+    expect_known_hash( z2b2@model , "ed1532565d" )
+    expect_equivalent( dim(precis(z2b2,2)) , c(5,6) )
+})
+
+# automated version
+# needs to build the merge code
+test_that("ulam continuous missing data 2",{
+    z2b_auto <- ulam(
+        alist(
+            admit ~ binomial(applications,p),
+            logit(p) <- a + b*male + bx*x,
+            x ~ normal( 0 , 1 ),
+            a ~ normal(0,4),
+            b ~ normal(0,1),
+            bx ~ normal(0,1)
+        ),
+        data=UCBadmit , sample=TRUE , log_lik=TRUE )
+    expect_known_hash( z2b_auto@model , "0656b0bec9" )
+    expect_equivalent( dim(precis(z2b_auto,2)) , c(5,6) )
+})
+
+UCBadmit$x2 <- rnorm(12)
+UCBadmit$x2[5:7] <- NA
+
+test_that("ulam continuous missing data 3",{
+    z2c <- ulam(
+        alist(
+            admit ~ binomial(applications,p),
+            logit(p) <- a + b*male + bx*x_merge + bx2*x2_merge,
+            x_merge ~ normal( 0 , 1 ),
+            x2_merge ~ normal( 0 , 1 ),
+            x_merge <- merge_missing( x , x_impute ),
+            x2_merge <- merge_missing( x2 , x2_impute ),
+            a ~ normal(0,4),
+            b ~ normal(0,1),
+            c(bx,bx2) ~ normal(0,1)
+        ), data=UCBadmit , sample=FALSE )
+    expect_known_hash( z2c$model , "db529f4ba3" )
+})
+
+
+# discrete missing data - need to mix over different terms
+data( UCBadmit )
+UCBadmit$male <- as.integer( ifelse( UCBadmit$applicant.gender=="male" , 1 , 0 ) )
+UCBadmit$dept <- rep( 1:6 , each=2 )
+UCBadmit$applicant.gender <- NULL
+UCBadmit$male2 <- UCBadmit$male
+UCBadmit$male2[1:2] <- (-9) # missingness code
+UCBadmit$male2 <- as.integer(UCBadmit$male2)
+
+test_that("ulam discrete missing 1",{
+    z <- ulam(
+        alist(
+            admit|male2==-9 ~ mixture(  
+                phi_male , 
+                binomial( applications , p_m1 ) , 
+                binomial( applications , p_m0 ) ),
+            admit|male2>-9 ~ binomial( applications , p ),
+            logit(p) <- a[dept] + b*male2,
+            logit(p_m1) <- a[dept] + b*1,
+            logit(p_m0) <- a[dept] + b*0,
+            male2|male2>-9 ~ bernoulli( phi_male ),
+            phi_male ~ beta(2,2),
+            a[dept] ~ normal(0,1),
+            b ~ normal(0,1)
+        ), data=UCBadmit , sample=TRUE , log_lik=FALSE )
+    expect_known_hash( z@model , "db97f3f2de" )
+})
+
+# same but with custom()
+test_that("ulam discrete missing 2",{
+    z2d <- ulam(
+        alist(
+            admit|male2==-1 ~ custom( log_mix( 
+                phi_male , 
+                binomial_lpmf(admit|applications,p_m1) , 
+                binomial_lpmf(admit|applications,p_m0) ) ),
+            admit|male2>-1 ~ binomial( applications , p ),
+            logit(p) <- a[dept] + b*male2,
+            logit(p_m1) <- a[dept] + b*1,
+            logit(p_m0) <- a[dept] + b*0,
+            male2|male2>-1 ~ bernoulli( phi_male ),
+            phi_male ~ beta(2,2),
+            a[dept] ~ normal(0,4),
+            b ~ normal(0,1)
+        ), data=UCBadmit , sample=FALSE , log_lik=FALSE )
+    expect_known_hash( z2d$model , "6ca98ae43e" )
+})
+
+# log_sum_exp form
+test_that("ulam discrete missing 3",{
+    z2e <- ulam(
+        alist(
+            admit|male2==-1 ~ custom( 
+                log_sum_exp( 
+                    log(phi_male) + binomial_lpmf(admit|applications,p_m1) , 
+                    log1m(phi_male) + binomial_lpmf(admit|applications,p_m0) ) ),
+            admit|male2>-1 ~ binomial( applications , p ),
+            logit(p) <- a + b*male2,
+            logit(p_m1) <- a + b*1,
+            logit(p_m0) <- a + b*0,
+            male2|male2>-1 ~ bernoulli( phi_male ),
+            phi_male ~ beta(2,2),
+            a ~ normal(0,4),
+            b ~ normal(0,1)
+        ), data=UCBadmit , sample=FALSE )
+    expect_known_hash( z2e$model , "eded9b0b6e" )
+})