---
title: "mathml: Translate R Expressions to MathML and LaTeX"
date: "2023-11-20"
abstract: >
  This R\ package translates R\ objects to suitable elements in MathML or LaTeX,
  thereby allowing for a pretty mathematical representation of R\ objects and
  functions in data analyses, scientific reports and interactive web content. In
  the R\ Markdown document rendering language, R\ code and mathematical content
  already exist side-by-side. The present package enables use of the same
  R\ objects for both data analysis and typesetting in documents or web content.
  This tightens the link between the statistical analysis and its verbal
  description or symbolic representation, which is another step towards
  reproducible science. User-defined hooks enable extension of the package by
  mapping specific variables or functions to new MathML and LaTeX entities. 
  Throughout the paper, examples are given for the functions of the package, and
  a case study illustrates its use in a scientific report.
draft: true
author:  
  - name: Matthias Gondan
    affiliation: Universität Innsbruck
    address:
    - Department of Psychology
    - Innsbruck, Austria
    url: https://www.uibk.ac.at/psychologie/mitarbeiter/gondan-rochon/index.html.en
    orcid: 0000-0001-9974-0057
    email:  Matthias.Gondan-Rochon@uibk.ac.at
  - name: Irene Alfarone
    url: https://www.uibk.ac.at/psychologie/mitarbeiter/alfarone/index.html.en
    email: Irene.Alfarone@uibk.ac.at
    orcid: 0000-0002-8409-8900
    affiliation: Universität Innsbruck
    address:
    - Department of Psychology
    - Innsbruck, Austria
type: package
output: 
  rjtools::rjournal_article:
    self_contained: yes
    toc: no
bibliography: bibliography.bib
preamble:
- \usepackage{cancel}
- \usepackage{amsmath}
---

```{r setup, include=FALSE}
library(mathml)
```

# Introduction

The R\ extension of the markdown language [@Xie2018;@rmarkdown] enables
reproducible statistical reports with nice typesetting in HTML, Microsoft Word,
and LaTeX. Moreover, since recently [@R, version 4.2], R\'s manual pages 
include support for mathematical expressions [@Sarkar2022;@Viechtbauer2022],
which is already a big improvement. However, except for special cases 
such as regression models [@equatiomatic] and R's own plotmath annotation, 
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 and pretty mathematical formulae 
such as \(P_{\mathrm{Bi}}(X \le k; N, p)\) require explicit LaTeX commands
like `P_{\mathrm{Bi}}\left(X \le k; N, p\right)`. Except for very basic use
cases, these commands are tedious to type 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 R\ Markdown documents [@Xie2020] and
Shiny Apps [@Chang2022]. The translation is done by an embedded Prolog
interpreter that maps nested expressions recursively to MathML and
LaTeX/MathJax, respectively. User-defined hooks enable extension of 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 potentially reduces mistakes, as will be illustrated below. Readers should
have basic knowledge of knitr and R\ Markdown to be able to follow this
article [@knitr;@rmarkdown], while to extend and customize the package, 
some basic knowledge of Prolog is needed.

The paper is organized as follows. We start with a description of the
technical background of the package, including the two main classes of rules for
translating R\ objects to mathematical expressions. The next section illustrates
the main features of the \CRANpkg{mathml} package, potential issues and their
workarounds using examples from the day-to-day perspective of a user. A case
study follows with a scientific report written with the help of the package. The
last section concludes with a discussion and ideas for further development.

# Background

Similar to other high-level programming languages, R is homoiconic, that is,
R\ commands (i.e., R\ "calls") are, themselves, symbolic data structures that
can be created, parsed and modified. Because the default response of the
R\ interpreter is to evaluate a call and return its result, this property is not
transparent to the general user. There exists, however, a number of built-in
R\ functions (e.g., `quote()`, `call()` etc.) that allow the user to create R\ calls
which can be stored in regular variables and then, for example, evaluated at
a later stage or in a specific environment [@Wickham2019]. The present package
includes a set of rules that translate such calls to a mathematical
representation in MathML and LaTeX. For a first illustration of the \pkg{mathml}
package, we consider the binomial probability.

