R/centilesPlot.R

In gamlss: Generalised Additive Models for Location Scale and Shape

# functions for centile plots
# last change on the 1 October 2019 
# to make sure that the function working even without xvar
# (Tim Cole sussestion)
# change MS Saturday, August 14, 2010
# centiles()  is modified by  Steve Ellison
# centiles.split() is  modified by Elaine Borghie 
# to include points at the end of the last interval 
#------------------------------------------------------------------------------------------
# centiles.split() is  modified by Elaine Borghie 
# to include points at the end of the last interval 
#------------------------------------------------------------------------------------------
centiles<-function (obj, 
                     xvar, 
                     cent = c(0.4, 2, 10, 25, 50, 75, 90, 98, 99.6),
                   legend = TRUE, 
                     ylab = "y", 
                     xlab = "x", 
                     main = NULL, 
                main.gsub = "@", 
                     xleg = min(xvar), #Note main and a substitution marker (an experiment!)
                     yleg = max(obj$y), 
                     xlim = range(xvar), 
                     ylim = range(obj$y), 
                     save = FALSE, 
                     plot = TRUE,
                   points = TRUE,  # this is new Saturday, August 14, 2010 MS
                      pch =  15,
                      cex = 0.5,
                      col = gray(0.7), 
             col.centiles = 1:length(cent)+2, 
             lty.centiles = 1, 
             lwd.centiles = 1,  #Handling for line appearance 
                     ...)        
{
if (!is.gamlss(obj))  stop(paste("This is not an gamlss object", "\n", ""))
# if(is.null(xvar)) stop(paste("The xvar argument is not specified", "\n", ""))
# attemp to take xvar argument out 1-10-2019
# new 2-10-19
if (missing(xvar))
{ 
    xvar <-  all.vars(obj$call$formula)[[2]]
  if (any(grepl("data", names(obj$call))))
    {
    DaTa <- eval(obj$call[["data"]]) #get(as.character(obj$call["data"])) get(as.character(obj$call["data"])) 
    xvar <- get(xvar, envir=as.environment(DaTa))
    } 
} 
xvarO <- deparse(substitute(xvar))  # get the name 
xvar <- try(xvar, silent = TRUE)    # get the vector
if  (any(class(xvar)%in%"try-error"))# if vector in DaTa not in the global Env 
{ # will fail therefore get it from DaTa
  DaTa <- eval(obj$call[["data"]])#get(as.character(obj$call["data"])) 
  xvar <- get(xvarO, envir=as.environment(DaTa))
}
# end of new 
  fname <- obj$family[1]
   qfun <- paste("q",fname,sep="")
  Title <- paste("Centile curves using", fname, sep = " ") #Note change to default title handling
   main <- if (is.null(main)) paste("Centile curves using", fname, sep=" ") 
else gsub(main.gsub,Title,main)    
  oxvar <- xvar[order(xvar)]
  oyvar <- obj$y[order(xvar)]  
if (is.matrix(obj$y)) # Monday, March 26, 2007 at 14:12
  {
  oyvar <-  obj$y[,1][order(xvar)] 
  ylim  <-  range(obj$y[,1])
  yleg = max(obj$y[,1])
  }
if (plot) 
  {
    lty.centiles <- rep(lty.centiles,length(cent))
    lwd.centiles <- rep(lwd.centiles,length(cent))
    col.centiles <- rep(col.centiles,length(cent))
    if (points==TRUE)
  {
      plot(oxvar, oyvar, type = "p", col = col, pch = pch, cex = cex,  
           xlab = xlab, ylab = ylab, xlim = xlim, ylim, ...) 
  }
  else
  {
      plot(oxvar, oyvar, type = "n", col = col, pch = pch, 
           xlab = xlab, ylab = ylab, xlim = xlim, ylim, ...) 
    }
    title(main)
  }      
   col <- 3 # set this to 1 if you do not want colour 
  lpar <- length(obj$parameters)
    ii <- 0
   per <- rep(0,length(cent))
  # per <- cent/100 
for(var in cent) 
  { 
    if(lpar==1) 
    {
      newcall <-call(qfun,var/100,
                     mu=fitted(obj,"mu")[order(xvar)]) 
    }
    else if(lpar==2)
    {
      newcall <-call(qfun,var/100,
                     mu=fitted(obj,"mu")[order(xvar)],
                     sigma=fitted(obj,"sigma")[order(xvar)]) 
    }
    else if(lpar==3)
    {
      newcall <-call(qfun,var/100,
                     mu=fitted(obj,"mu")[order(xvar)],
                     sigma=fitted(obj,"sigma")[order(xvar)],
                     nu=fitted(obj,"nu")[order(xvar)])
    }
    else 
    {
      newcall <-call(qfun,var/100,
                     mu=fitted(obj,"mu")[order(xvar)],
                     sigma=fitted(obj,"sigma")[order(xvar)],
                     nu=fitted(obj,"nu")[order(xvar)],
                     tau=fitted(obj,"tau")[order(xvar)]) 
    }    
    ii <- ii+1
    ll <- eval(newcall)
    if (plot)
    { 
      lines(oxvar, ll, col = col.centiles[ii], lty = lty.centiles[ii],lwd=lwd.centiles[ii], ...)
      # Modified by  Steve to include line types, colours etc
      # legend is moved below
      # if (legend==TRUE) legend(list(x=xleg,y=yleg), legend = cent, 
      #    col=c(3,4,5,6,7,8,9,10), lty=1, ncol=1, bg="white")# ,merge=TRUE)#, trace=TRUE)
    } 
    per[ii] <- (1-sum(oyvar>ll)/length(oyvar))*100
    if (!save) cat("% of cases below ", var,"centile is ", per[ii], "\n" )
  }
  if (plot) 
  {         #Legend moved outside plotting loop
    if (legend == TRUE) 
      legend(list(x = xleg, y = yleg), legend = cent, 
             col = col.centiles, lty = lty.centiles, lwd=lwd.centiles, 
             ncol = 1, ...)
  }     
  if (save) { return(cbind(cent,per))}
}
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------
# more than one centiles plot 
#------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------
# centiles.split will not work with less than 4 parameter distribution ??????

