Realizations in Biostatistics

1. I left off the UK. The number of signatures is over 3 million, and contains by far the largest percentage of signatories.
2. 5 of the 9 top countries are neighbors, including the top 2. The other 4 are Australia, the US, Canada, and New Zealand, who are all countries that have strong ties to the UK.
3. This assumes no petition fraud, which I can't guarantee. I saw at least one Twitter posting telling people to use her (if the profile pic is to be believed) residence code. There is a section of the petition data showing constituency, so I'm wondering if it would be possible to analyze the petition for fraud. I'm not as familiar with British census data as I am with US, but I imagine a mashup of the two would be useful.

I guess I've been living in a bubble for a bit, but apparently there are a lot of people who still mistrust R. I got asked this week why I used R (and, specifically, the package rpart) to generate classification and regression trees instead of SAS Enterprise Miner. Never mind the fact that rpart code has been around a very long time, and probably has been subject to more scrutiny than any other decision tree code. (And never mind the fact that I really don't like classification and regression trees in general because of their limitations.)

At any rate, if someone wants to pay the big bucks for me to use SAS Enterprise Miner just on their project, they can go right ahead. Otherwise, I have got a bit of convincing to do.

I work in an environment dominated by SAS, and I am looking to integrate R into our environment.

Why would I want to do such a thing? First, I do not want to get rid of SAS. That would not only take away most of our investment in SAS training and hiring good quality SAS programmers, but it would also remove the advantages of SAS from our environment. These advantages include the following:

• Many years of collective experience in pharmaceutical data management, analysis, and reporting
• Workflow that is second to none (with the exception of reproducible research, where R excels)
• Reporting tools based on ODS that are second to none
• SAS has much better validation tools than R, unless you get a commercial version of R (which makes IT folks happy)
• SAS automatically does parallel processing for several common functions

So, if SAS is so great, why do I want R?

• SAS’s pricing model makes it so that if I get a package that does everything I want, I pay thousands of dollars per year more than the basic package and end up with a system that does way more than I need. For example, if I want to do a CART analysis, I have to buy Enterprise Miner, which does way more than I would need.
• R is more agile and flexible than SAS
• R more easily integrates with Fortran and C++ than SAS (I’ve tried the SAS integration with DLLs, and it’s doable, but hard)
• R is better at custom algorithms than SAS, unless you delve into the world of IML (which is sometimes a good solution).

I’m still looking at ways to do it, although the integration with IML/IML studio is promising.

Let’s say you are in the middle of a salary negotiation, and you want to know whether you should be aggressive in your offering or conservative. One way to help with the decision is to make a decision tree. We’ll work with the following assumptions:

• You are at a job currently making 50ドルk
• You have the choices between asking 60ドルk (which will be accepted with probability 0.8) or 70ドルk (which will be accepted with probability 0.2).
• You get one shot. If your asking price is rejected, you stay at your current job and continue to make 50ドルk. (This is one of those simplifying assumptions that we might dispense with later.)

This simplification of reality can be represented with a decision tree:

blogpayoff

I went ahead and put in the expected payoff for each of these decisions. Because the more conservative approach has a higher expected payoff, this model suggests that you should take the conservative approach.

One shortcoming clearly is that this decision tree only shows two decisions, but really you have a range of decisions; you are not stuck with 60ドルk or 70ドルk for asking price. You might go with 65ドルk, or 62ドル.5k, or something else. So what would be the optimal asking price?

Again, we look at expected payoff, which is asking price*probability(offer accepted) + 50ドルk * probability(offer rejected). In this case, we need to model the probability that offer is accepted over the range of possible offers, not just the two points. The logistic model works very well for modeling probability, and that’s what I will use here to extend the two-point model. In fact, a logistic model with two parameters can be fit exactly to two points, and so that is what I will use here.

Here is my commented R code to implement this model:

my.offer <- function (x1=60,py.x1=.2,x2=70,py.x2=.8,ev.no=50,high=100,p.payoff=1) {
 # return the offer to maximize expected payoff
 # this assumes a game with one decision and one consequence
 # you give an offer, and it is taken or refused. If taken, you receive a salary of
 # (a function of) the offer. If refused, you stay at the old job and receive a
 # salary of ev.no (presumably a current salary, but set to 0 if you are
 # unemployed).
 # the probability of rejection is modeled with a logistic function defined by
 # two points (x1,py.x1) and (x2,py.x2)
 # for example, if you expected a 20% rej. prob. with an offer of 140k, then
 # x1,py.x1 = 140,.2. Similarly with x2,py.x2
 # the expected payoff is modeled as offer*P(Yes|offer) + ev.no*P(No|offer),
 # perhaps with modifications to account for benefits, negotiation, etc. This
 # is defined in payoff function below.
 # finally, high is defined as anything above what you would be expecting to offer
 # and is used to create the plot limits and set the bounds in the optimization
 # routine.  # model the probability of no given salary offer
 # here we have a logistic function defined by (x1,py.x1) and (x2,py.x2)
 # note that qlogis is the inverse logit function
 # also, matrices in R are defined in column-major form, not row-major form like
 # FORTRAN, so we have to use 1,1,x1,x2 rather than 1,x1,1,x2
 theta <- solve (matrix (c (1,1,x1,x2),nc=2),matrix (qlogis (c (py.x1,py.x2)),nc=1))  # for plot of probability function
 xseq <- seq (ev.no,high,length=100)
 yseq1 <- 1/(1+exp (-theta[1]-theta[2]*xseq))  # model the expected payoff of an offer
 # model negotiations, benefits, and other things here
 # (a simple way to model benefits though is just to change ev.no)
 payoff <- function (x) {
 tmp <- exp (-theta[1]-theta[2]*x)
 return ( (ev.no + ifelse (ev.no>x*p.payoff,ev.no,x*p.payoff)*tmp)/(1+tmp) )
 }   yseq <- payoff(xseq)  # plots
 par (mfrow=c (1,2))
 plot (xseq,yseq1,type='l',xlab='Offer',ylab='P(No|X)')
 plot (xseq,yseq,type='l',xlab='Offer',ylab='Expected salary')  # no sense in even discussing the matter if offer < ev.no
 return (optimize (payoff,interval=c (ev.no,high),maximum=TRUE))
}

