R/plotDiags.R

In ggplotFL: Using ggplot2 in FLR

# plotDiags.R - 
# ggplotFL/R/plotDiags.R

# Copyright 2003-2007 FLR Team. Distributed under the GPL 2 or later
# Maintainer: Iago Mosqueira, JRC, Laurie Kell, ICCAT
# $Id:  $

setGeneric('diags', function(object, ...)
  standardGeneric('diags'))
 

diagPlots=function() list(#### Specification of what is to be plotted in diagnostic plot

"FLSR"    =data.frame(name=c("Functional Form","Residuals By Year","AR(1) Residuals","Residuals by Predictor","Normal Q-Q Plot","Residuals by Predicted"),
                      x   =c("x",              "year",             "residual",       "x",                     "qqx",            "yHat"),
                      y   =c("y",              "residual",         "residualLag",    "residual",              "qqy",            "residual"),
                      hat =c("yHat",           "ref",              "ref",            "ref",                   "qqHat",          "ref")),
                  
"seine"   =data.frame(name=c("Functional Form","Residuals By Year","AR(1) Residuals","Residuals by Predictor","Normal Q-Q Plot","Residuals by Predicted"),
                      x   =c("x",              "year",             "residual",       "x",                     "qqx",            "yHat"),
                      y   =c("y",              "residual",         "residualLag",    "residual",              "qqy",            "residual"),
                      hat =c("yHat",           "ref",              "ref",            "ref",                   "qqHat",          "ref")),

"FLBioDym"=data.frame(name=c("Functional Form","Residuals By Year","AR(1) Residuals","Residuals by Predictor","Normal Q-Q Plot","Residuals by Predicted"),
                      x   =c("x",              "year",             "residual",       "x",                     "qqx",            "yHat"),
                      y   =c("y",              "residual",         "residualLag",    "residual",              "qqy",            "residual"),
                      hat =c("yHat",           "ref",              "ref",            "ref",                   "qqHat",          "ref")),

"FLSAM"   =data.frame(name=c("Functional Form","Residuals By Year","AR(1) Residuals","Residuals by Predictor","Normal Q-Q Plot","Residuals by Predicted"),
                      x   =c("x",              "year",             "residual",       "x",                     "qqx",            "yHat"),
                      y   =c("y",              "residual",         "residualLag",    "residual",              "qqy",            "residual"),
                      hat =c("yHat",           "ref",              "ref",            "ref",                   "qqHat",          "ref")),

"FLXSA"   =data.frame(name=c("Functional Form","Residuals By Year","AR(1) Residuals","Residuals by Predictor","Normal Q-Q Plot","Residuals by Predicted"),
                      x   =c("x",              "year",             "residual",       "x",                     "qqx",            "yHat"),
                      y   =c("y",              "residual",         "residualLag",    "residual",              "qqy",            "residual"),
                      hat =c("yHat",           "ref",              "ref",            "ref",                   "qqHat",          "ref")))



## local function to calculated expected QQ line
qqLine <- function(x,y){ 
  qtlx <- quantile(x, prob=c(0.25,0.75), na.rm=T)
  qtly <- quantile(y, prob=c(0.25,0.75), na.rm=T)
      
  a <- (qtly[1]- qtly[2]) / (qtlx[1] - qtlx[2])
  b <- qtly[1] - qtlx[1] * a
      
  res <- c(a,b)
        
  names(res) <- NULL
  names(res) <- c("a","b")

 return(res)}

## creates plotting df
reform=function(x,y,hat,name,object) 
             data.frame("name"=name,
                        "x"   =as.numeric(object[,ac(x)]),
                        "y"   =as.numeric(object[,ac(y)]),
                        "hat" =as.numeric(object[,ac(hat)]))
 
## function to create diags DF for plotting from residuals
getDiag=function(object,diagPlots=diagPlots()[[class(object)]]) {
    
    ## extract values  
    x        = object$x
    y        = object$y
    yHat     = object$hat
    residual = object$residual

    #residualLag = c(NA,rev(rev(residual)[-1]))
    residualLag = object$residualLag

    qq = qqnorm(residual,plot.it=FALSE)
       
    res = data.frame(x=x, y=y, yHat=yHat, residual=residual, residualLag=residualLag,
                                qqx=qq$x, qqy=qq$y, year=object$year)
    
    qqpar = qqLine(res$qqx,res$qqy)[c("a","b")]

    res =data.frame(res,qqHat=qqpar["a"]*res$qqx+qqpar["b"],"ref"=0)
    
    res=mdply(diagPlots, reform, object=res)
    
    return(res)}

## FLSR
# setMethod("diags", signature(object="FLSR"),
#   function(object,sctPlot=diagPlots()$FLSR){
#     #   function(object, col="#CCFF33", group=1, stuff=facet_wrap(~age,scale="free"),
#     
#       res=model.frame(FLQuants("x"       =ssb(      object),
#                                "y"       =rec(      object),
#                                "hat"     =predict(  object),
#                                "residual"=residuals(object)),drop=TRUE)
#  
#      res=getDiag(res,sctPlot)
#      
#      if ("Functional Form" %in% sctPlot$name)  
#         res=rbind.fill(res,cbind(name="Functional Form",  
#                 model.frame(FLQuants("hat" =predict(object,ssb=FLQuant(seq(0,min(ssb(object)),length.out=21))),
#                                      "x"   =FLQuant(seq(0,min(ssb(object)),length.out=21))),drop=TRUE)[,-1]))
#                 
#     # ggplot(diags(pSR,srPlot))                                  + 
#     #     geom_point(aes(x,y))+facet_wrap(~name,scale="free")    +
#     #     geom_line(aes(x,hat),colour="red")                     +
#     #     stat_smooth(aes(x,y))
# 
#     return(res)})

## function to create diags DF for plotting from residuals
getDiag=function(object,diagPlots=diagPlots()[[class(object)]]) {
         
    qqpar=qqLine(object$qqx,object$qqy)[c("a","b")]
    res  =data.frame(object,qqHat=qqpar["a"]*object$qqx+qqpar["b"],"ref"=0)
    res  =mdply(diagPlots, reform, object=res)
    
    return(res)}

setMethod("diags", signature(object="FLSAM"),
  function(object){
      
    res=FLSAM::diags(object)
    res=res[do.call(order, res[,c("fleet","age","x")]),]
 
    res2=ddply(res,.(fleet,age), getDiag,diagPlots()[["FLSAM"]])
    res2=res2[do.call(order, res2[,c("fleet","age","x")]),]
     
    p.=ggplot(res2) + 
        geom_point(aes(x,y,group=age:fleet))                           +
        facet_wrap(~name,scale="free")                 +
        geom_path(aes(x,hat,colour="red",group=factor(fleet):factor(age)))+geom_smooth(aes(x,y))
 
  print(p.)
  return(invisible(p.))})