---
title: Introduction to Numero
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 2
    number_sections: true
vignette: >
  %\VignetteIndexEntry{intro}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
body {
  text-align: justify;
}
p.caption {
  display: inline-block;
  padding-left: 0.5em;
  padding-right: 0.5em;
  border-color: #F4F4F4;
  border-width: 1px;
  border-style: solid;
  border-radius: 4px;
  background-color: #F6F6F6;
  text-align: justify;
}
div.caption {
  display: inline-block;
  padding-left: 0.5em;
  padding-right: 0.5em;
  border-color: #F4F4F4;
  border-width: 1px;
  border-style: solid;
  border-radius: 4px;
  background-color: #F6F6F6;
  text-align: justify;
}
h1 {
  margin-top: 1em;
  margin-bottom: 0.1em;
}
h2 {
  margin-top: 0.2em;
  margin-bottom: 0.1em;
}
</style>

**Author:** Ville-Petteri Mäkinen  
**Co-authors:** Song Gao, Stefan Mutter, Aaron E Casey  
**Version:** `r packageVersion("Numero")`

**Abstract:** In textbook examples, multivariable datasets are clustered into
distinct subgroups that can be clearly identified by a set of optimal
mathematical criteria. However, many real-world datasets arise from
synergistic consequences of multiple effects, noisy and partly redundant
measurements, and may represent a continuous spectrum of the different phases
of a phenomenon. In medicine, complex diseases associated with ageing are
typical examples. We postulate that population-based biomedical datasets (and
many other real-world examples) do not contain an intrinsic clustered
structure that would give rise to mathematically well-defined subgroups. From
a modeling point of view, the lack of intrinsic structure means that the data
points inhabit a contiguous cloud in high-dimensional space without abrupt
changes in density to indicate subgroup boundaries, hence a mathematical
criteria cannot segment the cloud reliably by its internal structure. Yet we
need data-driven classification and subgrouping to aid decision-making and to
facilitate the development of testable hypotheses. For this reason, we
developed the Numero package, a more flexible and transparent process that
allows human observers to create usable multivariable subgroups even when
conventional clustering frameworks struggle.

**Citation:** Gao S, Mutter S, Casey AE, Mäkinen V-P (2018) Numero: a
statistical framework to define multivariable subgroups in complex
population-based datasets, Int J Epidemiology,
https://doi.org/10.1093/ije/dyy113

***

