# [R] plot of Bernoulli data

Warnes, Gregory R gregory_r_warnes at groton.pfizer.com
Thu Oct 4 17:54:15 CEST 2001

```Hi Frank,

Actually the 'wapply' function I have does 'overlapping' windows.  (It would
have to if it uses 1/5 of the range at each of 10 points, n'est-pas?)  The
main reason I created it was to add a 'local variance' estimate to plots in
an attempt to more easily detect changes in variability. EG:

x <- rnorm(500,sd=10)
y <- rnorm(500,sd=sqrt(abs(x)))
plot(x,y)

sdplus <- function(x,sign=1)
{
if(length(x)>0)
mean(x) + sign * sqrt(var(x))
else
NA
}
lines(wapply(x,y,fun=sdplus,width=1/5),col="red")
lines(wapply(x,y,fun=sdplus,sign=-1,width=1/5),col="red")

Not very pretty, but does what I need.

I'm not familiar with the trellace/lattice functions.  Is there a better or
more elegant way to do this?

-Greg

>  -----Original Message-----
>  From: Frank E Harrell Jr [mailto:fharrell at virginia.edu]
>  Sent: Thursday, October 04, 2001 9:06 AM
>  To: Warnes, Gregory R
>  Cc: 'Bill Simpson'; 'r-help at stat.math.ethz.ch'
>  Subject: Re: [R] plot of Bernoulli data
>
>
>  Greg - My model for this type of application is to
>  either use
>
>  - shingle objects with Trellis/Lattice to allow overlapping
>    intervals (and for this application there is good reason
>    for intervals to overlap if you are not using the
>    better loess/lowess solution)
>  or
>  - the cut2 function in the soon (I hope) to be beta-released
>    Hmisc library, with tapply and all the other stratified
>    computation functions.  cut2 has the same arguments as
>    cut but also arguments g (number of quantile groups) and
>    m (number of observations per interval)
>
>  cut2 at present does not handle fixed windows as you
>  describe, but in general such windows are problematic.
>
>  The bigger point though is that there may be a slight
>  advantage in separating the stratification code (e.g.,
>  cut2) from the computation code (e.g., tapply and
>  summarize (in Hmisc)).
>
>  Frank Harrell
>
>
>  "Warnes, Gregory R" wrote:
>  >
>  > In the next release of the gregmisc library there is a
>  function "wapply"
>  > that applies a specified function over 'windows' of x-y
>  data (code below).
>  > You would use it in this case as:
>  >
>  >  x<-sort(runif(100,1,20))
>  >  p<-pnorm(x,10,3)
>  >  y<-as.numeric(runif(x)<p)
>  >  plot(x,y)
>  >  lines(x,p)
>  >
>  >  # draw a line with 10 points, each computed on regions
>  spanning 1/5th of
>  > the x range
>  >  lines(wapply(x,y,n=10,width=1/5),col="red")
>  >
>  > -Greg
>  >
>  > # \$Id: wapply.R,v 1.2 2001/08/31 23:45:45 warneg Exp \$
>  > #
>  > "wapply" _ function( x, y, fun=mean, method="range",
>  >                     width=1/10, n=50, ...)
>  > {
>  >   method <- match.arg(method,
>  c("width","range","nobs","fraction"))
>  >
>  >   if(method == "width" || method == "range" )
>  >     {
>  >       if(method=="range")
>  >         width <- width * range(x)
>  >
>  >
>  >       pts _ seq(min(x),max(x),length.out=n)
>  >
>  >       result _ sapply( pts, function(pts,y,width,fun,XX,...)
>  >                       {
>  >                         low _ min((pts-width/2),max(XX))
>  >                         high _ max((pts+width/2), min(XX))
>  >                         return (fun(y[(XX>= low) &
>  (XX<=high)],...))
>  >                       },
>  >                       y=y,
>  >                       width=width,
>  >                       fun=fun,
>  >                       XX = x,
>  >                       ...)
>  >
>  >       return(x=pts,y=result)
>  >     }
>  >   else # method=="nobs" || method=="fraction"
>  >     {
>  >       if( method=="fraction")
>  >         width <- floor(length(x) * width)
>  >
>  >       ord <- order(x)
>  >       x  <- x[ord]
>  >       y  <- y[ord]
>  >
>  >       n  <- length(x)
>  >       center  <- 1:n
>  >       below <- sapply(center - width/2, function(XX) max(1,XX) )
>  >       above <- sapply(center + width/2, function(XX) min(n,XX) )
>  >
>  >       retval  <- list()
>  >       retval\$x  <- x
>  >       retval\$y  <- apply(cbind(below,above), 1,
>  >                          function(x) fun(y[x[1]:x[2]],...) )
>  >
>  >       return(retval)
>  >     }
>  >
>  > }
>  >
>  >  >  -----Original Message-----
>  >  >  From: Bill Simpson [mailto:wsi at gcal.ac.uk]
>  >  >  Sent: Tuesday, October 02, 2001 7:57 AM
>  >  >  To: r-help at stat.math.ethz.ch
>  >  >  Subject: Re: [R] plot of Bernoulli data
>  >  >
>  >  >
>  >  >  >  x<-sort(runif(100,1,20))
>  >  >  >  p<-pnorm(x,10,3)
>  >  >  >  y<-as.numeric(runif(x)<p)
>  >  >  >  plot(x,y)
>  >  >  >  lines(x,p)
>  >  >  df<-data.frame(x,y)
>  >  >  aggregate(df,list(x=(x<5),(x>5)&(x<10),(x>10) &
>  >  >  (x<15),(x>15)), FUN=mean)
>  >  >  gives me what I want but if anyone has a better way
>  to collect the
>  >  >  observations into bins I'd like to hear it. It would
>  be nice to
>  >  >  pass along something like
>  >  >  breaks<-c(5,10,15,20)
>  >  >
>  >  >  Thanks
>  >  >
>  >  >  Bill Simpson
>  >  >
>  >  >  -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
>  >  >  -.-.-.-.-.-.-.-.-.-
>  >  >  r-help mailing list -- Read
>  >  >  http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
>  >  >  Send "info", "help", or "[un]subscribe"
>  >  >  (in the "body", not the subject !)  To:
>  >  >  r-help-request at stat.math.ethz.ch
>  >  >  _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
>  >  >  _._._._._._._._._._
>  >  >
>  >
>  > LEGAL NOTICE
>  > Unless expressly stated otherwise, this message is
>  confidential and may be privileged. It is intended for the
>  unauthorized. If you are not an addressee, any disclosure
>  or copying of the contents of this E-mail or any action
>  taken (or not taken) in reliance on it is unauthorized and
>  may be unlawful. If you are not an addressee, please inform
>  the sender immediately.
>  >
>  -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
>  -.-.-.-.-.-.-.-.-.-
>  > r-help mailing list -- Read
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._

--
Frank E Harrell Jr              Prof. of Biostatistics & Statistics
Div. of Biostatistics & Epidem. Dept. of Health Evaluation Sciences
U. Virginia School of Medicine  http://hesweb1.med.virginia.edu/biostat
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```