Code that works outside a custom function but not inside

["s b" <sbmailinglist@gmail.com>]

Dear S-plus users,

I am trying to plot some density histograms from a data.frame of simulated data. Since I want to plot several sets of simulated data, I have created this function myhisto which supposedly should do the work. Funny thing: the code within the function does not work when the data are passed to the function, BUT it does when it is used outside of a function (see the end of the code below). Can anyone explain what I am doing wrong in my function?

Thanks you in advance,
Sebastien

PS: if anyone can also explain why the attribute which.panel does not work in strip, I would also be grateful!


############# Code starts here

# Custom functions
# Random generation of 10000 individuals for each input sd

myranddf <- function(mysds){
  nsd <- length(mysds)
  tmp <- c()
  for (i in seq(nsd)){
    tmp <- c(tmp,rnorm(10000, mean = 0, sd = mysds[i]))
  }
  mydf <- data.frame(X=tmp, SD=rep(mysds, each = 10000))
  mydf
}


# Histogram plot

myhisto <- function(mydata,nBar){

  mysds <- levels(factor(mydata[,2]))
  nID <- length(mysds)
  
  myplot <- histogram(~mydata[,1]|mydata[,2],
                      data = "">                      as.table=TRUE,
                      type = "percent",
                      nint=nBar,
                      xlim=c(0,1),
                      xlab="X",
                      ylab="Probability density",
                      layout=c(2,2),
                      main=NULL,
                      strip = function(which.panel,...){
                        lines(c(0,1,1,0,0),c(0,0,1,1,0),col=1)
                        text(rep(0.5,nID),rep(0.5,nID),
                             paste("SD = ",mysds[which.panel],sep=""),
                             cex = 0.75)},
                      key=NULL)
  print(myplot)
}

# Create data set for typical logit transform with additive variability model
# X = exp(ln(TVX/(100-TVX))+eta)/(1+exp(ln(TVX/(100-TVX))+eta))

TVX = 50
mysds = c(0.05,0.1,0.25,0.75)*TVX

mydf <- myranddf(mysds)
mydf[,1] <- exp(log(TVX/(100-TVX))+mydf[,1])
mydf[,1] <- mydf[,1]/(1+mydf[,1])

# Plotting using the function

myhisto(mydf,50)

# Plotting outside a function

mydata <- mydf
nBar <- 50

  myplot <- histogram(~mydata[,1]|mydata[,2],
                      data = "">                      as.table=TRUE,
                      type = "percent",
                      nint=nBar,
                      xlim=c(0,1),
                      xlab="X",
                      ylab="Probability density",
                      layout=c(2,2),
                      main=NULL,
                      strip = function(which.panel,...){
                        lines(c(0,1,1,0,0),c(0,0,1,1,0),col=1)
                        text(rep(0.5,4),rep(0.5,4),
                             paste("SD = ",mysds[which.panel],sep=""),
                             cex = 0.75)},
                      key=NULL)
  print(myplot)

################ END of code