library('MASS') library('dplyr') y <- mcycle$accel x <- mcycle$times ## piecewise constant with 4 knots basis_expand_pc <- function(x) cbind(1, x>15, x>25, x>30, x>40) ## piecewise linear with 4 knots basis_expand_pl <- function(x) cbind(basis_expand_pc(x), x*basis_expand_pc(x)) ## piecewise quadratic with 4 knots basis_expand_pq <- function(x) cbind(basis_expand_pl(x), x^2*basis_expand_pc(x)) ## piecewise cubic with 4 knots basis_expand_cu <- function(x) cbind(basis_expand_pq(x), x^3*basis_expand_pc(x)) ## linear spline with 4 knots sp <- function(x,e,p=1) (x>e)*(x-e)^p basis_expand_sp <- function(x) cbind(1, x, sp(x,15,1), sp(x,25,1), sp(x,30,1), sp(x,40,1)) ## quadratic spline with 4 knots basis_expand_qs <- function(x) cbind(1, x, x^2, sp(x,15,2), sp(x,25,2), sp(x,30,2), sp(x,40,2)) ## cubic spline with 4 knots basis_expand_cs <- function(x) cbind(1, x, x^2, x^3, sp(x,15,3), sp(x,25,3), sp(x,30,3), sp(x,40,3)) ## natural cubic spline with 4 knots d <- function(x, e1, e2) (sp(x, e1, 3)-sp(x, e2, 3))/(e2-e1) basis_expand_ns <- function(x) cbind(1, x, d(x,15,40)-d(x,30,40), d(x,25,40)-d(x,30,40)) ## fit basis expansion and plot fit fit_basis <- function(y, x, basis_expand) { ## expand inputs x_e <- basis_expand(x) colnames(x_e) <- paste0('x',1:ncol(x_e)) ## fit linear model in expanded variables dat_e <- cbind(y, x_e) %>% as.data.frame fit_e <- lm(y ~ 0 + ., data=dat_e) ## get estimate of beta beta_e <- coef(fit_e) ## compute vcov for beta_e beta_e_vcov <- vcov(fit_e) ## plot data, fit, and pointwise 95% CB plot(x, y, xlab="Time (ms)", ylab="Acceleration (g)") abline(v=c(15,25,30,40), lty=3) x_plot <- seq(min(x),max(x),length.out=100) x_plot_e <- basis_expand(x_plot) colnames(x_plot_e) <- paste0('x',1:ncol(x_plot_e)) dat_plot_e <- as.data.frame(x_plot_e) pre_plot_e <- predict(fit_e, dat_plot_e, se.fit=T) y_hat_plot <- pre_plot_e$fit y_hat_plot_lo <- pre_plot_e$fit - pre_plot_e$se.fit * 1.96 y_hat_plot_hi <- pre_plot_e$fit + pre_plot_e$se.fit * 1.96 lines(x_plot, y_hat_plot, col="#882255", lwd=2) polygon(x=c(x_plot, rev(x_plot)), y=c(y_hat_plot_lo, rev(y_hat_plot_hi)), col="#882255AA", border=NA) } ## piecewise constant fit_basis(y, x, basis_expand_pc) ## piecewise linear fit_basis(y, x, basis_expand_pl) ## piecewise quadratic fit_basis(y, x, basis_expand_pq) ## piecewise cubic fit_basis(y, x, basis_expand_cu) ## linear spline fit_basis(y, x, basis_expand_sp) ## quadratic spline fit_basis(y, x, basis_expand_qs) ## cubic spline fit_basis(y, x, basis_expand_cs) ## natural cubic spline fit_basis(y, x, basis_expand_ns) ## fit basis expansion and plot fit with splines functions library('splines') fit_basis_splines <- function(y, x, form = form_dflt) { ## expand inputs ## fit linear model in expanded variables dat_e <- cbind(y, x) %>% as.data.frame fit_e <- lm(form, data=dat_e) ## get estimate of beta beta_e <- coef(fit_e) ## compute vcov for beta_e beta_e_vcov <- vcov(fit_e) ## plot data, fit, and pointwise 95% CB plot(x, y, xlab="Time (ms)", ylab="Acceleration (g)") abline(v=c(15,25,30,40), lty=3) x_plot_e <- seq(min(x),max(x),length.out=100) dat_plot_e <- data.frame(x=x_plot_e) pre_plot_e <- predict(fit_e, dat_plot_e, se.fit=T) y_hat_plot <- pre_plot_e$fit y_hat_plot_lo <- pre_plot_e$fit - pre_plot_e$se.fit * 1.96 y_hat_plot_hi <- pre_plot_e$fit + pre_plot_e$se.fit * 1.96 lines(x_plot_e, y_hat_plot, col="#882255", lwd=2) polygon(x=c(x_plot_e, rev(x_plot_e)), y=c(y_hat_plot_lo, rev(y_hat_plot_hi)), col="#882255AA", border=NA) } ## natural cubis spline (same as basis_expand_ns above) form_cs <- y ~ bs(x, knots=c(15,25,30,40), degree=3) fit_basis_splines(y, x, form_cs) fit_basis(y, x, basis_expand_cs)