### de Castro, Cancho, and Rodrigues (2010),
### A hands-on approach for fitting long-term survival models  
###   under the GAMLSS framework
### Computer Methods and Programs in Biomedicine 97(2), 168-177.
### doi:10.1016/j.cmpb.2009.08.002

library(gamlss)
library(gamlss.cens)
library(survival)
library(boot)

## Link function (log shifted to -1)
# tau >= -1
own.linkfun = function(mu) log(mu + 1 + 1e-05)
own.linkinv = function(eta) -1 + pmax(.Machine$double.eps, exp(eta))
own.mu.eta = function(eta) pmax(.Machine$double.eps, exp(eta))
own.valideta = function(eta) TRUE   
make.link.gamlss("own")

## Weibull distribution
# Chen, Ibrahim & Sinha (1999), JASA 94, 909-919
# mu = gamma_1, sigma = gamma_2, mu > 0, and sigma real
dWEI4 = function (y, mu = 1, sigma = 0, log = FALSE) 
{
    if (any(mu <= 0)) 
        stop(paste("mu must be greater than 0  ", "\n", ""))
    if (any(y <= 0)) 
        stop(paste("y must be greater than 0 ", "\n", ""))
    sigma2 = exp(-sigma / mu)
    fy = dweibull(y, shape = mu, scale = sigma2, log = log)
    fy
}

pWEI4 = function (q, mu = 1, sigma = 0, lower.tail = TRUE, log.p = FALSE) 
{
    if (any(mu <= 0)) 
        stop(paste("mu must be greater than 0 ", "\n", ""))
    if (any(q <= 0)) 
        stop(paste("y must be greater than 0 ", "\n", ""))
    sigma2 = exp(-sigma / mu)
    cdf = pweibull(q, shape = mu, scale = sigma2, lower.tail = lower.tail, 
        log.p = log.p)
    cdf
}

## Improper p.d.f.
# Competing causes ~ negative binomial
# Time to event ~ Weibull 4 (WEI4)
# mu = gamma_1, sigma = gamma_2, mu > 0, and sigma real
dNBWEI4 = function (y, mu = 2, sigma = 3, nu = 0.5, tau = 1, log = FALSE) 
{
    if (any(mu <= 0)) 
        stop(paste("mu must be greater than 0 ", "\n", ""))
    if (any(nu < 0) | any(nu > 1)) 
        stop(paste("nu must be a probability ", "\n", ""))
    if (any(tau < -1)) 
        stop(paste("tau must be >= -1 ", "\n", ""))       
    if (any(y <= 0)) 
        stop(paste("y must be greater than 0 ", "\n", ""))
    ly = max(length(y), length(mu), length(sigma), length(nu), length(tau))
    y = rep(y, length = ly)
    sigma = rep(sigma, length = ly)
    mu = rep(mu, length = ly)
    nu = rep(nu, length = ly)
    tau = rep(tau, length = ly)

    vf = dWEI4(y, mu = mu, sigma = sigma) # p.d.f.
    vF = pWEI4(y, mu = mu, sigma = sigma) # c.d.f.
    logfy = ifelse(tau == 0, log(-log(nu) * nu^vF * vf),
      log((1 + (nu^(-tau) - 1) * vF)^(-1 / tau - 1) *
          (nu^(-tau) - 1) / tau * vf))
    if (log == FALSE) 
        fy = exp(logfy)
    else fy = logfy
    fy   
}

