1## Summarize ssanova objects 2summary.ssanova <- function(object,diagnostics=FALSE,...) 3{ 4 y <- model.response(object$mf,"numeric") 5 w <- model.weights(object$mf) 6 offset <- model.offset(object$mf) 7 if (is.null(offset)) offset <- rep(0,length(y)) 8 ## Residuals 9 mf <- object$mf 10 if (!is.null(object$random)) mf$random <- object$random$z 11 res <- y - predict(object,mf) 12 ## Fitted values 13 fitted <- as.numeric(y-res) 14 ## (estimated) sigma 15 sigma <- sqrt(object$varht) 16 ## R^2 17 if (!is.null(w)) { 18 r.squared <- sum(w*(fitted-sum(w*fitted)/sum(w))^2) 19 r.squared <- r.squared/sum(w*(y-sum(w*y)/sum(w))^2) 20 } 21 else r.squared <- var(fitted)/var(y) 22 ## Residual sum of squares 23 if (is.null(w)) rss <- sum(res^2) 24 else rss <- sum(w*res^2) 25 ## Penalty associated with the fit 26 obj.wk <- object 27 obj.wk$d[] <- 0 28 if (!is.null(model.offset(obj.wk$mf))) obj.wk$mf[,"(offset)"] <- 0 29 penalty <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,])) 30 penalty <- as.vector(10^object$nlambda*penalty) 31 if (!is.null(object$random)) { 32 p.ran <- t(object$b)%*%object$random$sigma$fun(object$zeta, 33 object$random$sigma$env)%*%object$b 34 penalty <- penalty + p.ran 35 } 36 ## Calculate the diagnostics 37 if (is.null(object$partial)) labels.p <- NULL 38 else labels.p <- labels(object$partial$mt) 39 if (diagnostics) { 40 ## Obtain retrospective linear model 41 comp <- NULL 42 p.dec <- NULL 43 for (label in c(object$terms$labels,labels.p)) { 44 if (label=="1") next 45 if (label=="offset") next 46 comp <- cbind(comp,predict(object,object$mf,inc=label)) 47 jk <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,],inc=label)) 48 p.dec <- c(p.dec,10^object$nlambda*jk) 49 } 50 term.label <- object$terms$labels[object$terms$labels!="1"] 51 term.label <- term.label[term.label!="offset"] 52 term.label <- c(term.label,labels.p) 53 if (!is.null(object$random)) { 54 comp <- cbind(comp,predict(object,mf,inc=NULL)) 55 p.dec <- c(p.dec,p.ran) 56 term.label <- c(term.label,"random") 57 } 58 fitted.off <- fitted-offset 59 comp <- cbind(comp,yhat=fitted.off,y=fitted.off+res,e=res) 60 if (any(outer(term.label,c("yhat","y","e"),"=="))) 61 warning("gss warning in summary.ssanova: avoid using yhat, y, or e as variable names") 62 colnames(comp) <- c(term.label,"yhat","y","e") 63 ## Sweep out constant 64 if (!is.null(w)) 65 comp <- sqrt(w)*comp - outer(sqrt(w),apply(w*comp,2,sum))/sum(w) 66 else comp <- sweep(comp,2,apply(comp,2,mean)) 67 ## Obtain pi 68 comp1 <- comp[,c(term.label,"yhat")] 69 decom <- t(comp1) %*% comp1[,"yhat"] 70 names(decom) <- c(term.label,"yhat") 71 decom <- decom[term.label]/decom["yhat"] 72 ## Obtain kappa, norm, and cosines 73 corr <- t(comp)%*%comp 74 corr <- t(corr/sqrt(diag(corr)))/sqrt(diag(corr)) 75 norm <- apply(comp,2,function(x){sqrt(sum(x^2))}) 76 cosines <- rbind(corr[c("y","e"),],norm) 77 rownames(cosines) <- c("cos.y","cos.e","norm") 78 corr <- corr[term.label,term.label,drop=FALSE] 79 if (qr(corr)$rank<dim(corr)[2]) 80 kappa <- rep(Inf,len=dim(corr)[2]) 81 else kappa <- as.numeric(sqrt(diag(solve(corr)))) 82 ## Obtain decomposition of penalty 83 rough <- p.dec / penalty 84 names(kappa) <- names(rough) <- term.label 85 } 86 else decom <- kappa <- cosines <- rough <- NULL 87 ## Return the summaries 88 z <- list(call=object$call,method=object$method,fitted=fitted,residuals=res, 89 sigma=sigma,r.squared=r.squared,rss=rss,penalty=penalty, 90 pi=decom,kappa=kappa,cosines=cosines,roughness=rough) 91 class(z) <- "summary.ssanova" 92 z 93} 94