<CENTER>
<svg width="700" height="170" version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:cc="http://creativecommons.org/ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
 <defs>
  <marker id="Arrow2Mend" overflow="visible" orient="auto">
   <path transform="scale(-.6)" d="m8.7186 4.0337-10.926-4.0177 10.926-4.0177c-1.7455 2.3721-1.7354 5.6175-6e-7 8.0354z" fill-rule="evenodd" stroke-linejoin="round" stroke-width=".625"/>
  </marker>
 </defs>
 <metadata>
  <rdf:RDF>
   <cc:Work rdf:about="">
    <dc:format>image/svg+xml</dc:format>
    <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/>
    <dc:title/>
   </cc:Work>
  </rdf:RDF>
 </metadata>
 <g transform="translate(0 -882.36)">
  <rect y="882.36" width="700" height="170" fill="#fff"/>
 </g>
 <g transform="translate(0 -882.36)">
  <g transform="matrix(1.6924 0 0 1.6924 -7.5146 -727.96)">
   <g transform="matrix(.4203 .22499 -.22499 .4203 425.81 591.74)" stroke-linecap="round">
    <g transform="translate(253.07 859.01)" stroke-width=".41952px">
     <circle r="11.768" fill="#fffefc" stroke="#fffefc"/>
     <path d="m10.616-5.0796a11.768 11.768 0 0 0-7.9683-6.3871l4.4688-19.356a31.634 31.634 0 0 1 21.419 17.169z" fill="#ffd89f" stroke="#ffd89f"/>
     <path d="m10.59 5.1325a11.768 11.768 0 0 0 0.0255-10.212l17.92-8.5747a31.634 31.634 0 0 1-0.0685 27.451z" fill="#ff8f52" stroke="#ff8f52"/>
     <path d="m2.5901 11.48a11.768 11.768 0 0 0 8-6.3473l17.877 8.6639a31.634 31.634 0 0 1-21.504 17.062z" fill="#fff6e9" stroke="#fff6e9"/>
     <path d="m-7.3604 9.1825a11.768 11.768 0 0 0 9.9505 2.2973l4.3722 19.378a31.634 31.634 0 0 1-26.747-6.1751z" fill="#ace3ff" stroke="#ace3ff"/>
     <path d="m-11.768-0.0293a11.768 11.768 0 0 0 4.4079 9.2118l-12.425 15.501a31.634 31.634 0 0 1-11.849-24.762z" fill="#75ceff" stroke="#75ceff"/>
     <path d="m-7.3145-9.2191a11.768 11.768 0 0 0-4.4538 9.1898l-19.866-0.0496a31.634 31.634 0 0 1 11.972-24.703z" fill="#b2e5ff" stroke="#b2e5ff"/>
     <path d="m2.6473-11.467a11.768 11.768 0 0 0-9.9618 2.2476l-12.347-15.562a31.634 31.634 0 0 1 26.778-6.0416z" fill="#e9f8ff" stroke="#e9f8ff"/>
     <path d="m-13.523-28.598a31.634 31.634 0 0 0-11.741 9.5601l-15.865-11.955a51.499 51.499 0 0 1 19.114-15.564z" fill="#a5e1ff" stroke="#a5e1ff"/>
     <path d="m1.3162-31.606a31.634 31.634 0 0 0-14.839 3.0087l-8.4921-17.959a51.499 51.499 0 0 1 24.158-4.898z" fill="#a3e0ff" stroke="#a3e0ff"/>
     <path d="m15.854-27.374a31.634 31.634 0 0 0-14.538-4.232l0.8265-19.848a51.499 51.499 0 0 1 23.667 6.8896z" fill="#fffbf5" stroke="#fffbf5"/>
     <path d="m26.759-16.871a31.634 31.634 0 0 0-10.906-10.503l9.9559-17.191a51.499 51.499 0 0 1 17.754 17.099z" fill="#ff6817" stroke="#ff6817"/>
     <path d="m31.535-2.5031a31.634 31.634 0 0 0-4.7754-14.368l16.804-10.595a51.499 51.499 0 0 1 7.7742 23.391z" fill="#dc0f28" stroke="#dc0f28"/>
     <path d="m29.086 12.438a31.634 31.634 0 0 0 2.4489-14.942l19.803-1.572a51.499 51.499 0 0 1-3.9866 24.325z" fill="#ff6817" stroke="#ff6817"/>
     <path d="m19.974 24.531a31.634 31.634 0 0 0 9.112-12.092l18.266 7.8111a51.499 51.499 0 0 1-14.834 19.686z" fill="#ffedd4" stroke="#ffedd4"/>
     <path d="m6.286 31.003a31.634 31.634 0 0 0 13.688-6.4724l12.543 15.405a51.499 51.499 0 0 1-22.284 10.537z" fill="#cdeeff" stroke="#cdeeff"/>
     <path d="m-8.8419 30.373a31.634 31.634 0 0 0 15.128 0.6299l3.9474 19.469a51.499 51.499 0 0 1-24.628-1.0256z" fill="#6fc4ff" stroke="#6fc4ff"/>
     <path d="m-21.944 22.785a31.634 31.634 0 0 0 13.102 7.5881l-5.5525 19.074a51.499 51.499 0 0 1-21.33-12.353z" fill="#4891ff" stroke="#4891ff"/>
     <path d="m-30.019 9.9771a31.634 31.634 0 0 0 8.0751 12.808l-13.78 14.309a51.499 51.499 0 0 1-13.146-20.851z" fill="#1d76ff" stroke="#1d76ff"/>
     <path d="m-31.217-5.1164a31.634 31.634 0 0 0 1.1981 15.094l-18.852 6.2655a51.499 51.499 0 0 1-1.9505-24.572z" fill="#5096ff" stroke="#5096ff"/>
     <path d="m-25.264-19.038a31.634 31.634 0 0 0-5.9535 13.921l-19.604-3.2129a51.499 51.499 0 0 1 9.6921-22.664z" fill="#83d6ff" stroke="#83d6ff"/>
     <path d="m50.709-8.9894a51.499 51.499 0 0 0-5.6664-15.978l17.375-9.631a71.365 71.365 0 0 1 7.8522 22.141z" fill="#dc0f28" stroke="#dc0f28"/>
     <path d="m50.88 7.9628a51.499 51.499 0 0 0-0.1713-16.952l19.561-3.4676a71.365 71.365 0 0 1 0.2373 23.491z" fill="#fb5105" stroke="#fb5105"/>
     <path d="m45.538 24.052a51.499 51.499 0 0 0 5.3424-16.089l19.627 3.0716a71.365 71.365 0 0 1-7.4031 22.296z" fill="#ffb169" stroke="#ffb169"/>
     <path d="m35.261 37.535a51.499 51.499 0 0 0 10.277-13.483l17.566 9.278a71.365 71.365 0 0 1-14.241 18.684z" fill="#fff2de" stroke="#fff2de"/>
     <path d="m21.163 46.95a51.499 51.499 0 0 0 14.098-9.4154l13.602 14.479a71.365 71.365 0 0 1-19.536 13.047z" fill="#cfefff" stroke="#cfefff"/>
     <path d="m4.7712 51.278a51.499 51.499 0 0 0 16.391-4.3275l8.1633 18.111a71.365 71.365 0 0 1-22.714 5.9969z" fill="#74cdff" stroke="#74cdff"/>
     <path d="m-12.137 50.049a51.499 51.499 0 0 0 16.908 1.2292l1.8404 19.78a71.365 71.365 0 0 1-23.431-1.7033z" fill="#549cff" stroke="#549cff"/>
     <path d="m-27.73 43.396a51.499 51.499 0 0 0 15.593 6.6527l-4.6819 19.306a71.365 71.365 0 0 1-21.608-9.219z" fill="#106eff" stroke="#106eff"/>
     <path d="m-40.319 32.041a51.499 51.499 0 0 0 12.588 11.355l-10.697 16.74a71.365 71.365 0 0 1-17.444-15.736z" fill="#1437ee" stroke="#1437ee"/>
     <path d="m-48.538 17.213a51.499 51.499 0 0 0 8.219 14.828l-15.553 12.36a71.365 71.365 0 0 1-11.389-20.547z" fill="#2218e1" stroke="#2218e1"/>
     <path d="m-51.497 0.5204a51.499 51.499 0 0 0 2.9592 16.693l-18.723 6.6399a71.365 71.365 0 0 1-4.1007-23.132z" fill="#0a4ef6" stroke="#0a4ef6"/>
     <path d="m-48.876-16.229a51.499 51.499 0 0 0-2.6213 16.749l-19.864 0.2007a71.365 71.365 0 0 1 3.6324-23.21z" fill="#5da9ff" stroke="#5da9ff"/>
     <path d="m-40.958-31.219a51.499 51.499 0 0 0-7.9177 14.991l-18.853-6.2602a71.365 71.365 0 0 1 10.972-20.773z" fill="#c5ecff" stroke="#c5ecff"/>
     <path d="m-28.602-42.827a51.499 51.499 0 0 0-12.356 11.607l-15.799-12.043a71.365 71.365 0 0 1 17.122-16.085z" fill="#bce9ff" stroke="#bce9ff"/>
     <path d="m-13.146-49.793a51.499 51.499 0 0 0-15.456 6.9665l-11.033-16.52a71.365 71.365 0 0 1 21.417-9.6538z" fill="#73caff" stroke="#73caff"/>
     <g>
      <path d="m3.734-51.364a51.499 51.499 0 0 0-16.88 1.5706l-5.0711-19.207a71.365 71.365 0 0 1 23.392-2.1764z" fill="#e9f8ff" stroke="#e9f8ff"/>
      <path d="m20.21-47.368a51.499 51.499 0 0 0-16.476-3.9955l1.4404-19.813a71.365 71.365 0 0 1 22.831 5.5367z" fill="#fff6e9" stroke="#fff6e9"/>
      <path d="m34.495-38.24a51.499 51.499 0 0 0-14.286-9.1286l7.7957-18.272a71.365 71.365 0 0 1 19.796 12.65z" fill="#fb5105" stroke="#fb5105"/>
     </g>
     <path d="m45.042-24.967a51.499 51.499 0 0 0-10.547-13.272l13.306-14.751a71.365 71.365 0 0 1 14.616 18.392z" fill="#dc0f28" stroke="#dc0f28"/>
    </g>
    <circle transform="translate(253.07 859.01)" r="71.365" fill="none" stroke="#b2b2b2" stroke-width="1.2586px"/>
   </g>
   <g letter-spacing="0px" word-spacing="0px">
    <g font-family="'Nimbus Sans L'">
     <text x="337.21732" y="966.07263" fill="#000000" text-align="center" text-anchor="middle" style="line-height:0%" xml:space="preserve"><tspan x="337.21732" y="966.07263" font-family="'Nimbus Sans L'" font-size="10px" style="line-height:1.04999995">Map coloring</tspan></text>
     <text x="378.1792" y="992.26324" fill="#dc0f28" font-size="10px" style="line-height:0%" xml:space="preserve"><tspan x="378.1792" y="992.26324" style="line-height:1.04999995">High</tspan><tspan x="378.1792" y="1002.7632" style="line-height:1.04999995">regional</tspan><tspan x="378.1792" y="1013.2632" style="line-height:1.04999995">average</tspan><tspan x="378.1792" y="1023.7632" style="line-height:1.04999995">(red)</tspan></text>
     <text x="264.51056" y="998.60303" fill="#0a4ef6" font-size="10px" style="line-height:0%" xml:space="preserve"><tspan x="264.51056" y="998.60303" style="line-height:1.04999995">Low</tspan><tspan x="264.51056" y="1009.103" style="line-height:1.04999995">regional</tspan><tspan x="264.51056" y="1019.603" style="line-height:1.04999995">average</tspan><tspan x="264.51056" y="1030.103" style="line-height:1.04999995">(blue)</tspan></text>
    </g>
    <text x="105.95399" y="1026.8136" fill="#000000" font-family="sans-serif" style="line-height:0%" xml:space="preserve"><tspan x="105.95399" y="1026.8136" font-size="16px" style="line-height:1.25"> </tspan></text>
   </g>
   <path d="m100.92 1040-13.147-7.4044" fill="none" marker-end="url(#Arrow2Mend)" stroke="#000"/>
   <g font-family="'Nimbus Sans L'" letter-spacing="0px" word-spacing="0px">
    <text x="103.0471" y="1043.9761" fill="#000000" style="line-height:0%" xml:space="preserve"><tspan x="103.0471" y="1043.9761" font-size="10px" style="line-height:1.04999995">Data point</tspan></text>
    <text x="120.05804" y="990.17365" fill="#dc0f28" font-size="10px" style="line-height:0%" xml:space="preserve"><tspan x="120.05804" y="990.17365" style="line-height:1.04999995">High</tspan><tspan x="120.05804" y="1000.6736" style="line-height:1.04999995">value</tspan><tspan x="120.05804" y="1011.1736" style="line-height:1.04999995">(red)</tspan></text>
    <text x="7.8647575" y="1016.7698" fill="#0a4ef6" font-size="10px" style="line-height:0%" xml:space="preserve"><tspan x="7.8647575" y="1016.7698" style="line-height:1.04999995">Low</tspan><tspan x="7.8647575" y="1027.2698" style="line-height:1.04999995">value</tspan><tspan x="7.8647575" y="1037.7698" style="line-height:1.04999995">(blue)</tspan></text>
   </g>
   <g transform="translate(20.686 55.081)">
    <circle transform="rotate(27.522)" cx="462.66" cy="794.95" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="467.5" cy="825.41" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="471.76" cy="826.37" r="1.4011" fill="#f00"/>
    <g fill="#00f">
     <circle transform="rotate(27.522)" cx="464.16" cy="832.31" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="475.84" cy="828.12" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="471.27" cy="840.74" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="478.86" cy="844.81" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="465.82" cy="856.99" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="474.13" cy="850.43" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="499.76" cy="821.4" r="1.4011"/>
    </g>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="506.85" cy="819.6" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="498.3" cy="826.17" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="506.71" cy="825.25" r="1.4011"/>
    </g>
    <circle transform="rotate(27.522)" cx="491.89" cy="785.13" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="479.65" cy="803.67" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="483.4" cy="788.93" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="472.73" cy="813.52" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="449.58" cy="821.23" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="445.9" cy="804.81" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="484.4" cy="821.7" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="479.96" cy="823.51" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="479.93" cy="819.57" r="1.4011" fill="#00f"/>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="523.05" cy="796.02" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="515.7" cy="816.28" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="515.72" cy="804.78" r="1.4011"/>
    </g>
    <circle transform="rotate(27.522)" cx="475.99" cy="805.91" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="479.84" cy="809.65" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="465.42" cy="808.17" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="475.79" cy="816.76" r="1.4011" fill="#f00"/>
    <g fill="#00f">
     <circle transform="rotate(27.522)" cx="449.26" cy="831.42" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="459.51" cy="839.23" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="452.01" cy="847.12" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="464.98" cy="842.4" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="481.14" cy="836.13" r="1.4011"/>
    </g>
    <circle transform="rotate(27.522)" cx="486.19" cy="833.97" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="487.64" cy="843.11" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="495.76" cy="841.91" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="487.75" cy="829.24" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="497.32" cy="830" r="1.4011" fill="#00f"/>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="500.81" cy="807.12" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="495.69" cy="814.68" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="494.11" cy="800.07" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="501.76" cy="787.87" r="1.4011"/>
    </g>
    <circle transform="rotate(27.522)" cx="455.6" cy="802.04" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="456.5" cy="826.87" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="471.02" cy="821.41" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="481.83" cy="817.57" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="488.73" cy="818.49" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="491.81" cy="809.7" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="491.74" cy="827.51" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="485.27" cy="826.69" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="479.71" cy="829.42" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="482.4" cy="850.68" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="502.46" cy="847.7" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="494.83" cy="833.68" r="1.4011" fill="#00f"/>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="507.7" cy="835.64" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="515.1" cy="828.88" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="518.36" cy="824.89" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="502.5" cy="812.86" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="499.97" cy="802.91" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="509.8" cy="798.56" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="489.18" cy="799.67" r="1.4011"/>
    </g>
    <circle transform="rotate(27.522)" cx="485.45" cy="805.25" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="486.78" cy="809.75" r="1.4011" fill="#f00"/>
    <circle transform="rotate(27.522)" cx="476.25" cy="790.72" r="1.4011" fill="#00f"/>
    <circle transform="rotate(27.522)" cx="487.49" cy="822.74" r="1.4011" fill="#f00"/>
   </g>
   <text x="65.429413" y="965.96899" fill="#000000" font-family="'Nimbus Sans L'" letter-spacing="0px" text-align="center" text-anchor="middle" word-spacing="0px" style="line-height:0%" xml:space="preserve"><tspan x="65.429413" y="965.96899" font-family="'Nimbus Sans L'" font-size="10px" style="line-height:1.04999995">Scatter plot</tspan></text>
   <g transform="translate(171.28 -100.2)">
    <g transform="matrix(.4639 .24171 -.24171 .4639 1.0384 710.29)" fill="none" stroke="#999" stroke-linecap="round" stroke-linejoin="round" stroke-width="1.6727">
     <circle transform="translate(405.16 652.54)" r="11.768"/>
     <path d="m415.78 647.46a11.768 11.768 0 0 0-7.9683-6.3871l4.4688-19.356a31.634 31.634 0 0 1 21.419 17.169z"/>
     <path d="m415.75 657.68a11.768 11.768 0 0 0 0.0255-10.212l17.92-8.5747a31.634 31.634 0 0 1-0.0685 27.451z"/>
     <path d="m407.75 664.02a11.768 11.768 0 0 0 8-6.3473l17.877 8.6639a31.634 31.634 0 0 1-21.504 17.062z"/>
     <path d="m397.8 661.73a11.768 11.768 0 0 0 9.9505 2.2973l4.3722 19.378a31.634 31.634 0 0 1-26.747-6.1751z"/>
     <path d="m393.39 652.51a11.768 11.768 0 0 0 4.4079 9.2118l-12.425 15.501a31.634 31.634 0 0 1-11.849-24.762z"/>
     <path d="m397.85 643.32a11.768 11.768 0 0 0-4.4538 9.1898l-19.866-0.0496a31.634 31.634 0 0 1 11.972-24.703z"/>
     <path d="m407.81 641.08a11.768 11.768 0 0 0-9.9618 2.2476l-12.347-15.562a31.634 31.634 0 0 1 26.778-6.0416z"/>
     <path d="m391.64 623.94a31.634 31.634 0 0 0-11.741 9.5601l-15.865-11.955a51.499 51.499 0 0 1 19.114-15.564z"/>
     <path d="m406.48 620.94a31.634 31.634 0 0 0-14.839 3.0087l-8.4921-17.959a51.499 51.499 0 0 1 24.158-4.898z"/>
     <path d="m421.02 625.17a31.634 31.634 0 0 0-14.538-4.232l0.8265-19.848a51.499 51.499 0 0 1 23.667 6.8896z"/>
     <path d="m431.92 635.67a31.634 31.634 0 0 0-10.906-10.503l9.9559-17.191a51.499 51.499 0 0 1 17.754 17.099z"/>
     <path d="m436.7 650.04a31.634 31.634 0 0 0-4.7754-14.368l16.804-10.595a51.499 51.499 0 0 1 7.7742 23.391z"/>
     <path d="m434.25 664.98a31.634 31.634 0 0 0 2.4489-14.942l19.803-1.572a51.499 51.499 0 0 1-3.9866 24.325z"/>
     <path d="m425.14 677.07a31.634 31.634 0 0 0 9.112-12.092l18.266 7.8111a51.499 51.499 0 0 1-14.834 19.686z"/>
     <path d="m411.45 683.55a31.634 31.634 0 0 0 13.688-6.4724l12.543 15.405a51.499 51.499 0 0 1-22.284 10.537z"/>
     <path d="m396.32 682.92a31.634 31.634 0 0 0 15.128 0.6299l3.9474 19.469a51.499 51.499 0 0 1-24.628-1.0256z"/>
     <path d="m383.22 675.33a31.634 31.634 0 0 0 13.102 7.5881l-5.5525 19.074a51.499 51.499 0 0 1-21.33-12.353z"/>
     <path d="m375.14 662.52a31.634 31.634 0 0 0 8.0751 12.808l-13.78 14.309a51.499 51.499 0 0 1-13.146-20.851z"/>
     <path d="m373.95 647.43a31.634 31.634 0 0 0 1.1981 15.094l-18.852 6.2655a51.499 51.499 0 0 1-1.9505-24.572z"/>
     <path d="m379.9 633.5a31.634 31.634 0 0 0-5.9535 13.921l-19.604-3.2129a51.499 51.499 0 0 1 9.6921-22.664z"/>
     <path d="m455.87 643.55a51.499 51.499 0 0 0-5.6664-15.978l17.375-9.631a71.365 71.365 0 0 1 7.8522 22.141z"/>
     <path d="m456.04 660.51a51.499 51.499 0 0 0-0.1713-16.952l19.561-3.4676a71.365 71.365 0 0 1 0.2373 23.491z"/>
     <path d="m450.7 676.59a51.499 51.499 0 0 0 5.3424-16.089l19.627 3.0716a71.365 71.365 0 0 1-7.4031 22.296z"/>
     <path d="m440.42 690.08a51.499 51.499 0 0 0 10.277-13.483l17.566 9.278a71.365 71.365 0 0 1-14.241 18.684z"/>
     <path d="m426.33 699.49a51.499 51.499 0 0 0 14.098-9.4154l13.602 14.479a71.365 71.365 0 0 1-19.536 13.047z"/>
     <path d="m409.93 703.82a51.499 51.499 0 0 0 16.391-4.3275l8.1633 18.111a71.365 71.365 0 0 1-22.714 5.9969z"/>
     <path d="m393.03 702.59a51.499 51.499 0 0 0 16.908 1.2292l1.8404 19.78a71.365 71.365 0 0 1-23.431-1.7033z"/>
     <path d="m377.43 695.94a51.499 51.499 0 0 0 15.593 6.6527l-4.6819 19.306a71.365 71.365 0 0 1-21.608-9.219z"/>
     <path d="m364.84 684.58a51.499 51.499 0 0 0 12.588 11.355l-10.697 16.74a71.365 71.365 0 0 1-17.444-15.736z"/>
     <path d="m356.62 669.76a51.499 51.499 0 0 0 8.219 14.828l-15.553 12.36a71.365 71.365 0 0 1-11.389-20.547z"/>
     <path d="m353.67 653.06a51.499 51.499 0 0 0 2.9592 16.693l-18.723 6.6399a71.365 71.365 0 0 1-4.1007-23.132z"/>
     <path d="m356.29 636.31a51.499 51.499 0 0 0-2.6213 16.749l-19.864 0.2007a71.365 71.365 0 0 1 3.6324-23.21z"/>
     <path d="m364.2 621.32a51.499 51.499 0 0 0-7.9177 14.991l-18.853-6.2602a71.365 71.365 0 0 1 10.972-20.773z"/>
     <path d="m376.56 609.72a51.499 51.499 0 0 0-12.356 11.607l-15.799-12.043a71.365 71.365 0 0 1 17.122-16.085z"/>
     <path d="m392.02 602.75a51.499 51.499 0 0 0-15.456 6.9665l-11.033-16.52a71.365 71.365 0 0 1 21.417-9.6538z"/>
     <path d="m408.9 601.18a51.499 51.499 0 0 0-16.88 1.5706l-5.0711-19.207a71.365 71.365 0 0 1 23.392-2.1764z"/>
     <path d="m425.37 605.17a51.499 51.499 0 0 0-16.476-3.9955l1.4404-19.813a71.365 71.365 0 0 1 22.831 5.5367z"/>
     <path d="m439.66 614.3a51.499 51.499 0 0 0-14.286-9.1286l7.7957-18.272a71.365 71.365 0 0 1 19.796 12.65z"/>
     <path d="m450.2 627.58a51.499 51.499 0 0 0-10.547-13.272l13.306-14.751a71.365 71.365 0 0 1 14.616 18.392z"/>
     <circle transform="translate(405.16 652.54)" r="71.365"/>
    </g>
    <g>
     <circle transform="rotate(27.522)" cx="519.29" cy="947.04" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="520.78" cy="980.11" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="524.97" cy="980.11" r="1.4011" fill="#f00"/>
    </g>
    <g fill="#00f">
     <circle transform="rotate(27.522)" cx="520.78" cy="984.41" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="524.97" cy="984.41" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="526.57" cy="998.23" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="530.76" cy="998.23" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="526.57" cy="1002.5" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="530.76" cy="1002.5" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="560.2" cy="972.08" r="1.4011"/>
    </g>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="564.38" cy="972.08" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="560.2" cy="976.37" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="564.38" cy="976.37" r="1.4011"/>
    </g>
    <g>
     <circle transform="rotate(27.522)" cx="540.21" cy="937.14" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="538.34" cy="940.88" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="536.29" cy="937.14" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="511.25" cy="964.23" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="509.39" cy="967.97" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="507.33" cy="964.23" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="532.55" cy="973.94" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="530.68" cy="977.68" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="528.63" cy="973.94" r="1.4011" fill="#00f"/>
    </g>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="572.9" cy="958.62" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="571.04" cy="962.36" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="568.98" cy="958.62" r="1.4011"/>
    </g>
    <g>
     <circle transform="rotate(27.522)" cx="525.08" cy="953.39" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="529.26" cy="953.39" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="519.47" cy="962.17" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="523.66" cy="962.17" r="1.4011" fill="#f00"/>
    </g>
    <g fill="#00f">
     <circle transform="rotate(27.522)" cx="511.07" cy="986.09" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="515.25" cy="986.09" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="517.42" cy="994.49" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="521.6" cy="994.49" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="538.34" cy="992.63" r="1.4011"/>
    </g>
    <g>
     <circle transform="rotate(27.522)" cx="542.53" cy="992.63" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="547.31" cy="1001.6" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="551.49" cy="1001.6" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="548.43" cy="990.76" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="552.62" cy="990.76" r="1.4011" fill="#00f"/>
    </g>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="550.11" cy="970.58" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="554.3" cy="970.58" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="556.46" cy="943.49" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="560.65" cy="943.49" r="1.4011"/>
    </g>
    <g>
     <circle transform="rotate(27.522)" cx="512.81" cy="955.82" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="509.63" cy="976.37" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="519.72" cy="972.51" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="530.99" cy="965.72" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="538.96" cy="959.68" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="548.43" cy="961.8" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="548.37" cy="979.61" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="538.84" cy="981.67" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="530.31" cy="989.64" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="539.03" cy="1002.8" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="559.27" cy="997.42" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="558.71" cy="983.85" r="1.4011" fill="#00f"/>
    </g>
    <g fill="#f00">
     <circle transform="rotate(27.522)" cx="566.68" cy="990.07" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="571.72" cy="980.98" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="573.46" cy="970.52" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="561.76" cy="963.73" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="556.59" cy="955.01" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="566.43" cy="950.65" r="1.4011"/>
     <circle transform="rotate(27.522)" cx="548.99" cy="939.57" r="1.4011"/>
    </g>
    <g>
     <circle transform="rotate(27.522)" cx="547.37" cy="950.03" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="537.1" cy="949.47" r="1.4011" fill="#f00"/>
     <circle transform="rotate(27.522)" cx="527.88" cy="941.44" r="1.4011" fill="#00f"/>
     <circle transform="rotate(27.522)" cx="541.39" cy="970.7" r="1.4011" fill="#f00"/>
    </g>
   </g>
   <text x="201.5592" y="965.98694" fill="#000000" font-family="'Nimbus Sans L'" letter-spacing="0px" text-align="center" text-anchor="middle" word-spacing="0px" style="line-height:0%" xml:space="preserve"><tspan x="201.5592" y="965.98694" font-family="'Nimbus Sans L'" font-size="10px" style="line-height:1.04999995">Structured 2D layout</tspan></text>
  </g>
 </g>