Created by Pretty R at inside-R.org

And here are the graphs and result:

image

> my.offer()
$maximum
[1] 61.96761

$objective
[1] 58.36087

So this model suggests that the optimum offer is close to 62ドルk, with an expected payoff of around 58ドルk. As a side effect, a couple of graphs are produced: giving the probability of rejection as a function of the asking price, and the expected salary (payoff) as a function of asking price.

So a few comments are in order:

• The value in this model is in varying the inputs and seeing how that affects the optimum asking price.
  
• The function I provided is slightly more complicated than what I presented in this post. You can model things like negotiation (i.e. you may end up at a little less than your asking price if you are not turned down right away), differences in benefits, and so forth. Once you have a simple and reliable baseline model with which to work, you can easily modify it to account for other factors.
  
• Like all models, this is an oversimplification of the salary negotiation process, but a useful oversimplification. There are cases where you want to be more aggressive in your asking, and this model can point those out.
  
• I commented the code profusely, but the side effects are probably not the best programming practice. However, this really is a toy model, so feel free to rip off the code.
  
• This model of course extends to other areas where you have a continuous range of choices with payoffs and/or penalties.
  
• If you are able to gather data on the probability of rejection based on offer, so much the better. You can then, instead of fitting an exact probability model, perform a logistic regression and use that as the basis of the expected payoff calculation.

Many statistical algorithms are taught and implemented in terms of linear algebra. Statistical packages often borrow heavily from optimized linear algebra libraries such as LINPACK, LAPACK, or BLAS. When implementing these algorithms in systems such as Octave or MATLAB, it is up to you to translate the data from the use case terms (factors, categories, numerical variables) into matrices.

In R, much of the heavy lifting is done for you through the formula interface. Formulas resemble y ~ x1 + x2 + …, and are defined in relation to a data.frame. There are a few features that make this very powerful:

• You can specify transformations automatically. For example, you can do y ~ log(x1) + x2 + … just as easily.
• You can specify interactions and nesting.
• You can automatically create a numerical matrix for a formula using model.matrix(formula).
• Formulas can be updated through the update() function.

Recently, I wanted to create predictions via Bayesian model averaging method (bma library on CRAN), but did not see where the authors implemented it. However, it was very easy to create a function that does this:

predict.bic.glm <- function(bma.fit,new.data,inv.link=plogis) {
# predict.bic.glm
# Purpose: predict new values from a bma fit with values in a new dataframe
#
# Arguments:
# bma.fit - an object fit by bic.glm using the formula method
# new.data - a data frame, which must have variables with the same names as the independent
# variables as was specified in the formula of bma.fit
# (it does not need the dependent variable, and ignores it if present)
# inv.link - a vectorized function representing the inverse of the link function
#
# Returns:
# a vector of length nrow(new.data) with the conditional probabilities of the independent
# variable being 1 or TRUE
# TODO: make inv.link not be specified, but rather extracted from glm.family of bma.fit$call

form <- formula(bma.fit$call$f)[-2] # extract formula from the call of the fit.bma, drop dep var
des.matrix <- model.matrix(form,new.data)
pred <- des.matrix %*% matrix(bma.fit$postmean,nc=1)
pred <- inv.link(pred)
return(pred)
}

The first task of the function finds the formula that was used in the call of the bic.glm() call, and the [-2] subscripting removes the dependent variable. Then the model.matrix() function creates a matrix of predictors with the original function (minus dependent variable) and new data. The power here is that if I had interactions or transformations in the original call to bic.glm(), they are replicated automatically on the new data, without my having to parse it by hand. With a new design matrix and a vector of coefficients (in this case, the expectation of the coefficients over the posterior distributions of the models), it is easy to calculate the conditional probabilities.

In short, the formula interface makes it easy to translate from the use case terms (factors, variables, dependent variables, etc.) into linear algebra terms where algorithms are implemented. I’ve only scratched the surface here, but it is worth investing some time into formulas and their manipulation if you intend to implement any algorithms in R.

For a long time, I have relied on R-bloggers for new, interesting, arcane, and all around useful information related to R and statistics. Now my R-related material is appearing there.

If you use the R package at all, R-bloggers should be in your feed aggregator.

• Writing to CSV, then opening in MS Excel, then copying/pasting into Word (or even saving in another format)
• Going through HTML to get to Word (for example, the R2HTML package)
• Using a commercial solution such as Inference for R (has this been updated since 2009?!)
• Using Statconn, which seems broken in the later versions of Windows and is not cross platform in any case
• Going through LaTeX to RTF

• Use cat statements and the latex command from the Hmisc package (included with every R installation) or xtable package (downloadable from CRAN) to create a LaTeX document with the table.
• Call the l2r command included with the LaTeX2RTF installation. (This requires a bit of setup, see below.)
• Open up the result in Word and do any further manipulation.

I know that R is free and I am actually a Unix fan and think Open Source software is a great idea. However, for me personally and for most users, both individual and organizational, the much greater cost of software is the time it takes to install it, maintain it, learn it and document it. On that, R is an epic fail.