## ---- message=FALSE, warning=FALSE, echo=FALSE, fig.width=8.5, fig.height=6.5----

## Make grid of coordinates
makeCoords <- function(npts) {
  coords <- expand.grid(seq(1, 0, length.out = sqrt(npts)),
                        seq(0, 1, length.out = sqrt(npts)))[1:npts,]
  colnames(coords) <- c("y", "x")
  
  coords
}

## Define colors
colors <- c("#e19995", "#adaf64", "#4fbe9b", "#6eb3d9", "#d098d7")

## Create data.frame for points
set.seed(5)
df <- data.frame(color = factor(c(sample(colors, 16, replace = TRUE),
                                sample(colors, 120, replace = TRUE))),
                 size = c(abs(rnorm(16, 0.5, 0.25))+0.35,
                          abs(rnorm(120, 0.5, 0.25))+0.35),
                 set = c(rep('focal', 16),
                         rep('pool', 120)))

## Reorder factor level by colors
levels(df$color) <- colors

## Define focal and pool groups
focal <- df[df$set == 'focal',]
pool <- df[df$set != 'focal',]

## Sort by color
focal <- focal[order(focal$color, -focal$size),]

## Match ranges
library(nullranges)
set.seed(123)
x <- matchRanges(focal = focal,
                 pool = pool,
                 covar = ~color+size,
                 method = 'n', replace = TRUE)

## Sort by color
x <- x[order(x$color, -x$size),]

## Generate point grobs
library(grid)

## Focal set (sorted)
coords <- makeCoords(nrow(focal))
focalSet <- pointsGrob(x = unit(coords$x, 'native'),
                       y = unit(coords$y, 'native'),
                       pch = 21,
                       size = unit(focal$size, "char"),
                       gp = gpar(fill = as.character(focal$color), col = NA))

## Pool set (sorted)
coords <- makeCoords(nrow(pool))
poolSet <- pointsGrob(x = unit(coords$x, 'native'),
                      y = unit(coords$y, 'native'),
                      pch = 21,
                      size = unit(pool$size, "char"),
                      gp = gpar(fill = as.character(pool$color), col = NA))


## Matched set (sorted)
coords <- makeCoords(nrow(x))
matchedSet <- pointsGrob(x = unit(coords$x, 'native'),
                         y = unit(coords$y, 'native'),
                         pch = 21,
                         size = unit(x$size, "char"),
                         gp = gpar(fill = as.character(x$color), col = NA))


## Visualize sets
library(plotgardener)
pageCreate(width = 8.5, height = 6.5, showGuides = FALSE, xgrid = 0, ygrid =0)

## Pool set
plotGG(plot = poolSet, x = 1, y = 1, width = 2.5, height = 2.5)
plotText(label = "Pool Set",
            x = 2.25, y = 0.75,
            just = c("center", "bottom"),
            fontcolor = "#33A02C",
            fontface = "bold")

## Focal set
plotGG(plot = focalSet,
          x = 5.75, y = 1,
          width = (5/6)-(1/8),
          height = (5/6)-(1/8))
plotText(label = "Focal Set",
            x = 5.75 + ((5/6)-(1/8))/2,
            y = 0.75,
            just = c("center", "bottom"),
            fontcolor = "#1F78B4",
            fontface = "bold")

## Matched set
plotGG(plot = matchedSet,
          x = 5.75, y = 2.5,
          width = (5/6)-(1/8),
          height = (5/6)-(1/8))
plotText(label = "Matched Set",
            x = 5.75 + ((5/6)-(1/8))/2,
            y = 2.25,
            just = c("center", "bottom"),
            fontcolor = "#A6CEE3",
            fontface = "bold")

## Arrow and matchRanges label
plotSegments(x0 = 3.75, y0 = 2.5 + ((5/6)-(1/8))/2,
                x1 = 5.40, y1 = 2.5 + ((5/6)-(1/8))/2,
                arrow = arrow(type = "closed", length = unit(0.1, "inches")),
                fill = "black", lwd = 2)
plotText(label = "matchRanges()", fontfamily = 'mono',
            x = 4.625, y = 2.4 + ((5/6)-(1/8))/2,
            just = c("center", "bottom"))

## Matching plots
library(ggplot2)
smallText <- theme(legend.title = element_text(size=8),
                   legend.text=element_text(size=8),
                   title = element_text(size=8),
                   axis.title.x = element_text(size=8),
                   axis.title.y = element_text(size=8))