## Improper c.d.f.
# Competing causes ~ negative binomial
# Time to event ~ Weibull
# mu = gamma_1, sigma = gamma_2, mu > 0, and sigma real
pNBWEI4 = function (q, mu = 2, sigma = 3, nu = 0.5, tau = 1,
  lower.tail = TRUE,   log.p = FALSE)
{
    if (any(mu <= 0)) 
        stop(paste("mu must be greater than 0 ", "\n", ""))
    if (any(nu < 0) | any(nu > 1)) 
        stop(paste("nu must be a probability ", "\n", ""))
    if (any(tau < -1)) 
        stop(paste("tau must be >= -1 ", "\n", ""))       
    if (any(q <= 0)) 
        stop(paste("y must be greater than 0 ", "\n", ""))
    lq = max(length(q), length(mu), length(sigma), length(nu), length(tau))
    q = rep(q, length = lq)
    sigma = rep(sigma, length = lq)
    mu = rep(mu, length = lq)
    nu = rep(nu, length = lq)
    tau = rep(tau, length = lq)

    vF = pWEI4(q, mu = mu, sigma = sigma) # c.d.f.
    # cdf = 1 - Spop
    cdf = 1.0 - ifelse(tau == 0, nu^vF, (1 + (nu^(-tau) - 1) * vF)^(-1 / tau)) 

    if (lower.tail == TRUE) 
        cdf = cdf
    else cdf = 1 - cdf
    if (log.p == FALSE) 
        cdf = cdf
    else cdf = log(cdf)   
    cdf
}

