R/densityBy.R

In psych: Procedures for Psychological, Psychometric, and Personality Research

#added the ability to pass ylim   7/7/26
"violin"  <- function(x,data=NULL,var=NULL,grp=NULL,grp.name=NULL, xlab=NULL, ylab=NULL,main="Density plot",vertical=TRUE,
   dots=FALSE, rain=FALSE,jitter=.05,
	alpha=1,errors=FALSE,eyes=TRUE,adjust=1,
	restrict=TRUE,xlim=NULL,add=FALSE,col=NULL,pch=20,scale=NULL, ...) {
	#call Violinby 
	violinBy(x=x,var=var,grp=grp,data=data,grp.name=grp.name,xlab=xlab, ylab=ylab, main=main,vertical=vertical,
    dots=dots,rain=rain, jitter=jitter,
    alpha=alpha,errors=errors,eyes=eyes,adjust=adjust,
    restrict=restrict,xlim=xlim,add=add,col=col,pch=pch,scale=scale, 
    ylim=NULL,xlas=0, ylas=2,order=NULL,...)
	}

#switched from density = 50 to alpha =.5  to speed up the plotting 
"violinBy"  <- function(x,var=NULL,grp=NULL,data=NULL,grp.name=NULL,xlab=NULL, ylab=NULL, main="Density plot",vertical=TRUE,
    dots=FALSE,rain=FALSE, jitter=.05,alpha= 1,errors=FALSE,eyes=TRUE,adjust=1,
    restrict=TRUE,xlim=NULL,add=FALSE,col=NULL,pch=20,scale=NULL, 
    ylim=NULL,xlas=0, ylas=2,order=NULL,...) {
   if(is.null(scale)) SCALE=.3  #how wide are the plots?

count.valid <- function(x) {sum(!is.na(x)) }

if(is.null(col)) {col <- c("blue","red","grey","purple","green","yellow")}  #default colors

 formula <- FALSE
   if(inherits(x, "formula")) {  ps <- fparse(x)
   	formula <- TRUE
   if(is.null(data)) {data <- get(ps$y) 
      var <- NULL #we specified the entire dataframe} else { #stop("You must specify the data if you are using formula input") 
     x <- data} else {x <- data[,ps$y,drop=FALSE] }  #added the drop = FALSE 11/28/20
   grp <- data[ps$x]
   if(is.null(xlab)) xlab  <- ps$x 
   names <- ps$x
   grp.n <- ps$x
   if(is.null(ylab)) ylab <- ps$y
   }
  if(!formula) {ps <- list(x = grp, y = x)  #convert from x, grp   to y ~ x format 
     if(is.null(data))  {data <- x
       
       if(!is.null(var))  {data <- data[,var, drop=FALSE]}
        var <-  ps$y <- colnames(data)
        if(is.null(var))  var <- paste("V",1:NCOL(data))
    #    colnames(x) <- colnames(data) <- ps$y  <- var}
        colnames(data) <- ps$y  <- var}
        if(!is.null(grp)) {x <- cbind(data[,ps$y,drop=FALSE], x[,grp,drop=FALSE])} else {
     x <- data[,ps$y,drop=FALSE] }
     data <- x 
     if(is.null(ylab)) ylab = "Observed "
     if(is.null(xlab))  xlab = "" 
    }  
       if(!is.null(grp)) {
      
    if(!is.data.frame(grp) && !is.list(grp) && (length(grp) < NROW(x))) grp <- x[,grp,drop=FALSE]}

 if(!is.null(var)) {if(missing(ylab) & (length(var) ==1)) {ylab <- var}
          x <- x[,var , drop = FALSE]}
 nvar <- nvarg <- NCOL(x)
 x <- char2numeric(x)
 
 col <- adjustcolor(col,alpha.f =alpha)        
#if(!is.null(grp)) {
# if(!is.data.frame(grp) && !is.list(grp) && (length(grp) < NROW(x))) grp <- x[,grp]

 if(!is.null(grp)) {stats <- describeBy(data[,ps$y],group = data[,ps$x],quant=c(0,1,.5,.25,.75),mat=TRUE,skew=FALSE) } else { stats <- describe(data[,ps$y,drop=FALSE],quant=c(0,1,.5,.25,.75),skew=FALSE)}
 nvarg <- NROW(stats)
 if(is.null(order) ) order <- 1:NROW(stats)
  meanX <- stats[order,"mean"]
  minx <- stats[order,"min"]
  maxx <- stats[order,"max"]
  medx <- stats[order,"Q0.5"]
   Q25 <- stats[order,"Q0.25"]
  Q75  <- stats[order,"Q0.75"]
  
  n.col <- length(col)
  nX <- stats[,"n"]
  col <- rep(col,nvarg) #too many, but this will do

#  Qnt <-  apply(x,2,function(xx) by(xx,grp,quantile,prob=c(0,1,.5,.25,.75),na.rm=TRUE))
# meanX <- apply(x,2,function(xx) by(xx,grp,mean,na.rm=TRUE))
# nX <- apply(x,2,function(xx) by(xx,grp,count.valid))
# meanX <- matrix(unlist(meanX))
#    Qnt <- matrix(unlist(Qnt),nrow=5)
#    ngrp <- ncol(Qnt)/nvar
#    nvarg <- ncol(Qnt)
#    rangex <- matrix(c(Qnt[1,],Qnt[2,]),nrow=2,byrow=TRUE)
# } else {Qnt <- apply(x,2,quantile,prob=c(0,1,.5,.25,.75),na.rm=TRUE)
# meanX <- apply(x,2,mean,na.rm=TRUE)}
# minx <- Qnt[1,]
# maxx <- Qnt[2,]
# medx <- Qnt[3,]
# Q25 <-  Qnt[4,]
# Q75 <-  Qnt[5,]
# #rangex <- apply(x,2,range,na.rm=TRUE)
# rangex <- matrix(c(Qnt[1,],Qnt[2,]),nrow=2,byrow=TRUE)

tot.n.obs <- nrow(x)

 grp.name <- names(table(grp)) 
  grp.name <- grp.name[order]                                                     
names <- grp.name  
if(is.null(grp.name) & !is.null(grp)) { grp.name <- paste(ps$x[1],stats[,"group1"])} else {if(length(grp.name) != NROW(stats)) names <- rownames(stats)}

# names <- rep(grp.name, length(ps$x)) 
# #if(length(ps$x) > 1) ) ) {names <- paste(rep(names,each=ngrp),grp.name[1:ngrp],sep=" ")} else {names <- grp.name}
# if(length(ps$y) > 1) names <- rep(names,length(ps$x))
#      col <- rep(col,nvar* ngrp)
d <- list(nvar)
if(is.null(xlim)) xlim <- c(.5,nvarg+.5)



for (i in 1:nvar) {
 
#if(!is.null(grp)) { if(restrict) {d[[i]] <- by(x[,i], grp ,function(xx) density(xx,na.rm=TRUE,from=rangex[1,i],to=rangex[2,i]))}
if(!is.null(grp)) { if(restrict) {d[[i]] <- by(x[,i], grp ,function(xx) density(xx,na.rm=TRUE,adjust=adjust,from=min(xx,na.rm=TRUE),to=max(xx,na.rm=TRUE)))} else {
d[[i]] <- by(x[,i], grp ,function(xx) density(xx,na.rm=TRUE)) }} else {
if(restrict) {d[[i]] <- density(x[,i],na.rm=TRUE,adjust=adjust,from=minx[i],to=maxx[i])} else {
d[[i]] <- density(x[,i],na.rm=TRUE)} }
}

if(is.null(ylim)) ylim <- c(min(minx),max(maxx)) 
if(is.null(xlim))   xlim <- c(min(minx),max(maxx)) 
if(!vertical) {temp <- xlim
               xlim <- ylim
               ylim <- temp}
             
if(vertical){plot(1:nvarg,meanX,ylim=ylim,xlim=xlim,axes=FALSE,xlab=xlab,ylab=ylab,main=main,pch=pch,...)
  axis(1,1:nvarg,names,...)
  axis(2,...)} else {plot(meanX,1:nvarg,xlim=xlim ,ylim=c(.5,nvarg+.5),axes=FALSE,xlab=ylab,ylab=xlab,main=main,pch=pch,...)
      axis(1,...)
     axis(2, 1:nvarg,names,...)}
  box()
 
if(!is.null(grp)) d <- unlist(d,recursive=FALSE)
d <- d[order]
rev <- (length(d[[1]]$y):1) #this prevents a line down the middle


for(i in 1:nvarg) {
if(!is.null(scale)) {width <- scale*sqrt(nX[[i]]/tot.n.obs)/max(d[[i]]$y)} else {width <- SCALE/max(d[[i]]$y)}
#polygon(width*c(-d[[i]]$y,d[[i]]$y[rev])+i,c(d[[i]]$x,d[[i]]$x[rev]),density=density,col=col[i],...)
#if(vertical) {polygon(width*c(-d[[i]]$y,  d[[i]]$y[rev])+i,  c(d[[i]]$x,  d[[i]]$x[rev]),col=col[i ],...)} else {
 #        polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(-d[[i]]$y,  d[[i]]$y[rev])+i,col=col[i],...)}
zero <- rep(0,length(d[[i]]$x))


if(rain) {if(vertical) {polygon(width*c(d[[i]]$y,  zero)+i,c( d[[i]]$x,d[[i]]$x[rev]  ),col=col[i ],...)} else {
         polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(d[[i]]$y,  zero)+i,col=col[i],...)}} else {
         #not rain
     if(vertical) {polygon(width*c(-d[[i]]$y,  d[[i]]$y[rev])+i,  c(d[[i]]$x,  d[[i]]$x[rev]),col=col[i ],...)} else {
                   polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(-d[[i]]$y,  d[[i]]$y[rev])+i,col=col[i],...)}}
         
  # if(vertical) {polygon(width*c(-d[[i]]$y,  d[[i]]$y[rev])+i,  c(d[[i]]$x,  d[[i]]$x[rev]+i),col=col[i ],...)} else {
   #      polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(-d[[i]]$y,  d[[i]]$y[rev])+i,col=col[i],...)}}
         
