Introduction

The R extension of the markdown language (Xie, Dervieux, and Riederer 2020) enables reproducible statistical reports with nice typesetting in HTML, Microsoft Word, and Latex. Since recently (R Core Team 2022, version 4.2), R's manual pages include support for mathematical expressions (Sarkar and Hornik 2022; Viechtbauer 2022), which already is a big improvement. However, rules for the mapping of built-in language elements to their mathematical representation are still lacking. So far, R expressions such as pbinom(k, N, p) are printed as they are; pretty mathematical formulae such as \(P_{\mathrm{Bi}}(X \le k; N, p)\) require explicit Latex commands, that is, P_{\mathrm{Bi}}\left(X \le k; N, p\right). Except for minimalistic use cases, these commands are tedious to type in and their source code is hard to read.

The present R package defines a set of rules for the automatic translation of R expressions to mathematical output in RMarkdown documents (Xie, Dervieux, and Riederer 2020) and ShinyApp webpages (Chang et al. 2022). The translation is done by an embedded Prolog interpreter that maps nested expressions recursively to MathML and LaTeX/MathJax, respectively. User-defined hooks enable to extend the set of rules, for example, to represent specific R elements by custom mathematical signs.

The main feature of the package is that the same R expressions and equations can be used for both mathematical typesetting and calculations. This saves time and reduces mistakes, as will be illustrated below.

The paper is organized as follows. Section 2 presents the technical background of the package, including the two main classes of Prolog rules for translating R objects to mathematical expressions. Section 3 illustrates the main features of the mathml package, potential issues and workarounds using examples from the day-to-day perspective of a user. Section 4 present a case study for the use of the package in a scientific report. Section 5 concludes with a discussion and ideas for further development.

Background

The translation of R expressions to mathematical output is achieved through a Prolog interpreter. Prolog is a classical logic programming language with many applications in expert systems, computer linguistics and symbolic artificial intelligence. The main strength of Prolog is its concise representation of facts and rules for the representation of knowledge and grammar, as well as its efficient built-in search engine for closed world domains. As it is well-known, R is a statistical programming language for data analysis and statistical modeling which is widely used in academia and industry. Besides the core library, a lot of packages have been developed for all kinds of statistical problems, including statistics-based artificial intelligence tools such as neural networks for machine learning and deep learning. Whereas Prolog is weak in statistical computation, but strong in symbolic manipulation, the converse may be said for the R language. The rolog package (Gondan 2022) bridges this gap by providing an interface to a SWI-Prolog distribution in R. The communication between the two systems is mainly in the form of queries from R to Prolog, but two predicates allow Prolog to ring back and evaluate terms in R.

For a first illustration of the mathml package, we consider the binomial probability.

term <- quote(pbinom(k, N, p))
term

## pbinom(k, N, p)

The term is quoted to avoid its immediate evaluation (which would raise an error anyway since the variables k, N, p have not yet been defined). Experienced readers will remember that the quoted expression above is a short form for

term <- call("pbinom", as.name("k"), as.name("N"), as.name("p"))

As is seen from the output above, the variable term is not assigned the result of the calculation, but an R “call” (see, e.g., Wickham 2019, for details on “non-standard evaluation”). This call can eventually be evaluated with eval(),

k <- 10
N <- 22
p <- 0.4
eval(term)

## [1] 0.77195

The R package mathml can now be used to render the call in MathML, that is the dialect for mathematical elements on HTML webpages or in MathJax/LaTeX, as shown below (some of the curly braces are not really needed in this simple example, but are necessary in edge cases).

library(mathml)
mathjax(term)

## [1] "${P}_{\\mathrm{Bi}}{\\left({{X}{\\le}{k}}{{;}{{N}{{,}{p}}}}\\right)}$"

We can include the output in a RMarkdown document by specifying results='asis' in the R code chunk, as is done in the next example. The R function mathout() is a wrapper that invokes mathml() for HTML output and mathjax() for LaTeX output.

mathout(term)

$P_{\text{Bi}}\left(X\le k;N,p\right)$

