

# Consider an infinite population where the association
# y = x^2 + e 
# x ~ U(0, 15)
# e ~ N(0, 10)

# We seek a linear approximation of the relationship between
# hat(y) = a + b*x

# This function simulates "n" observations from the population.
simulate <- function(n=20) {
  x <- runif(n, 0, 15)

# relationship between x and y.
sam_fit <- function(dat)
  lm(y ~ x, data=dat)

# We can approximate the 'population' best linear approximation by
# statistics that converge to a population quantity.
pop_fit <- sam_fit(simulate(1e6))

## > pop_fit
## Call:
## Coefficients:
##   (Intercept)            x  
##        -37.53        15.00  


# This function creates a level confidence region for the intercept
# and slope of the best linear approximation, and tests whether
# the region includes the corresponding population values.
# confidence region larger or smaller.
  s <- sam_fit(dat)
  d <- coef(s) - val
  v <- vcov(s)
  c <- qchisq(level+adj, 2)

# drawing repeated samples from the population and checking whether
# the associated confidence regions include the population values:

## > coverage20
## [1] 0.855

# Because the true coverage is less than the nominal coverage,
# coverage is equal to the nominal value:

## > coverage20a
## [1] 0.9476



]]> 0



I will present examples with R code in a later post.

]]> 1


# create a random file name

# create temporary graphic file
png(pngfile, type = "cairo")
y <- x
z <- outer(x, y, f)
z[] <- 1
op <- par(bg = rgb(1,1,1,0))

# write headers
pngsize <-[["size"]]

# open pipe to stdout and pass image data

# remove intermediate graphic

]]> 0




]]> 0

]]> 0


> set.seed(42)
> x <- rbeta(5000, 1, 5)
> y <- rbinom(5000, 1, x)
+                      length.out=100))



# binary kernel smoothing, with pointwise confidence band
# y - vector of binary outcomes
# x - vector of probabilities
# bw - bandwidth
# conf - confidence level for pointwise confidence band
                sx = x, conf = 0.95) {
  # normal kernel estimate
  lsx   <- length(sx);
  lx    <- length(x)
  rsm   <- rowSums(wts)
  sms   <- wts %*% y
  est   <- sms / rsm
  # Clopper-Pearson intervals
  list(x = sx, est = est,
       clo = clo, chi = chi, conf = conf)

# Calibration curve
# y - vector of binary outcomes
# x - vector of probabilities (predictions)
# bw - bandwidth
# conf - confidence level for pointwise confidence band
  ox <- order(x)
  x  <- x[ox]
  y  <- y[ox]
  plot(hx$mids, hx$counts,
       type='h', ann=FALSE,
       xlim = c(0,1),
       ylim = c(0,1),
       xlab = "Model",
       ylab = "Empirical", ...)


> par(mfrow=c(2,2))
> replicate(4, { 
+   x <- rbeta(5000, 1, 5)
+   y <- rbinom(5000, 1, x)
+                        length.out=100))
+ })



  2. Eliminate the inner margin.
  3. Use mfrow/mfcol.
  4. For each plot:
    • Create a plot without titles, axes, or labels.
    • Add titles, axes, and labels manually.
> CalibPlot <- expression({
+   x <- rbeta(5000, 1, 5)
+   y <- rbinom(5000, 1, x)
+   calib(y, x,
+                  length.out=100),
+ })
> #1,1
> eval(CalibPlot)
> mtexti("Column I", 3)
> #1,2
> eval(CalibPlot)
> mtexti("Column II", 3)
> mtexti("Row I", 4)
> #2,1
> eval(CalibPlot)
> axis(1)
> mtexti("Model Probability", 1, 0.75)
> axis(2)
> mtexti("Empirical Probability", 2, 0.75)
> #2,2
> eval(CalibPlot)
> mtexti("Row II", 4)





]]> 4
Descending Text in Righthand Margin of R Graphics à la mtext Sat, 29 Jun 2013 04:59:27 +0000 There was an R-help thread in January regarding text in the righthand margin of an R graphic, where the text should be rendered in reading order from top to bottom. The base R function mtext is used to plot text in the margin. But, mtext is only able to render text from left to right and bottom to top. The mtext function does not honor the srt plotting parameter, which may be used to rotate text using the text function (see ?par and ?text). Indeed, R Core member Uwe Ligges responded that "You cannot do that with mtext, you rather need a dirty hack with text(), I believe." Below is my dirty hack.

The mtexti function defined below takes arguments that are similar to mtext, with one major exception. Rather than specifying the margin line on which to render the text, the offset (in inches) from the edge of the plotting region is specified instead. Hence, the "i" in mtexti is intended to remind the user of this distinction.