#draw the quartiles and  medians        
dmd <- max(which(d[[i]]$x <= medx[i]))
d25 <- max(which(d[[i]]$x <= Q25[i]))
d75 <- max(which(d[[i]]$x <= Q75[i]))

if(vertical) {
segments(x0=width*d[[i]]$y[dmd] +i ,y0=d[[i]]$x[dmd],x1=-(1-rain)*width*d[[i]]$y[dmd]+i,y1=d[[i]]$x[dmd],lwd=2)
segments(x0=width*d[[i]]$y[d25] +i ,y0=d[[i]]$x[d25],x1=-(1-rain)*width*d[[i]]$y[d25]+i,y1=d[[i]]$x[d25])
segments(x0=width*d[[i]]$y[d75] +i ,y0=d[[i]]$x[d75],x1=-(1-rain)*width*d[[i]]$y[d75]+i,y1=d[[i]]$x[d75])
  } else {
 
  segments(y0=width*d[[i]]$y[dmd] +i ,x0=d[[i]]$x[dmd],y1=(rain-1)*width*d[[i]]$y[dmd]+i,x1=d[[i]]$x[dmd],lwd=2)
segments(y0=width*d[[i]]$y[d25] +i ,x0=d[[i]]$x[d25],y1=(rain-1)*width*d[[i]]$y[d25]+i,x1=d[[i]]$x[d25])
segments(y0=width*d[[i]]$y[d75] +i ,x0=d[[i]]$x[d75],y1=(rain-1)*width*d[[i]]$y[d75]+i,x1=d[[i]]$x[d75])
  }
 } 

 
 if(errors) {   #currently only works for one DV and one IV

 if(vertical) {
    error.bars.by(x[,1],group=grp[,1],add=TRUE,lines=FALSE,eyes=eyes,v.labels=names,...) 

 } else {
           message("error bars are not implemented for horizontal plots")}
 }
 

