library('MASS') y <- mcycle$accel x <- mcycle$times ## smooth splines fit_ss <- smooth.spline(x, y, df=11) var_hat <- var(residuals(fit_ss)) x_plot <- seq(min(x),max(x),length.out=100) y_hat_plot <- predict(fit_ss, x_plot)$y ## extract smoothing spline matrix using special inputs ## remember S only depends on 'x' and 'df' smooth.matrix = function(x, df){ n = length(x); A = matrix(0, n, n); for(i in 1:n){ y = rep(0, n); y[i]=1; yi = predict(smooth.spline(x, y, df=df),x)$y; A[,i]= yi; } return(A) } S <- smooth.matrix(x, fit_ss$df) S_plot <- smooth.matrix(x_plot, fit_ss$df) y_hat_plot_vcov <- var_hat*S_plot%*%t(S_plot) plot(x, y, xlab="Time (ms)", ylab="Acceleration (g)", main="Smoothing Spline Fit") lines(x_plot, y_hat_plot, col="#882255", lwd=2) polygon(x=c(x_plot, rev(x_plot)), y=c(y_hat_plot + 1.96*sqrt(diag(y_hat_plot_vcov)), rev(y_hat_plot - 1.96*sqrt(diag(y_hat_plot_vcov)))), col="#882255AA", border=NA) legend('topleft',paste0("d.f. = ", round(fit_ss$df,1)), bty='n') image(S[ncol(S):1,], xaxt='n', yaxt='n', col=rainbow(10), main="Smoother Matrix")