R/qtlcharts is an R package to create interactive charts for quantitative trait locus (QTL) mapping data, for use with R/qtl.

The charts are saved to a temporary file and opened in a web browser, though you may also save them to a specific file, to be viewed later. We are focusing development on the Chrome browser. The graphs may also be viewed in other browsers, such as Safari, Opera and Firefox, but it can be hard to accommodate all possible browser differences.

Information on installing R/qtlcharts is available here; it requires installation of the packages R/qtl and htmlwidgets.

Basic use

You first need to load the package.

library(qtlcharts)

Let’s begin by considering the function iplotCorr, which creates a heatmap of a correlation matrix, linked to scatterplots of the underlying variables.

We’ll first load the geneExpr dataset, included with the R/qtlcharts package.

data(geneExpr)

This is a list with two components. The first component, geneExpr$expr, is a 491 × 100 matrix of gene expression data; the second component, geneExpr$genotype, is a vector of genotypes (of length 491) at a QTL that influences those 100 genes’ expression values. (The genes were selected from a larger expression genetics study, on that basis: that they are all influenced by this QTL.)

Let’s pull out those two components of geneExpr as separate objects, expr and geno.

expr <- geneExpr$expr
geno <- geneExpr$genotype

The simplest use of iplotCorr is with a numeric matrix, as with the expr dataset. For example:

iplotCorr(expr, reorder=TRUE)

This will open an interactive figure in a web browser, with a heat map of the correlation matrix of the genes on the left linked to the underlying scatterplots. With the argument reorder=TRUE, the genes are reordered (by hierarchical clustering with the R function hclust) to bring genes with similar expression patterns next to each other.

The following is a snapshot. (See a live example at the R/qtlcharts website.)

iplotCorr example

iplotCorr example

The left panel is a heat map of the correlation matrix. If you hover over a pixel in the panel, the individual correlations are shown. If you click on a pixel, the corresponding scatterplot is shown on the right.

We can have the scatterplot colored by the QTL genotypes contained in the vector geno, through the group argument.

iplotCorr(expr, geno, reorder=TRUE)

Here’s a snapshot of the result. (Again, see a live example at the R/qtlcharts website.)

iplotCorr example with scatterplot colored by genotype

iplotCorr example with scatterplot colored by genotype

Saving to a file

The interactive graphs produced by R/qtlcharts are, by default, saved to a temporary file and then opened in the default web browser. If you want to save a chart to a particular file, assign the result to some object and use the function saveWidget() in the htmlwidgets package, as follows:

corrplot <- iplotCorr(expr, reorder=TRUE)
htmlwidgets::saveWidget(corrplot, file="~/Desktop/iplotCorr_example.html")

The chart will not be opened for viewing. If you wish to view the chart, type the name of the object that you created, as so:

corrplot

Charts saved with saveWidget() will, by default, contain all necessary code/resources. A smaller file may be obtained by using selfcontained=FALSE in the call to saveWidget(), in which case the external resources will be placed in an adjacent directory.

Chart customization options

There are a number of ways in which one may customize the interactive graphs produced by R/qtlcharts, such as different colors or axis labels. The main method for such customization is through a function argument, chartOpts, which takes a named list.

For example, to change the dimensions of the image produced by iplotCorr, use the height option, which specifies the height of each panel in pixels; by default, each panel is square, with width = height:

iplotCorr(expr, reorder=TRUE, chartOpts=list(height=600))

If you want a different choice of color for the points in the scatterplot, use scatcolors.

iplotCorr(expr, reorder=TRUE,
          chartOpts=list(scatcolors="lightblue"))

You can modify both height and scatcolors at once.

iplotCorr(expr, reorder=TRUE,
          chartOpts=list(height=600, scatcolors="lightblue"))

If the points are to be colored by genotype, you need to specify the same number of colors as there are genotype groups.

iplotCorr(expr, geno, reorder=TRUE,
          chartOpts=list(scatcolors=c("lightblue", "lightgreen", "pink")))

If you want a caption below the chart, use caption. This is intended largely for charts saved to a file with htmlwidgets::saveWidget().

iplotCorr(expr, geno, reorder=TRUE,
          chartOpts=list(scatcolors=c("lightblue", "lightgreen", "pink"),
                         caption=paste("Click on the heatmap on the left to see",
                                 "corresponding scatterplots on the right.") ))

A full list of options is available in the chartOpts vignette.

Further control of chart dimensions

