rmcelreath / rethinking

Statistical Rethinking course and book package
2.16k stars 605 forks source link

compare "plot" method creates splom, not deviance plot (ref. p.228) #303

Open markdlindsay opened 3 years ago

markdlindsay commented 3 years ago

I cannot recreate the plot on page 228 of the text (2.ed) - appears the plot method in 'compare' is not being honoured as an splom is plotted instead.

image

HelenMylne commented 3 years ago

Not helpful for this particular plot, which I seem to have completely random as to whether it plots like the book or yours, but answering in case you have the same issues with the other plots of this type. If I just use plot(precis()) then I get one of these instead of the comparative plot. McElreath recommended I use precis_plot(precis()) and that works fine.

thomasajryan commented 3 years ago

in the code for the 'compare' function from rethinking, there's a 'setmethod' function for plotting (full code you need below). if you copy and paste that and run it individually, then re-run the plot(compare(m6.6, m6.7, m6.8) you get the full figure from the book that you're trying to replicate. This is a similar workaround I saw that someone used for a plot(coeftab.... issue they had in chapter 5. The precis_plot(precis( workaround mentioned above didn't work for me, so I tried this instead and it worked perfectly

setMethod("plot" , "compareIC" , function(x,y,xlim,SE=TRUE,dSE=TRUE,weights=FALSE,...) { dev_in <- x[[1]] - x[[5]]2 # criterion - penalty2 dev_out <- x[[1]] if ( !is.null(x[['SE']]) ) devSE <- x[['SE']] dev_out_lower <- dev_out - devSE dev_out_upper <- dev_out + devSE if ( weights==TRUE ) { dev_in <- ICweights(dev_in) dev_out <- ICweights(dev_out) dev_out_lower <- ICweights(dev_out_lower) dev_out_upper <- ICweights(dev_out_upper) } n <- length(dev_in) if ( missing(xlim) ) { xlim <- c(min(dev_in),max(dev_out)) if ( SE==TRUE & !is.null(x[['SE']]) ) { xlim[1] <- min(dev_in,dev_out_lower) xlim[2] <- max(dev_out_upper) } } main <- colnames(x)[1] set_nice_margins() dotchart( dev_in[n:1] , labels=rownames(x)[n:1] , xlab="deviance" , pch=16 , xlim=xlim , ... ) points( dev_out[n:1] , 1:n ) mtext(main)

standard errors

if ( !is.null(x[['SE']]) & SE==TRUE ) {
    for ( i in 1:n ) {
        lines( c(dev_out_lower[i],dev_out_upper[i]) , rep(n+1-i,2) , lwd=0.75 )
    }
}
if ( !all(is.na(x@dSE)) & dSE==TRUE ) {
    # plot differences and stderr of differences
    dcol <- col.alpha("black",0.5)
    abline( v=dev_out[1] , lwd=0.5 , col=dcol )
    diff_dev_lower <- dev_out - x$dSE
    diff_dev_upper <- dev_out + x$dSE
    if ( weights==TRUE ) {
        diff_dev_lower <- ICweights(diff_dev_lower)
        diff_dev_upper <- ICweights(diff_dev_upper)
    }
    for ( i in 2:n ) {
        points( dev_out[i] , n+2-i-0.5 , cex=0.5 , pch=2 , col=dcol )
        lines( c(diff_dev_lower[i],diff_dev_upper[i]) , rep(n+2-i-0.5,2) , lwd=0.5 , col=dcol )
    }
}

})