Here we will step through some of the new visualization tools built into RMINC 1.3. So, let’s start with getting a dataset ready - we’ll look at sex differences in the brain.

# load libraries
library(RMINC)
## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display
## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE
library(plotrix) # for legends

# some weirdness potentially particular to my mac. Images are flipped, so fix that:
options(RMINC_flip_image=TRUE)

# load an existing dataset
load("~/data/CREB/CREB-analyses/wtonly.RData")
table(gfwt$Sex)
## 
##  F  M 
## 38 42
mask <- "/Users/jason/data/CREB/mask.mnc"

Let’s start by running a linear model testing the effect of sex on local brain volume.

vs <- mincLm(relJac02 ~ Sex, gfwt, mask=mask)
## Method: lm
## Number of volumes: 80
## Volume sizes: 201 439 298
## N: 80 P: 2
## In slice 
##  0  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99  100  101  102  103  104  105  106  107  108  109  110  111  112  113  114  115  116  117  118  119  120  121  122  123  124  125  126  127  128  129  130  131  132  133  134  135  136  137  138  139  140  141  142  143  144  145  146  147  148  149  150  151  152  153  154  155  156  157  158  159  160  161  162  163  164  165  166  167  168  169  170  171  172  173  174  175  176  177  178  179  180  181  182  183  184  185  186  187  188  189  190  191  192  193  194  195  196  197  198  199  200 
## Done

And we take a look

# first step - load the background anatomy
anatVol <- mincGetVolume("/Users/jason/data/CREB/21aug15_est_conserv-nlin-3.mnc")
## Start: 0 0 0
## Count: 201 439 298
# and a series of slices through the brain
mincPlotSliceSeries(mincArray(anatVol),           # the anatomical volume
                    mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
                    anatLow=700, anatHigh=1400,   # set anatomy thresholds
                    low=2.5, high=10,             # set stats thresholds
                    symmetric=T)                  # show separate upper and lower

This introduces the first, and one of the most useful, of the visualization functions: mincPlotSliceSeries. It does what the name claims - shows a series of slices through the brain, in this case overlaying the mincLm results on the average anatomy background. Before proceeding, one note: mincArray is called for both the anatomical volume and the stats; for the moment it’s a shortcut to add dimension information to each volume. Future versions of RMINC will likely get rid of that and store the necessary information right after outputs are created by the likes of mincLm.

Anyway, let’s make the plot a bit prettier: there are too many slices at the beginning and end that show no useful info, so let’s get rid of them.

mincPlotSliceSeries(mincArray(anatVol),           # the anatomical volume
                    mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
                    anatLow=700, anatHigh=1400,   # set anatomy thresholds
                    low=2.5, high=10,             # set stats thresholds
                    symmetric=T,                  # show separate upper and lower
                    begin=50, end=-90)            # remove slices from both sides

Looks prettier already. But it would help to have an indication of where the slices are and what the colours represent. So let’s add that:

mincPlotSliceSeries(mincArray(anatVol),           # the anatomical volume
                    mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
                    anatLow=700, anatHigh=1400,   # set anatomy thresholds
                    low=2.5, high=10,             # set stats thresholds
                    symmetric=T,                  # show separate upper and lower
                    begin=50, end=-90,            # remove slices from both sides  
                    legend="t-statistics")

And let’s show the other dimensions:

mincPlotSliceSeries(mincArray(anatVol),           # the anatomical volume
                    mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
                    anatLow=700, anatHigh=1400,   # set anatomy thresholds
                    low=2.5, high=10,             # set stats thresholds
                    symmetric=T,                  # show separate upper and lower
                    begin=35, end=-35  ,          # remove slices from both sides  
                    legend="t-statistics",
                    dimension = 1)  

mincPlotSliceSeries(mincArray(anatVol),           # the anatomical volume
                    mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
                    anatLow=700, anatHigh=1400,   # set anatomy thresholds
                    low=2.5, high=10,             # set stats thresholds
                    symmetric=T,                  # show separate upper and lower
                    begin=25, end=-25  ,          # remove slices from both sides  
                    legend="t-statistics",
                    dimension = 3)

Let’s do something fancier - or at least requiring more customization. Let’s create a slice showing the estimated effect (i.e. the beta) with the significant regions contoured.

The core function underlying it all is mincImage. So we begin by putting together the plot with a simple rendition of the underlying anatomy.

mincImage(mincArray(anatVol), slice=250)

A good start, but the axes are quite unnecessary, so let’s get rid of them.

mincImage(mincArray(anatVol), slice=250, axes=F)

And let’s return to a cleaner setting of thresholds:

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)

Now we add the beta coefficients on top.

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=heat.colors(255), underTransparent = T)

Starting to look good, but we’ll want slighty different colours and both the positive and the negative effects shown:

poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)

This slice shows any effects above (or below) a log jacobians of 0.02. Is any of that significant? Let’s find out via FDR.

mincFDR(vs, mask=mask)
## 
## Computing FDR threshold for all columns
## Start: 0 0 0
## Count: 201 439 298
##   Computing threshold for  F-statistic 
##   Computing threshold for  tvalue-(Intercept) 
##   Computing threshold for  tvalue-SexM
## Multidimensional MINC volume
## Columns:       F-statistic tvalue-(Intercept) tvalue-SexM 
## NULL
## Degrees of Freedom: c(1, 78) 78 78 
## FDR Thresholds:
##      F-statistic tvalue-(Intercept) tvalue-SexM
## 0.01   12.720632           2.830987    3.569582
## 0.05    7.627336           2.159808    2.763280
## 0.1     5.579595           1.821244    2.363340
## 0.15    4.411686           1.601513    2.101459
## 0.2     3.597821           1.432507    1.897762

Sure enough, there is. So let’s add contours to where we have significance.

poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57), add=T)

We’ve added lines for 10% and 1%; let’s label them and make them easier to distinguish:

poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57), 
            labels=c("10%", "1%"), lty=c(2,1), add=T)

Fairly pretty, but I’d rather have the labels in a legend rather than on the image itself. Doable:

poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57), 
            drawlabels=F, lty=c(2,1), col="white", add=T)
legend(200, 200, c("10% FDR", "1% FDR"), lty=c(2,1), col="white",
       text.col="white", bg="black", bty="n")

Still missing the colour bars - let’s add those to the bottom. This will also need some mucking about with the par plot settings to get the whole background to be black and pretty.

poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)
graycolours = gray.colors(255)

opar <- par(bg=graycolours[1]) # set the background to be the same colour as the under colour

mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400, col=graycolours)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57), 
            drawlabels=F, lty=c(2,1), col="white", add=T)
legend(200, 200, c("10% FDR", "1% FDR"), lty=c(2,1), col="white",
       text.col="white", bg="black", bty="n")
color.legend(20, -20, 130, 0, c("-0.1", "-0.02"), rev(negcolours), col="white", align="rb")
color.legend(170, -20, 280, 0, c("0.02", "0.1"), poscolours, col="white", align="rb")

# restore the old settings
par(opar)