.packageName <- "mfp"
cox <- function() {
#
# Version 1.1.0         20Sep02
#
        structure(list(family = "Cox"), class = "family")
}
fp <- function(x, df = 4, select = NA, alpha = NA, scale=TRUE)
{
# Version 1.0.1	09 dec 03
#
    name <- deparse(substitute(x))
    attr(x, "df") <- df
    attr(x, "alpha") <- alpha
    attr(x, "select") <- select
    attr(x, "scale") <- scale
    attr(x, "name") <- name
    x
}
fp.fit <- function(X, Y, df, dfr, cox, gauss, shift, scale, ...)
{
#
# Version 1.1.1     09 dec 03
#
    X <- X[,apply(X, 2, function(x) !all(x==0)), drop=FALSE] # to avoid numerical problems
    x <- X[, 1] # x is first column of X                 # induced by fp.gen if df=0
    Xcov <- X[, -1, drop=FALSE]   # covariates
    ncov <- ncol(Xcov)
    nobs <- nrow(X)
    pwrs <- c(-2, -1, -0.5, 0, 0.5, 1, 2, 3)
    npwrs <- length(pwrs)
    dev2 <- dev4 <- 1000000000000000    
#
# Set up fitters
#
    if(cox) {
        if(exists("coxph.fit")) fitter <- get("coxph.fit")
        else fitter <- getFromNamespace("coxph.fit","survival")
        dv <- "loglik"
    }
    else {
        fitter <- get("glm.fit")
        dv <- "deviance"
    }
    n <- 1 + cox
#
# Null and linear models
#
    dev <- (-2)^cox * fitter(X, Y, ...)[[dv]]
    dev1 <- dev[n]
    if(!ncov)
        dev0 <- dev[1]
    else dev0 <- (-2)^cox * fitter(Xcov, Y, ...)[[dv]][n]
    if(df > 1) {
        for(i in 1:npwrs) {
#
# Find best single power transformation
#
            x.fp <- fp.gen(x, pwrs[i], shift, scale)
            dev <- (-2)^cox * fitter(cbind(x.fp, Xcov), Y, ...)[[dv
                ]][n]
            if(!is.null(dev) & dev < dev2) {
                dev2 <- dev
                pwr2 <- pwrs[i]
            }
            if(df == 4) {
#
# Find best two power transformation
#
                x.fp <- fp.gen(x, c(pwrs[i], pwrs[i]), shift, 
                  scale)
                dev <- (-2)^cox * fitter(cbind(x.fp, Xcov), Y, 
                  ...)[[dv]][n]
                if(!is.null(dev) & dev < dev4) {
                  dev4 <- dev
                  pwr4 <- c(pwrs[i], pwrs[i])
                }
                j <- i + 1
                while(j <= npwrs) {
                  x.fp <- fp.gen(x, c(pwrs[i], pwrs[j]), shift, scale)
                  dev <- (-2)^cox * fitter(cbind(x.fp, Xcov), Y,
                    ...)[[dv]][n]
                  if(!is.null(dev) & dev < dev4) {
                    dev4 <- dev
                    pwr4 <- c(pwrs[i], pwrs[j])
                  }
                  j <- j + 1
                }
            }
        }
    }
#
# Output
#   
    dev <- c(dev0, dev1, dev2, dev4)
    if(df < 4)
        dev[4] <- pwr4 <- NA
    if(df < 2)
        dev[3] <- pwr2 <- NA
    if(gauss)
        dev <- nobs * (1 + log((2 * pi * dev)/nobs))
    fit <- list(pwr4 = pwr4, pwr2 = pwr2, dev4 = dev[4], dev2 = dev[3], 
        dev1 = dev[2], dev0 = dev[1], nobs = nobs, dfr = dfr, df = df, 
        gauss = gauss)
    return(fit)
}
fp.gen <- function(x, pwrs, shift = NULL, scale = NULL)
{
#
#   Version 1.0.0   19 Oct 1999
#
    nobs <- length(x)
    pwr1 <- pwrs[1]
    pwr2 <- pwrs[2]
    X <- matrix(0, nrow = nobs, ncol = 2)   
#
# Deal with the first power & sort out rescaling
#
    if(!is.na(pwr1)) {
        if(is.null(scale) | is.null(shift)) {
            x.transform <- fp.scale(x, scaling=TRUE)
            shift <- x.transform$shift
            scale <- x.transform$scale
            cat("[Pre-transformation : x->(x+", shift, ")/", scale, 
                "]\n", sep = "")
            }
        x <- x + shift
        x <- x/scale
        x1 <- ifelse(pwr1 != rep(0, nobs), x^pwr1, log(x))
        X[, 1] <- x1
    }
#
# Other power
#
    if(!is.na(pwr2)) {
        if(pwr2 == pwr1)
            x2 <- log(x) * x1
        else x2 <- ifelse(pwr2 != rep(0, nobs), x^pwr2, log(x))
        X[, 2] <- x2
    }
    return(X)
}
fp.order <- function(x, y, cox, gauss, xnames, ...)
{
#
# Version 1.1.0     30.12.03
#
# Returns ordering of xvars (not intercept!)
# by Wald test ranking (Cox: by LR test)
#
    int <- as.numeric(!cox)
    nx <- ncol(x)
    if(cox) {
        if(exists("coxph.fit")) fitter <- get("coxph.fit")
        else fitter <- getFromNamespace("coxph.fit","survival")
        fit <- fitter(x, y, ...)
        se <- sqrt(diag(fit$var))
        dev <- -2 * fit$loglik
        p.value <- vector(length=length(xnames))
        for(i in unique(xnames)) {
            ld <- sum(xnames==i)
            ll <- fitter(x[,xnames!=i, drop=FALSE], y, ...)$loglik[2]
            p.value[xnames==i] <- pchisq(2*(fit$loglik[2]-ll), df=ld)
        }
        x.order <- order(p.value[(1 + int):nx])
    }
    else {
        fit <- glm.fit(x, y, ...)
        Iinv <- solve(t(fit$R) %*% fit$R)
        se <- sqrt(fit$dev/fit$df.residual * diag(Iinv))
        dev <- c(fit$null.deviance, fit$deviance)
t.value <- abs(fit$coef/se)
x.order <- rev(order(t.value[(1 + int):nx]))
    }
    if(gauss) {
        nobs <- nrow(x)
        dev <- nobs * (1 + log((2 * pi * dev)/nobs))
    }
    return(list(order = x.order, dev = dev))
}
fp.out <- function(..., pos)
{
    #
    # Version 1.0.1    28Aug01
    #
    x <- list(...)
    nx <- length(x)
    tabs <- 0
    cat("\n")
    for(i in 1:nx) {
        wordi <- x[[i]]
        if(pos[i] > tabs) {
            cat(rep("\t", (pos[i] - tabs)))
            tabs <- pos[i]
        }
        wordi <- wordi[!is.na(wordi)]
        # deal with powers 
        if(length(wordi) > 0) {
            cat(wordi)
            if(length(wordi) > 1) {
                nwordi <- 1 + sum(nchar(wordi))
            }
            else nwordi <- nchar(wordi)
            if(nwordi > 7) {
                tabs <- tabs + 1
            }
        }
        else cat("\t")
    }
    invisible()
}
fp.scale <- function(x, scaling)
{
#
#   Version 1.0.0   19 Oct 1999
#
scale <- 1; shift <- 0
  if(scaling) {
      if(min(x) <= 0) {
        z <- sort(x)[-1] - sort(x)[ - length(x)]
        shift <- min(z[z > 0]) - min(x)
    }
    else shift <- 0
    range <- max(x) - min(x)
    scale <- 10^(sign(log10(range)) * trunc(abs(log10(range))))
  }
    return(list(shift=shift, scale=scale))
}
fp.sel <- function(fit, alpha = 0.05, select = 1)
{
#
# Version 1.0.0     6Oct99
#
# Calculate deviance differences & p-values
#
    dfr <- fit$dfr - fit$df # residual df
    dd.null <- fit$dev0 - min(fit$dev1, fit$dev2, fit$dev4, na.rm=TRUE)
    if(fit$gauss) {
        f.null <- ((exp(dd.null/fit$nobs) - 1) * dfr)/fit$df
        p.null <- 1 - pf(f.null, fit$df, dfr)
    }
    else p.null <- 1 - pchisq(dd.null, fit$df)
    if(fit$df > 1) {
        dd.lin <- fit$dev1 - min(fit$dev2, fit$dev4, na.rm=TRUE)
        if(fit$gauss) {
            f.lin <- ((exp(dd.lin/fit$nobs) - 1) * dfr)/(fit$df - 1
                )
            p.lin <- 1 - pf(f.lin, fit$df - 1, dfr)
        }
        else p.lin <- 1 - pchisq(dd.lin, fit$df - 1)
        if(fit$df > 2) {
            dd.FP <- fit$dev2 - fit$dev4
            if(fit$gauss) {
                f.FP <- ((exp(dd.FP/fit$nobs) - 1) * dfr)/2
                p.FP <- 1 - pf(f.FP, 2, dfr)
            }
            else p.FP <- 1 - pchisq(dd.FP, 2)
        }
        else p.FP <- NA
    }
    else {
        p.lin <- NA
        p.FP <- NA
    }
    if(p.null > select) {
        df <- 0
        pwrs <- c(NA, NA)
        dev <- fit$dev0
    }
    else {
        if(fit$df > 1) {
            if(p.lin > alpha)
                df <- 1
            else {
                if(fit$df > 2) {
                  if(p.FP > alpha)
                    df <- 2
                  else {
                    df <- 4
                    pwrs <- fit$pwr4
                    dev <- fit$dev4
                  }
                }
                else df <- 2
            }
        }
        else df <- 1
    }
    if(df == 1) {
        pwrs <- c(1, NA)
        dev <- fit$dev1
    }
    if(df == 2) {
        pwrs <- c(fit$pwr2, NA)
        dev <- fit$dev2
    }
    results <- list(p.null = p.null, p.lin = p.lin, p.FP = p.FP, df = df, 
        pwrs = pwrs, dev = dev)
    return(list(results=results, fit=fit))
}
mfp <- function(formula = formula(data), data = parent.frame(), family = gaussian,
     subset, na.action, init, alpha = 0.05, select = 1, verbose = FALSE, x = TRUE, y = TRUE)
{
# Version 1.2.0     26072004
#
    call <- match.call()
    if (is.character(family))
        family <- get(family, mode="function")
    if (is.function(family))
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("`family' not recognized")
    }