To provide general control of the sizes of charts, when they are to be displayed in a web browser, use the function setScreenSize() to define the maximum height and width, in pixels. The default is 700 (height) × 1000 (width). For a smaller screen (600×900), use

setScreenSize("small")

For a larger screen (1200×1600), use

setScreenSize("large")

The default is “normal” (700×1000):

setScreenSize("normal")

You can also define a custom height and width:

setScreenSize(height=500, width=700)

Each charting function (e.g., iplotScanone and iplotPXG) has a defined default aspect ratio, and the height and width are chosen to fit within the maximum sizes defined by the call to setScreenSize.

You can also use chartOpts to define a specific height and width for a given plot.

iplotCorr(expr, chartOpts=list(height=400, width=800))

Including a chart in an R Markdown document

The interactive charts may be included within an R Markdown-based report. See the separate Rmarkdown vignette.

Examples of the charts’ use

iplotCorr

iplotCorr creates an heat map of a correlation matrix, linked to the underlying scatterplots. We’ve shown a few examples above, but there are several other features to consider.

By default, iplotCorr will calculate the correlation matrix to be displayed, using Pearson’s correlation; this can be sensitive to outliers. If you wish to use a different statistic, provide the matrix with the argument corr. For example, to use Spearman’s rank correlation:

data(geneExpr)
expr <- geneExpr$expr
geno <- geneExpr$genotype
spearman <- cor(expr, use="pairwise.complete.obs", method="spearman")

iplotCorr(expr, geno, corr=spearman)

If you provide your own correlation measure, the reorder argument is ignored. If you want to reorder the genes, to bring correlated genes together, you need to do it yourself. You can use hclust for this.

ord <- hclust(as.dist(-spearman))$order
iplotCorr(expr[,ord], geno, corr=spearman[ord,ord])

You can view just a portion of the correlation matrix by specifying rows and cols. For example, we could fit a regression line for each gene, relating the expression pattern to the eQTL genotype, in order to identify the genes with positive and negative associations with genotype.

beta <- apply(expr, 2, function(y,x) lm(y~x)$coef[2], geno)
pos <- which(beta > 0)
neg <- which(beta < 0)

We could then look at just the portion of the correlation matrix where the rows are genes with positive effects and the columns are genes with negative effects.

iplotCorr(expr, geno, rows=pos, cols=neg)

If you wish to use a different correlation statistic, you need to do the subsetting in advance, as the rows and cols arguments are ignored if corr is provided.

iplotCorr(expr, geno, corr=spearman[pos,neg])

The connection between expr and spearman[pos,neg] is taken from the column names of the former and both the row and column names of the latter.

iplot

iplot creates a simple, interactive scatterplot. Hover over a point to view an identifier. Points may be colored according to some provided group categories.

Let’s simulate a bit of data, as an example.

n <- 100
x <- rnorm(n)
grp <- sample(1:3, n, replace=TRUE)
y <- grp*x + rnorm(n)

We use iplot as follows.

iplot(x, y, grp)

The default is for the individual IDs to be simple sequential numbers. We can include more detail by providing character strings.

id <- paste(seq(along=x), " x =", round(x,1), "  y =", round(y, 1), " group =", grp)
iplot(x, y, grp, id)

If you want to change the x- and y-axis labels, you need to use the chartOpts argument. You can similarly change the point colors.

iplot(x, y, grp, id,
      chartOpts=list(xlab="x variable", ylab="y variable",
                     pointcolor=c("slateblue", "#999999", "violetred")))

iplotMap

iplotMap creates an interactive genetic map. Hover over a marker position to view the marker name. There’s also a search box, for searching for a particular marker.

Let’s first load the hyper data included with R/qtl and pull out the genetic map with pull.map. We need to load R/qtl, too.

library(qtl)
data(hyper)
map <- pull.map(hyper)

We use iplotMap as follows.

iplotMap(map)

If you want to shift each chromosome so that the initial marker is at 0, use shift=TRUE.

iplotMap(map, shift=TRUE)

iplotMap will also take a cross object, in which case it uses pull.map to extract and then plot the map.

iplotMap(hyper)

If you want to change the x- and y-axis labels, you need to use the chartOpts argument:

iplotMap(map, chartOpts=list(xlab="Linkage group", ylab="Position (Mbp)"))

If you want to plot a subset of the chromosomes, use the chr argument:

iplotMap(map, chr=c(5, 10, 15, "X"))

You can use the chart option title to put a title centered above the map panel.