# text - character, text to be plotted
# side - numeric, 1=bottom 2=left 3=top 4=right
# off  - numeric, offset in inches from the edge of the plotting region
# srt  - string rotation in degrees
# ...  - additional arguments passed to text()
mtexti <- function(text, side, off = 0.25,
                   srt = if(side == 2) 90  else
                         if(side == 4) 270 else 0, ...) {
    # dimensions of plotting region in user units
    usr <- par('usr')
    # dimensions of plotting region in inches
    pin <- par('pin')
    # user units per inch
    upi <- c(usr[2]-usr[1],
             usr[4]-usr[3]) / pin
    # default x and y positions
    xpos <- (usr[1] + usr[2])/2
    ypos <- (usr[3] + usr[4])/2
    if(1 == side)
        ypos <- usr[3] - upi[2] * off
    if(2 == side)
        xpos <- usr[1] - upi[1] * off
    if(3 == side)
        ypos <- usr[4] + upi[2] * off
    if(4 == side)
        xpos <- usr[2] + upi[1] * off
    text(x=xpos, y=ypos, text, xpd=NA, srt=srt, ...)

Here is an example:

plot(1, yaxt='n', xaxt='n', xlab='', ylab='', type='n')
mtexti("test", 1)
mtexti("test", 2)
mtexti("test", 3)
mtexti("test", 4)


]]> 4
Happy 2nd Birthday Luke! Tue, 23 Apr 2013 00:53:01 +0000 LukeBDay2013

]]> 2
Individual vs. Group Incentive for Weight Loss Thu, 11 Apr 2013 13:45:32 +0000 A new Annals of Internal Medicine article describes a study that compares two employer-sponsored financial incentive programs for promoting weight loss among obese employees. I first read about the article at the Pacific Standard. The study design is a randomized controlled prospective trial. The two programs are as follows:

  1. Program 1. Obese employees are given a monthly weight loss goal. If the goal is reached, the participant receives $100, otherwise the employer keeps the money. This is called the individual incentive.
  2. Program 2. Obese employees are organized into groups of five, and each participant is given a montly weight loss goal. A sum of $500 dollars is evenly split among those participants who achieve their monthly weight loss goal. In the event that no participant achieves their montly goal, the employer keeps the incentive money. This is called the group incentive.

The researchers found that the group incentive was associated with greater average weight loss than the individual incentive. This result is especially interesting from a psychological perspective, but I was most drawn to the issue of cost. I found it odd that the authors focused on the fact that "both designs used the same up-front allocation of resources". Presumably, this is to argue that the second program was more effective at no additional up-front cost. For example, the authors write: "Similar to that in the individual-incentive group, the up-front allocation of incentives for meeting weight-loss goals was $100 per participant per month (totaling $21 000)." But, the authors later write that, over a 24 week period: "Mean earnings were $514.70 (SD, $522.60) in the group-incentive group and $128.60 (SD, $165.50) in the individual-incentive group (mean between-group difference, $386.10 [CI, $201.00 to $571.30]; P < 0.001)." Hence, it's clear that the second program is more expensive, as one might expect. It's also a little odd that the study consisted mostly of women (89%). The allocation of race/ethnicity was also somewhat imbalanced.

I like that the authors used confidence intervals throughout to summarize the differences in average weight loss (and incentive earnings) between groups. They also used p-values, but I think this was unnecessary. The authors used multiple imputation for missing weights at 24 and 36 weeks. I've always had trouble accepting multiple imputation of outcomes, because the imputation depends so heavily on the method and model used for imputation. In the appendix, the authors write that weight was imputed "adjusting for incentive group, age, sex, race, education, household income, baseline weight, importance of controlling weight, and confidence in controlling weight". No additional details are given about the model, although the software used to implement the method is listed (SAS PROC MI and MIANALYZE). Finally, I felt this senctence was incomplete: "To maintain the type I error rate while testing the 3 hypotheses of primary interest, we used a Bonferroni correction to define an α of 0.0167 as our threshold for statistical significance." The authors neglected that this approach attempts to control the familywise type I error rate. This is an important omission.

]]> 0
Hurray! An R Connections API! Wed, 03 Apr 2013 13:00:34 +0000 I waited until April 3 to post this, so it wouldn't be taken as an April Fool's joke! A recent R News item announces that we now have a bona fide mechanism to create custom connections in R! This makes it possible to implement a custom connection in an R package. Until now, the only feasible alternative was to write a patch (but see also this trick), as I had done to implement a tty connection! Of course, this was a suboptimal solution that required an update for each new version of R. In fact, my last update was for R v2.13.1. This new development will make it much easier to share this code. I will also need to update my first draft of R Connection Internals, which is now three years old. Perhaps the document is better suited as a package vignette. More to follow later...

So that we know who to thank:

matt@deb6box$ svn log -r 62016
r62016 | urbaneks | 2013-02-21 14:29:44 -0500 (Thu, 21 Feb 2013) | 1 line

add API to create custom connections

Here is the header file itself, with credits. It looks like all of the Rconnection struct is made available. But, notice the warning!