if(dots) { if(!is.null(grp)) {
if(jitter >0) {if(length(ps$x) > 1 ) {stripchart(x[,ps$y]~ grp[,ps$x[1] ] +grp[,ps$x[2]] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...)} else {
        stripchart(x[,ps$y]~ grp[,ps$x[1] ] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...) }} else {
  if(length(ps$x) > 1 ) {stripchart(x[,ps$y]~ grp[,ps$x[1]]+ grp[,ps$x[,2]],vertical=vertical,add=TRUE,pch=pch,...)} else {
       {stripchart(x[,ps$y]~ grp[,ps$x[1]],vertical=vertical,add=TRUE,pch=pch,...)} 
       }}
       } else {
       if(jitter >0) {stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...)} else {
        stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...) }} 
 
                   }
                   
if(rain)  {  #does not work yet for reordered data

for(i in 1:nvar) {if(!is.null(grp)) {
  
 #tempat <- 1:length(table(grp)) -.1 + (i-1)*nvar
 tempat <- order[1:length(table(grp))] -.1 + (i-1)*nvar
if(jitter >0) {if(length(ps$x) > 1 ) {stripchart(x[,ps$y[i]]~ grp[,ps$x[1] ] +grp[,ps$x[2]] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch, at=tempat,...)} else {
      
        stripchart(x[,ps$y[i]]~ grp[,ps$x[1] ] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,at=tempat, ...) }} else {
  if(length(ps$x) > 1 ) {stripchart(x[,ps$y[i]]~ grp[,ps$x[1]]+ grp[,ps$x[,2]],vertical=vertical,add=TRUE,pch=pch,...)} else {
       {stripchart(x[,ps$y[i]]~ grp[,ps$x[1]],vertical=vertical,add=TRUE,pch=pch,...)} 
       }}
       }  else {
        tempat <- 1:nvar -.1 + (i-1)*nvar
       if(jitter >0) {stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,at=tempat,...)} else {
        stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,at=tempat,...) }} 
 
                   }
       
       }
       }
       
    
    #Improved 09/02/20

 "densityBy" <- function(x,var=NULL,grp=NULL,data=NULL,freq=FALSE,
   col=c("blue","red","black"),alpha=.5,adjust=1,
   ylim=NULL,xlim=NULL,xlab="Variable", 
    ylab="Density",main="Density Plot",legend=NULL) {
 formula <- FALSE
   if(inherits(x, "formula")) {  ps <- fparse(x)
   formula <- TRUE
   if(is.null(data)) stop("You must specify the data if you are using formula input") 
     var <- ps$y
   grp <- ps$x
    x <- data
   
   }
   

 n.grp <- length(table(x[grp]))
 if(length(col) < n.grp) col <- rainbow(n.grp)
if(!is.null(var)) x <- x[,c(var,grp),drop=FALSE]
 x <- char2numeric(x)
 col <- adjustcolor(col,alpha.f =alpha)
 if(missing(main) && freq) main="Frequency Plot"
 if(missing(ylab) && freq) ylab <- "N * density"
 if(missing(xlab) && (length(var) ==1)) xlab <- var
 if(!is.null(grp)) { 
    if(!is.data.frame(grp) && !is.list(grp) && (length(grp) < NROW(x))) grp <- x[,grp,drop=FALSE]
   d <- by(x[,var],grp,function(xx) density(xx, adjust=adjust,na.rm=TRUE))
   quants <-by(x[,var],grp,function(xx) quantile(xx,c(.25,.5,.75),na.rm=TRUE))
   maxiy <- rep(NA,length(d))
   rangex <- matrix(NA, ncol=2,nrow= length(d))
   for(i in 1:length(d) ) {if(freq){ maxiy[i] <- max(d[[i]]$n *d[[i]]$y)} else { 
                   maxiy[i] <- max(d[[i]]$y)}
                   rangex[i,] <- range(d[[i]]$x)
                   }
   maxy <- max(maxiy)
   ranges <- range(rangex)
   if(missing(ylim)) ylim <- c(0,maxy)
   if(missing(xlim)) xlim <- ranges
   plot(d[[1]],ylim=ylim,xlim=xlim,main=main,ylab=ylab,xlab=xlab ) 
   
   for (i in 1:length(d)) {
   if(freq) {scal <- d[[i]]$n
      d[[i]]$y <- d[[i]]$y*scal} 
   polygon(d[[i]] ,col=col[i])

   #find the y value corresponding to x median
   medy <- max(which(d[[i]]$x <= quants[[i]][2]))
   segments(x0 = quants[[i]][2],y0=0,x1=quants[[i]][2],y1=d[[i]]$y[medy],  lty= i)
   }
 } else {
 d <- density(x[,var],na.rm=TRUE,adjust=adjust)
   plot(d,main=main,ylab=ylab,xlab=xlab )
    polygon(d,col=col[1])
 }
 if(!missing(legend)) {
 location <- c("topleft","topright","top","left","right")
  grp.names <- paste(names(grp),names(table(grp)))
 leg.text <- grp.names
 legend(location[legend],legend=leg.text,col=c(1:n.grp),fill=c(1:n.grp), lty=c(1:n.grp))

 }
 invisible(d)
  }
