[R] plotmath vector problem; full program enclosed

[David Winsemius]

On Jul 6, 2010, at 1:41 PM, Duncan Murdoch wrote:

> On 06/07/2010 10:54 AM, Paul Johnson wrote:
>> Here's another example of my plotmath whipping boy, the Normal  
>> distribution.
>>
>> A colleague asks for a Normal plotted above a series of axes that
>> represent various other distributions (T, etc).
>>
>> I want to use vectors of equations in plotmath to do this, but have
>> run into trouble.  Now I've isolated the problem down to a relatively
>> small piece of working example code (below).  If you would please run
>> this, I think you will see the problem.  When plotmath meets one
>> vector of expressions, it converts all but one to math, so in the
>> figure output i get, in LaTeX speak
>>
>> b1         $\mu-1.0 \sigma$    $\mu$
>>
>> All values except the first come out correctly.
>>
>> This happens only when I try to use bquote or substitute to get
>> variables to fill in where the 1.96, 1.0, and so forth should be.  In
>> the figure output, you should see a second axis where all of the
>> symbols are resolved correctly.
>>
>> As usual, thanks in advance for your help, sorry if I've made an
>> obvious mistake or overlooked a manual.
>>
>> ### Filename: plotMathProblem.R
>> ### Paul Johnson July 5, 2010
>> ### email me <pauljohn at ku.edu>
>>
>>  sigma <- 10.0
>>  mu <- 4.0
>>
>>  myx <- seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)
>>
>>  myDensity <- dnorm(myx,mean=mu,sd=sigma)
>>
>>  ### xpd needed to allow writing outside strict box of graph
>>  ### Need big bottom margin to add several x axes
>>  par(xpd=TRUE, ps=10, mar=c(18,2,2,2))
>>
>>  plot(myx, myDensity, type="l", xlab="", ylab="Probability Density ",
>> main=myTitle1, axes=FALSE)
>>  axis(2, pos= mu - 3.6*sigma)
>>  axis(1, pos=0)
>>
>>  lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes
>>
>>
>>  addInteriorLine <- function(x, m, sd){
>>    for (i in 1:(length(x))){
>>      lines( c(x[i],x[i]), c(0, dnorm(x[i],m=m,sd=sd)), lty= 14,  
>> lwd=.2)
>>    }
>>  }
>>
>>
>>  dividers <- c(qnorm(0.025), -1, 0, 1, qnorm(0.975))
>>  addInteriorLine(mu+sigma*dividers, mu,sigma)
>>
>>  # bquote creates an expression that text plotters can use
>>  t1 <-  bquote( mu== .(mu))
>>  mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)
>>
>>
>>  addInteriorLabel <- function(pos1, pos2, m, s){
>>    area <- abs(100*( pnorm(m+pos1*s,m,s)-pnorm(m+pos2*s, m,s)))
>>    mid <- m+0.5*(pos1+pos2)*s
>>    text(mid, 0.5*dnorm(mid,m,s),label=paste(round(area,2),"%"))
>>  }
>>
>>
>>  addInteriorLabel(dividers[1],dividers[2],  mu, sigma)
>>  addInteriorLabel(dividers[2],dividers[3],  mu, sigma)
>>  addInteriorLabel(dividers[3],dividers[4],  mu, sigma)
>>  addInteriorLabel(dividers[4],dividers[5],  mu, sigma)
>>
>>
>>
>>
>> ### Following is problem point: axis will
>> ### end up with correct labels, except for first point,
>> ### where we end up with "b1" instead of "mu - 1.96*sigma".
>>   b1 <- substitute( mu - d*sigma, list(d=round(dividers[1],2)) )
>>   b2 <- substitute( mu - sigma )
>>   b3 <- substitute( mu )
>>   b4 <- substitute( mu + sigma )
>>   b5 <- substitute( mu + d*sigma, list(d=round(dividers[5],2)) )
>> ##   plot(-20:50,-20:50,type="n",axes=F)
>>   axis(1, line=4,at=mu+dividers*sigma,
>> labels=c(expression(b1),b2,b3,b4,b5), padj=-1)
>>
>
> You want "as.expression(b1)", not "expression(b1)".  The latter  
> means "the expression consisting of the symbol b1".  The former  
> means "take the object stored in b1, and convert it to an  
> expression.".
>
> It's not perfect, because you'll end up with "mu - -1.96sigma" (i.e.  
> two minus signs), but it's closer than what you had.

Easily addressed in this case with "~" instead of "-". The value of  
"d" provides the minus:

b1 <- substitute( mu ~ d*sigma, list(d=round(dividers[1],2)) )

>
> Duncan Murdoch
>
>>
>>
>>
>> ### This gets "right result" but have to hard code the dividers
>>  b1 <- expression( mu - 1.96*sigma )
>>  b2 <- expression( mu - sigma )
>>  b3 <- expression( mu )
>>  b4 <- expression( mu + sigma )
>>  b5 <- expression( mu + 1.96*sigma )
>>  axis(1, line=8,at=mu+dividers*sigma, labels=c(b1,b2,b3,b4,b5),  
>> padj=-1)
>>
>>
>>
>>
>>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius, MD
West Hartford, CT