matt@deb6box$ svn blame -r 62016
 11656     ripley /*
 11656     ripley  *  R : A Computer Language for Statistical Data Analysis
 62016   urbaneks  *  Copyright (C) 2000-2013   The R Core Team.
 11656     ripley  *
 11656     ripley  *  This program is free software; you can redistribute it and/or modify
 11656     ripley  *  it under the terms of the GNU General Public License as published by
 11656     ripley  *  the Free Software Foundation; either version 2 of the License, or
 11656     ripley  *  (at your option) any later version.
 11656     ripley  *
 11656     ripley  *  This program is distributed in the hope that it will be useful,
 11656     ripley  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 11656     ripley  *  GNU General Public License for more details.
 11656     ripley  *
 11656     ripley  *  You should have received a copy of the GNU General Public License
 42308     ripley  *  along with this program; if not, a copy is available at
 42308     ripley  *
 11656     ripley  */
 11656     ripley 
 62016   urbaneks #ifndef R_EXT_CONNECTIONS_H_
 62016   urbaneks #define R_EXT_CONNECTIONS_H_
 62016   urbaneks 
 11668     ripley #include 
 11668     ripley 
 62016   urbaneks #ifndef NO_C_HEADERS
 62016   urbaneks # include  /* for size_t */
 62016   urbaneks # include  /* for va_list */
 42677   urbaneks #endif
 42677   urbaneks 
 62016   urbaneks /* IMPORTANT: we do not expect future connection APIs to be
 62016   urbaneks    backward-compatible so if you use this, you *must* check the version
 62016   urbaneks    and proceed only if it matches what you expect
 62016   urbaneks 
 62016   urbaneks    We explicitly reserve the right to change the connection
 62016   urbaneks    implementation without a compatibility layer.
 62016   urbaneks  */
 62016   urbaneks #define R_CONNECTIONS_VERSION 1
 62016   urbaneks 
 45984     ripley /* this allows the opaque pointer definition to be made available 
 44013     ripley    in Rinternals.h */
 16472       luke #ifndef HAVE_RCONNECTION_TYPEDEF
 16472       luke typedef struct Rconn  *Rconnection;
 16472       luke #endif
 16472       luke struct Rconn {
 11656     ripley     char* class;
 11656     ripley     char* description;
 44013     ripley     int enc; /* the encoding of 'description' */
 11656     ripley     char mode[5];
 23228     ripley     Rboolean text, isopen, incomplete, canread, canwrite, canseek, blocking, 
 23228     ripley 	isGzcon;
 18583     ripley     Rboolean (*open)(struct Rconn *);
 11656     ripley     void (*close)(struct Rconn *); /* routine closing after auto open */
 11656     ripley     void (*destroy)(struct Rconn *); /* when closing connection */
 11656     ripley     int (*vfprintf)(struct Rconn *, const char *, va_list);
 11656     ripley     int (*fgetc)(struct Rconn *);
 32497     ripley     int (*fgetc_internal)(struct Rconn *);
 31166     ripley     double (*seek)(struct Rconn *, double, int, int);
 13305     ripley     void (*truncate)(struct Rconn *);
 11656     ripley     int (*fflush)(struct Rconn *);
 11656     ripley     size_t (*read)(void *, size_t, size_t, struct Rconn *);
 11656     ripley     size_t (*write)(const void *, size_t, size_t, struct Rconn *);
 59167     ripley     int nPushBack, posPushBack; /* number of lines, position on top line */
 11656     ripley     char **PushBack;
 12256         pd     int save, save2;
 32492     ripley     char encname[101];
 32492     ripley     /* will be iconv_t, which is a pointer. NULL if not in use */
 32492     ripley     void *inconv, *outconv;
 32497     ripley     /* The idea here is that no MBCS char will ever not fit */
 32497     ripley     char iconvbuff[25], oconvbuff[50], *next, init_out[25];
 32497     ripley     short navail, inavail;
 32497     ripley     Rboolean EOF_signalled;
 44101     ripley     Rboolean UTF8out;
 41765     ripley     void *id;
 41765     ripley     void *ex_ptr;
 11656     ripley     void *private;
 61527     ripley     int status; /* for pipes etc */
 16472       luke };
 11656     ripley 
 62016   urbaneks #ifdef  __cplusplus
 62016   urbaneks extern "C" {
 62016   urbaneks #endif
 11656     ripley 
 62016   urbaneks SEXP   R_new_custom_connection(const char *description, const char *mode, const char *class_name, Rconnection *ptr);
 62016   urbaneks size_t R_ReadConnection(Rconnection con, void *buf, size_t n);
 62016   urbaneks size_t R_WriteConnection(Rconnection con, void *buf, size_t n);
 13366     ripley 
 62016   urbaneks #ifdef  __cplusplus
 62016   urbaneks }
 25961     ripley #endif
 25961     ripley 
 62016   urbaneks #endif
]]> 0