library(faraway) #box-cox and profile likelihood require(MASS) data(savings, package="faraway") lmod <- lm(sr ~ pop15+pop75+dpi+ddpi,savings) boxcox(lmod, plotit=T) boxcox(lmod, plotit=T, lambda=seq(0.5,1.5,by=0.1)) #box-cox on a different dataset data(gala, package="faraway") lmod <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,gala) boxcox(lmod, lambda=seq(-0.25,0.75,by=0.05),plotit=T) #different transformation lmod <- lm(burntime ~ nitrogen+chlorine+potassium, leafburn) logtrans(lmod,plotit=TRUE, alpha=seq(-min(leafburn$burntime)+0.001,0,by=0.01)) #broken stick regression lmod1 <- lm(sr ~ pop15, savings, subset=(pop15 < 35)) lmod2 <- lm(sr ~ pop15, savings, subset=(pop15 > 35)) plot(sr ~ pop15,savings,xlab="Pop'n under 15", ylab="Savings Rate") abline(v=35,lty=5) segments(20,lmod1$coef[1]+lmod1$coef[2]*20,35, lmod1$coef[1]+lmod1$coef[2]*35) segments(48,lmod2$coef[1]+lmod2$coef[2]*48,35, lmod2$coef[1]+lmod2$coef[2]*35) lhs <- function(x) ifelse(x < 35,35-x,0) rhs <- function(x) ifelse(x < 35,0,x-35) lmod <- lm(sr ~ lhs(pop15) + rhs(pop15), savings) x <- seq(20,48,by=1) py <- lmod$coef[1] + lmod$coef[2]*lhs(x) + lmod$coef[3]*rhs(x) lines(x,py,lty=2) #polynomials summary(lm(sr ~ ddpi,savings)) summary(lm(sr ~ ddpi+I(ddpi^2),savings)) summary(lm(sr ~ ddpi+I(ddpi^2)+I(ddpi^3),savings)) #what does a shift do to parameter estimates? savings <- data.frame(savings,mddpi=savings$ddpi-10) summary(lm(sr ~ mddpi+I(mddpi^2),savings)) #ortogonal polynomials lmod <- lm(sr ~ poly(ddpi,4),savings) sumary(lmod) lmod <- lm(sr ~ polym(pop15,ddpi,degree=2),savings) #multidimensional polynomial pop15r <- seq(20, 50, len=10) ddpir <- seq(0, 20, len=10) pgrid <- expand.grid(pop15=pop15r, ddpi=ddpir) pv <- predict(lmod, pgrid) persp(pop15r, ddpir, matrix(pv, 10, 10), theta=45, xlab="Pop under 15", ylab="Growth", zlab = "Savings rate", ticktype="detailed", shade = 0.25) funky <- function(x) sin(2*pi*x^3)^3 x <- seq(0,1,by=0.01) y <- funky(x) + 0.1*rnorm(101) matplot(x,cbind(y,funky(x)),type="pl",ylab="y",pch=20, lty=1, col=1) #fit different polynomials - not too good g4 <- lm(y ~ poly(x,4)) g12 <- lm(y ~ poly(x,12)) matplot(x,cbind(y,g4$fit,g12$fit),type="pll", ylab="y",lty=c(1,2), pch=20, col=1) #splines require(splines) #regression spline knots <- c(0,0,0,0,0.2,0.4,0.5,0.6,0.7,0.8,0.85,0.9,1,1,1,1) bx <- splineDesign(knots,x) lmodb <- lm(y ~ bx -1) matplot(x, bx, type="l", col=1) matplot(x, cbind(y,lmodb$fit), type="pl", ylab="y", pch=20,lty=1, col=1) #smoothing spline ssf <- smooth.spline(x,y) matplot(x,cbind(y,ssf$y),type="pl",ylab="y", lty=1, pch=20, col=1) #additive models require(mgcv) gamod <- gam(sr ~ s(pop15) + s(pop75) + s(dpi) + s(ddpi), data=savings) plot(gamod)