library(MASS) attach(mcycle) par(mfrow = c(3,2)) plot(times, accel, main="Polynomial regression") lines(times, fitted(lm(accel ~ poly(times, 3)))) lines(times, fitted(lm(accel ~ poly(times, 6))), lty=3) legend(40, -100, c("degree=3", "degree=6"), lty=c(1,3),bty="n") plot(times, accel, main="Natural splines") lines(times, fitted(lm(accel ~ ns(times, df=5)))) lines(times, fitted(lm(accel ~ ns(times, df=10))), lty=3) lines(times, fitted(lm(accel ~ ns(times, df=20))), lty=4) legend(40, -100, c("df=5", "df=10", "df=20"), lty=c(1,3,4), bty="n") plot(times, accel, main="Smoothing splines") lines(smooth.spline(times, accel)) plot(times, accel, main="Lowess") lines(lowess(times, accel)) lines(lowess(times, accel, 0.2), lty=3) legend(40, -100, c("default span", "f = 0.2"), lty=c(1,3), bty="n") plot(times, accel, main ="ksmooth") lines(ksmooth(times, accel,"normal", bandwidth=5)) lines(ksmooth(times, accel,"normal", bandwidth=2), lty=3) legend(40, -100, c("bandwidth=5", "bandwidth=2"), lty=c(1,3), bty="n") plot(times, accel, main ="supsmu") lines(supsmu(times, accel)) lines(supsmu(times, accel, bass=3), lty=3) legend(40, -100, c("default", "bass=3"), lty=c(1,3), bty="n") detach() > attach(rock) > rock.lm <- lm(log(perm) ~ area + peri + shape) > summary(rock.lm) Call: lm(formula = log(perm) ~ area + peri + shape) Residuals: Min 1Q Median 3Q Max -1.809 -0.5413 0.1735 0.6493 1.479 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) 5.3331 0.5487 9.7200 0.0000 area 0.0005 0.0001 5.6021 0.0000 peri -0.0015 0.0002 -8.6228 0.0000 shape 1.7565 1.7559 1.0003 0.3226 Residual standard error: 0.8521 on 44 degrees of freedom Multiple R-Squared: 0.7483 F-statistic: 43.6 on 3 and 44 degrees of freedom, the p-value is 3.094e-13 Correlation of Coefficients: (Intercept) area peri area -0.1709 peri -0.2675 -0.8390 shape -0.7513 -0.3396 0.5066 > rock.gam <- gam(log(perm) ~ s(area) + s(peri) + s(shape), + control=gam.control(maxit=50, bf.maxit=50)) > summary(rock.gam) Call: gam(formula = log(perm) ~ s(area) + s(peri) + s(shape), control = gam.control( maxit = 50, bf.maxit = 50)) Deviance Residuals: Min 1Q Median 3Q Max -1.685537 -0.4696152 0.1253111 0.5424773 1.296309 (Dispersion Parameter for Gaussian family taken to be 0.7446037 ) Null Deviance: 126.9322 on 47 degrees of freedom Residual Deviance: 26.06537 on 35.00569 degrees of freedom Number of Local Scoring Iterations: 1 DF for Terms and F-values for Nonparametric Effects Df Npar Df Npar F Pr(F) (Intercept) 1 s(area) 1 3 0.341688 0.7952340 s(peri) 1 3 0.931302 0.4358301 s(shape) 1 3 1.433113 0.2496551 > anova(rock.lm, rock.gam) Analysis of Variance Table Response: log(perm) Terms Resid. Df RSS Test Df Sum of Sq 1 area + peri + shape 44.00000 31.94883 2 s(area) + s(peri) + s(shape) 35.00569 26.06537 1 vs. 2 8.994309 5.88346 F Value Pr(F) 1 2 0.878496 0.5531102 > anova(rock.lm, rock.gam) Analysis of Variance Table Response: log(perm) Terms Resid. Df RSS Test Df Sum of Sq 1 area + peri + shape 44.00000 31.94883 2 s(area) + s(peri) + s(shape) 35.00569 26.06537 1 vs. 2 8.994309 5.88346 F Value Pr(F) 1 2 0.878496 0.5531102 > par(mfrow=c(2,2)) > plot(rock.gam, se=T) > rock.gam1 <- gam(log(perm) ~ area + peri + s(shape)) > par(mfrow=c(2,2)) > plot(rock.gam1, se=T) > anova(rock.lm, rock.gam1, rock.gam) Analysis of Variance Table Response: log(perm) Terms Resid. Df RSS Test Df Sum of Sq 1 area + peri + shape 44.00000 31.94883 2 area + peri + s(shape) 41.00273 29.00181 1 vs. 2 2.997272 2.947017 3 s(area) + s(peri) + s(shape) 35.00569 26.06537 2 vs. 3 5.997037 2.936443 F Value Pr(F) 1 2 1.320478 0.2832617 3 0.657597 0.6838482 > rock.ns <- lm(log(perm) ~ area + peri + ns(shape, df=4)) > summary(rock.ns) Call: lm(formula = log(perm) ~ area + peri + ns(shape, df = 4)) Residuals: Min 1Q Median 3Q Max -1.643 -0.6276 0.1352 0.5468 1.397 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) 4.4102 0.7209 6.1175 0.0000 area 0.0004 0.0001 4.9336 0.0000 peri -0.0015 0.0002 -8.0285 0.0000 ns(shape, df = 4)1 1.2148 0.7003 1.7345 0.0903 ns(shape, df = 4)2 1.3893 0.7349 1.8905 0.0658 ns(shape, df = 4)3 3.2205 1.5357 2.0971 0.0422 ns(shape, df = 4)4 0.4817 0.7109 0.6777 0.5018 Residual standard error: 0.8415 on 41 degrees of freedom Multiple R-Squared: 0.7713 F-statistic: 23.04 on 6 and 41 degrees of freedom, the p-value is 1.088e-11 Correlation of Coefficients: (Intercept) area peri ns(shape, df = 4)1 area -0.0921 peri -0.1533 -0.8530 ns(shape, df = 4)1 -0.7124 -0.2882 0.2872 ns(shape, df = 4)2 -0.5355 -0.3417 0.3946 0.3549 ns(shape, df = 4)3 -0.8525 -0.2070 0.2431 0.6968 ns(shape, df = 4)4 -0.3711 -0.2033 0.3374 0.4592 ns(shape, df = 4)2 ns(shape, df = 4)3 area peri ns(shape, df = 4)1 ns(shape, df = 4)2 ns(shape, df = 4)3 0.6535 ns(shape, df = 4)4 0.0353 0.3855 > plot.gam(rock.ns, se=T) > attach(bwt) > age1 <- age*(ftv=="1"); age2 <- age*(ftv=="2+") > birthwt.gam <- gam(low ~ s(age) + s(lwt) + smoke + ptd + + ht + ui + ftv + s(age1) + s(age2) + smoke:ui, binomial, + bf.maxit=25) > summary(birthwt.gam) Call: gam(formula = low ~ s(age) + s(lwt) + smoke + ptd + ht + ui + ftv + s(age1) + s( age2) + smoke:ui, family = binomial, bf.maxit = 25) Deviance Residuals: Min 1Q Median 3Q Max -2.02435 -0.7179398 -0.4531256 0.7629086 2.210088 (Dispersion Parameter for Binomial family taken to be 1 ) Null Deviance: 234.672 on 188 degrees of freedom Residual Deviance: 170.3536 on 165.1782 degrees of freedom Number of Local Scoring Iterations: 7 DF for Terms and Chi-squares for Nonparametric Effects Df Npar Df Npar Chisq P(Chi) (Intercept) 1 s(age) 1 3.0 3.108912 0.3722996 s(lwt) 1 2.9 2.339239 0.4853169 smoke 1 ptd 1 ht 1 ui 1 ftv 2 s(age1) 1 3.0 3.250398 0.3465486 s(age2) 1 3.0 3.147244 0.3682947 smoke:ui 1 > table(low, predict(birthwt.gam) > 0) FALSE TRUE 0 115 15 1 28 31 > plot(birthwt.gam, ask=T, se=T) > detach()