iplotMap(map, chr=c(5, 10, 15, "X"),
         chartOpts=list(title="Selected chromosomes"))

iplotScanone

iplotScanone creates an interactive chart with LOD curves from a genome scan linked to estimated QTL effects.

Let’s first load the hyper data, included with R/qtl

data(hyper)
hyper <- calc.genoprob(hyper, step=1)
out <- scanone(hyper)

If you provide just the output from scanone, the only interactivity is that hovering over marker positions on the LOD curves will give information about the marker name, position, and LOD score.

iplotScanone(out)

You can use the chr to plot only selected chromosomes.

iplotScanone(out, chr=c(1, 4, 6, 15))

If you also provide the cross object, clicking on the markers will generate a QTL effect plot on the right. By default, means with ± 1 SE error bars will be shown.

iplotScanone(out, hyper)

By default, in calculating the means and standard errors, a single imputation is used to handle missing genotype information, using the R/qtl function fill.geno. Arguments may be passed to fill.geno as a list, with the argument fillgenoArgs. For example, to do the imputation by Viterbi, assuming a 1% genotyping error rate, do the following:

iplotScanone(out, hyper, fillgenoArgs=list(method="argmax", error.prob=0.01))

If you want to plot raw phenotypes vs genotypes, rather than mean ± SE, use pxgtype="raw", as follows:

iplotScanone(out, hyper, pxgtype="raw")

In the phenotype × genotype plot, pink points indicate that the genotype was imputed.

To change the color of the LOD curves, use the chartOpts option lod_linecolor. Colors are named as in html, or with hex RGB codes like "#842DCE" or #82C. (See examples here and a list of named colors here.)

iplotScanone(out, hyper, chartOpts=list(lod_linecolor="black"))

You can change the background colors with the chartOpts options darkrect and lightrect.

iplotScanone(out, hyper, chartOpts=list(darkrect="#CCC", lightrect="#EEE"))

Change the colors of the lines in the mean ± SE plot with the chartOpts option eff_linecolor.

iplotScanone(out, hyper, chartOpts=list(lod_linecolor="DarkViolet",
                                        eff_linecolor="DarkViolet"))

iplotMScanone

iplotMScanone creates a heatmap of LOD curves for a set of genome scans (for example, for each of a time course of phenotypes), linked to a plot of the individual genome scans, and also to the QTL effects.

We’ll consider the example dataset grav, included with R/qtlcharts. These are data on seedling gravitropism in Arabidopsis recombinant inbred lines; see Moore et al. (2013). The data are also available at the Mouse Phenome Database.

We first load the data.

data(grav)

We calculate QTL genotype probabilities across the genome, but reduce these probabilities to a 1 cM grid.

grav <- calc.genoprob(grav, step=1)
grav <- reduce2grid(grav)

We’ll perform a genome scan by Haley-Knott regression at each of the 241 time points.

phecol <- 1:nphe(grav)
out <- scanone(grav, phe=phecol, method="hk")

If we pass the scanone output to iplotMScanone, we get a basic plot.

iplotMScanone(out)

The upper-left panel is a heat map of the genome scans for each time point. Hover over this panel, and the LOD curves for horizontal slices (at specific time points) are shown below, and the LOD scores for vertical slices (at a given genomic position, across time points) are shown on the right.

The phenotypes here are measurements over time. The default result uses a qualitative y-axis for the LOD heatmap in the upper-right. We could instead have a quantitative y-axis scale by including the times at which the phenotypes were measured. Currently, the time points must be equally spaced.

The times at which the phenotypes were measured are included in the grav dataset as an attribute, "time".

times <- attr(grav, "time")

If we call iplotMScanone again with the times argument, we’ll get the quantitative y-axis scale. We can also include a y-axis label using the chart option lod_ylab.

iplotMScanone(out, times=times, chartOpts=list(lod_ylab="Time (hrs)"))

We can use the function estQTLeffects to get estimated QTL effects at each putative QTL position and for each time point.

eff <- estQTLeffects(grav, phe=phecol, what="effects")

If we pass both the genome scan results and the estimated QTL effects, the heat map will display signed LOD scores (blue for negative effects and red for positive effects) and QTL effects will be shown in the upper-right panel.

iplotMScanone(out, effects=eff, times=times)

As an alternative to first calling estQTLeffects, you can pass the genome-scan results and the cross object.

iplotMScanone(out, grav, times=times)

iplotScantwo

iplotScantwo creates an interactive plot of the output of the R/qtl function scantwo, of a two-dimensional, two-QTL genome scan.

