Empirical confidence intervals, redux

Edit:
Evidently I've re-invented the wheel (yet again!).

---

There's nothing quite like spending an afternoon trying to get R to Plot What I Mean. In hopes of getting a paper off my desk and out the door, I got to revisit the overplotting conundrum previously presented. I wanted to plot this:

but allow the viewer some concept of the relative density of points near 0% error . There's way too much overplotting here to tell by eye. With only a little reflection, it was obvious solution I offered before was jingus, chiefly because of the sparsity the data (and their non-rectangular hull) prevented the 2D kernel estimate from working very well. You can see this in action if you a command that almost Plots What I Mean--smoothScatter, which uses a 2D kernel estimate to plot most of the data, then overplots outliers:


A kernel density estimator is too fancy for this problem. All that you really need to do is divide the abscissa into a grid that has enough points at each slice that taking the pth percentile of all the points in the slice gives you a decent estimate of the true value of the percentile there. This is easy to do with quantcut (make the slices with equal number of observations) or cut2 (make the slices with at least n observations). Then you can connect the pth percentile at each slice with linear or loess interpolation. I even cooked up a panel function that you can drop straight into your favorite lattice xyplot:

> xyplot(y~x, data=comp3, cex=.5, xlim=c(-5, 1e3), ylab="Percentage error", xlab=expression(lambda), ylim=c(-100, 100), panel=panel.confbars, npt=15, conf=.05)
Make npt larger to get a denser interpolation grid, or smaller to get smoother interpolation. conf gives the 100-2*conf % confidence interval you wish to bound, so the default is a 95% confidence interval.

panel.confbars = function(x, y, npt=25, conf=.025, ...){
require(gtools)
panel.xyplot(x, y, ...)

xcut = quantcut(x, q=seq(0, 1, length.out=npt))
xval = quantile(x, p=seq(0, 1, length.out=npt))
yval = matrix(unlist(tapply(y, xcut, quantile, p=c(conf, 1-conf))), nrow=2)

lconf = approxfun(xval[-1], yval[1,])
uconf = approxfun(xval[-1], yval[2,])

panel.loess(xval, lconf(xval), type="l", col="red", lwd=2)
panel.loess(xval, uconf(xval), type="l", col="red", lwd=2)

}

Empirical confidence intervals on scatter plots

Update: see this instead for a better solution.

---

Here's a little R-trick for you nerds out there in cyberland, but mainly for me, so I won't have to rediscover this little bit of code.

>xyplot(results3gr[,1]~results3gr[,2], xlim=c(0, 5000), cex=.1, ylim=c(-200, 200))



So you've got a dense scatter plot. Damn, there's points EVERYWHERE on this stupid thing, but how many more are being overplotted along the x axis versus some other blue-colored spot?

Well, how about using a 2D-kernel density estimator to get the relative density of points across a regular grid? Then you could find the ordinates correspond to 95% of the mass at a given slice of the grid.

So, the package KernSmooth should do the trick:

> ks = bkde2D(results3gr, bandwidth=c(10, 30), range.x=list(c(-200, 200), c(0, 1e4)))

Then we need to find the marginal mass at each abscissa of the grid, and the cumulative mass:

> ksmarg = apply(ks$fhat, 2, sum)
> kscum = apply(ks$fhat, 2, cumsum)

Then divide the cumulative mass by the marginal mass, and get a list of linear approximations


> kscond = as.data.frame((t(t(kscum)/ksmarg)))
> kscdf = list()
> for (i in 1:ncol(kscond)){kscdf[[i]] = approxfun(kscond[[i]], ks$x1)}

to get a cumulative distribution function for each abscissa of the grid. Now we get a list of 1% and 99% ordinates for each value in the abscissa:

> confbars = matrix(nrow=2, ncol = ncol(kscond))
> for ( i in 1:ncol(kscond)){confbars[1, i] = kscdf[[i]](.025); confbars[2, i] = kscdf[[i]](.975)}


Now do some lattice witchcraft and plot:

> xyplot(results3gr[,1]~results3gr[,2], panel=function(x, y, ...){panel.xyplot(x, y, cex=.3); panel.xyplot(ks$x2, confbars[1,], type="l", col="black"); panel.xyplot(ks$x2, confbars[2, ], type="l", col="black")}, xlim=c(0, 5e3), ylim=c(-200, 200))



The bars appear to be over-conservative, probably because of sparsity...but oh well. It'll have to do.


Hat tip to Yihui to for the reference to KernSmooth.

The joy of awk: transforming group vCards

A system's analyst's life involves more data munging than I'd like to admit. Despite all of the advances in computing, sometimes the old tools are the best tools. Such is the case with awk and its even-more abstruse cousin, sed.

Awk really isn't all that bad. The way to envision it is a giant while-loop with branching if-blocks that use regular expressions to match parts of the current lines and then do stuff like print parts of lines or calculate things. Its syntax and grammar is more friendly to casual use than perl or bash or sed.

Anyway, here's an awk program I spent about 60 minutes writing that transforms group vCards into csv files, because Microsoft Outlook only reads the first entry from a group vCard and Address Book on Macintosh doesn't write anything except group vCards. Maybe it will be of use to you.

vcard2csv.awk

BEGIN {newcard = 0
  FS = ":"
  name = ""
  orgname = ""
  email = ""
  OFS = ","
  }
 $1 == "BEGIN" {
  newcard = 1
  next}
 newcard == 1 {
     if ($1 == "FN"){
      name = $2 
      }
     if ($1 == "ORG"){
      split($2, cleanedorg, ";")
      orgname = cleanedorg[1]
      }
     if ($1 == "EMAIL;type=INTERNET;type=WORK;type=pref"){
      email = $2
      }
    }
 $1 == "END" {
  print name, orgname, email
  name = ""
  orgname = ""
  email = ""
  newcard = 0
  }

You can run it by typing awk -f vcard2csv.awk < export.vcf > export.csv in a terminal. You might need to convert the line endings in the vcf to Unix (LF) first.