plot1 <-
  plotPropensity(x, sets=c('f','m','p')) +
  smallText +
  theme(legend.key.size = unit(0.5, 'lines'),
        title = element_blank())

plot2 <-
  plotCovariate(x=x, covar=covariates(x)[1], sets=c('f','m','p')) +
  smallText +
  theme(legend.text = element_blank(),
        legend.position = 'none')

plot3 <-
  plotCovariate(x=x, covar=covariates(x)[2], sets=c('f','m','p'))+
  smallText + 
  theme(legend.key.size = unit(0.5, 'lines'))


## Propensity scores
plotText(label = "plotPropensity()",
            x = 1.0, y = 4.24,
            just = c("left", "bottom"),
            fontface = "bold",
            fontfamily = 'mono')
plotText(label = "~color + size",
            x = 1.15, y = 4.5,
            just = c("left", "bottom"),
            fontsize = 10,
            fontfamily = "mono")
plotGG(plot = plot1,
          x = 0.9, y = 4.5, width = 2.5, height = 1.5,
          just = c("left", "top"))

## Covariate balance
plotText(label = "plotCovariate()",
            x = 3.65, y = 4.24,
            just = c("left", "bottom"),
            fontface = "bold",
            fontfamily = "mono")
plotText(label = covariates(x),
            x = c(3.9, 5.8), y = 4.5,
            just = c("left", "bottom"),
            fontsize = 10,
            fontfamily = "mono")
plotGG(plot = plot2,
          x = 3.40, y = 4.5, width = 1.8, height = 1.5,
          just = c("left", "top"))

plotGG(plot = plot3,
          x = 5.20, y = 4.5, width = 2.75, height = 1.5,
          just = c("left", "top"))

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(nullranges)
set.seed(123)
x <- makeExampleMatchedDataSet(type = 'GRanges')
x

## -----------------------------------------------------------------------------
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
                   pool = x[!x$feature1],
                   covar = ~feature2 + feature3)
mgr

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(GenomicRanges)
sort(mgr)

## ----paged.print=FALSE--------------------------------------------------------
overview(mgr)

## -----------------------------------------------------------------------------
plotPropensity(mgr)

## -----------------------------------------------------------------------------
plotCovariate(mgr)

## ---- message=FALSE, warning=FALSE, fig.height=6, fig.width=5-----------------
library(patchwork)
plots <- lapply(covariates(mgr), plotCovariate, x=mgr, sets = c('f', 'm', 'p'))
Reduce('/', plots)

## -----------------------------------------------------------------------------
matchedData(mgr)

## ---- results='hold'----------------------------------------------------------
covariates(mgr)
method(mgr)
withReplacement(mgr)

## ---- results='hold'----------------------------------------------------------
summary(focal(mgr))
summary(pool(mgr))
summary(matched(mgr))
summary(unmatched(mgr))

## -----------------------------------------------------------------------------
identical(matched(mgr), pool(mgr)[indices(mgr, set = 'matched')])

## -----------------------------------------------------------------------------
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
                   pool = x[!x$feature1],
                   covar = ~feature2 + feature3,
                   method = 'nearest',
                   replace = TRUE)
nn <- overview(mgr)
plotPropensity(mgr)

## ---- results='hold'----------------------------------------------------------
## Total number of duplicated indices
length(which(duplicated(indices(mgr))))

sum(table(indices(mgr)) > 1) # used more than once
sum(table(indices(mgr)) > 2) # used more than twice
sum(table(indices(mgr)) > 3) # used more than thrice

## -----------------------------------------------------------------------------
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
                   pool = x[!x$feature1],
                   covar = ~feature2 + feature3,
                   method = 'rejection',
                   replace = FALSE)
rs <- overview(mgr)
plotPropensity(mgr)

## -----------------------------------------------------------------------------
set.seed(123)
mgr <- matchRanges(focal = x[x$feature1],
                   pool = x[!x$feature1],
                   covar = ~feature2 + feature3,
                   method = 'stratified',
                   replace = FALSE)
ss <- overview(mgr)
plotPropensity(mgr)

## -----------------------------------------------------------------------------
## Extract difference in propensity scores
## between focal and matched sets
fmps <- sapply(c(nn, rs, ss), `[[`, "quality")
c('nearest', 'rejection', 'stratified')[which.min(fmps)]


## -----------------------------------------------------------------------------
sessionInfo()