R/nlAdditive.R

In gamlss.nl: Fitting non linear parametric GAMLSS models

# Non linear 
# This is an experimental non-linear additive function 
# created Tuesday, March 1, 2005 but finalized  December  2005 
# Mikis Stasinopoulos
#----------------------------------------------------------------------------------------
nl.obj <- function(formula, start, data)
   {
    nlform <- finterp(formula, .envir = substitute(data))
    # check that start and parameters in the formula have the same  length
    if (length(attr(nlform, "parameters"))!=length(start)) stop("the parameter or its starting values have the wrong length")
    attr(nlform,"start") <- start
    attr(nlform, "class") <- "nlobj"
    nlform
    }

#----------------------------------------------------------------------------------------
nl <-function(obj) #
{   
  scall <- deparse(sys.call()) 
    if (!is(obj, "nlobj")) stop("the object ", obj, " has to be generated by nl.obj()" )
      x <- eval(parse(text = attr(obj, "covariates")), envir=sys.frame(sys.parent(1)))
    len <- if(is.null(dim(x))) length(x) else dim(x)[1]
   xvar <- rep(0, len)
   attr(xvar, "nl.fun") <- obj
   real.call <- substitute(gamlss.nl(data[[scall]], z, w))
   attr(xvar, "call") <- real.call
   attr(xvar, "class") <- "smooth"
   xvar
}
#----------------------------------------------------------------------------------------
# the definition of the backfitting additive function
gamlss.nl <-function(x, y, w, xeval = NULL)
{
    nl.fun <- attr(x,"nl.fun")
     start <-  attr(attr(x,"nl.fun"),"start")
      optF <- function(p)  sum(w*(y-nl.fun(p))^2)
       fit <- nlm(optF,p=start, hessian=TRUE)
      coef <- fit$estimate
   #    eval(attr(attr(x,"nl.fun"),"start") <- coef, envir=sys.frame(sys.parent(1)))
   #   cat("s c",start, coef,"\n")
   hessian <- fit$hessian
     covar <- solve(hessian)
   staderr <- sqrt(diag(covar))
        fv <- nl.fun(coef)   
 residuals <- y-fv 
 list(fitted.values=fv, residuals=residuals,
     nl.df =length(coef), lambda=NA, 
     coefSmo = list(coef = coef, se = staderr, varcoeff = covar), var=fv )    # var=fv has to fixed
  if (is.null(xeval))
    {
   list(fitted.values=fv, residuals=residuals,
     nl.df =length(coef), lambda=NA, 
     coefSmo = list(coef = coef, se = staderr, varcoeff = covar), var=fv )    # var=fv has to fixed
    }
else 
    {
    stop("not prediction for the nl() function exist")
 #ndm <- as.matrix(attr(x,"design.matrix"))[seq(length(y)+1, 2),]
 #browser()
 #    ll <- dim(as.matrix(attr(x,"design.matrix")))[1]
 # nxvar <- as.matrix(attr(x,"design.matrix"))[seq(length(y)+1,ll),]
 #  pred <- drop(nxvar %*% coef(fit))
 #  pred
    }         

}