# #names <- grp.name  } else {if(length(grp.name) != NROW(stats)) names <- rownames(stats)}
# 
# # names <- rep(grp.name, length(ps$x)) 
# # #if(length(ps$x) > 1) ) ) {names <- paste(rep(names,each=ngrp),grp.name[1:ngrp],sep=" ")} else {names <- grp.name}
# # if(length(ps$y) > 1) names <- rep(names,length(ps$x))
# #      col <- rep(col,nvar* ngrp)
# d <- list(nvar)
# if(is.null(xlim)) xlim <- c(.5,nvarg+.5)
# 
# 
# 
# for (i in 1:nvar) {
#  
# #if(!is.null(grp)) { if(restrict) {d[[i]] <- by(x[,i], grp ,function(xx) density(xx,na.rm=TRUE,from=rangex[1,i],to=rangex[2,i]))}
# if(!is.null(grp)) { if(restrict) {d[[i]] <- by(x[,i], grp ,function(xx) density(xx,na.rm=TRUE,adjust=adjust,from=min(xx,na.rm=TRUE),to=max(xx,na.rm=TRUE)))} else {
# d[[i]] <- by(x[,i], grp ,function(xx) density(xx,na.rm=TRUE)) }} else {
# if(restrict) {d[[i]] <- density(x[,i],na.rm=TRUE,adjust=adjust,from=minx[i],to=maxx[i])} else {
# d[[i]] <- density(x[,i],na.rm=TRUE)} }
# }
# 
# if(!is.null(order)) {meanX <- meanX[order]
#                 }
# if(is.null(ylim)) ylim=c(min(minx),max(maxx)) 
# if(vertical){plot(1:nvarg,meanX,ylim=ylim,xlim=xlim,axes=FALSE,xlab=xlab,ylab=ylab,main=main,pch=pch,...)
#   axis(1,1:nvarg,names,...)
#   axis(2,...)} else {plot(meanX,1:nvarg,xlim=c(min(minx),max(maxx)),ylim=c(.5,nvarg+.5),axes=FALSE,xlab=ylab,ylab=xlab,main=main,pch=pch,...)
#       axis(1,...)
#      axis(2, 1:nvarg,names,...)}
#   box()
# if(!is.null(grp)) d <- unlist(d,recursive=FALSE)
# if(!is.null(order)) {d[] <- d[order]
#                 }
# rev <- (length(d[[1]]$y):1) #this prevents a line down the middle
# 
# 
# for(i in 1:nvarg) {
# if(!is.null(scale)) {width <- scale*sqrt(nX[[i]]/tot.n.obs)/max(d[[i]]$y)} else {width <- SCALE/max(d[[i]]$y)}
# #polygon(width*c(-d[[i]]$y,d[[i]]$y[rev])+i,c(d[[i]]$x,d[[i]]$x[rev]),density=density,col=col[i],...)
# #if(vertical) {polygon(width*c(-d[[i]]$y,  d[[i]]$y[rev])+i,  c(d[[i]]$x,  d[[i]]$x[rev]),col=col[i ],...)} else {
#  #        polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(-d[[i]]$y,  d[[i]]$y[rev])+i,col=col[i],...)}
# zero <- rep(0,length(d[[i]]$x))
# 
# 
# if(rain) {if(vertical) {polygon(width*c(d[[i]]$y,  zero)+i,c( d[[i]]$x,d[[i]]$x[rev]  ),col=col[i ],...)} else {
#          polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(d[[i]]$y,  zero)+i,col=col[i],...)}} else {
#          #not rain
#      if(vertical) {polygon(width*c(-d[[i]]$y,  d[[i]]$y[rev])+i,  c(d[[i]]$x,  d[[i]]$x[rev]),col=col[i ],...)} else {
#                    polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(-d[[i]]$y,  d[[i]]$y[rev])+i,col=col[i],...)}}
#          
#   # if(vertical) {polygon(width*c(-d[[i]]$y,  d[[i]]$y[rev])+i,  c(d[[i]]$x,  d[[i]]$x[rev]+i),col=col[i ],...)} else {
#    #      polygon(c(d[[i]]$x,  d[[i]]$x[rev]), width* c(-d[[i]]$y,  d[[i]]$y[rev])+i,col=col[i],...)}}
#          
# #draw the quartiles and  medians        
# dmd <- max(which(d[[i]]$x <= medx[i]))
# d25 <- max(which(d[[i]]$x <= Q25[i]))
# d75 <- max(which(d[[i]]$x <= Q75[i]))
# 
# if(vertical) {
# segments(x0=width*d[[i]]$y[dmd] +i ,y0=d[[i]]$x[dmd],x1=-(1-rain)*width*d[[i]]$y[dmd]+i,y1=d[[i]]$x[dmd],lwd=2)
# segments(x0=width*d[[i]]$y[d25] +i ,y0=d[[i]]$x[d25],x1=-(1-rain)*width*d[[i]]$y[d25]+i,y1=d[[i]]$x[d25])
# segments(x0=width*d[[i]]$y[d75] +i ,y0=d[[i]]$x[d75],x1=-(1-rain)*width*d[[i]]$y[d75]+i,y1=d[[i]]$x[d75])
#   } else {
#  
#   segments(y0=width*d[[i]]$y[dmd] +i ,x0=d[[i]]$x[dmd],y1=(rain-1)*width*d[[i]]$y[dmd]+i,x1=d[[i]]$x[dmd],lwd=2)
# segments(y0=width*d[[i]]$y[d25] +i ,x0=d[[i]]$x[d25],y1=(rain-1)*width*d[[i]]$y[d25]+i,x1=d[[i]]$x[d25])
# segments(y0=width*d[[i]]$y[d75] +i ,x0=d[[i]]$x[d75],y1=(rain-1)*width*d[[i]]$y[d75]+i,x1=d[[i]]$x[d75])
#   }
#  } 
# 
#  
#  if(errors) {   #currently only works for one DV and one IV
# 
#  if(vertical) {
#     error.bars.by(x[,1],group=grp[,1],add=TRUE,lines=FALSE,eyes=eyes,v.labels=names,...) 
# 
#  } else {
#            message("error bars are not implemented for horizontal plots")}
#  }
#  
# 
# if(dots) { if(!is.null(grp)) {
# if(jitter >0) {if(length(ps$x) > 1 ) {stripchart(x[,ps$y]~ grp[,ps$x[1] ] +grp[,ps$x[2]] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...)} else {
#         stripchart(x[,ps$y]~ grp[,ps$x[1] ] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...) }} else {
#   if(length(ps$x) > 1 ) {stripchart(x[,ps$y]~ grp[,ps$x[1]]+ grp[,ps$x[,2]],vertical=vertical,add=TRUE,pch=pch,...)} else {
#        {stripchart(x[,ps$y]~ grp[,ps$x[1]],vertical=vertical,add=TRUE,pch=pch,...)} 
#        }}
#        } else {
#        if(jitter >0) {stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...)} else {
#         stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,...) }} 
#  
#                    }
#                    
# if(rain)  {
# for(i in 1:nvar) {if(!is.null(grp)) {
#   
#  tempat <- 1:length(table(grp)) -.1 + (i-1)*nvar
# if(jitter >0) {if(length(ps$x) > 1 ) {stripchart(x[,ps$y[i]]~ grp[,ps$x[1] ] +grp[,ps$x[2]] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch, at=tempat,...)} else {
#       
#         stripchart(x[,ps$y[i]]~ grp[,ps$x[1] ] ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,at=tempat, ...) }} else {
#   if(length(ps$x) > 1 ) {stripchart(x[,ps$y[i]]~ grp[,ps$x[1]]+ grp[,ps$x[,2]],vertical=vertical,add=TRUE,pch=pch,...)} else {
#        {stripchart(x[,ps$y[i]]~ grp[,ps$x[1]],vertical=vertical,add=TRUE,pch=pch,...)} 
#        }}
#        }  else {
#         tempat <- 1:nvar -.1 + (i-1)*nvar
#        if(jitter >0) {stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,at=tempat,...)} else {
#         stripchart(x ,vertical=vertical,method="jitter",jitter=jitter,add=TRUE,pch=pch,at=tempat,...) }} 
#  
#                    }
#        
#        }
#        }