| # File src/library/stats/R/ppr.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright (C) 1998 B. D. Ripley |
| # Copyright (C) 2000-2016 The R Core Team |
| # |
| # This program is free software; you can redistribute it and/or modify |
| # it under the terms of the GNU General Public License as published by |
| # the Free Software Foundation; either version 2 of the License, or |
| # (at your option) any later version. |
| # |
| # This program is distributed in the hope that it will be useful, |
| # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| # GNU General Public License for more details. |
| # |
| # A copy of the GNU General Public License is available at |
| # https://www.R-project.org/Licenses/ |
| |
| ppr <- function(x, ...) UseMethod("ppr") |
| |
| ppr.formula <- |
| function(formula, data, weights, subset, |
| na.action, contrasts = NULL, ..., model = FALSE) |
| { |
| call <- match.call() |
| m <- match.call(expand.dots = FALSE) |
| m$contrasts <- m$... <- NULL |
| ## need stats:: for non-standard evaluation |
| m[[1L]] <- quote(stats::model.frame) |
| m <- eval(m, parent.frame()) |
| Terms <- attr(m, "terms") |
| attr(Terms, "intercept") <- 0L |
| X <- model.matrix(Terms, m, contrasts) |
| Y <- model.response(m) |
| w <- model.weights(m) |
| if(length(w) == 0L) w <- rep_len(1, nrow(X)) |
| fit <- ppr.default(X, Y, w, ...) |
| fit$na.action <- attr(m, "na.action") |
| fit$terms <- Terms |
| ## fix up call to refer to the generic, but leave arg name as `formula' |
| call[[1L]] <- as.name("ppr") |
| fit$call <- call |
| fit$contrasts <- attr(X, "contrasts") |
| fit$xlevels <- .getXlevels(Terms, m) |
| if(model) fit$model <- m |
| structure(fit, class=c("ppr.form", "ppr")) |
| } |
| |
| ppr.default <- |
| function(x, y, weights=rep(1,n), ww=rep(1,q), nterms, max.terms=nterms, |
| optlevel=2, sm.method=c("supsmu", "spline", "gcvspline"), |
| bass=0, span=0, df=5, gcvpen=1, trace = FALSE, ...) |
| { |
| call <- match.call() |
| call[[1L]] <- as.name("ppr") |
| sm.method <- match.arg(sm.method) |
| ism <- switch(sm.method, supsmu = 0L, spline = 1L, gcvspline = 2L) |
| if(trace) ism <- -(ism + 1L) |
| if(missing(nterms)) stop("'nterms' is missing with no default") |
| mu <- nterms; ml <- max.terms |
| x <- as.matrix(x) |
| y <- as.matrix(y) |
| if(!is.numeric(x) || !is.numeric(y)) |
| stop("'ppr' applies only to numerical variables") |
| n <- nrow(x) |
| if(nrow(y) != n) stop("mismatched 'x' and 'y'") |
| p <- ncol(x) |
| q <- ncol(y) |
| xnames <- if(!is.null(dimnames(x))) dimnames(x)[[2L]] else paste0("X", 1L:p) |
| ynames <- if(!is.null(dimnames(y))) dimnames(y)[[2L]] else paste0("Y", 1L:q) |
| msmod <- ml*(p+q+2*n)+q+7+ml+1 # for asr |
| nsp <- n*(q+15)+q+3*p |
| ndp <- p*(p+1)/2+6*p |
| .Fortran(C_setppr, |
| as.double(span), as.double(bass), as.integer(optlevel), |
| as.integer(ism), as.double(df), as.double(gcvpen) |
| ) |
| Z <- .Fortran(C_smart, |
| as.integer(ml), as.integer(mu), |
| as.integer(p), as.integer(q), as.integer(n), |
| as.double(weights), |
| as.double(t(x)), |
| as.double(t(y)), |
| as.double(ww), |
| smod=double(msmod), as.integer(msmod), |
| double(nsp), as.integer(nsp), |
| double(ndp), as.integer(ndp), |
| edf=double(ml) |
| ) |
| smod <- Z$smod |
| ys <- smod[q+6] |
| tnames <- paste("term", 1L:mu) |
| alpha <- matrix(smod[q+6L + 1L:(p*mu)],p, mu, |
| dimnames=list(xnames, tnames)) |
| beta <- matrix(smod[q+6L+p*ml + 1L:(q*mu)], q, mu, |
| dimnames=list(ynames, tnames)) |
| fitted <- drop(matrix(.Fortran(C_pppred, |
| as.integer(nrow(x)), |
| as.double(x), |
| as.double(smod), |
| y = double(nrow(x)*q), |
| double(2*smod[4L]))$y, |
| ncol=q, dimnames=dimnames(y))) |
| jt <- q + 7 + ml*(p+q+2*n) |
| gof <- smod[jt] * n * ys^2 |
| gofn <- smod[jt+1L:ml] * n * ys^2 |
| ## retain only terms for the size of model finally fitted |
| jf <- q+6+ml*(p+q) |
| smod <- smod[c(1L:(q+6+p*mu), q+6+p*ml + 1L:(q*mu), |
| jf + 1L:(mu*n), jf+ml*n + 1L:(mu*n))] |
| smod[1L] <- mu |
| structure(list(call=call, mu=mu, ml=ml, p=p, q=q, |
| gof=gof, gofn=gofn, |
| df=df, edf=Z$edf[1L:mu], |
| xnames=xnames, ynames=ynames, |
| alpha=drop(alpha), beta=ys*drop(beta), |
| yb=smod[5+1L:q], ys=ys, |
| fitted.values=fitted, residuals=drop(y-fitted), |
| smod=smod), |
| class="ppr") |
| } |
| |
| print.ppr <- function(x, ...) |
| { |
| if(!is.null(cl <- x$call)) { |
| cat("Call:\n") |
| dput(cl, control=NULL) |
| } |
| mu <- x$mu; ml <- x$ml |
| cat("\nGoodness of fit:\n") |
| gof <- setNames(x$gofn, paste(1L:ml, "terms")) |
| print(format(gof[mu:ml], ...), quote=FALSE) |
| invisible(x) |
| } |
| |
| summary.ppr <- function(object, ...) |
| { |
| class(object) <- "summary.ppr" |
| object |
| } |
| |
| print.summary.ppr <- function(x, ...) |
| { |
| print.ppr(x, ...) |
| mu <- x$mu |
| cat("\nProjection direction vectors ('alpha'):\n") |
| print(format(x$alpha, ...), quote=FALSE) |
| cat("\nCoefficients of ridge terms ('beta'):\n") |
| print(format(x$beta, ...), quote=FALSE) |
| if(any(x$edf >0)) { |
| cat("\nEquivalent df for ridge terms:\n") |
| edf <- setNames(x$edf, paste("term", 1L:mu)) |
| print(round(edf,2), ...) |
| } |
| invisible(x) |
| } |
| |
| plot.ppr <- function(x, ask, type = "o", cex = 1/2, |
| main = quote(bquote( |
| "term"[.(i)]*":" ~~ hat(beta[.(i)]) == .(bet.i))), |
| xlab = quote(bquote(bold(alpha)[.(i)]^T * bold(x))), |
| ylab = "", ...) |
| { |
| ppr.funs <- function(obj) |
| { |
| ## cols for each term |
| p <- obj$p; q <- obj$q |
| sm <- obj$smod |
| n <- sm[4L]; mu <- sm[5L]; m <- sm[1L] |
| jf <- q+6+m*(p+q) |
| jt <- jf+m*n |
| f <- matrix(sm[jf+1L:(mu*n)],n, mu) |
| t <- matrix(sm[jt+1L:(mu*n)],n, mu) |
| list(x=t, y=f) |
| } |
| obj <- ppr.funs(x) |
| if(!missing(ask)) { |
| oask <- devAskNewPage(ask) |
| on.exit(devAskNewPage(oask)) |
| } |
| for(i in 1L:x$mu) { |
| ord <- order(obj$x[ ,i]) |
| bet.i <- format(x$beta[[i]], digits = 3) |
| plot(obj$x[ord, i], obj$y[ord, i], type = type, cex = cex, |
| main = if(is.call(main)) eval(main) else main, |
| xlab = if(is.call(xlab)) eval(xlab) else xlab, |
| ylab = ylab, ...) |
| } |
| force(bet.i)# codetools |
| invisible() |
| } |
| |
| predict.ppr <- function(object, newdata, ...) |
| { |
| if(missing(newdata)) return(fitted(object)) |
| if(!is.null(object$terms)) { |
| newdata <- as.data.frame(newdata) |
| rn <- row.names(newdata) |
| # work hard to predict NA for rows with missing data |
| Terms <- delete.response(object$terms) |
| m <- model.frame(Terms, newdata, na.action = na.omit, |
| xlev = object$xlevels) |
| if(!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m) |
| keep <- match(row.names(m), rn) |
| x <- model.matrix(Terms, m, contrasts.arg = object$contrasts) |
| } else { |
| x <- as.matrix(newdata) |
| keep <- seq_len(nrow(x)) |
| rn <- dimnames(x)[[1L]] |
| } |
| if(ncol(x) != object$p) stop("wrong number of columns in 'x'") |
| res <- matrix(NA, length(keep), object$q, |
| dimnames = list(rn, object$ynames)) |
| res[keep, ] <- matrix(.Fortran(C_pppred, |
| as.integer(nrow(x)), |
| as.double(x), |
| as.double(object$smod), |
| y = double(nrow(x)*object$q), |
| double(2*object$smod[4L]) |
| )$y, ncol=object$q) |
| drop(res) |
| } |