</svg>
</CENTER>

``` {r echo=FALSE, results="hide"}
# Save space on screen printouts.
options(digits = 3)
```

***


---
###########################################################################
---

# Getting started {.tabset}

## Installation

You can install Numero using the standard procedure:

``` {r eval=FALSE}
# Install the package from a remote repository.
install.packages("Numero")
```

The functions come in two flavors: the ones starting with `numero` provide
high level pipeline components that will serve most users, whereas the `nro`
functions perform specific tasks and provide a more granular interface to
the package. In this introductory document, we will only use the `numero`
group of functions.

To activate and list the library functions, type

``` {r}
# Activate the library.
library("Numero")
packageVersion("Numero")
ls("package:Numero")
```

Each Numero function comes with a help page that contains a code example,
such as

``` {r eval=FALSE}
# Access function documentation (not shown in vignette).
? numero.create
```

You can run all the code examples in one go by typing

``` {r eval=FALSE}
# Run all code examples (not shown in vignette).
fn <- system.file("extcode", "examples.R", package = "Numero")
source(fn)
```

Please use the tabs at the top to navigate to the next page.

## Dataset

In the examples, we use data on diabetic kidney disease from a previous
publication (see the readme file below). While simplified, the data contain
enough information to replicate some of the findings from the original study.
The main clinical characteristics are summarized in Table 1, and you can
print the readme file on the screen by typing