At the Prolog end, a predicate math/2 translates the call pbinom(K, N, Pi) into a “function” fn/2 with the name P_Bi, one argument X =< K, and the two parameters N and Pi.

math(pbinom(K, N, Pi), M)
 => M = fn(subscript('P', "Bi"), (['X' =< K] ; [N, Pi])).

Thus, the predicate math/2 could be considered a “macro” that translates a mathematical element (here, pbinom(K, N, Pi)) to a different mathematical element, namely fn(Name, (Args ; Pars)). The low-level predicate ml/3 is used to convert these basic elements to MathML.

ml(Flags, fn(Name, (Args ; Pars)), M)
 => ml(Flags, Name, N),
    ml(Flags, paren(list(op(;), [list(op(','), Args), list(op(','), Pars)])), X),
    M = mrow([N, mo(&(af)), X]).

The relevant rule for ml/3 builds the MathML entity mrow([N, mo(&(af)), X]), with N representing the name of the function and X its arguments and parameters, enclosed in parentheses. A corresponding rule jax/3 does the same for MathJax/LaTeX. A list of flags can be used for context-sensitive translation (see, e.g., the errors below).

Package mathml in practice

mathml is an R package for pretty mathematical representation of R functions and objects in data analysis, scientific reports and interactive web content. The currently supported features are listed below, roughly following the order proposed by (Murrell and Ihaka 2000).

term <- quote(1 + -2L + a + abc + "a" + phi + Phi + varphi + hat(b)[i, j]^2L)
mathout(term)

term <- quote(NaN + NA + TRUE + FALSE + Inf + (-Inf))
mathout(term)

term <- quote(a - ((b + c)) - d*e + f*(g + h) + i/j + k^(l + m) + (n*o)^{p + q})
mathout(term)

term <- quote(dot(a, b) + frac(1L, nodot(c, d + e)) + dfrac(1L, times(g, h)))
mathout(term)

term <- quote(a^(b + c))
paste(term)

## [1] "^"       "a"       "(b + c)"

term <- quote(mean(X) %+-% 1.96 * s / sqrt(N))
mathout(term)

term <- quote('%+-%'(mean(X), 1.96 * s / sqrt(N)))
term <- quote(mean(X) %+-% {1.96 * s / sqrt(N)})   # the same
mathout(term)

term <- quote(sin(x) + sin(x)^2L + cos(pi/2L) + tan(2L*pi) * expm1(x))
mathout(term)

term <- quote(choose(N, k) + abs(x) + sqrt(x) + floor(x) + ceiling(x))
mathout(term)

sgn <- function(x)
{
  if(x == 0L) return(0L)
  if(x < 0L) return(-1L)
  if(x > 0L) return(1L)
}

mathout(sgn)

mathout(call("<-", quote(sgn(x)), sgn))

term <- quote(S[Y]^2L <- frac(1L, N) * sum(Y[i] - mean(Y))^2L)
mathout(term)

term <- quote(log(subsupscript(prod(L[i]), i==1L, N)) <- 
          subsupscript(sum(log(L[i])), i==1L, N))
mathout(term)

term <- quote(integrate(sin, 0L, 2L*pi))
mathout(term)

eval(term)

term <- quote(integrate(lower=0L, upper=2L*pi, sin))
canonical(term)

## integrate(sin, 0L, 2L * pi)

mathout(canonical(term))

v <- 1:3
mathout(call("t", v))

A <- matrix(data=11:16, nrow=2, ncol=3)
B <- matrix(data=21:26, nrow=2, ncol=3)
term <- call("+", A, B)
mathout(term)

hook(m_A, mean(X)["A"]) ; hook(s2_A, s["A"]^2L) ;
hook(n_A, n["A"])
hook(m_B, mean(X)["B"]) ; hook(s2_B, s["B"]^2L)
hook(n_B, n["B"]) ; hook(s2_p, s["pool"]^2L)

