R/summary.ssanova9.R

In gss: General Smoothing Splines

## Summarize ssanova objects
summary.ssanova9 <- function(object,diagnostics=FALSE,...)
{
    y <- model.response(object$mf,"numeric")
    nobs <- length(y)
    cov <- object$cov
    if (length(object$zeta)) ww <- cov$fun(object$zeta,cov$env)
    else ww <- cov$fun(cov$env)
    ww <- chol(ww)
    offset <- model.offset(object$mf)
    if (is.null(offset)) offset <- rep(0,length(y))
    ## Residuals
    mf <- object$mf
    res <- y - predict(object,mf)
    ## Fitted values
    fitted <- as.numeric(y-res)
    ## (estimated) sigma
    sigma <- sqrt(object$varht)
    ## R^2
    y.wk <- forwardsolve(t(ww),y)
    fitted.wk <- forwardsolve(t(ww),fitted)
    one.wk <- forwardsolve(t(ww),rep(1,nobs))
    mn.y <- sum(y.wk*one.wk)/sum(one.wk^2)
    mn.fitted <- sum(fitted.wk*one.wk)/sum(one.wk^2)
    r.squared <- sum((fitted.wk-mn.fitted*one.wk)^2)
    r.squared <- r.squared/sum((y.wk-mn.y*one.wk)^2)
    ## Residual sum of squares
    res.wk <- forwardsolve(t(ww),res)
    rss <- sum(res.wk^2)
    ## Penalty associated with the fit
    obj.wk <- object
    obj.wk$d[] <- 0
    if (!is.null(model.offset(obj.wk$mf))) obj.wk$mf[,"(offset)"] <- 0
    penalty <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,]))
    penalty <- as.vector(10^object$nlambda*penalty)
    ## Calculate the diagnostics
    if (is.null(object$partial)) labels.p <- NULL
    else labels.p <- labels(object$partial$mt)
    if (diagnostics) {
        ## Obtain retrospective linear model
        comp <- NULL
        p.dec <- NULL
        for (label in c(object$terms$labels,labels.p)) {
            if (label=="1") next
            if (label=="offset") next
            comp <- cbind(comp,predict(object,mf,inc=label))
            jk <- sum(obj.wk$c*predict(obj.wk,mf[object$id.basis,],inc=label))
            p.dec <- c(p.dec,10^object$nlambda*jk)
        }
        term.label <- object$terms$labels[object$terms$labels!="1"]
        term.label <- term.label[term.label!="offset"]
        term.label <- c(term.label,labels.p)
        fitted.off <- fitted-offset
        comp <- cbind(comp,yhat=fitted.off,y=fitted.off+res,e=res)
        comp <- forwardsolve(t(ww),comp)
        if (any(outer(term.label,c("yhat","y","e"),"==")))
            warning("gss warning in summary.ssanova: avoid using yhat, y, or e as variable names")
        colnames(comp) <- c(term.label,"yhat","y","e")
        ## Sweep out constant
        comp <- comp - outer(one.wk,apply(t(comp)%*%one.wk,1,sum))/sum(one.wk^2)
        ## Obtain pi
        comp1 <- comp[,c(term.label,"yhat")]
        decom <- t(comp1) %*% comp1[,"yhat"]
        names(decom) <- c(term.label,"yhat")
        decom <- decom[term.label]/decom["yhat"]
        ## Obtain kappa, norm, and cosines
        corr <- t(comp)%*%comp
        corr <- t(corr/sqrt(diag(corr)))/sqrt(diag(corr))
        norm <- apply(comp,2,function(x){sqrt(sum(x^2))})
        cosines <- rbind(corr[c("y","e"),],norm)
        rownames(cosines) <- c("cos.y","cos.e","norm")
        corr <- corr[term.label,term.label,drop=FALSE]
        if (qr(corr)$rank<dim(corr)[2])
            kappa <- rep(Inf,len=dim(corr)[2])
        else kappa <- as.numeric(sqrt(diag(solve(corr))))
        ## Obtain decomposition of penalty
        rough <- p.dec / penalty
        names(kappa) <- names(rough) <- term.label
    }
    else decom <- kappa <- cosines <- rough <- NULL
    ## Return the summaries
    z <- list(call=object$call,method=object$method,fitted=fitted,residuals=res,
              sigma=sigma,r.squared=r.squared,rss=rss,penalty=penalty,
              pi=decom,kappa=kappa,cosines=cosines,roughness=rough)
    class(z) <- "summary.ssanova"
    z
}