template.Rmd

---
title: "Sample Size Determination for Case-Control Microbiome Studies with Shiny"
author: "Federico Mattiello"
date: '`r format(Sys.time(), "%d/%m/%Y - %X")`'
output:
  html_document:
highlight: haddock
toc: no
---

This is the R Markdown report containing inputs and outputs coming from the Shiny 
application for determining power and sample sizes of case-control microbiome studies. 
Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. 
For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

*****


##INPUTS

```{r parameters, echo=FALSE}
aux1 <- c(
  'Stratification'     = input$strata, 
  'MonteCarlo Repl.s' = input$MC, 
  'Pi-One type'    = input$kindPiOne, 
  'Number of OTUs'      = input$numOTUs, 
  'Changed OTUs 1'    = input$diffOTUs1, 
  'Changed OTUs 2'    = input$diffOTUs2,
  'Most/Least Abundant 1' = input$mostLeastAb1, 
  'Most/Least Abundant 2' = input$mostLeastAb2)
aux2 <- c(
  'Sample 1 Size'           = input$n1, 
  'Sample 2 Size'           = input$n2, 
  'Min Sample Size'      = input$sampleSizes[1L], 
  'Max Sample Size'      = input$sampleSizes[2L],
  'Rel. Abund. Diff. 1'    = paste(input$relAbund1, '%'), 
  'Rel. Abund. Diff. 2'    = paste(input$relAbund2, '%'), 
  'Significance Level'       = input$alpha,
  ' ' = ' ')
auxMat <- cbind(names(aux1), aux1, names(aux2), aux2)
dimnames(auxMat) <- list(NULL, rep.int(c('Description', 'Value'), 2L))
kable(auxMat)#, format = 'markdown')
```


Theta parameter instead is.  
```{r thetaChunck, echo = FALSE}
theta <- isolate(piOne()$'theta')
theta <- round(as.numeric(theta), 3L)
names(theta) <- names(isolate(piOne()$'piOne'))
theta
```


##OUTPUTS


###Settings graphs  

The following graphs show desired relative abundance distributions with their respective 
estimates based on the two simulated samples of size `r input$n1` and `r input$n2`, 
respectively. 

```{r piPlotsSave, echo=FALSE, fig.height=4}
# counts <- isolate(generatedCounts()$piDirList)
counts <- NULL
one <- isolate(piOne()$'piOne')
two <- isolate(piTwo()$'piTwo')
theta <- isolate(piOne()$'theta')
for (strRun in seq_along(one))
{
drawPiPlot(
  countsData = counts[[strRun]], 
  piOneObj = one[[strRun]][seq_len(input$numOTUs)],
  piTwoObj = two[[strRun]][seq_len(input$numOTUs)], 
  main = names(one)[strRun], 
  ylab = 'Abundance Proportions', 
  theta = round(theta[[strRun]], 3L))
}
```


###Library Size Distribution
```{r libSizesPlot, echo=FALSE}
libSizeData <- isolate(totCountsGen()$nReads)
extremeLibSize <- sapply(libSizeData, quantile, probs = 0.975)
sample1 <- density(libSizeData[[1L]], from = 0, to = extremeLibSize[1L])
sample2 <- density(libSizeData[[2L]], from = 0)
plot(sample1, ylim = range(sample1$y, sample2$y), lwd = 2, 
     main = 'Library Size Distribution')
lines(sample2, lwd = 2, col = 'red')
rug(sample1$x, lwd = 0.1)
rug(sample2$x, col = 'red', lwd = 0.1)
legend(x = 'topright', col = c('black', 'red'), lwd = 2, 
       legend = c('Controls', 'Cases'))
```  


```{r tryPlot2, results='asis', eval=FALSE, include=FALSE}
plot.new()
code <- get('self', env = attr(powSimPlot, 'observable'))$.label
code <- strsplit(code, split = '\n', fixed = TRUE)[[1]]
code <- code[-c(1L, length(code))]
code <- paste(code, collapse = '\n ')
eval(parse(text = code))
# isolate(powSimPlot())
```


###Power simulation with current settings  

Power results for the current settings:
```{r powSimMC, echo=FALSE}
auxSingleWald()$"table2Display"
```


###Power and sample size graph  