```{r eval=FALSE}
# Show readme file on screen (not shown in vignette).
fn <- system.file("extdata", "finndiane.readme.txt", package = "Numero")
cat(readChar(fn, 1e5))
```

<DIV class="caption"> Table: Summary of the FinnDiane dataset. The mean and
standard deviation are reported for each continuous variable. Abbreviations:
urinary albumin excretion rate (uALB), triglycerides (TG), high density
lipoprotein subclass 2 (HDL2). P-values were estimated by the t-test for
continuous variables and by Fisher’s test for binary traits.</DIV>

Trait                              | No kidney disease | Diabetic kidney disease | P-value
---------------------------------- | ----------------- | ----------------------- | --------
Men / Women                        |     192 / 196     |     119 / 106           | 0.45
Age (years)                        |     38.8 ± 12.2   |     41.7 ± 9.7          | 0.0012
Type 1 diabetes duration (years)   |     25.3 ± 10.3   |     28.6 ± 7.8          | <0.001
Log10 uALB (mg/24h)                 |     1.20 ± 0.51   |     2.72 ± 0.59         | <0.001
Log10 TG (mmol/L)                  |     0.034 ± 0.201 |     0.159 ± 0.212       | <0.001 
Total cholesterol (mmol/L)         |     4.89 ± 0.77   |     5.35 ± 0.96         | <0.001 
HDL2 cholesterol (mmol/L)          |     0.54 ± 0.16   |     0.51 ± 0.18         | 0.027
Log10 serum creatinine (µmol/L)    |     1.94 ± 0.09   |     2.14 ± 0.24         | <0.001 
Metabolic syndrome                 |     90 (23.2%)    |     114 (50.7%)         | <0.001 
Macrovascular disease              |     16 (4.1%)     |     38 (16.9%)          | <0.001 
Diabetic retinopathy               |     133 (34.4%)   |     178 (79.1%)         | <0.001
Died during follow-up              |     13 (3.4%)     |     43 (19.1%)          | <0.001

The data are included in the package as a tab-delimited file:

``` {r}
# Import data.
fname <- system.file("extdata", "finndiane.txt", package = "Numero")
dataset <- read.delim(file = fname, stringsAsFactors = FALSE)
nrow(dataset)
colnames(dataset)
```

Please use the tabs at the top to navigate to the next page.

## Data integrity

Most datasets contain missing entries, duplicated rows and other unusable
features, therefore it is important to carefully inspect all data points
and their identification. In this case, the biochemical and clinical
information about individual participants is organized as rows and each
row is identified by the variable INDEX. To create a data frame where the
row names are set to INDEX, type

``` {r}
# Manage unusable entries and data identification.
dataset <- numero.clean(data = dataset, identity = "INDEX")
```

Almost all functions within Numero work only on numeric values, so if the
dataset contains factors or other categorical data types, please convert them
to integer labels before analyses.


---
###########################################################################
---

# Basic pipeline {.tabset}

## Prepare

A typical usage case of the Numero framework involves a set of training
variables that are considered to explain medical or other outcomes. Here,
we use the biochemical data as the training variables under the general
hypothesis that the metabolic profile (as measured by biochemistry)
predicts adverse clinical events. The kidney disease dataset includes five
biochemical measures:

``` {r}
# Select training variables.
trvars <- c("CHOL", "HDL2C", "TG", "CREAT", "uALB")
```

Centering and scaling are basic machine learning techniues to ensure that
the variables with large values due to the choice of the measurement unit do
not bias the results. To create a standardized training set, type

``` {r results="hide"}
# Center and scale the training data.
trdat.basic <- scale.default(dataset[,trvars])
```

All the training variables have now unit variance and we can expect their
impact on the SOM to be dependent only on the information content and
statistical associations:
``` {r}
# Calculate standard deviations.
apply(trdat.basic, 2, sd, na.rm = TRUE)
```

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Create

Conceptually, the SOM algorithm is analogous to a human observer who wants to
make sense of a set of objects. Suppose there is a pile of portrait photos of
ordinary people on a round table, and your task is to organize them on the
table in such a way that similar photos are placed next to each other and, by
consequence, dissimilar photos end up on the opposite sides of the table.
During the organization process, your brain interprets the multi-variable
visual data (age, gender, ethnicity, hair style, head shape, etc.) and
projects it onto two dimensions (vertical and horizontal positions).

The SOM algorithm, also known as Kohonen map, mimics the organization process
and was originally inspired by the plasticity of neuronal networks. As this
is a self-organizing process, the starting point may have an impact on the
outcomes. We have added an extra initialization step to ensure repeatable
results. The function `numero.create` applies K-means clustering to create
an initial map model, and then refines the model using the SOM algorithm:

``` {r}
# Create a new self-organizing map based on the training set.
modl.basic <- numero.create(data = trdat.basic)
summary(modl.basic)
```

The resulting data structure contains the SOM itself and additional
information that is needed for further statistical analyses, please see
the manual pages for details.

The output from `numero.create` contains the positions of the data points
on a virtual round "table"; we use the term *map* instead of table, and refer
to the positions as the *layout* of a dataset on the SOM. The results also
contain an internal model of the training dataset; the model comprises a set
of *district centroids* that summarize the salient features of the training
dataset, please see the section on terminology for further descriptions.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Quality

After the training is complete, it is prudent to examine the SOM for
potential problems with the data. The Numero package provides three
different quality tools:

**Histogram:** The distribution of data points on the map can be expressed as
a histogram that shows the counts of data points in the districts.
If any areas of the SOM are devoid, or disproportionally populated by data
points, it usually reflects unusual or discontinuous value distributions
within the data (it can also be caused by poorly chosen preprocessing
methods). The histogram is rarely flat and up to two-fold differences
between the least and most populated districts are common, but if the pattern
is more extreme, careful examination of the distributions of the input data
in relation to the preprocessing procedures is warranted.

**Coverage:** The SOM is instrinsically robust against missing values if the
missingness pattern is close to random. We define data coverage as the
mean ratio of usable values per a multi-variable data point. If the coverage
drops substantially in certain districts compared to others, it is possible
that the data point layout has been biased by a non-random pattern of missing
values.

**Model fit:** The location of each data point on the map depends on how
similar it is to district centroids; the district with the most similar
centroid is chosen as the optimal location. The similarity of a centroid and
a data point can be quantified using mathematical formulas, such as the
Euclidean distance. The Numero package provides two measures. The first is
the residual distance $d$ between a data point and a centroid. The second
version is indpependent of scale and is calculated as
$z = (d - d_{mean}) / d_{sd}$, where $d_{mean}$ and $d_{sd}$ are the mean
residual and standard deviation of residuals in the training set.

``` {r}
# Calculate map quality measures for the training data.
qc.basic <- numero.quality(model = modl.basic)
```

Outliers can be detected by plotting the distribution of quality measures:

``` {r dev="svg", results="hide", fig.width=7, fig.height=3, fig.align="center", fig.cap="Figure: Distribution of model residuals."}
# Plot frequencies of data points at different quality levels.
par(mar = c(5,4,1,0), mfrow = c(1,2))
hist(x = qc.basic$layout$RESIDUAL, breaks = 50,
     main = NULL, xlab = "RESIDUAL", ylab = "Number of data points",
     col = "#FFEFA0", cex = 0.8)
hist(x = qc.basic$layout$RESIDUAL.z, breaks = 50,
     main = NULL, xlab = "RESIDUAL.z", ylab = "Number of data points",
     col = "#FFEFA0", cex = 0.8)
```

In this example, the distribution of residuals is heavily skewed and there
are a few data points that are several standard deviations from the mean
(RESIDUAL.z). The skewness is often due to skewed distributions in the data;
urinary albumin and serum creatinine are highly skewed towards larger values.
The Numero framework includes automatic procedures based on log-transforms
or rank-based preprocessing to mitigate skewness. We will explore these in
separate examples later in the document.

The quality measures can be visualized on the map to reveal subgroups of
problematic data points:

``` {r dev="svg", results="hide", fig.width=9, fig.height=3, fig.align="center", fig.cap="Figure: Visualization of quality measures across map districts. "}
# Plot map quality measures.
numero.plot(results = qc.basic, subplot = c(1,4))
```

While coverage is uniform across the map (there were only a few missing
data), the model fit is noticeably worse on the top part of the map
(higher RESIDUAL). The poor model fit is also reflected in the sample
histogram that shows a substantial depletion of data points in the top
part of the map. In practice, these structural issues do not necessarily
present a problem (as we see in the next section), but they need to be
considered when interpreting the significance of map patterns.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Statistics

Organizing a pile of photos on a table is achievable for a human observer
-- remember the discussion from the previous section -- but a human would
struggle with rows of numerical data, especially if there was a large number
of data points. Also, it is of interest to examine each variable one at a
time (e.g. if there is a map region associated with high cholesterol), or
to see if the outcomes of interest (e.g. mortality) accumulate on a
particular section of the SOM.

Since we consider the SOM as a map, we can apply concepts from geography to
overcome these challenges. Just as a city map can be colored according to
the average regional house prices, we can color the SOM according to the
average regional cholesterol concentrations or mortality rates.

Before the SOM is visualized, however, we need to estimate descriptive
statistics so that the colors can be calibrated according to the strength
of the statistical evidence.

Any finite dataset will produce regional variation on the SOM. For example,
it is reasonable to expect that parts of the map will differ with respect to
average mortality simply by chance. But how to know when the regional
variation is within the random expectation or not? The Numero package
contains a permutation engine that provides the answer:

``` {r}
# Map statistics for the whole dataset.
stats.basic <- numero.evaluate(model = qc.basic, data = dataset)
summary(stats.basic)
```
The output contains data frames for the data point layout, district averages
for each variables (i.e. component planes) and z-scores and p-values for
the observed amount of regional variation. Of note, the function omits
P-values for columns that overlap with the training data as it does not make
sense to estimate statistical significance for the same data that were used
for constructing the model.

Overall, the results above comprise the minimal information needed to color
the SOM according to district averages.

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Statistically normalized colorings of the training variables in the kidney disease dataset. The color intensity depends on how likely the observed regional variation would arise by chance; intense reds and intense blues indicate that these extremes would be very unlikely if the data point layout was random. The numbers show the average values in original units for selected districts."}
# Plot map colorings of training variables.
numero.plot(results = stats.basic, variables = trvars, subplot = c(2,3))
```

As expected, all the training variables (including CHOL) produce strong
colors since they were the basis of data point locations and therefore
directly influenced the data point layout. On the other hand, we did not use
age or sex as training variables, yet both show significant regional variation:

``` {r}
stats.basic$statistics[c("CHOL","MALE","AGE","T1D_DURAT"),
                       c("TRAINING","Z","P.z")]