## Long-term model
# Competing causes ~ negative binomial
# Time to event ~ Weibull
NBWEI4 = function (mu.link = "log", sigma.link = "identity", nu.link = "logit",
  tau.link = "own")
{
    mstats = checklink("mu.link", "Neg. binom. WEI4 model", 
        substitute(mu.link), c("log","identity", "inverse", "sqrt"))
    dstats = checklink("sigma.link", "Neg. binom. WEI4 model", 
        substitute(sigma.link), c("log", "identity"))
    vstats = checklink("nu.link", "Neg. binom. WEI4 model",
        substitute(nu.link), c("logit", "probit", "cloglog"))
    tstats = checklink("tau.link", "Neg. binom. WEI4 model",
        substitute(tau.link), c("log", "identity", "own"))   
    structure(
        list(family = c("NBWEI4", "Neg. binom. WEI4 model"), 
        parameters = list(mu = TRUE, sigma = TRUE, nu = TRUE, tau = TRUE),
        nopar = 4, 
        type = "Continuous",
        mu.link = as.character(substitute(mu.link)),
        sigma.link = as.character(substitute(sigma.link)),
        nu.link = as.character(substitute(nu.link)),
        tau.link = as.character(substitute(tau.link)),             
        mu.linkfun = mstats$linkfun,  
        sigma.linkfun = dstats$linkfun,
        nu.linkfun = vstats$linkfun,
        tau.linkfun = tstats$linkfun,              
        mu.linkinv = mstats$linkinv, 
        sigma.linkinv = dstats$linkinv,
        nu.linkinv = vstats$linkinv,
        tau.linkinv = tstats$linkinv,              
        mu.dr = mstats$mu.eta, 
        sigma.dr = dstats$mu.eta,
        nu.dr = vstats$mu.eta,
        tau.dr = tstats$mu.eta,

        dldm = function(y, mu, sigma, nu, tau) {
           ymuesig = y^mu * exp(sigma)
           dfdm = y^(mu-1) * exp(sigma - ymuesig) * (1 + mu * log(y) *
             (1 - ymuesig))
           dFdm = log(y) * y^mu * exp(sigma - ymuesig)
           dldm = ifelse(tau == 0, log(nu) * dFdm + dfdm / dWEI4(y, mu, sigma),
             -(1 / tau + 1) * (nu^(-tau) - 1) / (1 + (nu^(-tau) - 1) *
             pWEI4(y, mu, sigma)) * dFdm + dfdm / dWEI4(y, mu, sigma))
           dldm },              
        
        dldd = function(y, mu, sigma, nu, tau) {
           ymuesig = y^mu * exp(sigma)
           dFdd = y^mu * exp(sigma - ymuesig)
           dfdd = dFdd * mu * (1 - ymuesig) / y
           dldd = ifelse(tau == 0, log(nu) * dFdd + dfdd / dWEI4(y, mu, sigma),
             -(1 / tau + 1) * (nu^(-tau) - 1) / (1 + (nu^(-tau) - 1) *
             pWEI4(y, mu, sigma)) * dFdd + dfdd / dWEI4(y, mu, sigma))
           dldd },            
             
        dldv = function(y, mu, sigma, nu, tau) {
           dldv = ifelse(tau == 0, (1 / log(nu) + pWEI4(y, mu, sigma)) / nu,
              tau * nu^(-tau - 1) * (-1 / (nu^(-tau) - 1) +
             (1 / tau + 1) * pWEI4(y, mu, sigma) / (1 + nu^(-tau - 1) *
             pWEI4(y, mu, sigma))))
           dldv },
             
        dldt = function(y, mu, sigma, nu, tau) {
           dldt = ifelse(tau == 0, 0,
             (-nu^(-tau) * log(nu)) / (nu^(-tau) - 1) - 1 / tau +
             (1 / tau + 1) * nu^(-tau) * log(nu) * pWEI4(y, mu, sigma) /
             (1 + (nu^(-tau) - 1) * pWEI4(y, mu, sigma)))
           dldt },            

        d2ldm2 = function(y, mu, sigma, nu, tau) {
            nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
                log = TRUE), "mu", delta = 1e-04)
            dldm = as.vector(attr(nd, "gradient"))
            d2ldm2 = -dldm * dldm
            d2ldm2 },

        d2ldmdd = function(y, mu, sigma, nu, tau) {
            nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
                log = TRUE), "mu", delta = 1e-04)
            dldm = as.vector(attr(nd, "gradient"))
            nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
                log = TRUE), "sigma", delta = 1e-04)
            dldd = as.vector(attr(nd, "gradient"))           
            d2ldmdd = -dldm * dldd
            d2ldmdd },

         d2ldmdv = function(y, mu, sigma, nu, tau) {
            nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
                log = TRUE), "mu", delta = 1e-04)
            dldm = as.vector(attr(nd, "gradient"))
            dldv = ifelse(tau == 0, (1 / log(nu) + pWEI4(y, mu, sigma)) / nu,
              tau * nu^(-tau - 1) * (-1 / (nu^(-tau) - 1) +
             (1 / tau + 1) * pWEI4(y, mu, sigma) / (1 + nu^(-tau - 1) *
             pWEI4(y, mu, sigma))))            
            d2ldmdv = -dldm * dldv
            d2ldmdv },

         d2ldmdt = function(y, mu, sigma, nu, tau) {
           nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
                log = TRUE), "mu", delta = 1e-04)
           dldm = as.vector(attr(nd, "gradient"))
           dldt = ifelse(tau == 0, 0,
             (-nu^(-tau) * log(nu)) / (nu^(-tau) - 1) - 1 / tau +
             (1 / tau + 1) * nu^(-tau) * log(nu) * pWEI4(y, mu, sigma) /
             (1 + (nu^(-tau) - 1) * pWEI4(y, mu, sigma)))           
           d2ldmdt = -dldm * dldt
           d2ldmdt },

         d2ldd2 = function(y, mu, sigma, nu, tau) {
            nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
               log = TRUE), "sigma", delta = 1e-04)
            dldd = as.vector(attr(nd, "gradient"))
            d2ldd2 = -dldd * dldd
            d2ldd2 },

         d2ldddv = function(y, mu, sigma, nu, tau) {
            nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
               log = TRUE), "sigma", delta = 1e-04)
            dldd = as.vector(attr(nd, "gradient"))
            dldv = ifelse(tau == 0, (1 / log(nu) + pWEI4(y, mu, sigma)) / nu,
              tau * nu^(-tau - 1) * (-1 / (nu^(-tau) - 1) +
             (1 / tau + 1) * pWEI4(y, mu, sigma) / (1 + nu^(-tau - 1) *
             pWEI4(y, mu, sigma))))                
            d2ldddv = -dldd * dldv
            d2ldddv },

         d2ldddt = function(y, mu, sigma, nu, tau) {
           nd = gamlss:::numeric.deriv(dNBWEI4(y, mu, sigma, nu, tau, 
              log = TRUE), "sigma", delta = 1e-04)
           dldd = as.vector(attr(nd, "gradient"))
           dldt = ifelse(tau == 0, 0,
             (-nu^(-tau) * log(nu)) / (nu^(-tau) - 1) - 1 / tau +
             (1 / tau + 1) * nu^(-tau) * log(nu) * pWEI4(y, mu, sigma) /
             (1 + (nu^(-tau) - 1) * pWEI4(y, mu, sigma)))           
           d2ldddt = -dldd * dldt
           d2ldddt },

        d2ldv2 = function(y, mu, sigma, nu, tau) {
          dldv = ifelse(tau == 0, (1 / log(nu) + pWEI4(y, mu, sigma)) / nu,
              tau * nu^(-tau - 1) * (-1 / (nu^(-tau) - 1) +
             (1 / tau + 1) * pWEI4(y, mu, sigma) / (1 + nu^(-tau - 1) *
             pWEI4(y, mu, sigma))))           
          d2ldv2 = -dldv * dldv
          d2ldv2 },

        d2ldvdt = function(y, mu, sigma, nu, tau) {
          dldv = ifelse(tau == 0, (1 / log(nu) + pWEI4(y, mu, sigma)) / nu,
              tau * nu^(-tau - 1) * (-1 / (nu^(-tau) - 1) +
             (1 / tau + 1) * pWEI4(y, mu, sigma) / (1 + nu^(-tau - 1) *
             pWEI4(y, mu, sigma))))             
          dldt = ifelse(tau == 0, 0,
             (-nu^(-tau) * log(nu)) / (nu^(-tau) - 1) - 1 / tau +
             (1 / tau + 1) * nu^(-tau) * log(nu) * pWEI4(y, mu, sigma) /
             (1 + (nu^(-tau) - 1) * pWEI4(y, mu, sigma)))                     
          d2ldvdt = -dldv * dldt
          d2ldvdt },           

        d2ldt2 = function(y, mu, sigma, nu, tau) {
          dldt = ifelse(tau == 0, 0,
             (-nu^(-tau) * log(nu)) / (nu^(-tau) - 1) - 1 / tau +
             (1 / tau + 1) * nu^(-tau) * log(nu) * pWEI4(y, mu, sigma) /
             (1 + (nu^(-tau) - 1) * pWEI4(y, mu, sigma)))                     
          d2ldt2 = -dldt * dldt
          d2ldt2 }, 

        G.dev.incr = function(y, mu, sigma, nu, tau,...) -2 * dNBWEI4(y,
          mu = mu, sigma = sigma, nu = nu, tau = tau, log = TRUE),
        rqres = expression(rqres(pfun = "pNBWEI4", type = "Continuous",
            ymin = 0, y = y, mu = mu, sigma = sigma, nu = nu, tau = tau)), 
        mu.initial = expression(mu = rep(1, length(y))),
        sigma.initial = expression(sigma = rep(-log(mean(y)), length(y))), 
        nu.initial = expression(nu = rep(0.5, length(y))),
        tau.initial = expression(tau = rep(0.5, length(y))),             
        mu.valid = function(mu) all(mu > 0),
        sigma.valid = function(sigma) TRUE,
        nu.valid = function(nu)  all(nu > 0) & all(nu < 1),
        tau.valid = function(tau) all(tau >= -1),             
        y.valid = function(y) all(y > 0)), class = c("gamlss.family", 
        "family"))
}

