R/ba.R

Defines functions gamlss.ba ba.control ba

Documented in ba ba.control gamlss.ba

# This is the interface to use the gam() and bam() function of Simon Wood 
# where fitting and initialization are done separately
# Authors: Daniil Kiose and Vlasiis Voudouris
# latest 03-OCT-16 MS
#--------------------------------------------------------------------------------------
ba <-function(formula, control = ba.control(...), ...)  
  {
  #------------------------------------------
  # function starts here
  #------------------------------------------    
  scall <- deparse(sys.call(), width.cutoff = 500L) 
  if (!is(formula, "formula")) 
    stop("formula argument in ba() needs a formula starting with ~")
  # get where "gamlss" is in system call, it can be in gamlss() or predict.gamlss()
  rexpr <- grepl("gamlss",sys.calls()) ## 
  for (i in length(rexpr):1) { 
    position <- i # get the position
    if (rexpr[i]==TRUE) break
  }
  gamlss.env <- sys.frame(position) #gamlss or predict.gamlss
  ## get the data
  ## this has been modified on the 12-12-14 to make sure that 
  ##  if model.frame.gamlss() is used as for example in term.plot() the
  ## function does not fail (It need to implemented to all smoother using formulea?)
  if (sys.call(position)[1]=="predict.gamlss()") { # if predict is used 
    Data <- get("data", envir=gamlss.env)
  } else if (sys.call(position)[1]=="gamlss()") { # if gamlss() is used
    if (is.null(get("gamlsscall", envir=gamlss.env)$data)) { # if no data argument but the formula can be interpreted
      Data <- model.frame(formula)	
    } else {# data argument in gamlss 
      Data <- get("gamlsscall", envir=gamlss.env)$data
    }
  } else {
    Data <- get("data", envir=gamlss.env)
  }
  Data <- data.frame(eval(substitute(Data)))
#-------------------------------------------------
  # new Daniil and Vlasis
  # Initialize bam
  formula <- as.formula(paste0("Y.var", deparse(formula, width.cutoff = 500L)))
  Data$Y.var = rep(0, nrow(Data))
    G = bam(formula, data=Data,
          offset=control$offset, method=control$method,
          control=control$control, select=control$select, 
          scale=control$scale, gamma=control$gamma, 
          knots=control$knots, sp=control$sp, min.sp=control$min.sp,
          paraPen=control$paraPen, chunk.size=control$chunk.size, 
          rho=control$rho, AR.start=control$AR.start, 
          discrete=control$discrete,
          cluster=control$cluster, nthreads=control$nthreads,
          gc.level=control$gc.level, use.chol=control$use.chol,
          samfrac=control$samfrac, 
          drop.unused.levels=control$bam$drop.unused.levels, 
          G=NULL, fit=FALSE, ...)
 #-------------------------------------------------- 
  xvar <- rep(0,  dim(Data)[1]) 
  attr(xvar,"formula")     <- formula
  attr(xvar,"control")     <- control
  attr(xvar, "gamlss.env") <- gamlss.env
  attr(xvar, "data")       <- as.data.frame(Data)
  attr(xvar, "call")       <- substitute(gamlss.ba(data[[scall]], z, w, ...)) 
  attr(xvar, "class")      <- "smooth"
  attr(xvar, "G")          <- G
  xvar
}

#--------------------------------------------------------------------------------------
ba.control = function(offset = NULL, 
                      method = "fREML", 
                 # optimizer = c("outer","newton"), 
                      select = FALSE,
                       scale = 0,
                       gamma = 1, 
                       knots = NULL,
                          sp = NULL, 
                      min.sp = NULL, 
                     paraPen = NULL,
                  chunk.size = 10000,
                         rho = 0,
                    AR.start = NULL,
                    discrete = FALSE,
                     cluster = NULL,
                    nthreads = NA,
                    gc.level = 1,
                    use.chol = FALSE,
                     samfrac = 1,
          drop.unused.levels = TRUE,
              drop.intercept = NULL, 
                               ...)  
{
  #gam()
  control <- gam.control(...)
  #ga()
  list( offset=offset, method=method, select=select, control=control,
        scale=scale, gamma=gamma, knots=knots, sp=sp, min.sp=min.sp,
        paraPen = paraPen, chunk.size = chunk.size,rho = rho,
        AR.start = AR.start,
        discrete=discrete,  cluster=cluster, nthreads=nthreads,
        gc.level=gc.level, use.chol=use.chol, samfrac=samfrac,
        drop.unused.levels = drop.unused.levels,
        drop.intercept=drop.intercept, ...)
}
#--------------------------------------------------------------------
#--------------------------------------------------------------------------------------
gamlss.ba <-function(x, y, w, xeval = NULL, ...) {
  if (is.null(xeval))
  {#fitting
             bam <- attr(x, "bam")
           Y.var <- y
           W.var <- w
               G <- attr(x,"G") 
             G$y <- Y.var
             G$w <- W.var
      G$mf$Y.var <- Y.var
G$mf$`(weights)` <- W.var
             fit <-  bam(G=G, fit=TRUE)
              df <- sum(fit$edf)-1 
              fv <- fitted(fit) 
       residuals <- y-fv
    list( fitted.values=fv, residuals=residuals,
         nl.df = df, lambda=fit$sp[1], #
         coefSmo = fit, var=NA)    # var=fv has to fixed
  } else { # predict 
      gamlss.env <- as.environment(attr(x, "gamlss.env"))
             obj <- get("object", envir=gamlss.env ) # get the object from predict
              TT <- get("TT", envir=gamlss.env ) # get wich position is now
              SL <- get("smooth.labels", envir=gamlss.env) # all the labels of the smoother
             fit <- eval(parse(text=paste("obj$", get("what", envir=gamlss.env), ".coefSmo[[",as.character(match(TT,SL)), "]]", sep="")))
           OData <- attr(x,"data") 
              ll <- dim(OData)[1]
            pred <- predict(fit,newdata = OData[seq(length(y)+1,ll),])
  }         
}

Try the gamlss.add package in your browser

Any scripts or data that you put into this service are public.

gamlss.add documentation built on May 30, 2017, 2:51 a.m.