[R] plot of Bernoulli data

Frank E Harrell Jr fharrell at virginia.edu
Thu Oct 4 18:19:55 CEST 2001


Thanks for elaborating Greg - I see this more clearly now.
Trellis/lattice handles overlapping intervals that are
used to make multiple panels (using shingles).  That
would be for a third variable.  I don't think shingles
will help with what you are doing, so I see the
utility of your function more better now.  Thanks -Frank


"Warnes, Gregory R" wrote:
> 
> 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
>  >  addressee(s) only. Access to this E-mail by anyone else is
>  >  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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
> _._

-- 
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list