term <- quote(t <- dfrac(m_A - m_B, 
    sqrt(denote(s2_p, frac((n_A - 1L)*s2_A + (n_B - 1L)*s2_B, n_A + n_B - 2L),
                "the pooled variance.") * (frac(1L, n_A) + frac(1L, n_B)))))
mathout(term)

m_A <- 1.5; s2_A <- 2.4^2; n_A <- 27; m_B <- 3.9; s2_B <- 2.8^2; n_B <- 20
print(eval(term))

## [1] -3.157427

t <- quote(dfrac(omit_right(mean(D) - mu[0L]), s / sqrt(N)))
mathout(t, flags=list(error='highlight'))

mathout(t, flags=list(error='fix'))

term <- quote(pbinom(successes, Ntotal, prob))
mathout(term)

hook(successes, k)
hook(quote(Ntotal), quote(N), quote=FALSE)
hook(prob, pi)
mathout(term)

A case study

This case study describes a model by Schwarz (1994) from mathematical Psychology using the features of package mathml. Schwarz presents a new explanation of redundancy gains that occur when observers respond to stimuli of different sources, and the same information is presented on two or more channels. In Schwarz’ model, decision-making builds on a process of noisy accumulation of information over time (e.g., Ratcliff et al. 2016). In redundant stimuli, the model assumes a superposition of channel-specific diffusion processes that eventually reach an absorbing barrier to elicit the response. For a detailed description the reader may refer to the original article.

Schwarz’ (1994) model refers to two stimuli A and B, presented either alone or in combination (AB, redundant stimuli), with the redundant stimuli being presented either simultaneously or with onset asynchrony \(\tau\). The channel activation is described as a two-dimensional Wiener process with drifts \(\mu_i\), variances \(\sigma^2_i\), and initial conditions \(X_i(t = 0) = 0, i = \mathrm{A, B}\). The channels may be correlated with \(\rho_{\mathrm{AB}}\), but we assume \(\rho_{\mathrm{AB}} = 0\) for simplicity.

A response is elicited when the process reaches an absorbing barrier \(c > 0\) for the first time. In single-target trials, the first passages at \(c\) are expected at

ED_single <- function(c, mu)
  dfrac(c, mu)

# display as E(D; mu), c is a scaling parameter
hook(ED_single(.C, .Mu), E(`;`(D, .Mu)))
mathout(call("=", quote(ED_single(c, mu)), ED_single))

$E\left(D;\mu \right)={\displaystyle \frac{c}{\mu }}$

In redundant stimuli, the activation of the channel-specific diffusion processes adds up, \(X_{\mathrm{AB}}(t) = X_{\mathrm A}(t) + X_{\mathrm B}(t)\), hence the name, superposition. \(X_{\mathrm{AB}}(t)\) is again a Wiener process with drift \(\mu_{\mathrm A} + \mu_{\mathrm B}\) and variance \(\sigma^2_{\mathrm A} + \sigma^2_{\mathrm B}\). For the expected first-passage time, we have

hook(mu_A, mu["A"])
hook(mu_B, mu["B"])
hook(sigma_A, sigma["A"])
hook(sigma_B, sigma["B"])
hook(mu_M, mu["M"])
hook(M, overline(X))

mathout(call("=", quote(E(D["AB"])), quote(ED_single(c, mu_A + mu_B))))

$E\left(D_{\text{AB}}\right)=E\left(D;{\mu }_{\text{A}}+{\mu }_{\text{B}}\right)$

For asynchronous stimuli, Schwarz [-Schwarz (1994); Eq. 10] derived the expected first-passage time as a function of the stimulus onset asyncrony \(\tau\),

ED_async <- function(tau, c, mu_A, sigma_A, mu_B)
{ dfrac(c, mu_A) + (dfrac(1L, mu_A) - dfrac(1L, mu_A + mu_B)) *
    ((mu_A*tau - c) * pnorm(dfrac(c - mu_A*tau, sqrt(sigma_A^2L*tau)))
      - (mu_A*tau + c) * exp(dfrac(2L*c*mu_A, sigma_A^2L))
        * pnorm(dfrac(-c - mu_A*tau, sqrt(sigma_A^2L*tau))))
}