````{r}
term <- quote(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
R\ users will recognize that the expression is a short form for

````{r}
term <- call("pbinom", as.name("k"), as.name("N"), as.name("p"))
term
````

As can be seen from the output, to the variable `term` is not assigned the 
result of the calculation, but instead an R\ call [see, e.g., @Wickham2019, 
for details on "non-standard evaluation"], which can eventually be evaluated 
with `eval()`,

````{r}
k <- 10
N <- 22
p <- 0.4
eval(term)
````

The R\ package \pkg{mathml} can now be used to render the call in MathML or in
MathJax/LaTeX. MathML is the dialect for mathematical elements on HTML webpages,
whereas LaTeX is typically used for typesetting printed documents, as shown
below.

````{r}
library(mathml)
substr(mathml(term), 1, 70)
mathjax(term)
````

Some of the curly braces are not really needed in the LaTeX output, but are
necessary in edge cases. The package also includes a function `mathout()` that
wraps a call to `mathml()` for HTML output and `mathjax()` for LaTeX output.
Moreover, the function `math(x)` adds the class `"math"` to its argument, such
that a special knitr printing function is
invoked [see the vignette on custom print methods in @knitr]. An R\ Markdown
code chunk with `mathout(term)` thus produces:

````{r, echo=FALSE}
math(term)
````

Similarly, `inline()` produces inline
output, `` `r "\x60r inline(term)\x60"` `` yields `r inline(term)`.

# Package \pkg{mathml} in practice

The currently supported R\ objects are listed below, roughly following the order
proposed by Murrell and Ihaka [-@murrell2000].

## Basic elements

\pkg{mathml} handles the basic elements of everyday mathematical expressions,
such as integers, floating-point numbers, Latin and Greek letters, multi-letter
identifiers, accents, subscripts, and superscripts.
 
```{r}
term <- quote(1 + -2L + a + abc + "a" + phi + Phi + varphi + roof(b)[i, j]^2L)
math(term)

term <- quote(round(3.1415, 3L) + NaN + NA + TRUE + FALSE + Inf + (-Inf))
math(term)
```

An expression such as `1 + -2` may be considered aesthetically unsatisfactory.
It is correct R\ syntax, though, and is reproduced accordingly, without the
parentheses. Parentheses around negative numbers or symbols can be added as
shown above for `+ (-Inf)`.

To avoid name clashes with package `stats`, `roof()` is used to put a hat
on a symbol (see next section for further decorations). Note that an 
R\ function `roof()` does not exist in base R\, it is provided by the package
for convenience and points to the identity function.

## Decorations

The package offers some support for different fonts as well as accents and
boxes etc. Internally, these decorations are implemented as identity functions,
so that they can be introduced into R\ expressions without side-effects.

```{r}
term <- quote(bold(b[x, 5L]) + bold(b[italic(x)]) + italic(ab) + italic(42L))
math(term)

term <- quote(tilde(a) + mean(X) + box(c) + cancel(d) + phantom(e) + prime(f))
math(term)
```

Note that the font styles only affect the display of identifiers, whereas
numbers, character strings etc. are left untouched.

## Operators and parentheses

Arithmetic operators and parentheses are translated as they are, as illustrated
below.

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

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

For multiplications involving only numbers and symbols, the multiplication sign
is omitted. This heuristic does not always produce the desired result;
therefore, \pkg{mathml} defines alternative R\ functions `dot()`, `nodot()`,
and `times()`. These functions calculate a product and produce the respective
multiplication signs. Similarly, `frac()` and `dfrac()` can be used for small
and large fractions.

For standard operators with known precedence, \pkg{mathml} is generally able to
detect if parentheses are needed; for example, parentheses are automatically
placed around `d + e` in the `nodot`-example. However, we note unnecessary
parentheses around `l + m` above. These parentheses are a consequence
of `quote(a^(b + c))` actually producing a nested R\ call of the
form `'^'(a, (b + c))` instead of `'^'(a, b + c)`:

```{r}
term <- quote(a^(b + c))
paste(term)
```

For the present purpose, this R\ feature is unfortunate because the extra
parentheses around `b + c` are not needed. The preferred result is obtained by
the functional form `quote('^'(k, l + m))` of the power, or curly braces as a
workaround (see `p + q` above).

## Custom operators

Whereas in standard infix operators, the parentheses typically follow the rules
for precedence, undesirable results may be obtained in custom operators.

```{r}
term <- quote(mean(X) %+-% 1.96 * s / sqrt(N))
math(term)
term <- quote('%+-%'(mean(X), 1.96 * s / sqrt(N))) # functional form of '%+-%'
term <- quote(mean(X) %+-% {1.96 * s / sqrt(N)})   # the same
math(term)
```

The example is a reminder that it is not possible to define the precedence of
custom operators in R\, and that expressions with such operators are evaluated
strictly from left to right. Again, the problem can be worked around by the
functional form of the operator or a curly brace to hide the parenthesis, and,
at the same time, enforce the correct operator precedence.

More operators are shown in Table \@ref(tab:custom-operators), including the
suggestions by Murrell and Ihaka [-@murrell2000] for graphical annotations and
arrows in R\ figures.

```{r custom-operators, echo=FALSE}
escape <- function(s)
{
  if(knitr::is_latex_output())
  {
    s = lapply(s, FUN=gsub, pattern="\\\\ ", replacement=" ")
    return(sapply(s, FUN=knitr:::escape_latex))
  }
  
  if(knitr::is_html_output())
    return(sapply(s, FUN=knitr:::escape_html))
  
  warning("escape: no knitr output specified")
  sapply(s, FUN=knitr:::escape_html)
}

op1 <- list(
  "A\\ %*%\\ B"=quote(A %*% B),
  "A\\ %.%\\ B"=quote(A %.% B),
  "A\\ %x%\\ B"=quote(A %x% B),
  "A\\ %/%\\ B"=quote(A %/% B),
  "A\\ %%\\ B"=quote(A %% B),
  "A\\ & B"=quote(A & B),
  "A\\ | B"=quote(A | B),
  "xor(A,\\ B)"=quote(xor(A, B)),
  "!A"=quote(!A),
  "A\\ ==\\ B"=quote(A == B),
  "A\\ <-\\ B"=quote(A <- B))

m1 <- lapply(op1, FUN=mathout, flags=list(cat=FALSE))

op1 <- names(op1)
op1 <- escape(op1)

op2 <- list(
  "A\\ !=\\ B"=quote(A != B),
  "A\\ ~\\ B"=quote(A ~ B),
  "A\\ %~~%\\ B"=quote(A %~~% B),
  "A\\ %==%\\ B"=quote(A %==% B),
  "A\\ %=~%\\ B"=quote(A %=~% B),
  "A\\ %prop%\\ B"=quote(A %prop% B),
  "A\\ %in%\\ B"=quote(A %in% B),
  "intersect(A,\\ B)"=quote(intersect(A, B)),
  "union(A,\\ B)"=quote(union(A, B)),
  "crossprod(A,\\ B)"=quote(crossprod(A, B)),
  "is.null(A)"=quote(is.null(A)))

m2 <- lapply(op2, FUN=mathout, flags=list(cat=FALSE))

op2 <- names(op2)
op2 <- escape(op2)

op3 <- list(
  "A\\ %<->%\\ B"=quote(A %<->% B),
  "A\\ %->%\\ B"=quote(A %->% B),
  "A\\ %<-%\\ B"=quote(A %<-% B),
  "A\\ %up%\\ B"=quote(A %up% B),
  "A\\ %down%\\ B"=quote(A %down% B),
  "A\\ %<=>%\\ B"=quote(A %<=>% B),
  "A\\ %=>%\\ B"=quote(A %=>% B),
  "A\\ %<=%\\ B"=quote(A %<=% B),
  "A\\ %dblup%\\ B" =quote(A %dblup% B),
  "A\\ %dbldown%\\ B"=quote(A %dbldown% B),
  " "="")

m3 <- lapply(op3, FUN=mathout, flags=list(cat=FALSE))

op3 <- names(op3)
op3 <- escape(op3)

t <- cbind(Operator=op1, Output=m1,
  Operator=op2, Output=m2,
  Operator=op3, Arrow=m3)

knitr::kable(t, caption="Custom operators in mathml",
  row.names=FALSE, escape=FALSE)
```



```{r base-stats, echo=FALSE}
op1 <- list(
  "sin(x)"=quote(sin(x)),
  "cosh(x)"=quote(cosh(x)),
  "tanpi(alpha)"=quote(tanpi(alpha)),
  "asinh(x)"=quote(asinh(x)),
  "log(p)"=quote(log(p)),
  "log1p(x)"=quote(log1p(x)),
  "logb(x,\\ e)"=quote(logb(x, e)),
  "exp(x)"=quote(exp(x)),
  "expm1(x)"=quote(expm1(x)),
  "choose(n,\\ k)"=quote(choose(n, k)),
  "lchoose(n,\\ k)"=quote(lchoose(n, k)),
  "factorial(n)"=quote(factorial(n)), 
  "lfactorial(n)"=quote(lfactorial(n)),
  "sqrt(x)"=quote(sqrt(x)),
  "mean(X)"=quote(mean(X)),
  "abs(x)"=quote(abs(x)))
m1 <- lapply(op1, FUN=mathout, flags=list(cat=FALSE))

op1 <- names(op1)
op1 <- escape(op1)

op2 <- list(
  "dbinom(k,\\ N,\\ pi)"=quote(dbinom(k, N, pi)),
  "pbinom(k,\\ N,\\ pi)"=quote(pbinom(k, N, pi)),
  "qbinom(p,\\ N,\\ pi)"=quote(qbinom(p, N, pi)),
  "dpois(k,\\ lambda)"=quote(dpois(k, lambda)),
  "ppois(k,\\ lambda)"=quote(ppois(k, lambda)),
  "qpois(p,\\ lambda)"=quote(qpois(p, lambda)),
  "dexp(x,\\ lambda)"=quote(dexp(x, lambda)),
  "pexp(x,\\ lambda)"=quote(pexp(x, lambda)),
  "qexp(p,\\ lambda)"=quote(qexp(p, lambda)),
  "dnorm(x,\\ mu,\\ sigma)"=quote(dnorm(x, mu, sigma)),
  "pnorm(x,\\ mu,\\ sigma)"=quote(pnorm(x, mu, sigma)),
  "qnorm(alpha/2L)"=quote(qnorm(alpha/2L)),
  "1L\\ -\\ pchisq(x,\\ 1L)"=quote(1L - pchisq(x, 1L)),
  "qchisq(1L\\ -\\ alpha,\\ 1L)"=quote(qchisq(1L-alpha, 1L)),
  "pt(t,\\ N\\ -\\ 1L)"=quote(pt(t, N-1L)),
  "qt(alpha/2L,\\ N\\ -\\ 1L)"=quote(qt(alpha/2L, N-1L)))
m2 <- lapply(op2, FUN=mathout, flags=list(cat=FALSE))

op2 <- names(op2)
op2 <- escape(op2)

t <- cbind(Function=op1, Output=m1, Function=op2, Output=m2)
knitr::kable(t, caption="R functions from base and stats",
  row.names=FALSE, escape=FALSE)
```

## Builtin functions

There is support for most functions from package `base`, with adequate use and
omission of parentheses.

```{r}
term <- quote(sin(x) + sin(x)^2L + cos(pi/2L) + tan(2L*pi) * expm1(x))
math(term)
term <- quote(choose(N, k) + abs(x) + sqrt(x) + floor(x) + exp(frac(x, y)))
math(term)
```

A few more examples are shown in Table \@ref(tab:base-stats), including
functions from `stats`. 

## Custom functions

For self-written functions, the matter is somewhat more complicated. For instance, 
if we consider a function such as `g <- function(...) ...`, the name _g_ is not 
transparent to R, because only the function body is represented. 
We can still display functions in the form `head(x) = body` if we embed the object to be shown into a call `"<-"(head, body)`.

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

math(sgn)
math(call("<-", quote(sgn(x)), sgn))
```

The function body is generally a nested R\ call of the form `'{'(L)`, with `L`
being a list of commands (the semicolon, not necessary in R, is translated to a
newline). The example also illustrates that \pkg{mathml} provides limited
support for control structures such as `if` that is internally represented
as `if(condition, action)`.

## Indices and powers

Indices in square brackets are rendered as subscripts, powers are rendered as
superscript. Moreover, \pkg{mathml} defines the
functions `sum_over(x, from, to)`, and `prod_over(x, from, to)` that simply
return their first argument. The other two arguments serve as
decorations (_to_ is optional), for example, for summation and product signs.

```{r}
term <- quote(S[Y]^2L <- frac(1L, N) * sum(Y[i] - mean(Y))^2L)
math(term)

term <- quote(log(prod_over(L[i], i==1L, N)) <- sum_over(log(L[i]), i==1L, N))
math(term)
```

## Ringing back to R

R\'s `integrate` function takes a number of arguments, the most important ones
being the function to integrate, and the lower and the upper bound of the
integration.

```{r}
term <- quote(integrate(sin, 0L, 2L*pi))
math(term)
eval(term)
```

For mathematical typesetting in the form
of \(\int f(x)\, dx\), \pkg{mathml} needs to find out the name of the
integration variable. For that purpose, the underlying Prolog bridge provides a
predicate `r_eval/2` that calls R\ from Prolog. This predicate is used to
evaluate `formalArgs(args(sin))` and returns the names of the arguments
of `sin()`, namely, _x_.

Above, the quoted term is an abbreviation
for `call("integrate", quote(sin), ...)`, with `sin` being an R\ symbol, not a
function. While the R\ function `integrate()` can handle both symbols and
functions, \pkg{mathml} needs the symbol because it is unable to determine the
function name of custom functions.

## Names and order of arguments

One of R's great features is the possibility to refer to function arguments by
their names, not only by their position in the list of arguments. At the other
end, the Prolog handlers for R\ calls are rather rigid, for
example, `integrate/3` accepts exactly three arguments in a particular order and
without names, that is, `integrate(lower=0L, upper=2L*pi, sin)`, would not print
the desired result.

To "canonicalize" function calls with named arguments and arguments in unusual
order, \pkg{mathml} provides an auxiliary R\ function `canonical(f, drop)` that
reorders the argument list of calls to known R\ functions and, 
if `drop=TRUE` (which is the default), also removes the names of the arguments.

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

```{r}
math(canonical(term))
```

This function can be used to feed mixtures of partially named and positional
arguments into the renderer. For details, see the R\ function `match.call()`.

## Matrices and Vectors

Of course, \pkg{mathml} also supports matrices and vectors.

```{r}
v <- 1:3
math(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)
math(term)
```

Note that the seemingly more convenient `term <- quote(A + B)` yields \(A + B\)
in the output---instead of the desired matrix representation. This behavior is
expected because quotation of R\ calls also quote the components of the
call (here, _A_ and _B_).

## Short mathematical names for R symbols

In typical R\ functions, variable names are typically longer than just single
letters, which may yield unsatisfactory results in the mathematical output.

```{r}
term <- quote(pbinom(successes, Ntotal, prob))
math(term)

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

To improve the situation, \pkg{mathml} provides a simple hook that can be used
to replace elements (e.g., verbose variable names) of the code by concise
mathematical symbols, as illustrated in the example. To simplify notation,
`hook()` uses non-standard evaluation of its arguments. If the `quote` flag
of `hook()` is set to `FALSE`, the user has to provide the quoted expressions.
Care should be taken to avoid recursive hooks such as `hook(s, s["A"])` that
endlessly replace the \(s\) from \(s_{\mathrm{A}}\) as
in \(s_{\mathrm{A}_{\mathrm{A}_{\mathrm{A}\cdots}}}\).

The hooks can also be used for more complex elements such as R\ calls, with
dotted symbols representing Prolog variables.

```{r}
hook(pbinom(.K, .N, .P), sum_over(dbinom(i, .N, .P), i=0L, .K))
math(term)
```

Further customization requires the assertion of new Prolog rules `math/2`,
`ml/3`, `jax/3`, as shown in the Appendix.

## Abbreviations

We consider the \(t\)-statistic for independent samples with equal variance. To
avoid clutter in the equation, the pooled variance \(s^2_{\mathrm{pool}}\) is
abbreviated, and a comment is given with the expression
for \(s^2_{\mathrm{pool}}\). For this purpose, \pkg{mathml} provides a
function `denote(abbr, expr, info)`, with `expr` actually being
evaluated, `abbr` being rendered, plus a comment of the
form "with `expr` denoting `info`".

```{r}
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)))))
math(term)
```

The term is evaluated below. `print()` is needed because the return value of an
assignment of the form `t <- dfrac(...)` is not visible in R.

```{r}
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))
```

## Context-dependent rendering

Consider an educational scenario in which we want to highlight a certain
element of a term, for example, that a student has forgotten to subtract the
null hypothesis in a \(t\)-ratio:

```{r}
t <- quote(dfrac(omit_right(mean(D) - mu[0L]), s / sqrt(N)))
math(t, flags=list(error="highlight"))
math(t, flags=list(error="fix"))
```

The R\ function `omit_right(a + b)` uses non-standard evaluation
techniques [e.g., @Wickham2019] to return only the left part an operation,
and cancels the right part. This may not always be desired, for example, when
illustrating how to fix the mistake.

For this purpose, the functions `mathml()`, `mathjax()`, `mathout()`
and `math()` have an optional argument `flags` which is a list with named
elements. In this example, we use this argument to tell \pkg{mathml} how to
render such erroneous expressions using the flag `error` which can be "asis",
"highlight", "fix", or "ignore". For more examples,
see Table \@ref(tab:mistakes).

```{r mistakes, echo=FALSE}
op1 <- list(
  "omit_left(a\\ +\\ b)"=quote(omit_left(a + b)),
  "omit_right(a\\ +\\ b)"=quote(omit_right(a + b)),
  "list(quote(a),\\ quote(omit(b)))"=list(quote(a), quote(omit(b))),
  "add_left(a\\ +\\ b)"=quote(add_left(a + b)),
  "add_right(a\\ +\\ b)"=quote(add_right(a + b)),
  "list(quote(a),\\ quote(add(b)))"=list(quote(a), quote(add(b))),
  "instead(a,\\ b)\\ +\\ c"=quote(instead(a, b) + c))

asis <- lapply(op1, FUN=mathout, flags=list(cat=FALSE, error="asis"))
high <- lapply(op1, FUN=mathout, flags=list(cat=FALSE, error="highlight"))
fix  <- lapply(op1, FUN=mathout, flags=list(cat=FALSE, error="fix"))
igno <- lapply(op1, FUN=mathout, flags=list(cat=FALSE, error="ignore"))

op1 <- names(op1)
op1 <- escape(op1)

t <- cbind(Operation=op1, 
  "error\\ =\\ asis"=asis, highlight=high, fix=fix, ignore=igno)
knitr::kable(t, caption="Highlighting elements of a term",
  row.names=FALSE, escape=FALSE)
```

# A case study

This case study describes a model by Schwarz [-@schwarz1994] from mathematical
psychology using the features of package \pkg{mathml}. Schwarz [-@schwarz1994]
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's [-@schwarz1994] model, decision-making builds on a
process of noisy accumulation of information over time [e.g., @ratcliff2016]. 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's [-@schwarz1994] 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 buildup of 
channel-specific activation 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 of \(c\) are
expected at

```{r}
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)))
math(call("=", quote(ED_single(c, mu)), ED_single))
```

One would typically use chunk option `echo=FALSE` to suppress the R\ code.

In redundant stimuli, the activation from 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

```{r}
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))

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

For asynchronous stimuli, Schwarz [-@schwarz1994] derived the expected
first-passage time as a function of the stimulus onset asynchrony \(\tau\),

```{r}
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))))
math(call("=", quote(E(D[tau])), ED_async))
```

For negative onset asynchrony (i.e., B before A), the parameters are simply
switched.

```{r}
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))))
math(call("=", quote(ED(tau, c, mu_A, sigma_A, mu_B, sigma_B)), ED))
```

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 as motor preparation and execution. Correspondingly, the
expected response time amounts to

```{r}
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))))
math(call("=", quote(E(T[tau])), ET))
```

Schwarz [-@schwarz1994] applied the model to data from a redundant signals 
task [@miller] with 13 onset
asynchronies \(0, \pm33, \pm67, \pm100, \pm133, \pm167, \pm\infty\) ms,
where \(\tau = 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 \@ref(tab:miller-data).

```{r miller-data, echo=FALSE}
tau  <- list(-Inf,-167L,-133L,-100L,-67L,-33L,0L,33L,67L,100L,133L,167L,Inf)
m    <- c(231, 234, 230, 227, 228, 221, 217, 238, 263, 277, 298, 316, 348)
s    <- c( 56,  58,  40,  40,  32,  28,  28,  28,  26,  30,  32,  34,  92)
n    <- c(400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400)

