--- a
+++ b/docs/getting-started-model-analysis.md
@@ -0,0 +1,647 @@
+# Getting Started with Model Analysis
+
+The following functions are provided to assist in the analysis and visualisation
+of a model fitted by `mixOmics`.  Focus has been given to sPLS-DA models where a
+single block is being fitted (single-omics) and to DIABLO models where multiple
+blocks are being fitted simulataneously (multi-omics).  The code below will
+generate known trained models and corresponding performance tests for each model
+type, which are then used to demonstrate the functions in the rest of this
+document.
+
+```R
+library(OmicsFold)
+library(mixOmics)
+library(dplyr)
+
+#######################################
+## sPLS-DA model and performance result
+#######################################
+
+# Prepare data
+data(srbct)
+splsda.X <- srbct$gene
+splsda.Y <- srbct$class
+
+# Model fit
+splsda.ncomp <- 3
+splsda.keepX <- c(9, 280, 30)
+splsda.model <- splsda(splsda.X, splsda.Y,
+                       ncomp = splsda.ncomp, keepX = splsda.keepX)
+
+# Preformance result -- takes about 1 min to run
+set.seed(1234) # for reproducibility
+splsda.perf <- perf(splsda.model, validation = "Mfold", folds = 5,
+                    dist = 'max.dist', nrepeat = 10,
+                    progressBar = FALSE)
+
+
+######################################
+## DIABLO model and performance result
+######################################
+
+# Prepare data
+data('breast.TCGA')
+diablo.data <- list(mRNA = breast.TCGA$data.train$mrna, 
+                    miRNA = breast.TCGA$data.train$mirna, 
+                    proteomics = breast.TCGA$data.train$protein)
+diablo.Y <- breast.TCGA$data.train$subtype
+diablo.design <- matrix(0.1, ncol = length(diablo.data),
+                        nrow = length(diablo.data), 
+                        dimnames = list(names(diablo.data), names(diablo.data)))
+diag(diablo.design) <- 0
+
+# Model fit
+diablo.ncomp <- 2
+diablo.keepX <- list(mRNA = c(6, 14), miRNA = c(5, 18), proteomics = c(6, 7))
+diablo.model <- block.splsda(X = diablo.data, Y = diablo.Y,
+                             ncomp = diablo.ncomp, keepX = diablo.keepX,
+                             design = diablo.design)
+
+# Preformance result -- takes about 1 min to run
+set.seed(4321) # for reproducibility
+diablo.perf <- perf(diablo.model, validation = 'Mfold', M = 10, nrepeat = 10, 
+                    dist = 'centroids.dist')
+```
+
+## Model variance analysis
+
+OmicsFold provides the [`get.block.centroids()`](#get-block-centroids) function
+to extract the centroids of variance across blocks in a DIABLO model and export
+them as a data frame and plot.  The
+[`get.model.variance()`](#get-model-variance) function is also provided which
+will take either an sPLS-DA model or a DIABLO model as input and creates a data
+frame showing the percentage contributions of each component to the model
+variance.
+
+### Get block centroids
+
+Applying the `get.block.centroids()` function to the DIABLO model from above
+gives the following:
+
+```R
+get.block.centroids(diablo.model)
+## A tibble: 150 x 4
+## Groups:   sample [150]
+##    sample label comp1   comp2
+##     <int> <fct> <dbl>   <dbl>
+##  1      1 Basal  2.38  2.22  
+##  2      2 Basal  1.74  2.03  
+##  3      3 Basal  2.28 -0.0206
+##  4      4 Basal  1.54  1.46  
+##  5      5 Basal  1.93 -1.17  
+## ... with 145 more rows
+```
+
+In addition to the tibble provided, a plot is produced showing the centroid
+across the blocks in the first and second component for each sample in the data.
+Each point in the plot represents a sample and is the mean of the three blocks
+loaded for the DIABLO model.  The circles show the distribution of all the
+points sharing the same label and gives a good indication of the separation of
+classifications in each of the first two components.
+
+![Block Centroids Plot](imgs/block-centroids-plot.png)
+
+### Get model variance
+
+The `get.model.variance()` function can be applied to the model to identify how
+much variation each component captures in the model for a given block:
+
+```R
+get.model.variance(diablo.model, "mRNA")
+##        Proportion Cumulative
+## comp1 0.794465823  0.7944658
+## comp2 0.008624793  0.8030906
+
+get.model.variance(diablo.model, "miRNA")
+##        Proportion Cumulative
+## comp1 0.669987181  0.6699872
+## comp2 0.001283394  0.6712706
+
+get.model.variance(diablo.model, "proteomics")
+##        Proportion Cumulative
+## comp1 0.770439163  0.7704392
+## comp2 0.002689493  0.7731287
+```
+
+Here we can see that the first component captures the vast majority of the
+variance for all three blocks but that slightly less of the variance in `miRNA`
+was captured than in `mRNA` or `proteomics`.
+
+## Feature analysis for sPLS-DA models
+
+OmicsFold provides functions for analysing the features selecting for fitting
+sPLS-DA models.
+
+- [`get.loadings.table()`](#get-feature-loadings-as-a-table) outputs a data
+  frame of loadings for features selected in each component of the model.
+- [`feature.selection.stability()`](#get-spls-da-feature-selection-stability)
+  gives the loadings from above for a single component, alongside stability
+  scores for the selected features.
+- [`merge.feature.stability()`](#merge-feature-stability-on-a-table-of-loadings)
+  takes a table such as that generated by `get.loadings.table()` and the output
+  from `feature.selection.stability()` merging them into a single data frame.
+- [`blockrank.splsda()`](#measure-variable-importance) takes an sPLS-DA model 
+  object and calculates a score measuring each features contribution to the variance in the response
+
+### Get feature loadings as a table
+
+By calling `get.loadings.table()` with the sPLS-DA model as a parameter, a data
+frame will be produced listing the loadings of each feature selected for the
+sparse feature set.  A separate column is included for each component and the
+features are arranged according to their loading value in the earliest component
+they are selected for.  Ordering is highest loading first.  Where a feature is
+selected in more than one component, it appears only once in the table, showing
+both loadings values.  Where features are not selected for a component, the
+loading value is `NA`.
+
+The code example below demonstrates typical output from this function.  The
+output has been truncated for brevity in this getting started guide.  The actual
+output is 320 lines long.
+
+```R
+get.loadings.table(splsda.model)
+##            comp1         comp2       comp3
+## g123  0.64922048            NA          NA
+## g846  0.44828969            NA          NA
+## g1606 0.30641602            NA          NA
+## ...
+## g1389         NA -2.312917e-01          NA
+## g246          NA -2.001212e-01          NA
+## g545          NA -1.986890e-01          NA
+## ...
+## g340          NA -4.226175e-02          NA
+## g742          NA  4.204118e-02  0.14003503
+## g1143         NA -4.124087e-02          NA
+## ...
+## g174          NA            NA -0.35418530
+## g1896         NA            NA -0.29939099
+## g603          NA            NA -0.27261886
+```
+
+### Get sPLS-DA feature selection stability
+
+Calling `feature.selection.stability()` with the model, performance test result
+and a component number as parameters generates a data frame showing the loadings
+values are were returned above, but for a single component only.  The second
+column in the data frame gives the stability value of the feature.  This is
+expressed as a fraction and indicates how many of the models in the performance
+test included this feature in the sparse set.  High stability indicates that a
+feature is very important to the fit of the model and is therefore probably
+important to separating the classes the model is trying to discriminate.
+
+The example below shows extraction of the stability values for the 9 top
+features included in component 1.  Note how all stabilities are an exact
+multiple of 0.02.  This is because our performance test used 10 repeats and 5
+folds, so 50 models were fitted, and each time a model selects a feature, this
+adds 0.02 to the stability of that feature.
+
+```R
+feature.selection.stability(splsda.model, splsda.perf, 1)
+##            value stability
+## g123  0.64922048      0.50
+## g846  0.44828969      0.50
+## g1606 0.30641602      0.40
+## g335  0.30031495      0.48
+## g836  0.26392532      0.46
+## g783  0.22022625      0.28
+## g758  0.22019689      0.38
+## g1386 0.16248889      0.38
+## g585  0.02058632      0.12
+```
+
+### Merge feature stability on a table of loadings
+
+After getting the loadings table for all the components of the model, it is
+probably desirable to create a table with stability values shown in a final
+column.  `merge.feature.stability()` can be called with the loadings table and
+stability table as parameters.  The returned data frame adds a column for the
+stability values matching the features in the loadings table.  `NA` is used
+where the stability value for a feature was not in the stability table.
+
+Output for the loadings table is truncated in the example below for brevity.
+You will see the full output in your own R instance.
+
+```R
+# Prepare the loadings table and stability value table we intend to merge on.
+loadings.table <- get.loadings.table(splsda.model)
+stability.comp3 <- feature.selection.stability(splsda.model, splsda.perf, 3)
+
+# Check the raw loadings table before merging any stability values.
+loadings.table
+##            comp1         comp2       comp3
+## g123  0.64922048            NA          NA
+## g846  0.44828969            NA          NA
+## g742          NA  4.204118e-02  0.14003503
+## g1143         NA -4.124087e-02          NA
+## g174          NA            NA -0.35418530
+## g1896         NA            NA -0.29939099
+
+# Merge on stability values for component 3.
+merge.feature.stability(loadings.table, stability.comp3)
+## feature      comp1         comp2       comp3 stability
+##    g123 0.64922048            NA          NA        NA
+##    g846 0.44828969            NA          NA        NA
+##    g742         NA  4.204118e-02  0.14003503      0.76
+##   g1143         NA -4.124087e-02          NA        NA
+##    g174         NA            NA -0.35418530      0.34
+##   g1896         NA            NA -0.29939099      0.32
+```
+
+### Measure variable importance
+
+With a trained model object `blockrank.splsda()` will create a single item list containing a vector of BlockRank scores for that model. For more detail see section [BlockRank](#blockrank) below.
+
+
+## Feature analysis for DIABLO models
+
+OmicsFold provides a number of functions for analysing the output from DIABLO
+models produced in `mixOmics`.
+
+- [`get.diablo.top.loadings()`](#get-diablo-top-loadings) outputs a data frame
+  of feature loadings for a specified block in the DIABLO model.
+- [`diablo.selection.stability()`](#get-diablo-feature-selection-stability)
+  extracts a data frame of stability scores from a DIABLO model performance test
+  for features in a particular block and component.
+- [`get.diablo.top.loadings.with.stability()`](#get-diablo-top-loadings-with-stability)
+  provides a concatenated data frame of feature loadings and stability scores
+  for a given block in the DIABLO model.
+- [`blockrank.diablo()`](#measure-feature-importance) takes a DIABLO model 
+  object and calculates a score measuring each features contribution to the variance in the response
+- [`plot.feature.stability()`](#plot-feature-stability) provides a simple
+  histogram plot of the highest stability features in a single component and
+  block, allowing visual comparison of their stability scores.
+- [`find.feature.associations()`](#find-diablo-feature-associations) outputs a
+  matrix of values indicating possible feature associations in the model based
+  on consistent feature measurement across samples.
+- [`export.matrix.as.network()`](#export-diablo-matrix-as-a-network) exports a
+  data frame showing relationships between features where feature association
+  strength exceeds a threshold.  The data frame is also written to a CSV file
+  which can be imported into a network plotting application such as
+  [Cytoscape](https://cytoscape.org/).
+
+### Get DIABLO top loadings
+
+`get.diablo.top.loadings()` allows the export of loadings from a DIABLO model
+into a data frame.  The resulting data frame has columns for each component of
+the model containing loadings for the features selected in that component.  Rows
+are sorted according to the absolute value of the loadings so that the highest
+impact features are at the top of the table.  All features for the first
+component are listed before listing begins for features of the second and
+subsequent components.  If a feature is selected in more than one component,
+only one row is given for that feature but multiple columns will contain a
+loading for that feature.  Where a feature was not selected in a component, the
+`NA` value is shown for the loading.
+
+In the example below, the loadings for the `proteomics` block are being arranged
+in the data frame.  Note that feature `AR` has been selected in both components.
+Because its loading value in component 2 is `-0.11183919`, this would normally
+place `AR` below `c-Kit`, but features are sorted by lower numbered components
+first.
+
+```R
+get.diablo.top.loadings(diablo.model$loadings$proteomics)
+##                   comp1       comp2
+## ER-alpha    -0.74995295          NA
+## GATA3       -0.62448770          NA
+## ASNS         0.16825721          NA
+## Cyclin_B1    0.12026614          NA
+## AR          -0.06842432 -0.11183919
+## JNK2        -0.01137347          NA
+## HER2                 NA -0.67182918
+## HER2_pY1248          NA -0.62130938
+## EGFR_pY1068          NA -0.33229248
+## c-Kit                NA  0.17791016
+## HER3_pY1289          NA -0.08706226
+## XRCC1                NA  0.02149536
+```
+
+### Get DIABLO feature selection stability
+
+`diablo.selection.stability()` is very similar to the [equivalent function for
+sPLS-DA models](#get-spls-da-feature-selection-stability).  By passing the
+performance test result, a selected component and a selected block as
+parameters, a data frame will be generated for the features in that block and
+component, showing the stability scores for those features.  Stability is a
+fraction indicating the number of models in the performance test which selected
+this feature when being fitted.  100 models were fitted when running the DIABLO
+performance test (10 repeats with 10 folds), so the stability scores shown are
+an exact multiple of `0.01`.  Features are presented in the data frame in
+reverse stability score order so that the most stable features are shown first.
+Those with a score of `1.0` are essential for the model to achieve the best fit
+and are therefore also likely to be relevant to the discrimination of classes
+the model is trying to achieve.
+
+The code below outputs the feature selection stability for the `proteomics`
+block in component 1.
+
+```R
+diablo.selection.stability(diablo.perf, comp = 1, block = 'proteomics')
+##   feature stability
+##  ER-alpha      1.00
+##     GATA3      1.00
+##      ASNS      0.99
+## Cyclin_B1      0.98
+##        AR      0.80
+##      JNK2      0.54
+##        PR      0.41
+## Cyclin_E1      0.22
+##    INPP4B      0.12
+```
+
+### Get DIABLO top loadings with stability
+
+`get.diablo.top.loadings.with.stability()` allows the combination of information
+from both the previous functions.  By passing in the model, performance test and
+block name, a data frame is generated with loadings of the selected features for
+each component alongside stability scores for those features.  Where a feature
+exists in more than one component, the stability score shown is for the lowest
+numbered component.
+
+The code below generates a data frame containing loading values and stability
+scores for the selected features in the `proteomics` block.  Note that, unless
+specified, the number of features returned will be up to `20` by default.
+
+```R
+get.diablo.top.loadings.with.stability(diablo.model, diablo.perf,
+                                       block = 'proteomics', feature.count = 50)
+##                   comp1       comp2 stability
+## ER-alpha    -0.74995295          NA      1.00
+## GATA3       -0.62448770          NA      1.00
+## ASNS         0.16825721          NA      0.99
+## Cyclin_B1    0.12026614          NA      0.98
+## AR          -0.06842432 -0.11183919      0.80
+## JNK2        -0.01137347          NA      0.54
+## HER2                 NA -0.67182918      1.00
+## HER2_pY1248          NA -0.62130938      1.00
+## EGFR_pY1068          NA -0.33229248      1.00
+## c-Kit                NA  0.17791016      1.00
+## HER3_pY1289          NA -0.08706226      0.95
+## XRCC1                NA  0.02149536      0.56
+```
+
+### Measure feature importance
+
+With a trained model object `blockrank.diablo()` will create a list of vectors of BlockRank scores for each feature. For more detail see section [BlockRank](#blockrank) below.
+
+### Plot feature stability
+
+`plot.feature.stability()` takes a single parameter with a table of stabilities,
+as obtained from
+[`diablo.selection.stability()`](#get-diablo-feature-selection-stability), and
+creates a bar-chart-style plot of stability for those features.  The shape of
+the plot will give and indication for how consistently the models choose the
+same features during the performance tests.  A plot that stays high for the
+first features and then drops sharply to a low level for the remainder shows
+that there is a clear preference for specific features and that those features
+should be of value to look at more closely in relation to the biological
+context.  If the curve is closer to a straight line and only a very small number
+of features achieve close to `1.0` stability, then the feature selection is less
+stable and it may be harder to draw conclusions about which are relevant to the
+biological context.
+
+The code below prepares a stability data frame for the DIABLO performance test
+for the `miRNA` block in component 2.  It then plots the stability scores as a
+bar-chart-style plot as shown below the code.  This plot shows a reasonably
+sharp drop off after about 10 features, which in turn means that the high
+stability features on the left are likely to have good relevance to the
+biological context.
+
+```R
+stability <- diablo.selection.stability(diablo.perf, comp = 2, block = 'miRNA')
+plot.feature.stability(stability)
+```
+
+![Feature Stability Plot](img/../imgs/plot-feature-stability.png)
+
+### Export DIABLO matrix as a network
+
+`export.matrix.as.network()` takes the output of the associations function above
+and both filters and converts the matrix to a flat format as a data frame as
+well as writing that to a CSV file.  The data frame has one line per association
+and does not show reverse associations.  Only associations which meet or exceed
+the given cutoff threshold are included.  The association can be positive or
+negative to satisfy inclusion by the cutoff.  The data frame shows both sides of
+the association as two feature columns and the strength of the association as a
+value column.  The block column provides the name of the block containing the
+left hand feature.
+
+The CSV file specified is created and contains the same information as the data
+frame.  The format of the CSV file is importable into programs for plotting
+networks, such as [Cytoscape](https://cytoscape.org/), which makes it easier to
+visualise how the features in the model might be relevant to each other.
+
+The code example below prepares the association matrix and then prepares a
+vector of strings where the string names are the names of the blocks each of the
+features are present in.  It is also possible to provide the parameter
+`block.feature.count` instead of the `block.association` in the unlikely event
+that each block contains the same number of features.  The output shown has been
+truncated to only show the associations included from the matrix above. The
+actual list of associations is much longer.  A CSV file named `network.csv` is
+created and populated with the information shown in the data frame output.
+
+```R
+associations <- find.feature.associations(diablo.model, block.count = 3)
+block.association <- character(0)
+for (block.name in names(diablo.model$keepX)) {
+  num.features.in.block <- diablo.model$keepX[[block.name]][1]
+  block.association <- append(block.association,
+                              rep(block.name, num.features.in.block))
+}
+export.matrix.as.network(associations, filename = "network.csv", cutoff = 0.7,
+                         block.association = block.association)
+##    feature.1    feature.2      value      block
+##        KDM4B       ZNF552  0.8247985       mRNA
+##        KDM4B        PREX1  0.7235040       mRNA
+##       ZNF552        LRIG1  0.7072751       mRNA
+##       ZNF552        PREX1  0.7308222       mRNA
+```
+
+## Network functions
+OmicsFold provides a set of functions to extract and process data from a DIABLO model 
+into a multi-omic interaction network visualization. The network functions identify associations between 
+top features in the model and reformat them into a network object 
+for direct visualization in R or in an external
+program such as Cytoscape.
+
+### Extract feature associations
+
+To extract feature associations from a multi-omic DIABLO model, use the 
+`find.feature.associations()` function.
+
+```R
+associations <- find.feature.associations(diablo.model, feature_number=50, score_type="blockrank")
+```
+
+The function requires the final trained DIABLO model. The `feature_number` parameter can be 
+adjusted to select the number of top discriminative features to be returned. The `score_type` function
+can be set to either "blockrank" or "loading" to determine which feature ranking metric to use.
+The function will return a matrix of associations between the most important features in the model according to either 
+blockrank or loading scores.
+
+### Filter and reformat network
+
+To filter and reformat the association matrix for network visualization, use the 
+`filter.network()` function.
+
+```R
+omicsfold_network <- filter.network(diablo.model=diablo.model, associations=associations, cutoff=0.7, 
+feature_list=c("phenol.sulfate", "Alistipes_putredinis"), remove_intrablock = FALSE)
+```
+
+The function requires the final trained DIABLO model and the association matrix from `find.feature.associations()`.
+A correlation cut off can be set with `cutoff` to remove associations between features that 
+do not meet the supplied threshold. A vector of feature names can be supplied to filter the network to only include
+associations to features of interest with `feature_list`, and associations between features from the same block
+can be removed with `remove_intrablock`. The function returns a dataframe describing associations 
+between features and the blocks that they belong to. The dataframe can
+be exported to Cytoscape or applied to the `plot.network()` function for visualization in R. 
+
+### Plot network
+
+To visualize the multi-omic network, use the `plot.network()` function. 
+
+```R
+plot.network(omicsfold_network=omicsfold_network, diablo.model=diablo.model)
+```
+
+![image](https://user-images.githubusercontent.com/49563541/182643757-33a50d8d-153c-496e-af17-dab19cddcdb8.png)
+
+The function requires the network dataframe from the `filter.network()` function and the final trained DIABLO model. Nodes in the network represent 
+features and edges represent correlations between nodes. Nodes are colored by block and edges are colored by correlation value, with darker red edges 
+indicating the higher end of the correlation values in the network and darker orange edges indicating the lower end of correlation values in the network.
+The function uses the Kamada-Kawai network layout algorithm to position nodes so that more strongly connected nodes will be pulled together,
+while more weakly connected nodes will pushed away from other nodes.
+
+## BlockRank
+
+OmicsFold implements a novel approach to feature analysis called BlockRank. 
+BlockRank allows the direct comparison of feature importance between different 
+data blocks, and can also be used to measure feature importance in a singleomic 
+setting.
+
+The BlockRank Score for each feature is a number between 0 and 1, with 0 
+indicating an insignificant feature and 1 indicating the most important feature 
+in the model.
+
+### BlockRank for sPLS-DA
+
+To measure the BlockRank score for a single-omic sPLS-DA model use the 
+`blockrank.splsda()` function.
+
+```R
+blockrank.scores <- blockrank.splsda(splsda.model)
+```
+
+`blockrank.scores` is a list of one numeric vector, the vector containing a 
+score for each feature in the model data.
+
+### BlockRank for DIABLO
+
+To measure the BlockRank score for a multi-omic DIABLO model use the 
+`blockrank.diablo()` function.
+
+```R
+blockrank.scores <- blockrank.diablo(diablo.model)
+```
+
+`blockrank.scores` is a list of numeric vectors. Each vector corresponds to one 
+block of data, and the vector contains a score for each feature in the block.
+
+### Plotting BlockRank Scores
+
+The function `plot.blockrank.scores()` returns a bar chart of the top 'n' BlockRank scores.
+
+```R
+plot.blockrank.scores(blockrank.scores,
+             nscores = 20,
+             feature.font.size = 8,
+             model = "DIABLO",
+             data_source = "Breast TCGA")
+```
+![Example Blockrank Plot](imgs/blockrank_diablo.png)
+
+The function requires a list of BlockRank scores to plot, and optionally has parameters for how many scores are displayed, `nscores`; the font size of the feature labels, `feature.font.size`; the name of type of model plotted for the title of the graph, `model`, and the name of the data source for the graph title, `data_source`.
+
+## Model predictivity
+
+`plot.predicted.projection()` is useful when an sPLS-DA single-omics model is
+used to predict the classifications of unseen data and you want to visualise the
+accuracy of that prediction.  The built in `predict()` function from `mixOmics`
+will take a trained model and some unseen observations of the same set of
+features, returning an object which is complicated to parse, but contains the
+predicted classifications of the new data.
+
+The projection plot generated by our function contains a large centroid point in
+the first two components for each of the classifications surrounded by smaller
+points for each observation in the unseen set.  The smaller points are coloured
+according to their actual classification.  A predictive model will render the
+smaller points nearer to the large points of the same colour.  Smaller points of
+the wrong colour compared to their closest large point were predicted
+incorrectly.
+
+The code example below prepares a new sPLS-DA model from a smaller proportion of
+the full dataset, described as the training data set.  This training data set is
+used to fit the predictive model and then that model is made to predict the
+outcome of the held back observations.
+
+Printing the projection plot generated by our function renders the plot shown
+below the code.  Here we can see that the predictions are mostly accurate with
+the exception of one `EWS` point being too close to the `RMS` centroid.  There
+are also some other points that are situated approximately halfway between two
+centroids, but they still appear to be closest to the correct centroid.  The
+model is therefore able to reasonably accurately predict the classification of
+the unseen set with at least 1 error shown.
+
+```R
+# Prepare data
+data(srbct)
+all.X <- srbct$gene
+all.Y <- srbct$class
+all.index <- 1:length(all.Y)
+
+# Subset the data for a training set
+set.seed(1234) # for reproducibility of subsetting
+train.index <- sort(sample(all.index, length(all.index) * 0.8))
+train.X <- all.X[train.index,]
+train.Y <- all.Y[train.index]
+
+# Model fit with the training set
+predict.model <- splsda(train.X, train.Y, ncomp = 3, keepX = c(9, 280, 30))
+
+# Prepare the prediction set
+predict.index <- setdiff(all.index, train.index)
+predict.X <- all.X[predict.index,]
+predict.Y <- all.Y[predict.index]
+
+# Create a prediction for the predict set and plot using plot.predicted.projection()
+prediction <- predict(predict.model, predict.X)
+projection.plot <- plot.predicted.projection(prediction, predict.Y)
+print(projection.plot)
+```
+
+![Prediction Projection Plot](imgs/plot-predicted-projection.png)
+
+## Utility functions
+
+`center.truncate()` is a helper function for plot labels, though it may have
+other uses.  Where labels are long, as is typical for features of some -omics
+data, this can cause plots to be dominated by axis labels.  By passing a string
+to this function, it will be truncated to exactly 43 characters or fewer.
+Characters in the middle of the string will be replaced by `...` instead.
+
+The examples below show both a single string being truncated, as well as how to
+apply the truncation to a vector of strings using the built in R function
+`sapply()`.
+
+```R
+center.truncate("When a string is particularly long it will be truncated back to 43 characters")
+## "When a string is par...back to 43 characters"
+
+labels = c("Short name",
+           "Very long name would need truncating for plotting",
+           "Plots don't work well when the axis is forced over by long labels")
+unname(sapply(labels, center.truncate))
+## [1] "Short name"
+## [2] "Very long name would...uncating for plotting"
+## [3] "Plots don't work wel...d over by long labels"
+```