hook(ED_async(.Tau, .C, .MA, .SA, .MB), E(`;`(D[.Tau], `,`(.MA, .SA, .MB))))
mathout(call("=", quote(ED_async(tau, c, mu_A, sigma_A, mu_B)), ED_async))

$E\left(D_{\tau };{\mu }_{\text{A}},{\sigma }_{\text{A}},{\mu }_{\text{B}}\right)={\displaystyle \frac{c}{{\mu }_{\text{A}}}}+\left({\displaystyle \frac{1}{{\mu }_{\text{A}}}}-{\displaystyle \frac{1}{{\mu }_{\text{A}}+{\mu }_{\text{B}}}}\right)\cdot \left[\left({\mu }_{\text{A}}\tau -c\right)\cdot \Phi \left({\displaystyle \frac{c-{\mu }_{\text{A}}\tau }{\sqrt{{\sigma }_{\text{A}}^2\tau }}}\right)-\left({\mu }_{\text{A}}\tau +c\right)\cdot \exp {\displaystyle \frac{2c{\mu }_{\text{A}}}{{\sigma }_{\text{A}}^2}}\cdot \Phi \left({\displaystyle \frac{-c-{\mu }_{\text{A}}\tau }{\sqrt{{\sigma }_{\text{A}}^2\tau }}}\right)\right]$

For negative onset asynchrony (i.e., B before A), the parameters need to be switched.

ED <- function(tau, c, mu_A, sigma_A, mu_B, sigma_B)
{
  if(tau == Inf) return(ED_single(c, mu_A)) ;
  if(tau == -Inf) return(ED_single(c, mu_B)) ;
  if(tau == 0L) return(ED_single(c, mu_A + mu_B)) ;
  if(tau > 0L) return(ED_async(tau, c, mu_A, sigma_A, mu_B)) ;
  if(tau < 0L) return(ED_async(abs(tau), c, mu_B, sigma_B, mu_A))
}

hook(ED(.Tau, .C, .MA, .SA, .MB, .SB), E(`;`(D[.Tau], `,`(.MA, .SA, .MB, .SB))))
mathout(call("=", quote(ED(tau, c, mu_A, sigma_A, mu_B, sigma_B)), ED))