ttau <- lapply(tau, FUN=mathout, flags=list(cat=FALSE))
tm   <- lapply(as.integer(m), FUN=mathout, flags=list(cat=FALSE))
ts   <- lapply(as.integer(s), FUN=mathout, flags=list(cat=FALSE))
tn   <- lapply(as.integer(n), FUN=mathout, flags=list(cat=FALSE))

t <- cbind('\\(\\tau\\)'=ttau, '\\(m\\)'=tm, '\\(s\\)'=ts, '\\(n\\)'=tn)
knitr::kable(t, caption="Miller (1986) data", row.names=FALSE, escape=FALSE)
```

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

```{r}
z <- function(m, s, n, tau, c, mu_A, sigma_A, mu_B, sigma_B, mu_M)
  dfrac(m - denote(mu[tau], ET(tau, c, mu_A, sigma_A, mu_B, sigma_B, mu_M),
              "the expected mean response time"),
    s / sqrt(n))

math(call("=", quote(z[tau]), z))
```

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.

```{r}
zv <- Vectorize(z, vectorize.args = c('m', 's', 'n', 'tau'))
hook(zv(.M, .S, .N, .Tau, .C, .MA, .SA, .MB, .SB, .MM), z[.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)

math(call("=", quote(X["8 df"]^2L), gof))
```

with the degrees of freedom given by the difference between the number of
observations (`r length(tau)`) and the number of free model
parameters \(\theta = \langle\mu_{\mathrm A}, \sigma_{\mathrm A},
\mu_{\mathrm B}, \sigma_{\mathrm B}, \mu_{\mathrm M}\rangle\); the barrier \(c\)
is only a scaling parameter.

```{r, echo=FALSE}
lsq <- function(tau, m, s, n)
{
  invisible(theta <- c(mu_A=0.5, sigma_A=5, mu_B=0.5, sigma_B=5, mu_M=150))
  optim(theta, gof, tau=tau, m=m, s=s, n=n, hessian=TRUE)
}

math(call("=", quote(hat(theta)), lsq))
```

The best fitting parameter values and their confidence intervals are given in
Table \@ref(tab:params).

```{r params, echo=FALSE}
fit <- lsq(tau, m, s, n)
theta <- fit$par
se <- sqrt(diag(solve(fit$hessian)))
ci <- mapply(FUN=qnorm, theta, se, MoreArgs=list(p=c(0.025, 0.975)),
    SIMPLIFY=FALSE)

ntheta <- lapply(lapply(FUN=as.symbol, names(theta)), FUN=mathout,
    flags=list(cat=FALSE))
ttheta <- lapply(theta, FUN=mathout, flags=list(cat=FALSE))
tci <- lapply(ci, FUN=mathout, flags=list(cat=FALSE, sep=","))

t <- cbind(Parameter=ntheta, Estimate=ttheta, CI=tci)
knitr::kable(t, caption="Model fit", row.names=FALSE, escape=FALSE)
```

The goodness-of-fit statistic indicates some lack of
fit, \(X^2(8\ \mathrm{df}) = `r sprintf("%.2f", fit$value)`, 
  p = `r sprintf("%.4f", 1-pchisq(fit$value, df=8))`\). Given the large trial
numbers in the original study, this is not an unexpected result. For more
detail, especially on fitting the observed standard deviations, the reader is
referred to the original paper [@schwarz1994].

# Conclusion

This package allows R\ to render its terms in pretty mathematical equations. It
extends the current features of R\ and existing packages for displaying
mathematical formulas in R\ [@murrell2000], but most
importantly, \pkg{mathml} bridges the gap between computational needs,
presentation of results, and their reproducibility. The package supports both
MathML and LaTeX/MathJax for use in R\ Markdown documents, presentations and
Shiny App webpages.

Researchers or teachers can already use R\ Markdown to conduct analyses and show
results, and \pkg{mathml} smoothes this process and allows for integrated
calculations and output. As shown in the case study of the previous
section, \pkg{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 [@green1977],
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 [-@schwarz1994].

````{r}
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)))))
}

math(f1)
````

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

````{r}
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))))
}

math(f2)
````

As the reader may know from their 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.

One limitation of the package is the lack of a convenient way to insert line
breaks. This is mostly due to lacking support by MathML and LaTeX renderers. 
For example, in its current stage, the LaTeX package breqn [@breqn] is mostly
a proof of concept. Moreover, \pkg{mathml} only works in one direction,
that is, it is not possible to translate from LaTeX or HTML back to R\ 
[see @latex2r, for an example].

The package \pkg{mathml} is available for R\ version 4.2 and later, and can be
easily installed using the usual `install.packages("mathml")`. At its present
stage, it supports output in HTML, LaTeX, and Microsoft 
Word [via pandoc, @pandoc]. The source code of the package is found
at https://github.com/mgondan/mathml.

# Appendix: Customizing the package

## Implementation details

For convenience, the translation of the R\ expressions is achieved through a
Prolog interpreter provided by another R\ package \CRANpkg{rolog} [@rolog]. If a
version of SWI-Prolog [@swipl] is found on the system, \CRANpkg{rolog} connects
to it. Alternatively, the SWI-Prolog runtime libraries can be conveniently
accessed by installing the R\ package \CRANpkg{rswipl} [@rswipl]. Prolog is a
classical logic programming language with many applications in expert systems,
computer linguistics and symbolic artificial intelligence. The strength of
Prolog lies in its concise representation of facts and rules for knowledge and
grammar, as well as its efficient built-in search engine for closed world
domains. Whereas Prolog is weak in statistical computation, but strong in
symbolic manipulation, the converse may be said for the R\ language. \pkg{rolog}
bridges this gap by providing an interface to a SWI-Prolog
distribution [@swipl] in R. The communication between the two systems is mainly
in the form of queries from R\ to Prolog, but two Prolog functions allow ring
back and evaluation of terms in R. 

The proper term for a Prolog "function" is predicate, and it is typically
written with name and arity (i.e., number of arguments), separated by a forward
slash. Thus, at the Prolog end, the predicate `math/2` translates the
representation of the R\ call `pbinom(K, N, Pi)` into a more general 
representation of an R\ function `fn/2` with the name `P_Bi`, one 
argument `X =< K`, and the two parameters `N` and `Pi`, as shown below.

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

`math/2` operates like a macro that translates one mathematical
element (here, `pbinom(K, N, Pi)`) to another mathematical element,
namely `fn(Name, (Args ; Pars))`. The low-level predicate `ml/3` is used to
convert these basic elements to MathML.

````prolog
ml(fn(Name, (Args ; Pars)), M, Flags)
 => ml(Name, N, Flags),
    ml(paren(list(op(;), [list(op(','), Args), list(op(','), Pars)])), X, Flags),
    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 section on errors above).

Several ways exist for translating new R\ terms to their mathematical
representation. We have already seen above how to use "hooks" to translate long
variable names from R\ to compact mathematical signs, as well as functions such
as cumulative probabilities \(P(X \le k)\) to different representations
like \(\sum_{i=0}^k P(X = i)\). Obviously, the hooks require that there already
exists a rule to translate the target representation into MathML and MathJax.

In this appendix we describe a few more ways to extend the set of translations
according to a user's needs. As stated in the background section, the Prolog end
provides two classes of rules for translation, macros `math/2,3,4` mirroring the
R\ hooks mentioned above, and the low-level predicates `ml/3` and `jax/3` that
create proper MathML and LaTeX terms.

## Linear models

To render the model equation of a linear model such
as `lm(EOT ~ T0 + Therapy, data=d)` in mathematical
form [see also @equatiomatic], it is sufficient to map the `Formula`
in `lm(Formula, Data)` to its respective equation. This can be done in two ways,
using either the hooks described above, or a new `math/2` macro at the Prolog
end.

````r
hook(lm(.Formula, .Data), .Formula)
````

The hook is simple, but is a bit limited because only R's tilde-form
of linear models is shown, and it only works for a call with exactly two
arguments.

Below is an example of how to build a linear equation of the
form \(Y = b_0 + b_1X_1 + ...\) using the Prolog macros from \pkg{mathml}.

````prolog
math_hook(LM, M) :-
    compound(LM),
    LM =.. [lm, ~(Y, Sum) | _Tail],
    summands(Sum, Predictors),
    findall(subscript(b, X) * X, member(X, Predictors), Terms),
    summands(Model, Terms),
    M = (Y == subscript(b, 0) + Model + epsilon).
````

The predicate `summands/2` unpacks an expression `A + B + C` to a
list `[C, B, A]` and vice-versa (see the file `lm.pl` for details).

````{r}
rolog::consult(system.file(file.path("pl", "lm.pl"), package="mathml"))

term <- quote(lm(EOT ~ T0 + Therapy, data=d, na.action=na.fail))
math(term)
````

See Section 3.11 for replacing lengthy R\ labels by short mathematical symbols.

## \(n\)-th root

Base R\ does not provide a function like `cuberoot(x)` or `nthroot(x, n)`, and
the present package does not support the respective representation. To obtain a
cube root, a programmer would typically type `x^(1/3)` or better `x^{1/3}` (see
the practice section why the curly brace is preferred in an exponent), resulting
in \(x^{1/3}\) which may still not match everyone's taste. Here we describe the
steps needed to represent the \(n\)-th root as \(\sqrt[n]x\).

We assume that `nthroot(x, n)` is available in the current
namespace [manually defined, or from R\ package \CRANpkg{pracma}, @pracma], so that
the names of the arguments and their order are accessible to `canonical()` if
needed. As we can see below, \pkg{mathml} uses a default
representation `name(arguments)` for such unknown functions.

```{r}
nthroot <- function(x, n)
  x^{1L/n}

term <- canonical(quote(nthroot(n=3L, 2L)))
math(term)
```

A proper MathML term is obtained by `mlx/3` (the x in mlx indicates that it is
an extension and is prioritized over the default ml/3 rules). `mlx/3`
recursively invokes `ml/3` for translating the function arguments _X_ and _N_,
and then constructs the correct MathML entity `<mroot>...</mroot>`.

```prolog
mlx(nthroot(X, N), M, Flags) :-
    ml(X, X1, Flags),
    ml(N, N1, Flags),
    M = mroot([X1, N1]).
```

The explicit unification `M = ...` in the last line serves to avoid clutter in
the head of `mlx/3`. The Prolog file `nthroot.pl` also includes the respective
rule for LaTeX and can be consulted from the package folder via the underlying
package \pkg{rolog}.

```{r}
rolog::consult(system.file(file.path("pl", "nthroot.pl"), package="mathml"))

term <- quote(nthroot(a * (b + c), 3L)^2L)
math(term)
term <- quote(a^(1L/3L) + a^{1L/3L} + a^(1.0/3L))
math(term)
```

The file `nthroot.pl` includes three more statements `precx/3` and `parenx/3`,
as well as a `math_hook/2` macro. The first sets the operator precedence of the
cubic root above the power, thereby putting a parentheses around nthroot
in \((\sqrt[3]{\ldots})^2\). The second tells the system to increase the counter
of the parentheses below the root, such that the outer parenthesis becomes a
square bracket.

The last rule maps powers like `a^(1L/3L)` to `nthroot/3`, as shown in the
first summand. Of course, \pkg{mathml} is not a proper computer algebra system.
As is illustrated by the other terms in the sum, such macros are limited to
purely syntactical matching, and terms like `a^{1L/3L}` with the curly brace
or `a^(1.0/3L)` with a floating point number in the numerator are not detected.

# Acknowledgment

Supported by the Erasmus+ program of the European
Commission (2019-1-EE01-KA203-051708).
