inst/extdata/fth/plot_fth_diggle.R

In sje30/eglen2015: Bivariate spatial point patterns in the visual system: a reproducible review

## Generate Diggle plot of TH data.
## Wed 12 Mar 2003
## This file can be run in batch:
## R BATCH plot_fth_diggle.R
##
## This generates figure 5 of the 2003 paper.

## Here we use data that was generated by dmin_diggle.r.  Data is stored
## in a separate directory for safe-keeping.
load("./thdiggle.Rda")

res <- f0103.dres

kl <- function(k, show.l=TRUE) {
  ## Return K or L function for plotting.
  if (show.l)
    sqrt(k/pi)
  else
    k
}


diggle.th.fig <- function() {
  ## Generate the figure.  Call this function.
  ## Try the Arial fonts.
  ##postscript(file="fthdiggle.ps", width=7, height=8, onefile=F,
  ##           pointsize=12, horizontal=F)

  postscript(file="fthdiggle.ps", width=5, height=7, onefile=F,
             family = c("~/langs/R/afm/arialbd.afm", #note bold = default.
               "~/langs/R/afm/arial.afm",
               "~/langs/R/afm/ariali.afm",
               "~/langs/R/afm/arialbi.afm"),
             pointsize=12, horizontal=F)

  on.exit(dev.off())

  ## cex can change on order. 12pt text should be default.
  ## 12pt for (a,b,c), 10pt for axis labels and 8pt for numbers.

  ##par(mar=c(5,4,5,1), las=1)          #rug above plot.
  par(mar=c(6.5,3,2,0.3), las=1)            #rug below plot.

  par(mfrow=c(2,2), cex.axis=0.8, bty="n")
  ##par(mgp=c(2,0.5, 0))
  ##par(lwd=0.2)

  diggle.plot(res$csr.1, 1, "a  Spatial randomness", "(INL)", res$steps)
  diggle.plot(res$csr.2, 2, "b  Spatial randomness", "(GCL)", res$steps)
  ##diggle.plot(res$csr.0, 0, "c  INL+GCL", "CSR", res$steps)
  diggle.plot(res$csr.0, 0,
              "c  Spatial randomness ",
              expression((INL+GCL)),
              res$steps)
  diggle.plot(res$tor,  12,
              "d  Spatial independence",
              expression( (INL %*% GCL) ),
              res$steps)
}

diggle.plot <- function(l, my.ylab, label, minor, steps) {
  ## Do one of the sub plots.
  if(my.ylab == 1) {
    the.ylab <- expression(paste(L[i], " (", mu, "m)"))
  } else if (my.ylab == 2) {
    the.ylab <- expression(paste(L[g], " (", mu, "m)"))
  } else if (my.ylab == 0) {
    the.ylab <- expression(paste(L[i+g], " (", mu, "m)"))
  } else {
    the.ylab <- expression(paste(L[i%*%g], " (", mu, "m)"))
  }

  my.lwd <- 0.5
  par(mgp=c(2,0.5, 0))
  ## Plot out only to 400 um, to focus on differences.
  
  xmax <- ymax <- 400
  plot(steps, kl(l$real), type="l",
       col="red",
       xlab=expression(paste("Distance (", mu, "m)")),
       yaxt="n",
       xaxt="n",
       lwd=my.lwd*3,
       xlim=c(0, xmax), ylim=c(0, ymax),
       asp=1,
       ylab=the.ylab)
  
  par(mgp=c(2, 0.7, 0)); #need more space for the y axis when las=1
  axis(2);
  par(mgp=c(1, 0.7, 0)); 
  axis(1);

  segments(0, 0, xmax, xmax, col="black",
           lty=2, lwd=my.lwd) #leading diagonal.
  lines(steps, kl(l$lower), col="blue", lwd=my.lwd);
  lines(steps, kl(l$upper), col="blue", lwd=my.lwd) 

  ## for rug above plot, use line=4; use line=0 when rug below plot.
  mtext(label, side=3, adj=0, xpd=NA, cex=1, at=-20, line=0.1) #12 point text.
  mtext(minor, side=3, adj=0, xpd=NA, cex=1, at=20, line=-1.1) #12 point text.
  ## for rug above plot.
  ##draw.rank.rug(400, l$u, 500)
  draw.rank.rug(xmax, l$u, -180)         #draw rug below plot.
}


draw.rank.rug <- function(xhi, u, yht) {
  ## Draw the ranks of the experiments over the top of the plot.
  ## XHI is the x-coord for the maximum end of the rug plot.
  ## U is the vector of u values, with u[1] representing the
  ## experimental field.
  ## YHT is the y coordinate of where the rug should be positioned.
  
  xlo <- 0;
  ## normalise the rank to the maximum.
  rank2 <- u/max(u) * xhi               #absolute value not important.
  ##yht <- par()$usr[2] + 60
  nsim <- length(u) - 1
  axis(1, pos=yht-25, at=c(0, xhi), labels=c("0", "max"))
  axis(1, pos=yht-35, at=c(xhi/2), labels=c("u"), tick=F, cex.axis=1.0)
  
  xs <- rank2
  ys <- rep(yht, length(rank2));
  dy <- c(20, rep(10, nsim))
  segments(xs, ys-dy, xs, ys+dy, xpd=NA,
           col=c("red", rep("blue", nsim)),
           lwd=c(1.2, rep(0.4, nsim) ))
}

diggle.th.fig()