$E\left(D_{\tau };{\mu }_{\text{A}},{\sigma }_{\text{A}},{\mu }_{\text{B}},{\sigma }_{\text{B}}\right)=\{\begin{array}{}E\left(D;{\mu }_{\text{A}}\right)\text{,}\text{if}\tau =\infty E\left(D;{\mu }_{\text{B}}\right)\text{,}\text{if}\tau =-\infty E\left(D;{\mu }_{\text{A}}+{\mu }_{\text{B}}\right)\text{,}\text{if}\tau =0E\left(D_{\tau };{\mu }_{\text{A}},{\sigma }_{\text{A}},{\mu }_{\text{B}}\right)\text{,}\text{if}\tau >0E\left(D_{\left|\tau \right|};{\mu }_{\text{B}},{\sigma }_{\text{B}},{\mu }_{\text{A}}\right)\text{,}\text{if}\tau <0\end{array}$

The observable response time is assumed to be the sum of \(D\), the time employed to reach the threshold for the decision, and a residual \(M\) denoting other processes such a motor preparation and execution. Correspondingly, the expected response time amounts to

ET <- function(tau, c, mu_A, sigma_A, mu_B, sigma_B, mu_M)
  ED(tau, c, mu_A, sigma_A, mu_B, sigma_B) + mu_M

hook(ET(.Tau, .C, .MA, .SA, .MB, .SB, .MM), 
     E(`;`(T[.Tau], `,`(.MA, .SA, .MB, .SB, .MM))))
mathout(call("=", quote(E(T[tau])), ET))

$E\left(T_{\tau }\right)=E\left(D_{\tau };{\mu }_{\text{A}},{\sigma }_{\text{A}},{\mu }_{\text{B}},{\sigma }_{\text{B}}\right)+{\mu }_{\text{M}}$

Schwarz (1994) applied the model to data from a redundant signals task with 13 onset asynchronies \(0, \pm33, \pm67, \pm100, \pm133, \pm167, \pm\infty\) ms, where \(0\) refers to the synchronous condition, and \(\pm\infty\) to the single-target presentations. Each condition was replicated 400 times. The observed mean response times and their standard deviations are given in Table 4 below.

Assuming that the model is correct, the observable mean reaction times follow an approximate Normal distribution around the model prediction \(E(T)\). We can, therefore, use a standard goodness-of-fit measure by \(z\)-standardisation.

z <- function(m, s, n, tau, c, mu_A, sigma_A, mu_B, sigma_B, mu_M)
  dfrac(m - denote(E, ET(tau, c, mu_A, sigma_A, mu_B, sigma_B, mu_M),
              "the expected mean response time"),
    s / sqrt(n))

mathout(call("=", quote(z[tau]), z))

$z_{\tau }={\displaystyle \frac{m-E}{s/\sqrt{n}}}$, with $E=E\left(T_{\tau };{\mu }_{\text{A}},{\sigma }_{\text{A}},{\mu }_{\text{B}},{\sigma }_{\text{B}},{\mu }_{\text{M}}\right)$ denoting the expected mean response time

The overall goodness-of-fit is the sum of the squared \(z\)-statistics for each onset asynchrony. Assuming again that the architecture of the model is correct, but the parameters are adjusted to the data, it follows a \(\chi^2(8~\mathrm{df})\)-distribution.

zv <- Vectorize(z, vectorize.args = c('m', 's', 'n', 'tau'))

gof <- function(par, tau, m, s, n)
  sum(zv(m, s, n, tau, c=100L, mu_A=par["mu_A"], sigma_A=par["sigma_A"], 
         mu_B=par["mu_B"], sigma_B=par["sigma_B"], mu_M=par["mu_M"])^2L)

hook(zv(.M, .S, .N, .Tau, .C, .MA, .SA, .MB, .SB, .MM), z[.Tau])
mathout(call("=", quote(X["8 df"]^2L), gof))

$X_{\text{8 df}}^2=\sum z_{\tau }^2$

The degrees of freedom result from the 13 predicted mean response times, minus 5 free parameters \(\theta = \langle\mu_{\mathrm A}, \sigma_{\mathrm A}, \mu_{\mathrm B}, \sigma_{\mathrm B}, \mu_{\mathrm M}\rangle\) (\(c\) is a redundant scaling parameter).

$\widehat{\theta }=\mathrm{argmin}\mathrm{gof}\left(\theta \right)$

The best fitting parameter values and their confidence intervals are given in Table 5 below.

The statistic goodness-of-fit value is \(X^2(8) = 30.54\), \(p = .0002\), indicating some lack of fit. Given the large trial numbers in the original study, some lack of fit is expected. For more detail, especially on fitting the observed standard deviations, the reader is referred to the original paper (Schwarz 1994).

Conclusions

This package expands the possibilities of R by allowing calculations and graphical rendering of the same terms. Building on current features of R and existing packages for displaying formulas in R (Murrell and Ihaka 2000; Allaire et al. 2018)), the mathml package bridges the gap between computational needs, presentation of results, and their reproducibility. Two possible modes for typesetting are available to export R Markdown documents in different formats or ShinyApp webpages: MathML and MathJax. A further, not-yet-mentioned, use of mathml is related to presentations made using RMarkdown. By adopting mathml researchers or teachers can use the RMarkdown to both conduct analysis and show results, smoothing the process and uniforming the graphical output.

The case study presented above shows that mathml can help to improve data analyses and statistical reports from an aesthetical perspective, as well as regarding reproducibility of research. Furthermore, the package may also allow for a better detection of possible mistakes in R programs.

Similar to most programming languages (Green 1977), R code is notoriously hard to read, and the poor legibility of the language is one of the main sources of mistakes. For illustration, we consider again Equation 10 in Schwarz (1994).