cox <- (family$family == "Cox")
if(cox) require(survival)
m <- match.call(expand = FALSE)
temp <- c("", "formula", "data", "weights", "subset", "na.action")
m <- m[match(temp, names(m), nomatch = 0)]
special <- c("strata", "fp")
Terms <- if (missing(data))
        terms(formula, special)
else terms(formula, special, data = data)
    m$formula <- Terms
    m$alpha <- m$select <- m$scale <- m$family <- m$verbose <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval(m, sys.parent())
    Y <- model.extract(m, "response") #
    weights <- model.extract(m, "weights")
#
    offset <- attr(Terms, "offset")
    tt <- length(offset)
    offset <- if (tt == 0 & cox)
        rep(0, nrow(Y))
    else if (tt == 0 & !cox)
        rep(0, length(Y))
    else if (tt == 1)
        m[[offset]]
    else {
        ff <- m[[offset[1]]]
        for (i in 2:tt) ff <- ff + m[[offset[i]]]
        ff
    }
    attr(Terms, "intercept") <- 1
    dropx <- NULL
    strats <- attr(Terms, "specials")$strata
    if (length(strats)) {
        temp <- untangle.specials(Terms, "strata", 1)
        dropx <- temp$terms
        if (length(temp$vars) == 1)
            strata.keep <- m[[temp$vars]]
        else strata.keep <- strata(m[, temp$vars], shortlabel=TRUE)
        strats <- as.numeric(strata.keep)
    }
    if (length(dropx)) # dropx is number of stratification variables
        newTerms <- Terms[-dropx]  # newTerms are Terms w/o stratif. var
    else newTerms <- Terms
    X <- model.matrix(newTerms, m)
    if (missing(init))
        init <- NULL