```

These signals indicate that the biochemical profiles are associated
with age and sex, which is not a surprise considering the effects of ageing
and sex dimorphism in human populations. Of note, P-values are only reported
for the non-biochemical variables that were not used for training.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Subgroups

To recap, we hypothesized that the metabolic profile indicates and/or predicts
adverse health outcomes. For this reason, we chose to use only the biochemical
data to construct the SOM, which allows us to use the clinical traits as
validation outcomes from a modelling perspective. The next step is to divide
the dataset into a limited number of subgroups that summarize the typical
metabolic profiles we can expect too see in the data.

It is important to use only the training variables in this step so that we
remain agnostic to the clinical traits and preserve the original unsupervised
study design.

The Numero framework includes an interactive tool to conduct the
subgrouping based on SOM colorings:

``` {r eval=FALSE}
# Interactive subgrouping based on the training set.
subgr.basic <- numero.subgroup(results = stats.basic, variables = trvars)
```

In addition, the colorings can be saved as an interactive web page, which is
useful when the standard R plotting gets too slow or unusable due to multiple
colorings or large map size (see the chapter on large datasets and maps).

Unfortunately, the R vignette format is not compatible with the interactive
procedure, therefore we have to use a naive workaround (i.e. districts with
the highest and lowest urinary albumin) to demonstrate the features of the
pipeline.

**THE WORKAROUND IS ONLY FOR THE VIGNETTE AND WE ADVICE AGAINST USING IT IN
REAL STUDIES.**

``` {r echo=FALSE, results=FALSE}
# Workaround for the vignette document, DO NOT USE IN REAL STUDIES!
x <- (stats.basic$planes$DIAB_KIDNEY)*(stats.basic$planes$uALB)
tops <- which(x >= quantile(x, 0.75, na.rm=TRUE))
bottoms <- which(x <= quantile(x, 0.25, na.rm=TRUE))
workaround <- as.data.frame(stats.basic$map$topology)
workaround$REGION.label <- " "
workaround$REGION.color <- ""
workaround$REGION <- "not_selected"
workaround$REGION.label[tops] <- "A"
workaround$REGION.color[tops] <- "#FF000050"
workaround$REGION[tops] <- "HighALB"
workaround$REGION.label[bottoms] <- "B"
workaround$REGION.color[bottoms] <- "#89AB0070"
workaround$REGION[bottoms] <- "LowALB"
```

``` {r}
# Interactive selection not available in vignette.
if(!exists("subgr.basic")) subgr.basic <- workaround
```

``` {r origgrp1, dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Statistically normalized colorings of the training variables in the kidney disease dataset. The labels show the results from the subgrouping procedure."}
# Plot results from the subgrouping procedure.
numero.plot(results = stats.basic, variables = trvars,
            topology = subgr.basic, subplot = c(2,3))
```

It is evident from Figure 2 that the subgrouping workaround could be improved
to cover the most important patterns on the map. Serum creatinine and urinary
albumin are biomarkers of kidney disease, so they are highly concordant,
whereas cholesterol and triglycerides show patterns that do not fit well
with the kidney biomarkers. Therefore, additional subgroups and manual
adjustments of the district labels is the preferred course of action.

With the interactive tool, users can choose their own subgroups based on
the SOM patterns. This is the greatest strength of the pipeline: multiple
people will be able to participate in the process of defining the subgroups,
and different subgroups can be created for different settings. For example,
cardiovascular disease is the most common cause of death in people with type
1 diabetes, and therefore a subgrouping that takes cholesterol into account
may be warranted.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Compare

We omitted the clinical end-points intentionally from the previous analyses to
minimize their influence on the subgrouping step, as required by the
study design. Since we have now defined the subgroups, it is safe to plot the
colorings for all the variables. The subplot option was used here to ensure a
pleasing look for the vignette, but it can be omitted in most circumstances.

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Statistically normalized colorings of selected variables in the kidney disease dataset."}
# Plot results from subgrouping procedure for non-biochemical variables.
numero.plot(results = stats.basic,
            variables = c("AGE",
                          "MALE",
                          "DIAB_KIDNEY",
                          "DIAB_RETINO",
                          "MACROVASC",
                          "DECEASED"),
            topology = stats.basic$map$topology, subplot = c(2,3))
```

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of selected variables with subgroup labels."}
# Plot results from subgrouping procedure for non-biochemical variables.
numero.plot(results = stats.basic,
            variables = c("AGE",
                          "MALE",
                          "DIAB_KIDNEY",
                          "DIAB_RETINO",
                          "MACROVASC",
                          "DECEASED"),
            topology = subgr.basic, subplot = c(2,3))
```

The subgroup with the highest urinary albumin (A: HighALB) overlaps with the
clinical classification of diabetic kidney disease (top-right coloring).
This is expected, since albuminuria is one of the key diagnostic markers
for kidney disease. Accordingly, the same subgroup overlaps with high
prevalence of retinopathy (bottom-left) and high mortality rate
(bottom-right).

In research papers, it is often easier to report conventional comparison
statistics between the subgroups you have defined rather than the SOM
visualizations (however, the colorings should always be included as
supplements). The Numero package includes the function `numero.summary()`
that combines the visual subgroupings with the original data:

``` {r}
# Compare subgroups.
report.basic <- numero.summary(results = stats.basic, topology = subgr.basic)
colnames(report.basic)
```

For instance, the report contains the statistical signals for mortality
across the subgroups:

``` {r}
# Show results for mortality rate.
rows <- which(report.basic$VARIABLE == "DECEASED")
report.basic[rows,c("SUBGROUP","N","MEAN","P.chisq")]
```

The first row of the output shows the chi2 test for the full contingency
table. The largest subgroup is always chosen as the reference, hence
P = 1 for the individuals who were not selected to either of the albumin
subgroups. Although the LowALB subgroup (B) shows a lower mortality rate
compared to the reference, the difference is not as strong statistically.
Instead, it seems that the HighALB group is responsible for most of the
signal.

Lastly, it is often useful to identify the members of these subgroups.
The output contains that information as attributes, please see the
documentation on `numero.summary()` and `nroSummary()` for details.

---
###########################################################################
---

# Preprocessing and adjustment {.tabset}

## Prepare

The introductory analyses revealed that both sex and age were associated with
the metabolic profiles via the SOM. Moreover, men and women display anatomical
and metabolic differences, which usually complicate the interpretation of the
SOM. For this reason, we recommend using a sex-specific standardization
procedure that eliminates the differences. If necessary, separate
visualizations can be made afterwards for men and women using the same map.

Ageing is another factor that is often considered a confounder, especially
within cross-sectional study designs. For this reason, we will reduce
the co-variation with AGE and T1D_DURAT in the training set prior to
creating the SOM.

To prepare the training data, we can use the tool provided with the
Numero package:

``` {r}
# Mitigate stratification and confounding factors.
trdata.adj <- numero.prepare(data = dataset, variables = trvars, batch = "MALE",
                             confounders = c("AGE", "T1D_DURAT"))
colnames(trdata.adj)
```

In addition to adjusting for the confounders, the function also applies
the logarithm to variables that are highly skewed. While binary data can
be used for training, we recommend using continuous variables only for
building a SOM since they are easier to adjust for confounders and less
prone to numerical instability.

We also need to know which identities belong to men or women so that
the sexes can be separated in later analysis steps. As a by-product of
the batch correction for MALE, the identities of batch members are included
in the output according to the values in the MALE column
(1 = man, 0 = woman):

``` {r}
subsets <- attr(trdata.adj, "subsets")
women <- subsets[["0"]]
men <- subsets[["1"]]
c(length(men), length(women))
```

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Create

The SOM is created the same way as within the basic pipeline, but this time
using the sex-adjusted training variables:

``` {r}
# Create a new self-organizing map based on sex-adjusted data.
modl.adj <- numero.create(data = trdata.adj)
summary(modl.adj)
```

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Quality

After the training is complete, the structure of the SOM can be inspected
as we did in the previous example:

``` {r}
# Calculate map quality measures for sex-adjusted data.
qc.adj <- numero.quality(model = modl.adj)
```

``` {r dev="svg", results="hide", fig.width=7, fig.height=3, fig.align="center", fig.cap="Figure: Distribution of model residuals. The training data were adjusted for age and sex."}
# Plot frequencies of data points at different quality levels.
par(mar = c(5,4,1,0), mfrow = c(1,2))
hist(x = qc.adj$layout$RESIDUAL, breaks = 20,
     main = NULL, xlab = "RESIDUAL", ylab = "Number of data points",
     col = "#FFEFA0", cex = 0.8)
hist(x = qc.adj$layout$RESIDUAL.z, breaks = 20,
     main = NULL, xlab = "RESIDUAL.z", ylab = "Number of data points",
     col = "#FFEFA0", cex = 0.8)
```

The preprocessing resulted in a substantial improvement of the residual
distributions. Although the residuals are still skewed, the most extreme
residuals have been removed:

``` {r}
# Maximum residuals from unadjusted and adjusted analyses.
c(max(qc.basic$layout$RESIDUAL.z, na.rm=TRUE),
  max(qc.adj$layout$RESIDUAL.z, na.rm=TRUE))
```

Furthermore, the spread of data points across the map is now more
balanced

``` {r}
# Variation in point density in unadjusted and adjusted analyses.
c(sd(qc.basic$planes$HISTOGRAM, na.rm=TRUE),
  sd(qc.adj$planes$HISTOGRAM, na.rm=TRUE))
```

and the depletion of data points in the poor quality area is attenuated:

``` {r dev="svg", results="hide", fig.width=9, fig.height=3, fig.align="center", fig.cap="Figure: Visualization of data point quality measures across map districts from age and sex adjusted analysis."}
# Plot map quality measures.
numero.plot(results = qc.adj, subplot = c(1,4))
```

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Statistics

We can estimate the map statistics using the whole dataset as before, except
that the model is now age and sex adjusted. Note how we use the original
unadjusted data as input `data = dataset` but the evaluation
is done via the map layout created from adjusted data:

``` {r}
# Map statistics for the whole dataset.
stats.adj <- numero.evaluate(model = qc.adj, data = dataset)
```

Furthermore, it is also possible to conduct sex-specific analyses after we
identified men and women in the dataset:

``` {r results="hide"}
# Map statistics for women.
stats.adjW <- numero.evaluate(model = qc.adj, data = dataset[women,])
```

``` {r results="hide"}
# Map statistics for men.
stats.adjM <- numero.evaluate(model = qc.adj, data = dataset[men,])
```

We could have used `data = trdata.adj` instead of `data = dataset`,
however, using the original unadjusted data to describe the sex-adjusted
subgroups is easier to interpret in relation to what you might see when
meeting the participants in the real world.

If the adjustments of the training data are successful, we expect
to see non-significant P-values for MALE and AGE when using the original
data:

``` {r}
stats.adj$statistics[c("MALE","AGE","T1D_DURAT"), c("Z","P.z")]
```

Indeed, the regional variation in the proportion of men to women is now
within the range of random sampling fluctuation. That said, the sex
adjustment only affects the individuals' positions on the map, thus the
well known sex differences in the average biomarker concentrations are
not affected:

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of training variables using the data from female participants."}
numero.plot(results = stats.adjW, variables = trvars, subplot = c(2,3))
```

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of training variables using the data from male participants."}
numero.plot(results = stats.adjM, variables = trvars, subplot = c(2,3))
```

The colorings look very similar, but the regional averages are systematically
different between the sexes. For instance, HDL2C is higher in women across
all areas of the map, while the color patterns are almost indistinguishable
between the sexes.

As we estimated the statistics separately for men and women, the color
scales are likely to differ, which may make it more difficult to compare
the data patterns. Another option is to calibrate the colors between the
sexes against a common reference, that is, the same shade of blue will
always correspond to the same concentration of cholesterol, for example,
regardless of which subset of individuals are being visualized.
This is accomplished by using the results of the full dataset as the
reference point `reference = stats.adj`:

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of training variables using the data from female participants. Color scales were derived from the full dataset."}
numero.plot(results = stats.adjW, variables = trvars,
            subplot = c(2,3), reference = stats.adj)
```

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of training variables using the data from male participants. Color scales were derived from the full dataset."}
numero.plot(results = stats.adjM, variables = trvars,
            subplot = c(2,3), reference = stats.adj)
```

It is now easier to see from the colors that women tend to have higher
HDLC2 (bright red areas on the right side and less blue on the left side),
whereas they tend to have lower serum creatinine (more blue for women
compared to men).

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Subgroups

To recap, we followed the procedures in the basic pipeline, except that the
training data were adjusted for age and sex differences, and we also analyzed
men and women separately. The subgrouping step proceeds as before.

The Numero framework includes an interactive tool to conduct the
subgrouping based on SOM colorings

``` {r eval=FALSE}
# Interactive subgrouping based on the training set.
subgr.adj <- numero.subgroup(results = stats.adj, variables = trvars)
```

As mentioned previously, the vignette format is not compatible with the
interactive tool, therefore we have to use a naive workaround (i.e.
districts with the highest and lowest urinary albumin) to demonstrate the
features of the pipeline.

**THE WORKAROUND IS ONLY FOR THE VIGNETTE AND WE ADVICE AGAINST USING IT IN
REAL STUDIES.**

``` {r echo=FALSE, results=FALSE}
# Workaround for the vignette document, DO NOT USE IN REAL STUDIES!
x <- (stats.adj$planes$DIAB_KIDNEY)*(stats.adj$planes$uALB)
tops <- which(x >= quantile(x, 0.75, na.rm=TRUE))
bottoms <- which(x <= quantile(x, 0.25, na.rm=TRUE))
workaround <- as.data.frame(stats.adj$map$topology)
workaround$REGION.label <- " "
workaround$REGION.color <- ""
workaround$REGION <- "not_selected"
workaround$REGION.label[tops] <- "A"
workaround$REGION.color[tops] <- "#FF000050"
workaround$REGION[tops] <- "HighALB"
workaround$REGION.label[bottoms] <- "B"
workaround$REGION.color[bottoms] <- "#89AB0070"
workaround$REGION[bottoms] <- "LowALB"
```

``` {r}
# Interactive selection not available in vignette.
if(!exists("subgr.adj")) subgr.adj <- workaround
```

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Statistically normalized colorings of the training variables in the kidney disease dataset. The labels show the results from the subgrouping procedure. The map was created from sex-adjusted data."}
# Plot results from subgrouping procedure.
numero.plot(results = stats.adj, variables = trvars,
            topology = subgr.adj, subplot = c(2,3))
```

Although the sex adjustment had a minor effect on the map colorings, we can
observe some changes in the separation of low and high urinary albumin.
As was discussed earlier, the subgroups could be made more reasonable, which
highlights the importance of being able to set the subgroup labels yourself
for the most useful or parsimonious results.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Compare

In the last step, the subgroups are compared using conventional statistics
as we did earlier. It is always useful to plot the subgroup labels in
relation to the variables that were not used for training:


``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of selected variables with subgroup labels. The map was created from age and sex adjusted data."}
# Plot results from subgrouping procedure for non-biochemical variables.
numero.plot(results = stats.adj,
            variables = c("AGE",
                          "MALE",
                          "DIAB_KIDNEY",
                          "DIAB_RETINO",
                          "MACROVASC",
                          "DECEASED"),
            topology = subgr.adj, subplot = c(2,3))
```

The effect of the sex adjustment on the variable MALE is immediately
apparent: the regional variation in sex ratio is so close to the random
expectation that the color scale has been squashed flat. In the unadjusted
model, the difference between the predominantly male and the
predominantly female districts was substantial, whereas
here the gap has been reduced:

``` {r}
# District averages from unadjusted analysis.
summary(stats.basic$planes$MALE)
```

``` {r}
# District averages from adjusted analysis.
summary(stats.adj$planes$MALE)
```

The results for mortality have changed as a consequence of the changes in
subgrouping, but not substantially:

``` {r}
# Compare subgroups.
report.adj <- numero.summary(results = stats.adj, topology = subgr.adj)
```

``` {r}
# Show results for mortality rate.
rows <- which(report.adj$VARIABLE == "DECEASED")
report.adj[rows,c("SUBGROUP","N","MEAN","P.chisq")]
```

Again, the statistical signal is mostly coming from the high albuminuria
subgroup. It is important to note that while the preprocessing affected
the map quality and the results in general, the overall interpretations
of the results did not change dramatically.


---
###########################################################################
---

# Multiple datasets {.tabset}

## Prepare

Men and women are biologically different, and also have different
biochemical profiles as we observed earlier. In this example, we
use the sexes as if they are two different epidemiological cohorts.
Indeed, differences in data collection between epidemiological studies
lead to differences in the mean concentrations of biomarkers. In this
respect, the male-female comparison is an illustrative and easy
proxy for discussing multi-cohort applications of the SOM.

A typical multi-cohort setting involves a discovery cohort and one or
more replication cohorts with additional meta-analyses across all
available data. We will focus on two aspects: 1) defining subgroups
in the discovery cohort and 2) assigning participants from the replication
cohort into these subgroups. For the meta-analysis part, established
statistical methods can be applied after the subgroups have been defined
(not discussed here).

We also use this example as an opportunity to show how rank-based
preprocessing of the data affects the quality measures by setting
`method = "tapered"`. The SOM framework includes two rank-based methods:
"uniform" for normalized ranks between -1 and +1 and "tapered" for a version
that puts more samples around zero and less around -1 or +1.

To facilitate more intuitive discussion, we refer to the sexes as two
different cohorts from now on. To simulate incompatible datasets, we
add a dataset of men with the metabolic syndrome as a third cohort:

``` {r}
ds.discov <- dataset[women,]
ds.replic <- dataset[men,]
ds.mets <- ds.replic[which(ds.replic$METAB_SYNDR == 1),]
```

We create four preprocessed datasets (one for discovery/training and three
for replication) to demonstrate how data can be used across cohorts. Summary
statistics are shown for urinary albumin as it is the most important
clinical biomarker:

``` {r results="hide"}
# Discovery cohort and training set.
trdata.discov <- numero.prepare(data = ds.discov, variables = trvars,
                                confounders = c("AGE", "T1D_DURAT"),
                                method = "tapered")
```

``` {r}
summary(trdata.discov$uALB)
```

``` {r results="hide"}
# Replication cohort, version A.
param <- attr(trdata.discov, "pipeline")
trdata.replicA <- numero.prepare(data = ds.replic, pipeline = param)
```

``` {r}
summary(trdata.replicA$uALB)
```

The first replication dataset A is analogous to a situation where two
cohorts/batches are compared without batch correction, that is, we assume
that there is no measurement or other bias between the descriptive cohort
statistics. This is the safe starting point, since the quality control
will reveal if the assumption is violated, and we can adjust the second
round of analyses accordingly.


``` {r results="hide"}
# Replication cohort, version B.
trdata.replicB <- numero.prepare(data = ds.replic, variables = trvars,
                             confounders = c("AGE", "T1D_DURAT"),
                             method = "tapered")
```

``` {r}
summary(trdata.replicB$uALB)
```

The second dataset B is adjusted for discrepancies in average values.
We still assume that both cohorts represent the same underlying population
that is sampled the same way. This situation is similar to having two
separate birth cohorts of approximately the same age from the same broader
population, but with biochemical data measured in different labs.

``` {r results="hide"}
# Replication cohort, MetS version.
trdata.mets <- numero.prepare(data = ds.mets, variables = trvars,
                            confounders = c("AGE", "T1D_DURAT"),
                            method = "tapered")
```

``` {r}
summary(trdata.mets$uALB)
```

The third dataset reflects a situation where the study cohort have been
recruited differently or are otherwise demographically incompatible.
This represents a challenge to the subgrouping analysis and can lead to
erroneous conclusions. Bias of this type is difficult to distinguish from
bias from measurement techniques, and therefore needs to be accounted for
by re-sampling or calibration techniques before using the SOM framework.

Please use the tabs at the top to navigate to the next page.

## Create

The SOM is created the usual way and trained with the discovery dataset,
exept that we set the map radius to be the same as before to make it
easier to compare results across examples:

``` {r}
# Create a new self-organizing map based on sex-adjusted data.
radius.basic <- attr(modl.basic$map$topology, "radius")
modl.discov <- numero.create(data = trdata.discov, radius = radius.basic)
summary(modl.discov)
```

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Quality

In this section, we inspect the structure of the SOM for the discovery cohort
and the versions of the replication cohort. Notice how the model is the same,
but new data are used for the replication cohorts:

``` {r results="hide"}
# Calculate map quality measures.
qc.discov <- numero.quality(model = modl.discov)
qc.replicA <- numero.quality(model = modl.discov, data = trdata.replicA)
qc.replicB <- numero.quality(model = modl.discov, data = trdata.replicB)
qc.mets <- numero.quality(model = modl.discov, data = trdata.mets)
```

The training data were transformed into ranks during the preprocessing step.
In many cases, this is the safest approach since it mitigates the effect of
outliers or highly skewed distributions. On the other hand, the rank transform
discards the original differences in magnitudes between data points, which
may be undesirable in some circumstances.

In this example, the rank preprocessing resulted in less skewed residuals
compared to the earlier examples:

``` {r}
# Define comparable histogram bins.
rz <- c(qc.adj$layout$RESIDUAL.z, qc.discov$layout$RESIDUAL.z)
rz.breaks <- seq(min(rz, na.rm=TRUE), max(rz, na.rm=TRUE), length.out=20)
``` 

``` {r dev="svg", results="hide", fig.width=7, fig.height=3, fig.align="center", fig.cap="Figure: Distribution of model residuals when the training data were preprocessed by scaling & centering or by tapered ranking."}
# Plot frequencies of data points at different quality levels.
par(mar = c(5,4,1,0), mfrow = c(1,2))
hist(x = qc.adj$layout$RESIDUAL.z, breaks = rz.breaks,
     main = NULL, xlab = "RESIDUAL.z (scale & center)",
     ylab = "Number of data points", col = "#FFEFA0", cex = 0.8)
hist(x = qc.discov$layout$RESIDUAL.z, breaks = rz.breaks,
     main = NULL, xlab = "RESIDUAL.z (rank)",
     ylab = "Number of data points", col = "#FFEFA0", cex = 0.8)
```

Residuals for the replication datasets are expected to be higher than those
for the training set (the training set always overfits to itself).
Nevertheless, all the rank-transformed replication sets had lower maximum
residuals than the previous centered and scaled training set:

``` {r}
# Comparison of maximum residuals across examples.
r <- c(max(qc.basic$layout$RESIDUAL.z, na.rm=TRUE),
       max(qc.adj$layout$RESIDUAL.z, na.rm=TRUE),
       max(qc.discov$layout$RESIDUAL.z, na.rm=TRUE),
       max(qc.replicA$layout$RESIDUAL.z, na.rm=TRUE),
       max(qc.replicB$layout$RESIDUAL.z, na.rm=TRUE),
       max(qc.mets$layout$RESIDUAL.z, na.rm=TRUE))
names(r) <- c("basic", "adj", "discov", "replicA", "replicB", "mets")
print(r)
```

More balanced data point density is another potential benefit of the
rank-transform:

``` {r dev="svg", results="hide", fig.width=10, fig.height=3, fig.align="center", fig.cap="Figure: Visualization of quality measures across map districts for the discovery cohort. "}
# Plot map quality measures.
numero.plot(results = qc.discov, subplot = c(1,4))
```

On the other hand, the centered and scale analysis in the previous chapter
revealed a connection between weak model fit and high burden of disease;
this aspect of the data was lost due to the rank transform. The concentric
pattern of good fit at the center and weaker fit on the edges (see RESIDUAL)
is often seen when using rank-based preprocessing, which suggests that the
pattern should not be interpreted as being a characteristic of any
particular dataset.

Comparisons between the replication sets revealed important clues on how
well the data fit together:

``` {r dev="svg", results="hide", fig.width=10, fig.height=3, fig.align="center", fig.cap="Figure: Visualization of quality measures across map districts for the replication dataset A. "}
# Plot map quality measures.
numero.plot(results = qc.replicA, subplot = c(1,4))
```

For replication set A, the data point histogram is very uneven with the
majority of data points clustered in one quadrant of the map. To recap,
replication A was preprocessed using the same parameters as the discovery
set under the assumption that the two datasets were samples from the same
underlying population with identical measurement protocols. Evidently, this
assumption is wrong since the multi-variable data did not overlap well.

``` {r dev="svg", results="hide", fig.width=10, fig.height=3, fig.align="center", fig.cap="Figure: Visualization of quality measures across map districts for the replication dataset B. "}
# Plot map quality measures.
numero.plot(results = qc.replicB, subplot = c(1,4))
```

Replication B was assumed to come from the same population, but with
potential bias in the measured values. As the training data were preprocessed
independently of the discovery set, resulting in identical standardized
distributions, the data point histogram and residuals are similar to the
discovery set. This does not mean that the initial assumption is correct,
but it is possible to investigate demographic factors such as age, sex
or clinical disease profiles to check that the overall characteristics are
compatible across the different regions on the map.

``` {r dev="svg", results="hide", fig.width=10, fig.height=3, fig.align="center", fig.cap="Figure: Visualization of quality measures across map districts for the replication dataset with the metabolic syndrome. "}
# Plot map quality measures.
numero.plot(results = qc.mets, subplot = c(1,4))
```

Lastly, the metabolic syndrome subset was created as an example with biased
sampling and where the replication is not compatible with the discovery
dataset. However, this is hard to detect from the quality plots. As mentioned
above, the independent preprocessing will reproduce the even data point
histogram and residuals (in most settings), therefore the quality measures
alone cannot identify problems arising from these types of problems in
study design.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Statistics

To create colorings of the SOM for the different datasets, we combine
several features of the Numero pipeline. The first task is to identify
where the data points go on the map. For the discovery cohort (that
includes the training set), this information is contained in the model
itself. However, for the other cohorts, the district assignments
are obtained as part of the quality control; the output from
`numero.quality()` contains the layout of data points that were used
for calculating the histograms and residuals in the previous section.

The estimation of map statistics proceeds as before, except that the
quality control results are used for the replication datasets:

``` {r results="hide"}
# Map statistics for discovery and replication datasets.
stats.discov <- numero.evaluate(model = modl.discov, data = ds.discov)
stats.replicA <- numero.evaluate(model = qc.replicA, data = ds.replic)
stats.replicB <- numero.evaluate(model = qc.replicB, data = ds.replic)
stats.mets <- numero.evaluate(model = qc.mets, data = ds.mets)
```

The second task is to choose the coloring scheme. In the previous chapter,
we experimented with sex-specific colorings and here we will apply the same
technique of a common reference `reference = stats.discov`. This will
ensure that colors between the datasets are directly comparable, that is,
the same data value will always correspond to the same color.

In addition, we reduce the "gain" or global amplification of the color range.
This has the effect of making the available value-color mapping wider so that
values that exceed the reference range wash out less:

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of training variables for the discovery dataset."}
numero.plot(results = stats.discov, variables = trvars,
            gain = 0.8, subplot = c(2,3))
```

Notice how the strongest reds and blues are not used for the discovery set
due to the reduced (<1.0) gain. This means that they are now available to use
to distinguish more extreme values on the map.

The results for replication A are very similar to the discovery set both
for the color patterns and also the district averages. Yet we know from the
quality analysis that the two datasets did not overlap well on the SOM and
the assumptions on study design were invalid -- this is one of the reasons
it is important to check the quality of the SOM.

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of selected variables for replication A."}
numero.plot(results = stats.replicA, variables = trvars,
            gain = 0.8, subplot = c(2,3), reference = stats.discov)
```

Replication B reproduces the patterns of the discovery set in relative
terms (areas of lower and higher concentrations match), but the absolute
values are different. This is expected since we assumed that the data
patterns would be the same but there was systematic bias in the measurement
data.

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of selected variables for replication B."}
numero.plot(results = stats.replicB, variables = trvars,
            gain = 0.8, subplot = c(2,3), reference = stats.discov)
```

The metabolic syndrome dataset produced colorings that were substantially
different from the discovery set:

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Colorings of selected variables for MetS dataset."}
numero.plot(results = stats.mets, variables = trvars,
            gain = 0.8, subplot = c(2,3), reference = stats.discov)
```

While the relative patterns appear somewhat similar, the district averages
reveal that the SOM of the discovery set does not replicate well. The main
observation that demonstrates the poor fit is biological. For example, the
lowest serum creatinine district average in the discovery set is within the
range that you would expect to see in healthy adults (60 - 110 umol/L),
however, the district averages for the MetS dataset are higher in those map
areas, which may indicate a substantial difference in kidney health.

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Subgroups

In previous chapters, the subrgouping was done based on the training
variables only so that the other traits could be used as methodologically
independent data. Here, the discovery set enables us to define the subgroups
based on all variables (if we so choose) since we did not use the replication
datasets for training. For this example, a good option is to use clinically
interesting variables for the subgrouping:

``` {r}
# Selection of clinically interesting variables.
clinvars <- c("uALB", "AGE", "DIAB_KIDNEY", "DIAB_RETINO",
              "MACROVASC", "DECEASED")
```

As mentioned earlier, the Numero framework includes an interactive tool to
conduct the subgrouping based on SOM colorings

``` {r eval=FALSE}
# Interactive subgrouping based on the training set.
subgr.discov <- numero.subgroup(results = stats.discov, variables = clinvars)
```

However, the vignette format is not compatible with the interactive tool,
therefore we have to use a naive workaround (i.e. districts with the highest
and lowest prevalence of kidney disease) to demonstrate the features of the
pipeline.

**THE WORKAROUND IS ONLY FOR THE VIGNETTE AND WE ADVICE AGAINST USING IT IN
REAL STUDIES.**

``` {r echo=FALSE, results=FALSE}
# Workaround for the vignette document, DO NOT USE IN REAL STUDIES!
x <- (stats.discov$planes$DIAB_KIDNEY)*(stats.discov$planes$uALB)
tops <- which(x >= quantile(x, 0.75, na.rm=TRUE))
bottoms <- which(x <= quantile(x, 0.25, na.rm=TRUE))
workaround <- as.data.frame(stats.discov$map$topology)
workaround$REGION.label <- " "
workaround$REGION.color <- ""
workaround$REGION <- "not_selected"
workaround$REGION.label[tops] <- "A"
workaround$REGION.color[tops] <- "#FF000050"
workaround$REGION[tops] <- "HighDiabKD"
workaround$REGION.label[bottoms] <- "B"
workaround$REGION.color[bottoms] <- "#89AB0070"
workaround$REGION[bottoms] <- "LowDiabKD"
```

``` {r}
# Interactive selection not available in vignette.
if(!exists("subgr.discov")) subgr.discov <- workaround
```

``` {r dev="svg", results="hide", fig.width=6, fig.height=4, fig.align="center", fig.cap="Figure: Statistically normalized colorings of the training variables in the kidney disease dataset. The labels show the results from the subgrouping procedure."}
# Plot results from subgrouping procedure.
numero.plot(results = stats.discov, variables = clinvars,
            topology = subgr.discov, subplot = c(2,3))
```

Please use the tabs at the top to navigate to the next page.

---
#---------------------------------------------------------------------------
---

## Compare

In the last step, the subgroups are compared using conventional statistics.
We use the dataset-specific map statistics to characterize the
subgroups in the discovery and replication sets:

``` {r eval=FALSE}
# Compare subgroups.
report.discov <- numero.summary(results = stats.discov,
                                topology = subgr.discov)
report.replicA <- numero.summary(results = stats.replicA,
                                 topology = subgr.discov)
report.replicB <- numero.summary(results = stats.replicB,
                                 topology = subgr.discov)
report.mets <- numero.summary(results = stats.mets,
                              topology = subgr.discov)
```

``` {r echo=FALSE, results="hide"}
suppressWarnings({
report.discov <- numero.summary(results = stats.discov,
                                topology = subgr.discov)
report.replicA <- numero.summary(results = stats.replicA,
                                 topology = subgr.discov)
report.replicB <- numero.summary(results = stats.replicB,
                                 topology = subgr.discov)
report.mets <- numero.summary(results = stats.mets,
                              topology = subgr.discov)})
```

The mortality rates in the discovery set are as expected with the
highest mortality observed for the districts that had the highest
kidney disease prevalence:

``` {r}
# Show results for mortality rate in the discovery set.
rows <- which(report.discov$VARIABLE == "DECEASED")
report.discov[rows,c("SUBGROUP","N","MEAN","P.chisq")]
```

The same pattern is observable also for replication A, however, the
subgroup sizes are poorly balanced and the statistics weaker as
a consequence:


``` {r}
# Show results for mortality rate in replication A.
rows <- which(report.replicA$VARIABLE == "DECEASED")
report.replicA[rows,c("SUBGROUP","N","MEAN","P.chisq")]
```

The outcomes in the discovery set are replicated well in dataset B
and the subgroups are also well balanced in their size:

``` {r}
# Show results for mortality rate in replication B.
rows <- which(report.replicB$VARIABLE == "DECEASED")
report.replicB[rows,c("SUBGROUP","N","MEAN","P.chisq")]
```

The mortality outcomes for the metabolic syndrome dataset are
not statistically significant. This is partly explained by the smaller
size of the dataset (= less power) but the overall mortality risk was
higher already due to the biased sampling, and therefore direct comparisons
with the discovery set are problematic.

``` {r}
# Show results for mortality rate in metabolic syndrome subset.
rows <- which(report.mets$VARIABLE == "DECEASED")
report.mets[rows,c("SUBGROUP","N","MEAN","P.chisq")]
```

In all studies, we recommend careful verification of the design and the
validity of assumptions using traditional statistical tests, SOM quality
measures, SOM colorings and also the final subgroups comparisons.

---
###########################################################################
---

# Large datasets and maps

The Numero package is also suitable for analyzing larger datasets (>100k
rows of data). Given the increase in statistical power, the map radius
can be safely increased, but this usually causes difficulties in plotting
the colorings on screen due to the inefficient way the default R graphics
are implemented. In particular, it may not be feasible to do the
interactive part using the standard R plotting.

It is possible to use web browsers for a better user experience but, to view
the map colorings, the plots need to be saved in an HTML file first. Here,
we are using the results from the basic example that was introduced in the
beginning of the document:

``` {r eval=FALSE}
# Save map colorings of training variables.
numero.plot(results = stats.basic, variables = trvars,
            folder = "/tmp/Results")
```

``` {r echo=FALSE}
s <- paste("\n", 
    "*** numero.plot ***\n", 
    "Thu Sep  5 15:44:02 2019\n", 
    "\n", 
    "Resources:\n", 
    "5 column(s) included\n", 
    "destination folder '/tmp/Results'\n", 
    "\n", 
    "Figure 1:\n", 
    "5 subplot(s)\n", 
    "file name '/tmp/Results/figure01.svg'\n", 
    "97622 bytes saved in '/tmp/Results/figure01.svg'\n", 
    "99982 bytes saved in '/tmp/Results/figure01.html'\n", 
    "\n", 
    "Summary:\n", 
    "1 figure(s) -> '/tmp/Results'\n", sep="")
cat(s)
```

You can replace the folder to something suitable in your computer. Here we
used the temporary folder of a Linux system as an example. The command
creates both a Scalable Vector Graphics (SVG) file that you can open in
editors such as Inkscape, and a web page (HTML) that you can open in
a browser. Importantly, the web page is interactive, and you can select
regions of the SOM and click on empty space between plots to save
the topology and regions in a tab-delimited text file.

Assuming we assigned districts to different regions and downloaded the
results in the file 'Downloads/regions.txt', we can then import the
spreadsheet into R by typing

``` {r eval=FALSE}
# Import topology and region assignments.
subgr <- read.delim(file = "Downloads/regions.txt", stringsAsFactors=FALSE)
```

and calculate subgroup statistics with


``` {r eval=FALSE}
# Calculate subgroup statistics.
report <- numero.summary(results = stats.basic, topology = subgr)
```

as before.

**Useful hints**: If it takes a long time to train a SOM for a large dataset,
try setting the `subsample` parameter to something less than the number of
rows in the training set. It may also take a long time to evaluate the map
statistics; to reduce waiting time, you can assign data columns into smaller
batches and then run each batch on a different processor in parallel.

---
###########################################################################
---

# Collinear inputs

High-dimensional datasets are usually characterized by strong correlations
between groups of variables. The correlations can arise from multiple factors
and the SOM is in general a safe approach to use even if some unwanted
collinearity is present. In an extreme case, data collection can be biased
to an extent that the SOM patterns start to reflect the limitations of
the collection protocol rather than the phenomena under study.

Consider the following example: a mass spectrometry platform is optimized to
capture a large number of lipid molecules from a biological sample, but
only a few non-lipid metabolites are available. If the goal of the SOM is
to model the biochemical profiles, then using such data as such may end up
reducing the SOM to a model that is too heavily influenced by the amount of
lipid in the sample at the expence of other molecular information.

The Numero offers a network-based tool to merge collinear variables in a
dataset. This is beneficial particularly when training the SOM as it allows
for a more balanced model fit, but it is also of general utility as a means
to reduce the dimensionality of a given dataset in a way that is also easy
to interpret.

To demonstrate how the procedure works, we start with a random dataset:

``` {r eval=TRUE}
# Random data matrix.
x <- matrix(rnorm(100000), ncol=100)
```

Next, correlation structure is introduced by adding values together across
columns:

``` {r eval=TRUE}
# Create correlation modules.
x[,12:20] <- (x[,12:20] + x[,11])
x[,32:40] <- (x[,32:40] + x[,31])
x[,62:90] <- (x[,62:90] + x[,61])
x[,50] <- (x[,20] + x[,90]) # connecting node
```

Note the addition of the connecting node. This will produce an edge between
two network modules to illustrate a scenario where correlation structure is
too complex to be divided into fully separate modules.

To keept track of things, the columns are named in a way that reflects the
modular correlations. This helps us to verify that the results are correct.
``` {r eval=TRUE}
# Set column names.
cnames <- paste0("X", 1:ncol(x))
cnames[11:20] <- paste0("M1.", cnames[11:20])
cnames[31:40] <- paste0("M2.", cnames[31:40])
cnames[61:90] <- paste0("M3.", cnames[61:90])
colnames(x) <- cnames
```
To merge redundant columns into a single column of aggregated score, we
invoke the appropriate option in `numero.prepare()`:
``` {r eval=TRUE}
# Merge collinear modules.
y <- numero.prepare(x, coalesce=TRUE, method="")
```

For a more granular control of the algorithm, refer to the help page of
`nroCoalesce()`.

We can now compare the columns before and after merging:
``` {r eval=TRUE}
# Columns that were removed or added.
print(setdiff(colnames(x), colnames(y)))
print(setdiff(colnames(y), colnames(x)))
```

The output from `numero.prepare()` also contains detailed information
of the modules:
``` {r eval=TRUE}
# Show module members.
modules <- attr(y,"pipeline")$modules
mnames <- names(modules)
lapply(mnames, function(k, x) {
    cat("\n", k, "\n", sep="")
    print(x[[k]]$weights)
}, x=modules)
```

Importantly, the same pipeline can be applied to new datasets (via
`numero.prepare()`) to allow the same preprocessing procedures to be used
across multiple datasets.

---
###########################################################################
---

# Terminology

**Best-matching centroid**: The BMC is the district centroid that is the most
similar to a data point. The concept is closely related to the data point
layout: the location of a data point on the map is determined by finding the
BMC for that data point.

**Coloring**: The Numero framework always creates a single map model. However, the map districts can be painted with different colors. This enables the user to create multiple colorings of the map to visualize regional differences. These colorings can be made for each variable, which helps to identify which parts of the map are particularly important for a specific phenomenon. Again, this is similar to a real city map where the districts are colored according to the income level of the local residents, or according to the mean age, smoking rates, obesity etc.

**Data point**: We define the term data point as a single uniquely identifiable row in a spreadsheet of data (with variables as columns). For instance, in the diabetic kidney disease dataset in the basic pipeline example, a data point refers to a patient (and vice versa) as there is only one row per patient. On the
other hand, if a patient dataset included multiple examination visits, and the
visits were organized on separate rows, a patient could be linked to multiple
data points.

**District**: A district refers to a pre-defined division of the map into uniformly sized areas. The districts are created mainly for technical reasons: using districts speeds up calculations and enables the estimation of map-related statistics. This is analogous to a real city being divided into districts to estimate regional demographics, for instance.

**District centroid**: The SOM algorithm works through the districts during the optimization of the data point layout on the map. The computational process eventually converges to a stable configuration that is stored as a set of district centroids. From a practical point of view, a district centroid represents the typical averaged profile that captures the characteristics of the data points within the district. In technical terms, the district centroid (also known as the prototype) contains the mean weighted data values across all the data points, where the weights are determined by the neighborhood function used in the SOM algorithm.

**Layout**: We make a distinction between what is a map, and what is the layout of data points on it. The layout is a vector of data point locations as coordinates, whereas the map is a more integrated concept that also includes information that is necessary to find the locations of new previously unseen data points, and to draw and paint the map in visual form.

**Map**: A map is a general term to describe the two-dimensional canvas onto which the multivariable data points are projected. The concept is analogous to a geographic map that indicates where people live, except that the location is not based on geography (i.e. physical distances), but comes from the data (i.e. distances = data-based similarities).

**Subgroup**: We expect that most uses of Numero will result in the subgrouping of a complex dataset. Visually, we define a subgroup via a contiguous set of adjacent districts on the map. Consequently, all the data points that are located within the set of districts are the subgroup members. Of note, our concept of a subgroup is different from subgroup discovery in data mining; we use the term as a replacement to the word 'cluster' to emphasize the lack of intrinsic structure in the original data space.

# Build information

```{r echo=FALSE}
sessionInfo()
Sys.time()
```