Let’s first run scantwo using a coarse step size and focusing on selected chromosomes. For the sake of speed, we’ll use Haley-Knott regression, even though it does not work well for selectively genotyped data.

data(hyper)
hyper <- hyper[c(1,4,6,15),]
hyper <- calc.genoprob(hyper, step=5)
out2 <- scantwo(hyper, method="hk", verbose=FALSE)

The simplest use of iplotScantwo is to just provide the scantwo output, but it’s best to also provide the QTL cross data, so that QTL effect plots may be produced.

iplotScantwo(out2, hyper)

There are drop-down menus to select which LOD scores to display in the heat map. Hover over the heatmap to view the LOD scores, and click to view cross-sectional slices below and the QTL effect plot to the right.

You can control the size of the heatmap through the chart option pixelPerCell — the number of pixels, in each direction, for each small rectangle within the heatmap.

iplotScantwo(out2, hyper, chartOpts=list(pixelPerCell=12))

You can orient the plot with the diagonal going from the top-left to the lower right (so that chromosome 1 is at the top rather than at the bottom) using the chart option oneAtTop, as follows.

iplotScantwo(out2, hyper, chartOpts=list(oneAtTop=TRUE))

The estimated QTL effects in the right-hand panels are derived following a single imputation to fill in missing data, and so are a bit crude.

iplotCurves

iplotCurves creates a plot of a set of curves linked to one or two scatterplots.

As an example, we’ll again consider the phenotypes in the dataset grav, included with R/qtlcharts.

data(grav)
times <- attr(grav, "time")
phe <- grav$pheno

The object phe is an 162×241 matrix: curves for 162 individuals at 241 time points, with the time points specified by times.

The simplest use of iplotCurves is the following.

iplotCurves(phe)

This plots the rows of phe as a set of curves. If you hover over a curve, it will be highlighted, and the individual ID (taken from the rownames) will be displayed on the right. If the input matrix has no row names, numeric indices are used.

The values on the x-axis are simply column indices. To use the actual values in times, include it as a second argument.

iplotCurves(phe, times)

To change the axis labels, pass curves_xlab and curves_ylab using chartOpts.

iplotCurves(phe, times, chartOpts=list(curves_xlab="Time (hrs)",
                                       curves_ylab="Response"))

The curves_ bit is to distinguish this panel from the optional scatterplot panels; more on this below. You can use just xlab and ylab, as follows.

iplotCurves(phe, times, chartOpts=list(xlab="Time (hrs)",
                                       ylab="Response"))

More interesting is to have the curves linked to scatterplots, of the response at particular times. For example, 2 hrs vs 4 hrs and 4 hrs vs 6 hrs. This is done by providing two-column matrices whose rows correspond to the main matrix.

iplotCurves(phe, times, phe[,times==2 | times==4], phe[,times==4 | times==6],
            chartOpts=list(curves_xlab="Time (hours)", curves_ylab="Response",
                           scat1_xlab="Angle at 2 hrs", scat1_ylab="Angle at 4 hrs",
                           scat2_xlab="Angle at 4 hrs", scat2_ylab="Angle at 6 hrs"))

You can also have just one linked scatterplot.

iplotCurves(phe, times, phe[,times==2 | times==4],
            chartOpts=list(curves_xlab="Time (hours)", curves_ylab="Response",
                           scat1_xlab="Angle at 2 hrs", scat1_ylab="Angle at 4 hrs"))

You can also include a numeric vector of groups assignments for the curves, to have groups of curves and points with common colors. For example, we can pull out the genotype for a marker and color the curves by their genotype. We’ll use the fill.geno function in R/qtl to impute any missing genotypes.

g <- pull.geno(fill.geno(grav))[,"BF.206L-Col"]
iplotCurves(phe, times, phe[,times==2 | times==4], phe[,times==4 | times==6],
            group=g,
            chartOpts=list(curves_xlab="Time (hours)", curves_ylab="Response",
                           scat1_xlab="Angle at 2 hrs", scat1_ylab="Angle at 4 hrs",
                           scat2_xlab="Angle at 4 hrs", scat2_ylab="Angle at 6 hrs"))

To change the colors that are used in the plot, you need to provide a vector of colors with the chartOpts options strokecolor, strokecolorhilit, pointcolor, and pointcolorhilit. If you want the points and curves to be the same, you can just use color and colorhilit.