#
# Set up lists for fitter
#
    nx <- ncol(X) - 1
    nobs <- nrow(X)
    df.list <- rep(1, nx)
    scale.list <- rep(FALSE, nx)
    alpha.list <- rep(alpha, nx)
    select.list <- rep(select, nx)
#
# Deal with FP terms
#
    fp.pos <- grep("fp", dimnames(X)[[2]])
    fp.mpos <- attributes(Terms)$specials$fp
    fp.xpos <- unlist(attributes(Terms)$specials) - 1
    if(length(fp.pos) > 0) {
        fp.pos <- fp.pos - 1    # without intercept
        fp.data <- m[, fp.mpos, drop=FALSE]
        df.list[fp.pos] <- unlist(lapply(fp.data, attr, "df"))
        scale.list[fp.pos] <- unlist(lapply(fp.data, attr, "scale"))
        alpha.list[fp.pos] <- unlist(lapply(fp.data, attr, "alpha"))
        alpha.list[sapply(alpha.list, is.na)] <- alpha
        select.list[fp.pos] <- unlist(lapply(fp.data, attr, "select"))
        select.list[sapply(select.list, is.na)] <- select
        names <- dimnames(X)[[2]]
        names[fp.pos + 1] <- unlist(lapply(fp.data, attr, "name"))
        xnames <- names[-1]
        xnames[-fp.pos] <- attr(Terms,"term.labels")[-fp.xpos]
        dimnames(X)[[2]] <- names
    }
    unlist(lapply(m, attr, "name"))