## Improper p.d.f.
# Competing causes ~ Poisson
# Time to event ~ Weibull 4 (WEI4)
# mu = gamma_1, sigma = gamma_2, mu > 0, and sigma real
dPOWEI4 = function (y, mu = 1, sigma = 3, nu = 0.5, log = FALSE) 
{
    if (any(mu <= 0)) 
        stop(paste("mu must be greater than 0 ", "\n", ""))
    if (any(nu < 0) | any(nu > 1)) 
        stop(paste("nu must be a probability ", "\n", ""))
    if (any(y <= 0)) 
        stop(paste("y must be greater than 0 ", "\n", ""))
    ly = max(length(y), length(mu), length(sigma), length(nu))
    y = rep(y, length = ly)
    sigma = rep(sigma, length = ly)
    mu = rep(mu, length = ly)
    nu = rep(nu, length = ly)

    vf = dWEI4(y, mu = mu, sigma = sigma) # p.d.f.
    vF = pWEI4(y, mu = mu, sigma = sigma) # c.d.f.   
    logfy = log(-log(nu) * nu^vF * vf)
    if (log == FALSE) 
        fy = exp(logfy)
    else fy = logfy
    fy   
}

## Improper c.d.f.
# Competing causes ~ negative binomial
# Time to event ~ Weibull
# mu = gamma_1, sigma = gamma_2, mu > 0, and sigma real
pPOWEI4 = function (q, mu = 2, sigma = 3, nu = 0.5, lower.tail = TRUE,
  log.p = FALSE) 
{
    if (any(mu <= 0)) 
        stop(paste("mu must be greater than 0 ", "\n", ""))
    if (any(nu < 0) | any(nu > 1)) 
        stop(paste("nu must be a probability ", "\n", ""))
    if (any(q <= 0)) 
        stop(paste("y must be greater than 0 ", "\n", ""))
    lq = max(length(q), length(mu), length(sigma), length(nu))
    q = rep(q, length = lq)
    sigma = rep(sigma, length = lq)
    mu = rep(mu, length = lq)
    nu = rep(nu, length = lq)

    vF = pWEI4(q, mu = mu, sigma = sigma) # c.d.f.       
    cdf = 1 - nu^vF # = 1 - Spop

    if (lower.tail == TRUE) 
        cdf = cdf
    else cdf = 1 - cdf
    if (log.p == FALSE) 
        cdf = cdf
    else cdf = log(cdf)   
    cdf
}

