lms calibrartion is called before reconstructiong the object

Author: mstasinopoulos

R/lms.R

@@ -1,15 +1,20 @@
 # Authors Mikis Stasinopoulos Bob Rigby and Vlasios Voudouris
 # created 11-04-12
-# upadated  8--6-14
-#-------------------------------------------------------------------------------
+# updated  8--6-14
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 # what is new June 2014
 # i)  the power transformation (x^p) function ptrans() now is defined to include 
 #      zero as log(x)
 # ii) the power transformation uses GAIC instead GD which was not reliable 
 #      since different transformation will use differt effective df's 
 # iii) a prediction function for an lms object is created
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 #This function is design to help the user to construct centile estimation.
 #It is only applicable when only "one" explanatory variable is available 
 #   (usually age).
@@ -23,13 +28,18 @@
 # If the response variable contains negative values and/or zeros then use 
 # the argument theSHASH theSHASH <-  c("NO", "SHASHo") or add any other 
 #  distribution which you think is appropriate
-#-------------------------------------------------------------------------------
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 #  the LMS familily of distributions
       LMS <- c("BCCGo",  "BCPEo", "BCTo")
 # the SHASH
  theSHASH <-  c("NO", "SHASHo")
-#-------------------------------------------------------------------------------
-#------------------------------------------------------------------------------- 
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 lms <- function(y, x,
         families = LMS,
             data = NULL, 
@@ -51,7 +61,7 @@ lms <- function(y, x,
               ... 
                 )  
 {
-#-------------------------------------------------------------------------------
+################################################################################
 # local function 
 findPower <- function(y, x, data = NULL,  lim.trans = c(0, 1.5), prof=FALSE, k=2,  c.crit = 0.01, step=0.1)  
 {
@@ -79,27 +89,30 @@ findPower <- function(y, x, data = NULL,  lim.trans = c(0, 1.5), prof=FALSE, k=2
   }  
   par
 }
-#-------------------------------------------------------------------------------
+################################################################################
 ptrans<- function(x, p) if (p==0) log(x) else I(x^p)
-#-------------------------------------------------------------------------------
+################################################################################
 # end of local function
-#-------------------------------------------------------------------------------
+################################################################################
+################################################################################
 ## the families to fit
-        FAM <- families      
+        FAM <- families
+################################################################################        
 ## which method
   method.pb <- match.arg(method.pb)
 ## get the variables  
        ylab <- deparse(substitute(y))
        xlab <- deparse(substitute(x))
           y <- if (!is.null(data)) get(deparse(substitute(y)), envir=as.environment(data)) else y
           x <- if (!is.null(data)) get(deparse(substitute(x)), envir=as.environment(data)) else x
-## ----------------------------------------------------------------------------- 
+################################################################################ 
 ## if need to check for transformation in x    
   if (is.null(fix.power))
     {
     if (trans.x) # if x^p
     {
-      par <- findPower(y, x,   lim.trans = lim.trans, prof=prof, k=k,  c.crit = c.crit, step=0.1)    
+      par <- findPower(y, x,   lim.trans = lim.trans, prof=prof, k=k,  
+                  c.crit = c.crit, step=0.1)    
       ox <- x
       x <-  ptrans(x,par)
     } 
@@ -112,6 +125,7 @@ ptrans<- function(x, p) if (p==0) log(x) else I(x^p)
 ##  starting  model for fitted values for mu (we assuming that this will work).
 ##  Note no sigma is fitted here
 ##  fit the model -------------------------------------------------------------- 
+################################################################################  
     cat('*** Initial  fit***'," \n")   
     switch(method.pb, 
         "ML"= {m0 <- gamlss(y~pb(x), sigma.formula=~1, data=data, c.crit = 0.01)},
@@ -122,7 +136,8 @@ ptrans<- function(x, p) if (p==0) log(x) else I(x^p)
         aic <- AIC(m0, k=k)
        fits <- c(fits, aic) 
   whichdist <- 0
-## fitting the diferent models in FAM ------------------------------------------       
+################################################################################  
+## fitting the different models in FAM 
   for (i in 1:length(FAM)) 
   {
     cat('*** Fitting', FAM[i], "***","\n")  
@@ -167,7 +182,21 @@ if(whichdist==0)
                               sigma.formula=~pb(x), 
                         data  =data, c.crit = 0.01)}) ## initial fit  finish
 }
-## changing the call t look better in the output -------------------------------
+################################################################################
+################################################################################
+## calibration -----------------------------------------------------------------
+  if (calibration)
+  {
+    calibration(m0, xvar=x, cent=cent, pch = 15, cex = 0.5, col = gray(0.7), ylab=ylab, xlab=xlab, legend=legend)	
+  } 
+  else 
+  {
+    centiles(m0, xvar=x, cent=cent, pch = 15, cex = 0.5, 
+             col = gray(0.7), ylab=ylab, xlab=xlab, legend=legend)		
+  }
+################################################################################
+################################################################################  
+## changing the call to look better in the output ------------------------------
 m0$call$mu.start <- NULL # this works OK
     m0$call$data <- substitute(data) # this is OK
   m0$call$family <- if(whichdist==0) "NO" else FAM[whichdist] # this is OK
@@ -194,7 +223,6 @@ m0$call$mu.start <- NULL # this works OK
 # FaM
 # m0$call$family <- substitute(Fam)
 # m0$call
-# 
 # sub("FAM[i]", as.character(FAM[i]), m0$call,)
 ## transformation needed -------------------------------------------------------        
      if (trans.x||!is.null(fix.power))   
@@ -211,16 +239,6 @@ m0$call$mu.start <- NULL # this works OK
     m0$ylab <- ylab
     m0$xlab <- xlab
   if (!is.null(data)) m0$call$data  <- substitute(data)
-## calibration -----------------------------------------------------------------
-  if (calibration)
-  {
-    calibration(m0, xvar=x, cent=cent, pch = 15, cex = 0.5, col = gray(0.7), ylab=ylab, xlab=xlab, legend=legend)	
-  } 
-  else 
-  {
-    centiles(m0, xvar=x, cent=cent, pch = 15, cex = 0.5, 
-             col = gray(0.7), ylab=ylab, xlab=xlab, legend=legend)		
-  }
  ## saving the fitted functions for mu sigma nu and tau  for prediction --------
 if ("mu"%in%m0$par)
 {
@@ -242,7 +260,7 @@ if ("tau"%in%m0$par)
   tauFun <- splinefun(x, fitted(m0,"tau"), method="natural")
   m0$tauFun <- tauFun
 }
- #------------------------------------------------------------------------------
+################################################################################
 #   deparse((m0$call))   
 #  toString(m0$call)  
 #  toString(format(m0$call))
@@ -253,9 +271,10 @@ if ("tau"%in%m0$par)
 class(m0) <- c("lms", class(m0))
   m0  # save the last model
 }
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 # se ??
 predict.lms <- function(object, 
                         what = c("all","mu", "sigma", "nu", "tau"),
@@ -321,9 +340,10 @@ predict.lms <- function(object,
   }  
 out    
 }
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 # this function is appropriate to used when fitted model fails to c
 calibration <- function(object, xvar, cent=c(0.4, 2, 10, 25, 50, 75, 90, 98, 99.6),#100*pnorm((-4:4)*2/3), 
               legend=FALSE, fan=FALSE,  ...)
@@ -337,18 +357,15 @@ calibration <- function(object, xvar, cent=c(0.4, 2, 10, 25, 50, 75, 90, 98, 99.
  }
  else
  {
-  cc <- centiles(object,  cent=percent, legend=legend, save=TRUE, ...)
+  cc <- centiles(object, xvar=xvar, cent=percent, legend=legend, save=TRUE, ...)
   cpp<-cbind(target=cent, calib.=cc[,1], sample=cc[,2])
   print(cpp, digits=3)
  }
 }
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
-#------------------------------------------------------------------------------- 
-#-------------------------------------------------------------------------------
-#-------------------------------------------------------------------------------
+################################################################################
+################################################################################
+################################################################################
+################################################################################
 z.scores <- function(object, y,x)
 {
   if (!is(object,"lms"))  stop(paste("This is not an lms object", "\n", "")) 
@@ -371,6 +388,9 @@ z.scores <- function(object, y,x)
   rqres <- qnorm(cdf)
   rqres  
 }
-
+ ################################################################################
+ ################################################################################
+ ################################################################################
+ ################################################################################