iplotCurves(phe, times, phe[,times==2 | times==4], phe[,times==4 | times==6],
            group=g,
            chartOpts=list(curves_xlab="Time (hours)", curves_ylab="Response",
                           scat1_xlab="Angle at 2 hrs", scat1_ylab="Angle at 4 hrs",
                           scat2_xlab="Angle at 4 hrs", scat2_ylab="Angle at 6 hrs",
                           color=c("pink", "lightblue"), colorhilit=c("red", "blue")))

iboxplot

iboxplot creates an interactive graph for a large set of box plots (rendered as lines connecting the quantiles), linked to underlying histograms.

Let’s first simulate some data. These data are not particularly interesting: 10,000 measurements from 500 normal distributions with a common SD but shifts in the mean.

# simulated data
n.ind <- 500
n.gene <- 10000
mat <- matrix(rnorm(n.ind * n.gene, rnorm(n.ind, 5)), ncol=n.gene)
dimnames(mat) <- list(paste0("ind", 1:n.ind),
                      paste0("gene", 1:n.gene))

iboxplot expects a matrix, with the rows corresponding to individuals and the columns corresponding to measurements; the plotted boxplots will characterize the distributions in each row.

iboxplot(mat)

In the upper panel, each curve represents a particular quantile and connects the values for the different individuals, who are arranged along the x-axis. (By default, the individuals are sorted by their median; use orderByMedian=FALSE to prevent this reordering.) Hover over a different individuals and the corresponding histograms are shown below. Click to have a particular histogram retained; click again (if you can do so precisely!) to remove.

If you want to order the individuals in a different way, you need to rearrange the rows of the matrix in advance and then use orderByMedian=FALSE. For example, to order the individuals by their 99.9th percentiles, you could do:

iboxplot(mat[order(apply(mat, 1, quantile, 0.99)),], orderByMedian=FALSE)

By default, the function plots the 0.1, 1, 10, 25, and 50, 75, 90, 99, and 99.9th percentiles. To plot a different set of percentiles, use the argument qu, giving a vector of numeric values between 0 and 0.5. The median will always be included, and for every quantile \(q\), the corresponding \(1-q\) quantile will be included as well.

To leave off the 0.1 and 99.9 percentiles, do the following.

iboxplot(mat, qu=c(0.25, 0.05, 0.01))

To use different colors for the quantile curves, include the option qucolors through the chartOpts argument. The vector should have length length(qu)+1 (so, by default, 5). For each quantile \(q\), the same color is used for the corresponding quantile \(1-q\). Here’s an example. The first color will correspond to the median, and then they move from the inside out.

iboxplot(mat, qu=c(0.25, 0.05, 0.01),
         chartOpts=list(qucolors=c("black", "green", "red", "blue")))

You can change the bins in the histograms in the lower panel through the argument breaks. If a single numeric value, it indicates the number of bins; if a vector, it indicates the endpoints for the bins. The vector of breakpoints do not need to span the range of the data, but a warning will be issued if they don’t.

Here’s an example specifying the number of bins.

iboxplot(mat, breaks=151)

Here’s an example specifying actual bins. (This will issue a warning.)

iboxplot(mat, breaks=seq(-2, 12.1, by=0.1))

To change the axis labels, use the chartOpts options xlab and ylab, for the x- and y-axes, respectively, in the top panel. The y-axis label in the top panel will be used as the x-axis label in the lower panel.

iboxplot(mat, chartOpts=list(xlab="Mice", ylab="Expression level"))

iheatmap

iheatmap creates an interactive heatmap of a numeric matrix, with pixels linked to plots of horizontal and vertical slices.

We’ll first generate a set of data to plot: \(z(x,y) = x y \exp(-x^2 - y^2)\), for \(x,y \in [-2, 2]\)

n <- 101
x <- y <- seq(-2, 2, len=n)
z <- matrix(ncol=n, nrow=n)
for(i in seq(along=x))
  for(j in seq(along=y))
    z[i,j] <- x[i]*y[j]*exp(-x[i]^2 - y[j]^2)

We can now call iheatmap with these data. The first argument is a matrix with the \(z\)-values. The next two optional arguments are vectors that indicate the \(x\) and \(y\) values.

iheatmap(z, x, y)

By default, the iheatmap takes as the range for the z-axis, \([-M, M]\) where \(M = \max |z|\). If you want a different range, you can use the option zlim with chartOpts; you need to provide three values: the minimum, the center, and the maximum. The color of the heatmap can be controlled with the option colors; again this is a vector of three values: the color for the minimum, the color for the middle value, and the color for the maximum value.