## Long-term model
# Competing causes ~ Poisson
# Time to event ~ Weibull
POWEI4 = function (mu.link = "log", sigma.link = "identity", nu.link = "logit")
{
    mstats = checklink("mu.link", "Promotion time cure model", 
        substitute(mu.link), c("log","identity", "inverse", "sqrt"))
    dstats = checklink("sigma.link", "Promotion time cure model", 
        substitute(sigma.link), c("log", "identity"))
    vstats = checklink("nu.link", "Promotion time cure model",
        substitute(nu.link), c("logit", "probit", "cloglog"))
    structure(
        list(family = c("POWEI4", "Promotion time cure model"), 
        parameters = list(mu = TRUE, sigma = TRUE, nu = TRUE),
        nopar = 3, 
        type = "Continuous",
        mu.link = as.character(substitute(mu.link)),
        sigma.link = as.character(substitute(sigma.link)),
        nu.link = as.character(substitute(nu.link)),
        mu.linkfun = mstats$linkfun,  
        sigma.linkfun = dstats$linkfun,
        nu.linkfun = vstats$linkfun,
        mu.linkinv = mstats$linkinv, 
        sigma.linkinv = dstats$linkinv,
        nu.linkinv = vstats$linkinv,
        mu.dr = mstats$mu.eta, 
        sigma.dr = dstats$mu.eta,
        nu.dr = vstats$mu.eta,
             
        dldm = function(y, mu, sigma, nu) {
           ymuesig = y^mu * exp(sigma)
           dfdm = y^(mu-1) * exp(sigma - ymuesig) * (1 + mu * log(y) *
             (1 - ymuesig))
           dFdm = log(y) * y^mu * exp(sigma - ymuesig)
           dldm = log(nu) * dFdm + dfdm / dWEI4(y, mu, sigma)           
           dldm },     
        
        dldd = function(y, mu, sigma, nu) {
           ymuesig = y^mu * exp(sigma)
           dFdd = y^mu * exp(sigma - ymuesig)
           dfdd = dFdd * mu * (1 - ymuesig) / y
           dldd = log(nu) * dFdd + dfdd / dWEI4(y, mu, sigma)             
           dldd },     
             
        dldv = function(y, mu, sigma, nu) {
           dldv = (1 / log(nu) + pWEI4(y, mu, sigma)) / nu
           dldv },
      
        d2ldm2 = function(y, mu, sigma, nu) {
            nd = gamlss:::numeric.deriv(dPOWEI4(y, mu, sigma, nu,
                log = TRUE), "mu", delta = 1e-04)
            dldm = as.vector(attr(nd, "gradient"))
            d2ldm2 = -dldm * dldm
            d2ldm2 },
             
        d2ldmdd = function(y, mu, sigma, nu) {
            nd = gamlss:::numeric.deriv(dPOWEI4(y, mu, sigma, nu, 
                log = TRUE), "mu", delta = 1e-04)
            dldm = as.vector(attr(nd, "gradient"))
            nd = gamlss:::numeric.deriv(dPOWEI4(y, mu, sigma, nu, 
                log = TRUE), "sigma", delta = 1e-04)
            dldd = as.vector(attr(nd, "gradient"))           
            d2ldmdd = -dldm * dldd
            d2ldmdd },

        d2ldmdv = function(y, mu, sigma, nu) {
            nd = gamlss:::numeric.deriv(dPOWEI4(y, mu, sigma, nu, 
                log = TRUE), "mu", delta = 1e-04)
            dldm = as.vector(attr(nd, "gradient"))
            dldv = (1 / log(nu) + pWEI4(y, mu, sigma)) / nu            
            d2ldmdv = -dldm * dldv
          d2ldmdv },

        d2ldd2 = function(y, mu, sigma, nu) {
            nd = gamlss:::numeric.deriv(dPOWEI4(y, mu, sigma, nu, 
               log = TRUE), "sigma", delta = 1e-04)
            dldd = as.vector(attr(nd, "gradient"))
            d2ldd2 = -dldd * dldd
            d2ldd2 },

        d2ldddv = function(y, mu, sigma, nu) {
            nd = gamlss:::numeric.deriv(dPOWEI4(y, mu, sigma, nu, 
               log = TRUE), "sigma", delta = 1e-04)
            dldd = as.vector(attr(nd, "gradient"))
            dldv = (1 / log(nu) + pWEI4(y, mu, sigma)) / nu            
            d2ldddv = -dldd * dldv
            d2ldddv },

        d2ldv2 = function(y, mu, sigma, nu) {
            dldv = (1 / log(nu) + pWEI4(y, mu, sigma)) / nu          
            d2ldv2 = -dldv * dldv
            d2ldv2 },
        G.dev.incr = function(y, mu, sigma, nu, ...)
          -2 * dPOWEI4(y, mu = mu, sigma = sigma, nu = nu, log = TRUE),
        rqres = expression(rqres(pfun = "pPOWEI4", type = "Continuous",
            ymin = 0, y = y, mu = mu, sigma = sigma, nu = nu)), 
        mu.initial = expression(mu = rep(1, length(y))),
        sigma.initial = expression(sigma = rep(-log(mean(y)), length(y))),    
        nu.initial = expression(nu = rep(0.5, length(y))),
        mu.valid = function(mu) all(mu > 0),
        sigma.valid = function(sigma) TRUE,
        nu.valid = function(nu)  all(nu > 0) & all(nu < 1),
        y.valid = function(y) all(y > 0)), class = c("gamlss.family", 
        "family"))
}