#
# Cox fit
#
    if(cox) {
        if(!inherits(Y, "Surv"))
            stop("Response must be a survival object")
        type <- attr(Y, "type")
        if(type != "right")
            stop("The data must be right censored")
        X <- X[, -1, drop=FALSE]  # remove intercept
        control <- coxph.control()
        method <- "efron"
        fit <- mfp.fit(X, Y, TRUE, FALSE, df.list, scale.list, alpha.list, select.list,
                  verbose = verbose, strata=strats, offset=offset, init, control,
                  weights = weights, method = method, rownames = row.names(m),
                  xnames = xnames)
        if(is.character(fit)) {
            fit <- list(fail = fit)
            attr(fit, "class") <- c("mfp", "coxph")
        }
        else attr(fit, "class") <- c("mfp", fit$method)
        fit$n <- nobs
        fit$family <- family
    }
    else {
#
# GLM fit
#
        gauss <- (family$family == "gaussian")
        fit <- mfp.fit(X, Y, FALSE, gauss, df.list, scale.list, alpha.list, select.list,
                  verbose = verbose, family = family, xnames = xnames)
        attr(fit, "class") <- c("mfp", "glm", "lm")
    }
# glm add-on (to derive var)
   if(!cox) {
        dispersion <-
        if (any(fit$family$family == c("poisson","binomial")))
            1
        else if (fit$df.residual > 0) {
            if (any(fit$weights == 0))
                warning("observations with zero weight ", "not used for calculating dispersion")
            sum(fit$weights * fit$residuals^2)/fit$df.residual
        } else Inf
   p <- fit$rank
    if (p > 0) {
        p1 <- 1:p
        Qr <- fit$qr
        aliased <- is.na(coef(fit))
        coef.p <- fit$coefficients[Qr$pivot[p1]]
        covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
        dimnames(covmat.unscaled) <- list(names(coef.p), names(coef.p))
        fit$var <- dispersion * covmat.unscaled
    }
   }