f1 <- function(tau)
{ dfrac(c, mu["A"]) + (dfrac(1L, mu["A"]) - dfrac(1L, mu["A"] + mu["B"]) * 
    ((mu["A"]*tau - c) * pnorm(dfrac(c - mu["A"]*tau, sqrt(sigma["A"]^2L*tau)))
      - (mu["A"]*tau + c) * exp(dfrac(2L*mu["A"]*tau, sigma["A"]^2L))
        * pnorm(dfrac(-c - mu["A"]*tau, sqrt(sigma["A"]^2L*tau)))))
}

mathout(f1)

${\displaystyle \frac{c}{{\mu }_{\text{A}}}}+\left\{{\displaystyle \frac{1}{{\mu }_{\text{A}}}}-{\displaystyle \frac{1}{{\mu }_{\text{A}}+{\mu }_{\text{B}}}}\cdot \left[\left({\mu }_{\text{A}}\tau -c\right)\cdot \Phi \left({\displaystyle \frac{c-{\mu }_{\text{A}}\tau }{\sqrt{{\sigma }_{\text{A}}^2\tau }}}\right)-\left({\mu }_{\text{A}}\tau +c\right)\cdot \exp {\displaystyle \frac{2{\mu }_{\text{A}}\tau }{{\sigma }_{\text{A}}^2}}\cdot \Phi \left({\displaystyle \frac{-c-{\mu }_{\text{A}}\tau }{\sqrt{{\sigma }_{\text{A}}^2\tau }}}\right)\right]\right\}$

The first version has a wrong parenthesis, which is barely visible in the code, whereas the wrong curly brace in the mathematical representation is immediately obvious. The correct version is shown below for comparison.

f2 <- function(tau)
{ dfrac(c, mu["A"]) + (dfrac(1L, mu["A"]) - dfrac(1L, mu["A"] + mu["B"])) * 
    ((mu["A"]*tau - c) * pnorm(dfrac(c - mu["A"]*tau, sqrt(sigma["A"]^2L*tau)))
      - (mu["A"]*tau + c) * exp(dfrac(2L*mu["A"]*tau, sigma["A"]^2L))
        * pnorm(dfrac(-c - mu["A"]*tau, sqrt(sigma["A"]^2L*tau))))
}

mathout(f2)

${\displaystyle \frac{c}{{\mu }_{\text{A}}}}+\left({\displaystyle \frac{1}{{\mu }_{\text{A}}}}-{\displaystyle \frac{1}{{\mu }_{\text{A}}+{\mu }_{\text{B}}}}\right)\cdot \left[\left({\mu }_{\text{A}}\tau -c\right)\cdot \Phi \left({\displaystyle \frac{c-{\mu }_{\text{A}}\tau }{\sqrt{{\sigma }_{\text{A}}^2\tau }}}\right)-\left({\mu }_{\text{A}}\tau +c\right)\cdot \exp {\displaystyle \frac{2{\mu }_{\text{A}}\tau }{{\sigma }_{\text{A}}^2}}\cdot \Phi \left({\displaystyle \frac{-c-{\mu }_{\text{A}}\tau }{\sqrt{{\sigma }_{\text{A}}^2\tau }}}\right)\right]$

As the reader may know from own experience, missed parentheses are frequent causes of wrong results and errors that are hard to locate in programming code. This particular example shows that mathematical rendering can help to substantially reduce the amount of careless errors in programming.

In its current version mathml has some limitations. For example, it is currently not possible to use the functions of mathml for writing intratextual formulas, here, the user has to adopt the usual LateX notation. However, we do not consider this to be a serious limitation, since the most productive use of mathml still comes from more complicated calculations.

Also: In long equations, we did not yet find a convenient way to insert linebreaks. This is mostly because it is not properly supported by MathML and LaTeX renderers. In its current stage, the LaTeX package breqn (Robertson et al. 2021) is mostly a proof of concept.

mathml is available for R version (…) and later, and can be easily installed using the usual install.packages("mathml"). The source code of the package is found at https://github.com/mgondan/mathml.