"plot.lm" <- function(reg.object) { res<-residuals(reg.object) pred<-fitted(reg.object) plot(pred,res,ylab="Residuals",xlab="Fitted values") lines(lowess(pred,res)) abline(h=0,lty=2) title("Plot of residuals versus fitted values") qqnorm(residuals(reg.object),ylab="Residuals") abline(0,sqrt(var(res)),lty=2) title("Normal plot of residuals") } "diag.plots" <- function(reg.object) { layout(2,2) plot(reg.object) resids <- residuals(reg.object) pred <- fitted(reg.object) plot(pred, resids, type = "n",main="Labelled plot") text(pred, resids, 1:length(pred)) lrplot(reg.object) layout(1, 1) invisible() } "is.influential" <- function(infmat) { # Argument is result of using influence.measures # Returns a matrix of logicals structured like the argument n <- nrow(infmat) k <- ncol(infmat) - 4 if(n <= k) stop("Too few cases, n 1, absmat[, k + 1] > 3 * sqrt(k/(n - k)), abs(1 - infmat[, k + 2]) > (3 * k)/(n - k), qf(infmat[, k + 3], k, n - k) > 0.9, infmat[, k + 4] > (3 * k)/n) dimnames(result) <- dimnames(infmat) result } "influence.measures" <- function(lm.object) { # Calculates DFBETAS, DFFITS, CovRatio, Cooks D. # Argument is lm object. # Returns the results as a matrix. inf <- lm.influence(lm.object) hats <- inf$hat si <- inf$sigma s <- summary(lm.object)$sigma e <- residuals(lm.object) beta <- coefficients(lm.object) betaij <- inf$coefficients k <- ncol(betaij) # dfbetas ses <- outer(si, sqrt(diag(summary(lm.object)$cov.unscaled))) dfbetas <- betaij/ses #dffits dffits <- (e * sqrt(hats))/(si * (1 - hats)) #covratio cov.ratio <- (si/s)^(2 * k)/(1 - hats) # Cook's D cooks.d <- ((e/(s * (1 - hats)))^2 * hats)/k infmat <- cbind(dfbetas, dffits, cov.ratio, cooks.d, hats) is.star <- apply(is.influential(infmat), 1, any) temp.df <- data.frame(round(infmat, 3), inf = ifelse(is.star, "*", " ")) print(temp.df) invisible(infmat) } "funnel" <- function(reg.object) { # plots log sds vs log means # and then squared residuals versus fitted values # returns a smoothed version of the plot as an # estimate of the variance function par(mfrow=c(1,2)) pred<-fitted(reg.object) res<-residuals(reg.object) cut.points<-quantile(pred,c(0.,.2,.4,.6,.8,1.0)) + c(-0.01,0,0,0,0,1.01) group<-cut(pred,cut.points) log.means<-log(tapply(pred+res,group,mean)) log.sds<-log(sqrt(tapply(pred+res,group,var))) plot(log.means,log.sds) res.sq<-res^2 plot(pred,res.sq) low.stuff<-lowess(res.sq~pred) lines(sort(pred),fitted(low.stuff)[order(pred)],lty=2) cat("Slope:",lm(log.sds~log.means)$coef[2],"\n") layout(1,1) invisible(fitted(low.stuff)) } "boxcoxplot" <- function (X,y,p = seq(-1, 2, length = 20), ...) { # X matrix of explanatory variables # y response variable l <- length(p) boxcox <- seq(l) n <- length(y) sumlog <- sum(log(y)) for (i in seq(l)) { y.p<- (y^p[i] -1)/p[i] trans.res<-residuals(lm(y.p~X)) ResSS.p<-sum(trans.res^2) boxcox[i] <- n * log(ResSS.p)/2 - (p[i] - 1) * sumlog } plot(p, boxcox, type = "l", ylab = "Profile liklihood", main = "Box-Cox plot", ...) } "partial.plots" <- function(lm.obj,intercept=T) { X <- model.matrix(lm.obj$terms, model.frame(lm.obj)) y<-X%*%coefficients(lm.obj)+residuals(lm.obj) X<-X[,-1] k<-dim(X)[2] for(i in (1:k)) { z1 <- lsfit(X[, - i], X[, i], intercept = intercept) z2 <- lsfit(X[, - i], y, intercept = intercept) plotname <- dimnames(X)[[2]][i] plot(z1$residuals, z2$residuals, xlab = plotname, ylab = "Residuals", main = paste("Partial plot of", plotname)) abline(h=0,lty=2) lines(lowess(residuals(z1),residuals(z2))) } invisible() } "all.poss.regs" <- function(X,y,nbest=3) { selection.stuff<-stepwise(X,y,method="ex",nbest=nbest) n<-length(y) rssp<-selection.stuff$rss p<-selection.stuff$size Rsq<-1-rssp/sum((y-mean(y))^2) sigma2<-rssp/(n-p-1) adjRsq<- 1-(1-Rsq)*(n-1)/(n-p-1) Cp<-rssp/sigma2[length(p)] + 2*(p+1)-n round(cbind(rssp,sigma2,adjRsq,Cp,selection.stuff$which),3) } "wireframe.plot" <- function(x, y, z, ...) { lowess.fit <- loess(z ~ x * y, span = 1, degree = 2) x.marg <- seq(min(x), max(x), length = 50) y.marg <- seq(min(y), max(y), length = 50) grid <- expand.grid(list(x = x.marg, y = y.marg)) fit <- predict(lowess.fit, grid) wireframe(fit ~ grid$x * grid$y, ...) } "do.prediction" <- function(reg.object, new.df) { pred.object <- predict(reg.object, new.df, se.fit = T) yhat <- pred.object$fit StdErrA <- pred.object$se.fit sigma.hat <- pred.object$residual.scale StdErrB <- sqrt(StdErrA^2 + sigma.hat^2) conf.int.lower <- yhat - qt(0.975, pred.object$df) * StdErrA conf.int.upper <- yhat + qt(0.975, pred.object$df) * StdErrA pred.int.lower <- yhat - qt(0.975, pred.object$df) * StdErrB pred.int.upper <- yhat + qt(0.975, pred.object$df) * StdErrB result <- cbind(yhat, conf.int.lower, conf.int.upper, pred.int.lower, pred.int.upper) dimnames(result) <- list(dimnames(new.df)[[1]], c("prediction", "c i lower", "c i upper", "pred int lower", "pred int lower")) result } "stepwise.reg" <- function(X,y) { # function to calculate stepwise regression # and format the output nicely stepstuff<-stepwise(X,y) inoutmat<-ifelse(stepstuff$which==T,"in","out") labels<-dimnames(stepstuff$which) var.added<-substring(labels[[1]],first=2) k<-length(var.added) result<-cbind(1:k,stepstuff$size,var.added,inoutmat) dimnames(result)<-list(rep("",k), c("Step", "No. in model","Var in/out", labels[[2]])) cat("Final model:\n", labels[[2]][stepstuff$which[k,]],"\nHistory:\n") print(result,quote=F) invisible(result) } lrplot<-function(lm.object) { hats<-lm.influence(lm.object)$hat plot(residuals(lm.object)^2, hats, xlab= "Squared residuals", ylab = "influence",type="n") text(residuals(lm.object)^2, hats, 1:length(hats)) title("Leverage-residual plot") invisible() } ######################################## layout<-function(x,y) { par(mfrow=c(x,y)) invisible() } ######################################## is.influential<-function(infmat) { # Argument is result of using inluence.measures # Returns a matrix of logicals structured like the argument n<- nrow(infmat) k<- ncol(infmat)-4 if (n<=k) stop("too few cases, n 1, absmat[, k+1] > 3* sqrt(k/(n-k)), abs(1-infmat[, k+2]) > (3*k)/(n-k), pf(infmat[, k+3], k, n-k) > 0.9, infmat[, k+4] > (3*k)/n) dimnames(result) <- dimnames(infmat) result } ######################################## par.residual.plots<-function(glm.object) { X<-model.matrix(glm.object$terms, model.frame(glm.object)) X<-X[,-1] k<-dim(X)[2] old.row<-par("mfrow") mfr<-if(k>2) c(2,2) else if(k==1) c(1,1) else c(1,2) par(mfrow=mfr) p.res<-residuals(glm.object,type="pearson") weights<-if (is.null(glm.object$prior.weights))1 else(glm.object$prior.weights) partial.res<-p.res/sqrt((weights*glm.object$fitted.values*(1- glm.object$fitted.values))) for(i in (1:k)){ plotname<-dimnames(X)[[2]][i] plot(X[, i],partial.res+glm.object$coefficient[i+1]*X[, i], ylab="Partial residuals", xlab=plotname,main=paste("Partial residual plot of", plotname)) lines(lowess(X[, i],partial.res+glm.object$coefficient[i+1]*X[, i])) if((i%%4)==0) pause() } par(mfrow=old.row) invisible() } ######################################## glm.diag.plots<-function(glm.model) { layout(2,2) p.res<-residuals(glm.model,type="pearson") d.res<-residuals(glm.model,type="deviance") plot(d.res,type="n",xlab="Observation number", ylab="Deviance Residuals",main="Index plot of deviance residuals") text(d.res) hats<-lm.influence(glm.model)$hat plot(hats,type="n",ylim=c(0,max(hats)),xlab="Observation Number" ,ylab="Leverage", main="Leverage plot") text(hats) p<-glm.model$rank cooks<-(p.res^2*hats)/(p*(1-hats)^2) plot(cooks,type="n",ylim=c(0,max(cooks)),xlab="Observation number", ylab="Cook's Distance", main="Cook's Distance Plot") text(cooks) dev.changes<-d.res^2+p.res^2*(hats/(1-hats)) plot(dev.changes,type="n",ylim=c(0,max(dev.changes)), xlab="Observation number",ylab="Deviance changes", main="Deviance Changes Plot") text(dev.changes) invisible() }