# back-transformation of coefs and vars
  if(length(fit$df.final)>1) {
   i2 <- if(cox) 0 else 1
   tp <- if(cox) double() else 1
   for(i in seq(fit$df.final)) {
    if(fit$df.final[i]!=0) {
    p <- 1; if(fit$df.final[i]==4) p <- 2
    i1 <- i2+1; i2 <- i2+p
    tp <- c(tp, fit$scale[i,2]^fit$powers[i,1:p])
    fit$coefficients[i1:i2] <- fit$coefficients[i1:i2]/(fit$scale[i,2]^fit$powers[i,1:p])  # coefs
    }
   }
  if(length(tp)>0) fit$var <- fit$var/(tp%*%t(tp))  # vars
  } else {
    if(fit$df.final!=0) {
      p <- 1; if(fit$df.final==4) p <- 2
      if(cox) {
       tp <- as.vector(fit$scale[2]^fit$powers[1:p])
       fit$coefficients <- fit$coefficients/(fit$scale[2]^fit$powers[1:p])  # coefs
       fit$var <- fit$var/(tp%*%t(tp))  # vars
      } else {
       tp <- c(1,as.vector(fit$scale[2]^fit$powers[1:p]))
       fit$coefficients[-1] <- fit$coefficients[-1]/(fit$scale[2]^fit$powers[1:p])  # coefs
       fit$var <- fit$var/(tp%*%t(tp))  # vars
      }
    }
  }
# Cox add-on (needed for summary.coxph)
    if(cox) {
      if(length(fit$coefficients) && is.null(fit$wald.test)) {
        nabeta <- !is.na(fit$coefficients)
        if(is.null(init))
          temp <- fit$coefficients[nabeta]
        else temp <- (fit$coefficients - init)[nabeta]
      }
    if(exists("coxph.wtest")) tester <- get("coxph.wtest")
    else tester <- getFromNamespace("coxph.wtest","survival")

    fit$wald.test <- tester(fit$var[nabeta, nabeta],
                temp, .Machine$double.eps^0.75)$test
    }
