######################################################################### ### Illustration of smoothing cardiovascular mortality series from ### Shumway & Stoffer (2006) [http://www.stat.pitt.edu/stoffer/tsa2/textRcode.htm] ### The data are average weekly cardiovascular mortality in Los Angeles ### County over the 10-year period 1970-1979. ######################################################################### # read data and convert to a time series object mort = ts(scan("Data/cmort.dat"),start=1970, frequency=52) # # Moving Average smooothers: ma5 = filter(mort, sides=2, rep(1,5)/5) ma53 = filter(mort, sides=2, rep(1,53)/53) plot(mort, type="p", xlab="week", ylab="mortality") lines(ma5) # yearly smooth (trend) lines(ma53) # monthly smooth (trend + seasonality) # # Polynomial and Harmonic Regression smoothers: wk = time(mort) - mean(time(mort)) # week (centered) - wk is basically t/52 ... wk2 = wk^2 # ... because frequency=52 for mort wk3 = wk^3 c = cos(2*pi*wk) s = sin(2*pi*wk) reg1 = lm(mort~wk + wk2 + wk3, na.action=NULL) reg2 = lm(mort~wk + wk2 + wk3 + c + s, na.action=NULL) plot(mort, type="p", xlab="week", ylab="mortality") lines(fitted(reg1)) # trend lines(fitted(reg2)) # trend + seasonality # # Kernel Smoothing: plot(mort, type="p", xlab="week", ylab="mortality") lines(ksmooth(time(mort), mort, "normal", bandwidth=2)) # trend lines(ksmooth(time(mort), mort, "normal", bandwidth=10/52)) # trend + seasonality # # Lowess smoothing: plot(mort, type="p", xlab="week", ylab="mortality", main="lowess") lines(lowess(mort, f=.02)) # trend + seasonality lines(lowess(mort, f=2/3)) # trend # # Splines smoothing: plot(mort, type="p", xlab="week", ylab="mortality") lines(smooth.spline(time(mort), mort, spar=0)) # trend + seasonality lines(smooth.spline(time(mort), mort, spar=1)) # trend