Here’s an example with expanded z-axis limits and different colors.

iheatmap(z, x, y,
         chartOpts=list(zlim=c(-1, 0, 1),
                        colors=c("purple", "white", "orangered")))

The labels for the axes can be controlled with the options xlab, ylab, and zlab. Here’s an example.

iheatmap(z, x, y, chartOpts=list(xlab="theta", ylab="psi", zlab="response"))

iplotRF

iplotRF is an interactive version of the R/qtl function plot.rf, a diagnostic plot of pairwise linkage between markers.

We’ll look at the badorder dataset included with R/qtl. We first run est.rf to estimate the recombination fraction between each pair of markers and to further calculate a LOD score indicating evidence for linkage for each marker pair.

data(badorder)
badorder <- est.rf(badorder)

Making the plot is simple.

iplotRF(badorder)

The main plot is a heatmap of the LOD scores indicating pairwise linkage; blue indicates large LOD score with recombination fraction < 1/2; red indicates large LOD score with recombination fraction > 1/2. (Think: red is bad; blue is good.) Hover over the heat map to view the LOD and recombination fraction values. Click to generate additional views: to the right, a cross-tabulation of two-locus genotypes for the pair of markers; below, plots of the LOD scores for linkage, for the selected markers versus every other marker across the genome.

Hover over the row and column headings in the cross-tab to view conditional distributions (genotypes as percentages of the non-missing observations; missing value as percentage of the total individuals). Click on the points in the lower panels to refresh the cross-tab and lower panels with the selected marker.

By default, the LOD scores in the heatmap are truncated at 12. This can be changed with the chart option lodlim, which takes an interval. Marker pairs with LOD below the lower limit will be omitted from the heat map, which can give greater performance in the case of a very large number of markers.

iplotRF(badorder, chartOpts=list(lodlim=c(2, 15)))

ipleiotropy

ipleiotropy is for exploring potential pleiotropy for a pair of traits: a double-slider can be used to select two QTL positions on a chromosome which are then used to color the points in a scatterplot of two phenotypes.

Let’s look at the gravitropism phenotype for a pair of times in the grav data (Moore et al. 2013).

data(grav)
grav <- calc.genoprob(grav, step=1)
out <- scanone(grav, pheno.col=c(156, 241), method="hk")

We then call ipleiotropy with the cross, the scanone output, and indications of the phenotype columns, LOD score columns, and chromosome (a single one chromosome must be selected).

ipleiotropy(grav, out, pheno.col=c(156, 241), lodcolumn=1:2, chr=1)

scat2scat

scat2scat produces a pair of linked scatterplots, where the points in the first scatterplot drive the contents in the second scatterplot.

For example, consider a large set of scatterplots, with a pair of statistics summarizing each. The first plot in scat2scat will be a scatterplot of those summary statistics, and if you click on a point, the corresponding scatterplot will be shown on the right.

The input data is a \(n \times 2\) matrix that will be shown in first scatterplot, and a list (of length \(n\)) of 2-column matrices that will be shown shown in the second scatterplot when clicking on the corresponding point in the first scatterplot. (The matrices can have three columns, in which case the third column is used as a grouping factor for coloring the points in the scatterplot.)

To illustrate the plot, we’ll first simulate some data: random pairs of SDs and correlations, and multivariate normal datasets with those SDs and correlations.

p <- 500
n <- 300
SD <- runif(p, 1, 5)
r <- runif(p, -1, 1)
scat2 <- vector("list", p)
V <- SD[i]^2*matrix(c(1, r[i], r[i], 1), ncol=2)
for(i in 1:p)
   scat2[[i]] <- matrix(rnorm(2*n), ncol=2) %*% chol(V)
scat1 <- cbind(SD=SD, r=r)

We then create the plot as follows:

scat2scat(scat1, scat2)

itriplot

itriplot produces a scatterplot on an equilateral triangle, with each point representing a trinomial probability, (\(p_0\), \(p_1\), \(p_2\)).

To illustrate the plot, we’ll first simulate some data.

n <- 26
p <- matrix(runif(3*n), ncol=3)
p <- p / colSums(p)
ind <- LETTERS
g <- sample(1:3, n, replace=TRUE)

We then create the plot as follows. The indID argument is a vector of individual identifiers for tool tip labels, and the group argument is a vector of categories for coloring the points.

itriplot(p, indID=ind, group=g)