The next graph shows the **power** _vs._ the **sample size** behaviour where the latter 
goes from `r input$sampleSizes[1L]` to `r input$sampleSizes[2L]`.

```{r powPlot, echo=FALSE, fig.cap='Power vs. Sample Size plot'}
plot.new()
seqSizes <- isolate(mcHmpWaldResults()$'seqSizes')
powerData <- isolate(mcHmpWaldResults()$'pow')

### interpolate points with natural splines
## interpolated Wald data
waldLine <- spline(
  x = seqSizes, 
  y = powerData[, 1L], method = 'natural')
waldLine$y[waldLine$y > 1] <- 1

if (input$'wmwTest')
{
  ## interpolated WMW data, total power: at least one DA detected
  wmwLineTot <- spline(
      x = seqSizes, 
      y = powerData[, "WMW"], method = "natural")
  wmwLineTot$y[wmwLineTot$y > 1] <- 1
  ## interpolated WMW data, single OTU power average across DA OTUs
  wmwLineAvg <- spline(
      x = seqSizes, 
      y = powerData[, "WMW avg"], method = "natural")
  wmwLineAvg$y[wmwLineAvg$y > 1] <- 1
} else {}


par(mar = c(4, 4, 1, 1))
if(input$powerSimStart == 0 && is.null(powerData))
  {
  plot(0, 0 , type = 'n', xlim = input$sampleSizes, ylim = c(0, 1), 
       xlab = 'sample size', ylab = 'power')
  text(mean(input$sampleSizes), .5, labels = 'NOT YET STARTED', 
       cex = 2)
  } else {}

plot(seqSizes, powerData[, "Wald"], col = "black", bg = "red2",
    xlim = input$sampleSizes + c(-2L, 2L), 
    ylim = c(0, 1.1), pch = 21, lwd = 2, 
    type = "p", #main = "Power vs. Sample Size",
    xlab = "sample size", ylab = "power", cex = 2)
lines(waldLine, col = "red2", lty = 1, lwd = 2)

if (input$"wmwTest")
{
  points(seqSizes, powerData[, "WMW"], 
      col = "black", bg = "blue2", pch = 22, cex = 2, lwd = 2)
  lines(wmwLineTot, col = "blue2", lty = 1, lwd = 2)
  points(seqSizes, powerData[, "WMW avg"], 
      col = "black", bg = "blue2", pch = 23, cex = 2, lwd = 2)
  lines(wmwLineAvg, col = "blue2", lty = 2, lwd = 2)
  legend(x = "left", pch = 21:23, 
      pt.bg = c("red2", "blue2", "blue2"), col = "black", 
      lty = c(1, 1, 2), lwd = 2, pt.cex = 1.5,
      legend = c("Wald", "WMW", "WMW avg"))
} else {}


abline(h = c(0, input$alpha, 1), lty = 4, col = 'gray70', lwd = 2)
#text(x =  min(input$sampleSizes) + .1 * diff(input$sampleSizes), 
text(x =  diff(input$sampleSizes), 
     y = input$alpha, pos = 3, 
     labels = paste0('alpha = ', input$alpha), cex = 1)



### draw point and lines for reading power
clickCoords <- list(
  'x' = input$sampleSizes[1L] + .75 * diff(input$sampleSizes), 'y' = 0.4)
tmpCoords <- isolate(powPlotCoords)

if (!is.null(tmpCoords$x))
{
  clickCoords <- list('x' = tmpCoords$x, 'y' = tmpCoords$y)
} else {}

## function creating horizontal or vertical line
lineCoords <- function(dat, direction = c('horizontal', 'vertical'))
  {
  if (direction == 'horizontal')
    {
    list('x' = c(0, dat['x']), 'y' = rep.int(dat['y'], 2L))
    } else
      {
        list('x' = rep.int(dat['x'], 2L), 'y' = c(0, dat['y']))
        }
  }# END: function - lineCoords

## function that interpolates data
specificPowFunWald <- splinefun(
  x = seqSizes, 
  y = powerData[, 1L])

if (input$'wmwTest')
  {
  specificPowFunWMW <- splinefun(
    x = seqSizes, 
    y = powerData[, 2L])
  } else {}

coordsPointWald <- c(
  'x' = clickCoords$x, 
  'y' = min(1, specificPowFunWald(clickCoords$x)))

if (input$'wmwTest')
  {
  coordsPointWMW <- c(
    'x' = clickCoords$x, 
    'y' = min(1, specificPowFunWMW(clickCoords$x)))
  } else {}


## cross corresponding to click
points(clickCoords, pch = 3, lwd = 2, cex = 2, col = 'black')

### Wald power line
points(as.list(coordsPointWald), pch = 3, lwd = 2, cex = 1, col = 'red2')
## horizontal line
lines(lineCoords(coordsPointWald, direction = 'horizontal'), 
      lty = 4, col = 'red2', lwd = 2)
## vertical line
#              lines(x = c(resCoordsWald['x'], resCoordsWald['x']), 
#                  y = c(0, resCoordsWald['y']), 
#                  lty = 4, col = 'red2', lwd = 2)
lines(lineCoords(coordsPointWald, direction = 'vertical'), 
      lty = 4, col = 'red2', lwd = 2)
## text showing coordinates values
#              text(x = clickCoords$x, y = clickCoords$y, 
text(as.list(coordsPointWald), 
     labels = paste0(
       'Size=', round(clickCoords$x, 2L),
       '\n Power=', 
       round(coordsPointWald['y'], 3L)), 
     pos = 3, offset = 1, cex = 1)

### Wilcoxon-Mann-Whitney power line
if (input$'wmwTest')
  {
  points(as.list(coordsPointWMW), pch = 3, lwd = 2, cex = 1, col = 'blue2')
  ## horizontal line
  lines(lineCoords(coordsPointWMW, direction = 'horizontal'), 
        lty = 4, col = 'blue2', lwd = 2)
  ## vertical line
  lines(lineCoords(coordsPointWMW, direction = 'vertical'), 
        lty = 4, col = 'blue2', lwd = 2)
  ## text showing coordinates values
  #                text(x = clickCoords$x, y = clickCoords$y, 
  text(as.list(coordsPointWMW), 
       labels = paste0(
         'Size=', round(clickCoords$x, 2L),
         '\n Power=', 
         round(coordsPointWMW['y'], 3L)), 
       pos = 3, offset = 1, cex = 1)
  }# END: if - WMW test
```