#
    if (x) fit$X <- X
    if (y) fit$y <- Y
    fit$terms <- Terms
    fit$call <- call
    fit
}
mfp.fit <- function(x, y, cox, gauss, df, scaling, alpha, select, verbose = TRUE, xnames = NULL, ...)
{
#
# Version 1.2.0     27072004
#
    int <- as.numeric(!cox) # intercept
    nx <- ncol(x) - int
    nobs <- nrow(x)
    x.names <- dimnames(x)[[2]][int + seq(nx)]
    if(sum((df == 1 | df == 2 | df == 4), na.rm=TRUE) != nx)
        stop("df is invalid")
    if(sum((alpha > 0 & alpha <= 1), na.rm=TRUE) != nx)
        stop("alpha is invalid")
    if(sum((select > 0 & select <= 1), na.rm=TRUE) != nx) stop(
            "select is invalid")    #
#
# Order variables
#
    x.order <- fp.order(x, y, cox, gauss, xnames, ...)
    x <- x[, c(int, int + x.order$order), drop=FALSE]
    x.names <- x.names[x.order$order]
    df <- df[x.order$order]
    alpha <- alpha[x.order$order]
    select <- select[x.order$order]
    scaling <- scaling[x.order$order]
#
# Set up powers & working matrix
#
    pwrs.mx <- matrix(c(rep(1, nx), rep(NA, nx)), ncol = 2, byrow=FALSE, 
        dimnames = list(x.names, c("power1", "power2")))    
    # NA = no power
    pvals.mx <- matrix(rep(NA, nx * 6), ncol = 6, dimnames = list(x.names,
        c("p.null", "p.lin", "p.FP", "power2", "power4.1", "power4.2")))
    scale.mx <- matrix(NA, ncol = 2, nrow = nx, dimnames = list(x.names, c(
        "shift", "scale")))
    df.work <- rep(1, nx)
    if(nx > 1)
        pwrs.comp <- matrix(ncol = 2, nrow = nx * (nx - 1))
    x.work <- matrix(0, nrow = nobs, ncol = int + 2 * nx)
    if(int)
        x.work[, 1] <- x[, 1]
    x.work[, 2 * seq(nx) - 1 + int] <- x[, seq(nx) + int]
    x.work.names <- c(paste(rep(x.names, each = 2), ".", rep(1:2, nx), sep
         = ""))
    if(cox)
        dimnames(x.work) <- list(1:nobs, x.work.names)
    else dimnames(x.work) <- list(1:nobs, c("Intercept", x.work.names))
    pwrs.stable <- 0    
#
# Backfitting loop
#
    its <- 0
        if(verbose) {
      pos <- c(1, 3, 5) # output formatting
      fp.out("Variable", "Deviance", "Power(s)", pos = pos)
      cat("\n", rep("-", 48), sep = "")
        }
    while(!pwrs.stable) {
        its <- its + 1
        j <- 0
        while(j < nx & !pwrs.stable) {
            j <- j + 1  #
#
# Check convergence at start of loop
#
            if(nx > 1) {
                num <- (j - 1) * (nx - 1) + seq(nx - 1)
                if(its > 1) {
                  mx.old <- pwrs.comp[num,  , drop=FALSE]
                  mx.new <- pwrs.mx[ - j,  , drop=FALSE]
                  pwrs.stable <- mx.com(mx.old, mx.new)
                }
                pwrs.comp[num,  ] <- pwrs.mx[ - j,  , drop=FALSE]
                dfr <- nobs - sum(df.work[ - j]) - int
            }
            else dfr <- nobs - int
            if(!pwrs.stable) {
#
# Set up matrices and fit and find best FP
#
                num <- 2 * (j - 1) + int + seq(2)
                xj <- x[, j + int]
                if(its == 1) {
                          xj.transform <- fp.scale(xj, scaling[j])
                  scale.mx[j, 1] <- xj.transform$shift
                  scale.mx[j, 2] <- xj.transform$scale
                }
# if(its==1 & j==2) browser()
                fitj <- fp.fit(cbind(xj, x.work[,  - num, drop=FALSE]),   
                  y, df[j], dfr, cox, gauss, scale.mx[j, 1], 
                  scale.mx[j, 2], ...)
                res <- fp.sel(fitj, alpha[j], select[j])
                best.fitj <- res$results
                fit.fitj <- res$fit
                pwrsj <- best.fitj$pwrs
                pwrs.mx[j,  ] <- pwrsj
                df.work[j] <- best.fitj$df
                x.work[, num] <- fp.gen(xj, pwrsj, scale.mx[j, 1], 
                 scale.mx[j, 2])   
#
# Verbose output
#
        if(verbose) {
                if(j==1) cat("\nCycle", its)
            pos <- c(1, 3, 5)   
                namej <- x.names[[j]]
            fp.out(namej, round(fit.fitj$dev0, 3), " ", pos = pos)
            fp.out("", round(fit.fitj$dev1, 3), "1", pos = pos)
            fp.out("", round(fit.fitj$dev2, 3), fit.fitj$pwr2, pos = pos)
            fp.out("", round(fit.fitj$dev4, 3), fit.fitj$pwr4, pos = pos)
# best fit
#     pos <- c(0, 3, 5) 
#     fp.out("selected", round(best.fitj$dev, 3), pwrsj, pos = pos)
            cat("\n")   #
        }
#
            if(df[j] < 4)
                pvals.mx[j,  ] <- c(best.fitj$p.null,
                   best.fitj$p.lin, best.fitj$p.FP, fitj$pwr2, NA, NA)
            else pvals.mx[j,  ] <- c(best.fitj$p.null,
                best.fitj$p.lin, best.fitj$p.FP, fitj$pwr2, fitj$pwr4)
            } # end of "if(!powers.stable)"
        }   # end of "while(j < nx & !pwrs.stable)"
        if(nx == 1)
            pwrs.stable <- 1
    } # end of "while(!pwrs.stable)"
    if(verbose) cat("\n")   # end of cycle its
#
# Set up final matrix
#
    num <- int + seq(2 * nx)[is.na(as.vector(t(pwrs.mx)))]
    if(sum(num, na.rm=TRUE) > 0)
        x.work <- x.work[,  - num, drop=FALSE]
    if(cox) {
        if(exists("coxph.fit")) fitter <- get("coxph.fit")
        else fitter <- getFromNamespace("coxph.fit","survival")
    }
    else {
        fitter <- get("glm.fit")
    }
    fit <- fitter(x.work, y, ...)
    fit$x <- x.work
    fit$powers <- pwrs.mx
    fit$pvalues <- pvals.mx
    fit$scale <- scale.mx
    fit$df.initial <- matrix(df, dimnames = list(x.names, "df.initial"))
    fit$df.final <- matrix(df.work, dimnames = list(x.names, "df.final"))
    fit$dev <- best.fitj$dev
    fit$dev.lin <- x.order$dev[2]
    fit$dev.null <- x.order$dev[1]
#
# Final output
#
    power1 <- pwrs.mx[, 1]
    power2 <- pwrs.mx[, 2]
    power1[is.na(power1)] <- "."
    power2[is.na(power2)] <- "."
    fit$fptable <- data.frame(df.initial = df, select, alpha, df.final = df.work,
        power1, power2, row.names = x.names)
if(verbose) {
    cat("\nTansformation\n")
    print(fit$scale)
    cat("\nFractional polynomials\n")
    print(fit$fptable)
    cat("\n")
    cat("\nNull model: "); cat(fit$dev.null)
       cat("\nLinear model: "); cat(fit$dev.lin)
          cat("\nFinal model: "); cat(fit$dev)
             cat("\n")
}
fit
}
mx.com <- function(x, y)
{
#
# Version 1.0.0     8 May 1998
#
    if(is.matrix(x) & is.matrix(y)) {
        if(sum(dim(x) == dim(y)) != 2)
            stop("matrices have different dimensions")
    }
    x <- as.vector(x)
    y <- as.vector(y)
    nx <- length(x)
    if(nx != length(y))
        stop("vector lengths unequal")
    na.same <- is.na(x) == is.na(y) # NA status of x & y
    na.diff <- sum(!na.same)
    if(!na.diff) {
        xy.same <- x[!is.na(x)] == y[!is.na(y)] 
    # number status of x & y
        xy.diff <- sum(!xy.same)
    }
    else xy.diff <- 1
    return(!xy.diff)
}
plot.mfp <- function (x, var=NULL, ...) 
{
    if (!inherits(x, "mfp")) 
        stop("This is not an mfp object")
    name <- dimnames(x$powers)[[1]]
    choices <- name
    if(is.null(var)) {
      pick <- seq(name)[-which(is.na(x$powers[,1]))]
    } else { pick <- which(name %in% var) }
    int <- as.numeric(x$family[["family"]] != "Cox")
    for(ip in pick) {
        npwrsx <- sum(!is.na(x$powers[ip, ]))
        if (npwrsx > 0) {
            if (ip > 1) 
                posx <- int + sum(!is.na(x$powers[seq(ip-1), ])) + seq(npwrsx)
            else posx <- int + seq(npwrsx)
            fx <- x$x[, posx, drop = FALSE] %*% x$coef[posx]
        }
        namex <- name[ip]
        if (is.null(x$X)) 
            stop("you did not specify x=T in the fit")
        if (any(dimnames(x$X)[[2]] == namex, na.rm = TRUE)) {
            tmpx <- x$X[, namex]
            ix <- which(dimnames(x$X)[[2]] == namex)
        }
        else {
            tmpx <- eval(as.name(namex))
        }
        ord <- order(tmpx)
        if (int) {
            if (npwrsx > 0) {
                plot(tmpx[ord], fx[ord], xlab = namex, ylab = paste("Linear predictor", 
                  sep = ""), type = "l", ...)
                pres <- x$residuals + fx
                plot(tmpx, pres, xlab = namex, ylab = "Partial residuals", 
                  ...)
                fl <- lowess(tmpx[ord], pres[ord], iter = 0)
                lines(fl$x, fl$y, lwd = 1, col = "red")
            }
        }
        else {
            require(survival)
            x0 <- coxph(x$y ~ 1)
            res0 <- resid(x0, type = "mart")
            plot(tmpx, res0, xlab = namex, ylab = "Martingale residuals", 
                type = "p", ...)
            fl <- lowess(tmpx[ord], res0[ord], iter = 0)
            lines(fl$x, fl$y, lwd = 1, col = "red")
            if (npwrsx > 0) {
                plot(tmpx[ord], fx[ord], xlab = namex, ylab = "Linear predictor", 
                  type = "l", ...)
                pres <- x$residuals + fx
                plot(tmpx, pres, xlab = namex, ylab = "Partial residuals", 
                  ...)
                fl <- lowess(tmpx[ord], pres[ord], iter = 0)
                lines(fl$x, fl$y, lwd = 1, col = "red")
            }
        }
    }
}