centiles.split<- function(obj, 
                          xvar, 
                          xcut.points = NULL, 
                          n.inter = 4 , 
                          cent = c(.4,2,10,25,50,75,90,98,99.6), 
                          legend = FALSE, 
                          main = NULL,
                          main.gsub = "@",   
                          ylab = "y",
                          xlab = "x",
                          ylim = NULL, # ms Friday, February 6, 2004 at 19:40
                          overlap = 0, 
                          save = TRUE,
                          plot = TRUE,
                          ...
)
{
  # function 1 -------------------------
  check.overlap <- function(interval)
  {
    if (!is.matrix(interval) ) {stop(paste("The interval specified is not a matrix."))}
    if (dim(interval)[2] !=2) 
    {
      stop(paste("The interval specified is not a valid matrix.\nThe number of columns should be equal to 2."))
    }
    crows = dim(interval)[1]
    for (i in 1:(crows-1))
    {
      #if (interval[i,2] != interval[i+1,1]) {interval[i+1,1]=interval[i,2]}
      if (!(abs(interval[i,2]-interval[i+1,1])<0.0001)) {interval[i+1,1]=interval[i,2]}
    }
    return(interval)
  }
  #------------------------------------------
  # function 2 
  get.intervals <- function (xvar, xcut.points ) 
  {
    if (!is.vector(xcut.points))  {stop(paste("The interval is not a vector."))}
    if ( any((xcut.points < min(xvar)) | any(xcut.points > max(xvar))))
    {stop(paste("The specified `xcut.points' are not within the range of the x: (", min(xvar),
                " , ", max(xvar), ")"))}
    int <- c(min(xvar), xcut.points, max(xvar))
    ii <- 1:(length(int)-1)
    r <- 2:length(int)
    x1 <- int[ii]
    xr <- int[r]
    if (any(x1>xr)) {stop(paste("The interval is are not in a increasing order."))}
    cbind(x1,xr)
  }
  #------------------------------------------
  # the main function starts here
  #------------------------------------------
  if (!is.gamlss(obj))  stop(paste("This is not an gamlss object", "\n", ""))
  # if(is.null(xvar)) 
  # { xvar <- seq(1,length(obj$y),1)
  # warning(paste("The xvar has been replace by an index 1:n", "\n", ""))  
  # }    # do we need this for centiles() 
  # new 2-10-19
  # new 2-10-19
  if (missing(xvar))
  { 
    xvar <-  all.vars(obj$call$formula)[[2]]
    if (any(grepl("data", names(obj$call))))
    {
      DaTa <- eval(obj$call[["data"]]) #get(as.character(obj$call["data"])) 
      xvar <- get(xvar, envir=as.environment(DaTa))
    } 
  } 
  xvarO <- deparse(substitute(xvar))  # get the name 
  xvar <- try(xvar, silent = TRUE)    # get the vector
  if  (any(class(xvar)%in%"try-error"))# if vector in DaTa not in the global Env 
  { # will fail therefore get it from DaTa
    DaTa <-eval(obj$call[["data"]])# get(as.character(obj$call["data"])) 
    xvar <- get(xvarO, envir=as.environment(DaTa))
  }
  #
  if(is.null(xcut.points)) 
  { # getting the intervals automatic
    g.in <- co.intervals(xvar, number=n.inter, overlap=overlap )
    if (overlap==0) g.in <- check.overlap(g.in) 
  }                 
  else
  { # if xcut.points is set
    g.in <- get.intervals(xvar, xcut.points ) 
  }
  howmany <- dim(g.in)[1]
  fname <- obj$family[1]
  qfun <- paste("q",fname,sep="")
  Title <-vector(length = howmany)
  for (i in 1:howmany)  Title[i] <- paste("(",paste(g.in[i,1],g.in[i,2], sep=", "), ")" )  
  if (is.null(main)) main <- Title
  else for (i in 1:howmany) main[i] <-gsub(main.gsub,Title[i],main[i])         
  #  if (is.null(main)) 
  #      {#rep(paste("Centile curves using", fname, sep=" "),howmany )
  #       main<-vector(length = howmany)
  #       for (i in 1:howmany)  main[i] <- paste("(",paste(g.in[i,1],g.in[i,2], sep=", "), ")" )
  #      }               
  #  else main #gsub(main.gsub,Title,main)    
  if (length(main)!=howmany) rep(main, lenght=howmany) 
  ncent <- length(cent)
  X <- matrix(0, nrow = ncent, ncol = howmany, 
              dimnames=list(as.character(seq(1,ncent)),
                            as.character(seq(1,howmany))))
  dimnames(X)[[1]] <- as.character(cent)
  dimnames(X)[[2]] <-  paste(substr(as.character(g.in[1:howmany,1]),1,7),"to",
                             substr(as.character(g.in[1:howmany,2]),1,7))
  var1 <- ceiling(howmany/2)
  
  #   op <- par(mfrow=c(var1,2),  col.axis="blue4", col.main="blue4", 
  #              col.lab="blue4",  col="darkgreen", bg="beige", ... )#
  if (howmany==2)
    layout(matrix(c(0,1,1,1,1,1,1,0,0,2,2,2,2,2,2,0), 8, 2))
  else if (howmany==3)
    layout(matrix(c(1,1,1,1,0,0,0,0,1,1,1,1,3,3,3,3,2,2,2,2,3,3,3,3,2,2,2,2,0,0,0,0), 8, 4))
  else if (howmany==4)
    layout(matrix(c(1,1,1,1,3,3,3,3,1,1,1,1,3,3,3,3,2,2,2,2,4,4,4,4,2,2,2,2,4,4,4,4), 8, 4))
  else {
    op <- par(mfrow=c(var1,2))#
  }   
  oxvar <- xvar[order(xvar)]
  oyvar <- obj$y[order(xvar)] 
  for (i in 1:howmany)
  {
    if(i==howmany) 
    {   ##### Include points at the end of the last interval HB
      yvar1 <- subset(oyvar, oxvar>=g.in[i,1]&oxvar<=g.in[i,2])
      xvar1 <- subset(oxvar, oxvar>=g.in[i,1]&oxvar<=g.in[i,2])
      mu1 <- subset(obj$mu.fv[order(xvar)],oxvar>=g.in[i,1]&oxvar<=g.in[i,2])
      sigma1 <- subset(obj$sigma.fv[order(xvar)],oxvar>=g.in[i,1]&oxvar<=g.in[i,2])
      nu1 <- subset(obj$nu.fv[order(xvar)], oxvar>=g.in[i,1]&oxvar<=g.in[i,2])
      tau1 <- subset(obj$tau.fv[order(xvar)],oxvar>=g.in[i,1]&oxvar<=g.in[i,2])
    }
    else
    {
      yvar1 <- subset(oyvar, oxvar>=g.in[i,1]&oxvar<g.in[i,2])
      xvar1 <- subset(oxvar, oxvar>=g.in[i,1]&oxvar<g.in[i,2]) 
      mu1 <- subset(obj$mu.fv[order(xvar)],oxvar>=g.in[i,1]&oxvar<g.in[i,2])
      sigma1 <- subset(obj$sigma.fv[order(xvar)],oxvar>=g.in[i,1]&oxvar<g.in[i,2])
      nu1 <- subset(obj$nu.fv[order(xvar)], oxvar>=g.in[i,1]&oxvar<g.in[i,2])
      tau1 <- subset(obj$tau.fv[order(xvar)],oxvar>=g.in[i,1]&oxvar<g.in[i,2])
    }   
    xlim1 <- c(g.in[i,1], g.in[i,2])
    ylim1 <- if (is.null(ylim)) range(yvar1)
    else ylim      
    xleg  <- min(xvar1) 
    yleg  <- max(yvar1) 
    obj1 <- obj
    obj1$y <- yvar1
    obj1$mu.fv <- mu1
    obj1$sigma.fv <- sigma1
    obj1$nu.fv <- nu1
    obj1$tau.fv <- tau1
    main1 <- main[i]
    sampleval <- centiles(obj1, xvar=xvar1, cent, legend=legend, main=main1, main.gsub = main.gsub, 
                          ylab, xlab ,xleg ,yleg, xlim=xlim1, ylim=ylim1, 
                          save=TRUE, plot=plot, ...)
    X[,i] <- sampleval[,2] 
  }
  if (howmany>4) par(op)
  layout(mat=1)             
  if (save) { return(X)}
}