### Power of individual OTUs
The following table appears only in case WMW test was selected.
Column __Type__ indicates the type of the OTU: 

- type2DA: second type of Differentially Abundant OTUs, the ones with the bigger 
  default effect size (percentage increase);
- type1DA: first type of Differentially Abundant OTUs, the ones with the smaller 
  default effect size.


\small
```{r chunkForOtuPowTable, echo=FALSE, results='asis'}
# tmp <- auxSingleWald$"otuPowTable"
tmp <- auxSingleWald()$"otuPowTable"
# otuType <- el(piTwo$"otuType")
otuType <- el(piTwo()$"otuType")
otuType <- otuType[otuType != 0]

if (length(unique(otuType)) == 2L)
{
  auxOtuType <- c("0", "type1DA", "type2DA")[otuType]
} else
{
#   auxOtuType <- factor(otuType, labels = c("compens", "type1DA", "type2DA"))
  auxOtuType <- c("compens", "0", "type1DA", "type2DA")[otuType]
}# END - ifelse: all unrequested OTUs are used for compensation

if (NCOL(tmp) > 1L)
{
  if (NCOL(tmp) > 2L) {
    auxNames <- rep(names(piOne()$"piOne"), each = 2L)
    auxLabels <- rep.int(c("ID", "Power"), length(piOne()$"piOne"))
    dimnames(tmp)[[2L]] <- paste(auxNames, auxLabels, sep = ".")
    count <- 1L
    for (stRun in 1L:(ncol(tmp)/2))
    {
      aux <- cbind(tmp[, count:(count + 1L)], "Type" = auxOtuType)
      print(kable(aux))
      count <- count + 2L
    }
  } else
  {
    tmp <- cbind(tmp, "Type" = auxOtuType)
    print(kable(tmp))
  }
} else {}# END - if: _wmwTest_ individual OTU powers
```


*****
\newpage

##R Session Information  

The R session information (including the OS info, R version and all 
packages used):

```{r session-info}
sessionInfo()
Sys.time()
```