fitplotR - quickly assessing model recovery

fitplotR

This ggplot2-based function allows quickly assessing paramter recovery. It shows the point estimates, the corresponding confidence/credibility intervals and the true value. The deviation of the point estimates is highlighted. Different models (e.g. using different estimation algorithms or different prior specifications) can be compared.

The function is a straightforward call of ggplot2:

library(ggplot2)

fitplotr = function(est, lw, up, true) {

    if (is.null(dim(est))) {
        nvar = length(est)
        if (is.null(names(est))) {
            varnames = paste0("v", 1:nvar)
        } else {
            varnames = names(est)
        }
        nmodel = 1
        model = factor(1)
    } else {
        nvar = nrow(est)
        if (is.null(rownames(est))) {
            varnames = paste0("v", 1:nvar)
        } else {
            varnames = rownames(est)
        }
        nmodel = ncol(est)

        if (is.null(colnames(est))) {
            model = factor(1:nmodel)
        } else {
            model = colnames(est)
        }
    }


    d = data.frame(x = rep(varnames, nmodel + 1), model = c(rep(model, each = nvar), 
        rep("true", nvar)), y = c(as.vector(est), true), Lo = c(as.vector(lw), 
        rep(NA, nvar)), Hi = c(as.vector(up), rep(NA, nvar)))

    # computing the y-coordinates for the geom_segment part
    starts <- ends <- drop(matrix(d$y[d$model == "true"], 1, sum(d$model != 
        "true")))
    foo0 = d$y[d$model == "true"] > d$y[d$model != "true"]
    starts[foo0] = (d$y[d$model != "true"])[foo0]
    ends[!foo0] = (d$y[d$model != "true"])[!foo0]
    d$Lo2[d$model != "true"] = starts
    d$Hi2[d$model != "true"] = ends

    # plot the main part
    p = ggplot(d[d$model != "true", ], aes(x, y, color = model)) + geom_point(size = 4, 
        position = position_dodge(width = 0.3)) + geom_errorbar(aes(ymin = Lo, 
        ymax = Hi), width = 0.15, position = position_dodge(width = 0.3)) + 
        geom_linerange(aes(ymin = Lo2, ymax = Hi2), size = 2, position = position_dodge(width = 0.3))
    # plot the true value line
    p + geom_crossbar(aes(ymin = y, ymax = y), data = d[d$model == "true", ])

}

Example using simulated data:

estimates = cbind(c(3, 4, 7), c(3, 5, 7) + 1)
true = c(4, 5, 6)
lowbound = estimates - runif(6) * 4
upbound = estimates + runif(6) * 4
fitplotr(estimates, lowbound, upbound, true)

## ymax not defined: adjusting position using y instead

hview - quick preview of html output in the rstudio viewer pane

Sometimes, you want to quickly check what the html output of xtable or similar function will look like in your markdown code, but you don't want to open up a new markdown document. Right from the console, you can create a temp html file, put the output in and view it in the viewer pane. This way you can quickly look at model output in that pane, as opposed to on the console.

require(xtable)

## Loading required package: xtable

hview = function(htmllines) {
    htmlFile <- tempfile(fileext = ".html")
    writeLines(htmllines, htmlFile)
    rstudio::viewer(htmlFile)
}

data(tli)
fm2 <- lm(tlimth ~ sex * ethnicty, data = tli)
fm2.table <- xtable(fm2)

print(fm2.table, type = "html")

	Estimate	Std. Error	t value	Pr(&gt |t|)
(Intercept)	73.6364	4.2502	17.33	0.0000
sexM	-1.6364	5.8842	-0.28	0.7816
ethnictyHISPANIC	-9.7614	6.5501	-1.49	0.1395
ethnictyOTHER	15.8636	10.8360	1.46	0.1466
ethnictyWHITE	4.7970	4.9687	0.97	0.3368
sexM:ethnictyHISPANIC	10.6780	8.7190	1.22	0.2238
sexM:ethnictyWHITE	5.1230	7.0140	0.73	0.4670

hview(print(fm2.table, type = "html"))```

GHK Algorithm for Simulation of Multivariate Normal Rectangle Probabilities in RCPP

When starting new projects in RCPP, the biggest hurdle is to re-implement existing code snippets for standard tasks that are not (yet) part of RCPP. I have implementations for certain full-conditional draws coded up efficiently in vectorized R. I have implementations in plain C and C++ that were previously accessed in compiled form throug dyn.load or .C. I started to slowly transition more and more snippets to RCPP C++.
Here's a reimplementation of the GHK simulator from bayesm. Some minor mofidications were necessary. Now it is ready for implementation within larger RCPP projects.

Have fun writing MCMC samplers or doing simulated ML.

brackets [[ for RCPP

One cool feature of R is the ease of data handling. Lists, data frames and matrices can easily be accessed and indexed, making vector-oriented programming easy. When switching over to C/RCPP, things are a little different. Personally, I like the Armadillo library, since it gets me as close to R or Matlab as possible. One thing I have been missing though is quick access to matrix elements through brackets and logical statements. I have written some functions for illustration:

Try it out

Truncated Normal [RcppArmadillo]

Rcpp and especially RcppArmadillo with its nice syntax have developed into an extremely powerful tool for R. The concept of vectors and matrices is transported into C, and the extra effort in coding isn't that big anymore. RStudio even syntax highlights the C++ code and offers a 'source' button which runs the sourceCpp comand from the Rcpp package.
More tutorials are being written, but I just want to provide some simple examples. Let's start with a simple inverse CDF approach to drawing from a univariate truncated Normal distribution:

Then it is straightforward to source the code in R.