##------------------------------------------------------------------------## ## Script for web appendix on nonparametric regression ## ## An R and S-PLUS Companion to Applied Regression ## ## John Fox ## ## Sage Publications, 2002 ## ##------------------------------------------------------------------------## # preliminaries options(width=65) options(digits=5) library(car) # for data sets ## Local polynomial regression # How it works tricube <- function(z) { ifelse (abs(z) < 1, (1 - (abs(z))^3)^3, 0) } data(Prestige) attach(Prestige) par(mfrow=c(2,2)) plot(income, prestige, xlab="Average Income", ylab="Prestige", type='n', main="(a)") ord <- order(income) inc <- income[ord] pre <- prestige[ord] x0 <- inc[80] diffs <- abs(inc - x0) which.diff <- sort(diffs)[50] abline(v=c(x0-which.diff, x0+which.diff), lty=2) abline(v=x0) points(inc[diffs > which.diff], pre[diffs > which.diff], pch=16, col=gray(.75)) points(inc[diffs <= which.diff], pre[diffs <= which.diff]) x.n <- inc[diffs <= which.diff] y.n <- pre[diffs <= which.diff] plot(range(income), c(0,1), xlab="Average Income", ylab="Tricube Weight", type='n', main="(b)") abline(v=c(x0-which.diff, x0+which.diff), lty=2) abline(v=x0) xwts <- seq(x0-which.diff, x0+which.diff, len=250) lines(xwts, tricube((xwts-x0)/which.diff), lty=1, lwd=2) points(x.n, tricube((x.n - x0)/which.diff)) plot(income, prestige, xlab="Average Income", ylab="Prestige", type='n', main="(c)") abline(v=c(x0-which.diff, x0+which.diff), lty=2) abline(v=x0) points(x.n, y.n) mod <- lm(y.n ~ x.n, weights=tricube((x.n-x0)/which.diff)) reg.line(mod, lwd=2, col=1) points(x0, predict(mod, data.frame(x.n=x0)), pch=16, cex=1.8) plot(income, prestige, xlab="Average Income", ylab="Prestige", main="(d)") lines(lowess(income, prestige, f=0.5, iter=0), lwd=2) # multiple regression library(modreg) mod.lo <- loess(prestige ~ income + education, span=.5, degree=1) summary(mod.lo) inc <- seq(min(income), max(income), len=25) ed <- seq(min(education), max(education), len=25) newdata <- expand.grid(income=inc,education=ed) fit.prestige <- matrix(predict(mod.lo, newdata), 25, 25) persp(inc, ed, fit.prestige, theta=45, phi=30, ticktype='detailed', xlab='Income', ylab='Education', zlab='Prestige', expand=2/3, shade=0.5) mod.lo.inc <- loess(prestige ~ income, span=.7, degree=1) mod.lo.ed <- loess(prestige ~ education, span=.7, degree=1) anova(mod.lo.inc, mod.lo) anova(mod.lo.ed, mod.lo) # smoothing splines plot(income, prestige) inc.100 <- seq(min(income), max(income), len=100) pres <- predict(mod.lo.inc, data.frame(income=inc.100)) lines(inc.100, pres, lty=2, lwd=2) mod.lo.inc lines(smooth.spline(income, prestige, df=3.85), lwd=2) # additive regression models library(mgcv) mod.gam <- gam(prestige ~ s(income) + s(education)) mod.gam fit.prestige <- matrix(predict(mod.gam, newdata), 25, 25) persp(inc, ed, fit.prestige, theta=45, phi=30, ticktype='detailed', xlab='Income', ylab='Education', zlab='Prestige', expand=2/3, shade=0.5) par(mfrow=c(1,2)) plot(mod.gam) fit.prestige <- matrix(predict(mod.gam, newdata), 25, 25) persp(inc, ed, fit.prestige, theta=45, phi=30, ticktype='detailed', xlab='Income', ylab='Education', zlab='Prestige', expand=2/3, shade=0.5) # generalized nonparametric regression detach(Prestige) remove(list=objects()) # clean up everything data(Mroz) attach(Mroz) k5f <- factor(k5) k618f <- factor(k618) mod.1 <- gam(lfp ~ s(age) + s(inc) + k5f + k618f + wc + hc, family=binomial) summary(mod.1) par(mfrow=c(1,2)) plot(mod.1) mod.2 <- gam(lfp ~ age + s(inc) + k5f + k618f + wc + hc, family=binomial) anova(mod.2, mod.1, test='Chisq') mod.3 <- gam(lfp ~ s(age) + inc + k5f + k618f + wc + hc, family=binomial) anova(mod.3, mod.1, test='Chisq') mod.4 <- gam(lfp ~ s(inc) + k5f + k618f + wc + hc, # omits age family=binomial) anova(mod.4, mod.1, test='Chisq') ## Some methods for gam objects in the R mgcv package summary.gam <- function(object, smooth.coef=FALSE, ...){ result <- list() result\$call <- object\$call result\$family <- object\$family terms <- terms(formula(object)) trms <- attr(terms, "term.labels") vars <- attr(terms, "variables") sm.posn <- grep("^s\\(.*\\)", trms) nonsm.posn <- setdiff(1:length(trms), sm.posn) environ <- if (is.null(object\$call\$data)) attr(terms(object\$formula), ".Environment") else eval(object\$call\$data) vars.eval <- eval(vars, envir=environ)[-1] factors <- sapply(vars.eval, is.factor) levs <- sapply(sapply(vars.eval[factors], levels), length) if (object\$nsdf > 0){ coefs <- coefficients(object)[1:object\$nsdf] se <- sqrt(diag(object\$Vp)[1:length(coefs)]) z <- coefs/se table <- cbind(coefs, se, z, 2*(1-pnorm(abs(z)))) rownames(table) <- names(coefs) colnames(table) <- c('Estimate', 'Std. Error', 'z value', 'Pr(>|z|)') result\$table <- table } if (smooth.coef){ sel <- (object\$nsdf + 1):length(coefficients(object)) s.coefs <- coefficients(object)[sel] s.se <- sqrt(diag(object\$Vp)[sel]) z.s <- s.coefs/s.se s.table <- cbind(s.coefs, s.se, z.s, 2*(1-pnorm(abs(z.s)))) rownames(s.table) <- names(s.coefs) colnames(s.table) <- c('Estimate', 'Std. Error', 'z value', 'Pr(>|z|)') result\$s.table <- s.table } df.terms <- rep(1, length(trms)) df.terms[sm.posn] <- object\$edf if (length(levs) > 0) df.terms[factors] <- levs - 1 if (attr(terms, "intercept") == 1) { df.terms <- c(1, df.terms) names(df.terms) <- c('(Intercept)', trms) } else names(df.terms) <- trms result\$df.terms <- df.terms n <- object\$df.null result\$n <- n df.residual <- n - sum(object\$edf) - object\$nsdf result\$df.residual <- df.residual result\$gcv <- n * object\$sig2/df.residual result\$deviance <- object\$deviance result\$null.deviance <- object\$null.deviance result\$dispersion <- object\$sig2 class(result) <- 'summary.gam' result } print.summary.gam <- function (x, digits = max(3, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"), ...){ x\$dispersion <- if (is.element(x\$family[[1]], c('binomial', 'Poisson'))) "taken as 1" else x\$dispersion cat('Call:\n') print(x\$call) print(x\$family) if (!is.null(x\$table)){ cat('Parametric Coefficients\n') print.coefmat(x\$table, digits = digits, signif.stars = signif.stars, ...) cat('\n') } if (!is.null(x\$s.table)){ cat('Smooth Coefficients\n') print.coefmat(x\$s.table, digits = digits, signif.stars = signif.stars, ...) cat('\n') } cat('\nDegrees of freedom for terms in the model:\n') print(x\$df.terms) cat("\nNull Deviance: ", x\$null.deviance) cat("\nDeviance: ", x\$deviance) cat("\nGCV score: ", x\$gcv) cat("\nResidual Degrees of Freedom :", x\$df.residual) cat("\nDispersion: ", x\$dispersion) cat("\nNumber of Observations: ", x\$n, "\n") invisible(x) } anova.gam <- function (object, ..., dispersion = NULL, test = NULL) # adapted from anova.glm { object\$df.residual <- object\$df.null - sum(object\$edf) - object\$nsdf dotargs <- list(...) named <- if (is.null(names(dotargs))) rep(FALSE, length(dotargs)) else (names(dotargs) != "") if (any(named)) warning(paste("The following arguments to anova.glm(..)", "are invalid and dropped:", paste(deparse(dotargs[named]), collapse = ", "))) dotargs <- dotargs[!named] is.gam <- unlist(lapply(dotargs, function(x) inherits(x, "gam"))) dotargs <- dotargs[is.gam] if (length(dotargs) == 0) stop('need two or more models to compare') for (i in 1:length(dotargs)) { class(dotargs[[i]]) <- c(class(dotargs[[i]]), 'glm') dotargs[[i]]\$df.residual <- dotargs[[i]]\$df.null - sum(dotargs[[i]]\$edf) - dotargs[[i]]\$nsdf } anova.glmlist(c(list(object), dotargs), test = test) }