## Profile deviance
# Function kindly provided by an anonymous referee
prof.term <-
function (model = NULL, criterion = "GD", penalty = 2.5, other = NULL, 
    min = NULL, max = NULL, step = NULL, type = "o", xlabel = NULL, 
    plot = TRUE, term = TRUE, perc = 95, ...) 
{
    if (is.null(model)) 
        stop("you have not defined the model")
    if (is.null(min)) 
        stop("you have not defined the minimum")
    if (is.null(max)) 
        stop("you have not defined the maximum")
    if (is.null(step)) 
        stop("you have not defined the step")
    if (!criterion %in% c("IC", "GD")) 
        stop("criterion should be IC or GD")
    interval <- seq(from = min, to = max, by = step)
    I.C <- rep(0, length(interval))
    call <- model
    if (!is.null(model$data)) {
        a <- model$data
        attach(eval(substitute(a)))
        on.exit(detach(eval(substitute(a))))
    }
    for (i in length(interval):1) {
        this <<- this <- interval[i]
        if (!is.null(other)) 
            eval(other)
        mod.1 <- eval(call)
        call <- mod.1$call
        call$mu.start <- fitted(mod.1, "mu")
        if ("sigma" %in% mod.1$parameters) 
            call$sigma.start <- fitted(mod.1, "sigma")
        if ("nu" %in% mod.1$parameters) 
            call$nu.start <- fitted(mod.1, "nu")
        if ("tau" %in% mod.1$parameters) 
            call$tau.start <- fitted(mod.1, "tau")
        I.C[i] <- if (criterion == "GD") 
            deviance(mod.1)
        else IC(mod.1, penalty)
    }
    xlab <- if (!is.null(xlabel)) 
        xlabel
    else "parameter"
    ylab <- if (criterion == "GD") 
        "Global Deviances"
    else paste("GAIC pen=", penalty)
    main <- if (criterion == "GD") 
        "Profile Global Deviance"
    else "Profile GAIC"
    prof.out <- cbind(interval, I.C)
    if (criterion == "GD") {
        Gmin <- min(I.C)
        mx <- which.min(I.C)
        min.parameter <- interval[mx]
        lim <- Gmin + qchisq((perc/100), 1)
        xl <- as.vector(interval)
        m <- length(interval)
    }
    if (plot) {
        plot(interval, I.C, xlab = xlab, ylab = ylab, main = main, 
            col = "darkgreen", frame.plot = TRUE, type = type)
        plims <- par("usr")
        if (criterion == "GD" & term == TRUE) {
            if (lim < max(I.C)) {
                abline(h = lim, lty = 3)
                y0 <- plims[3]
                scal <- (1/10 * (plims[4] - y0))/par("pin")[2]
                scx <- (2/10 * (plims[2] - plims[1]))/par("pin")[1]
                text(xl[1] + scx, lim + scal, paste(perc, "%"))
                la <- xl[mx]
                if (mx > 1 && mx < m) 
                  segments(la, y0, la, Gmin, lty = 3)
            }
            ind <- range((1:m)[I.C < lim])
            if (I.C[1] > lim) {
                i <- ind[1]
                x <- xl[i - 1] + ((lim - I.C[i - 1]) * (xl[i] - 
                  xl[i - 1]))/(I.C[i] - I.C[i - 1])
                min.ci.par <- x
                segments(x, y0, x, lim, lty = 3)
            }
            if (I.C[m] > lim) {
                i <- ind[2] + 1
                x <- xl[i - 1] + ((lim - I.C[i - 1]) * (xl[i] - 
                  xl[i - 1]))/(I.C[i] - I.C[i - 1])
                max.ci.par <- x
                segments(x, y0, x, lim, lty = 3)
            }
            cat("*******************************************************************", 
                "\n")
            cat("Best estimate of the fixed parameter is ", min.parameter, 
                "\n")
            cat("with a Global Deviance equal to ", Gmin, " at position ", 
                mx, "\n")
            if ((I.C[1] > lim) && (I.C[m] > lim)) {
                cat("A ", perc, "% Confidence interval is: (", 
                  min.ci.par, ",", max.ci.par, ") \n")
            }
            cat("*******************************************************************", 
                "\n")
        }
        return(list(prof.out = invisible(prof.out), ci = c(min.ci.par, max.ci.par)))
    }
    else {
        return(list(prof.out = prof.out, ci = c(min.ci.par, max.ci.par)))
    }
}