require(rms)
getHdata(support)
totcst <- support$totcst
totcst <- ifelse(totcst == 0., NA, totcst)

doplot <- function(predictor, type=c('spline','quadratic')) {
  type <- match.arg(type)
  r <- range(predictor, na.rm=T)
  xs <- seq(r[1], r[2], length=150)
  f <- switch(type,
              spline   = ols(log(totcst) ~ rcs(predictor, 5)),
              quadratic= ols(log(totcst) ~ pol(predictor, 2)))
  print(f)
  print(anova(f))
  p <- Predict(f, predictor=xs)
  pan <- function(...) {
    plsmo(predictor, log(totcst), add=T, trim=0, col='green', lwd=3, grid=TRUE)
    scat1d(predictor, grid=TRUE)
  }
  print(plot(p, predictor=xs, xlab=label(predictor), addpanel=pan,
             main=paste('n=', f$stats['n']), adj=0))
  invisible()
}

with(support,
     {
       doplot(pafi)
       doplot(scoma, 'quadratic')
     } )