Solved – Repeated loess smoothing for time series data

I have 3 time series and I would like to determine if there are any common patterns in the data. Zuur mentioned in one of his books that a reasonably simple approach to test this is based on repeated loess smoothing. Zuur states that to visualize the general patterns of the time series, we fit a loess smoother through each series and then plot them all on the same graph. This generally captures the long term variations, to capture the short term variations (which I am interested in) we have to apply loess smoothing again, but now on the residuals.

Below, is an example data set where we have hourly values of ozone concentrations for 3 countries, I would like to perform what Zuur has explained, but to this example data set. Could anyone offer some suggestions? Please note that my knowledge of R and statistics in general is low, I apologize if this is a simplistic question.

require(plyr) require(lattice)  TopFolder <- list(""                   ,""                   ,"" ) ## create variable for the data  data = ldply(TopFolder, header = TRUE, read.table, sep = "", skip = 3) ## define Ozone levels Ozone <- data$Value     Ozone[Ozone==-999] <- NA     Ozone <- data.frame(Ozone)     ## define Datetime -  need to concatenate arrays     DateTime <- paste(data$Date,data$Hour, sep = " ")     Date <- as.POSIXct(DateTime, format = "%d.%m.%Y %H:%M")     ## define countries      Countries <- c("Country1","Country2","Country3")     Country <- data.frame(Country = rep(Countries, each = 8760))     ## bind together     Dat <- cbind(Ozone, Country = Country)     Dat <- transform(Dat, Date=Date, Doy = as.numeric(format(Date,format = "%j")),                      Tod = as.numeric(format(Date,format = "%H")),                      DecTime = rep(seq(1,365, length = 8760),by = 3))     ## Select data for a month during the summer      NewDat <- subset(Dat, as.Date(Dat$Date) >= '2009-07-01 00:00:00' & as.Date(Dat$Date) <= '2009-08-01 00:00:00')   ## plot the data xyplot(Ozone~DecTime | Country, data = NewDat, type  = "l", col = 1,        strip = function(bg = 'white',...)strip.default(bg = 'white',...)) 

So, here I have imported the data and taken a subset (for one month), and then plotted the values. I would now like to apply loess smoothing to represent the short term variation i.e. hourly.

Try this:

PlotData <- lapply(Countries,function(country) {   fit1 <- loess(Ozone~DecTime ,data = NewDat[NewDat$Country==country,],na.action=na.exclude)    fit2 <- loess(residuals(fit1)~DecTime ,data = NewDat[NewDat$Country==country,],na.action=na.exclude)    data.frame(DecTime=NewDat[NewDat$Country==country,"DecTime"],SmoothResiduals=predict(fit2),Country=country)   })  PlotData <-'rbind',PlotData)  xyplot(SmoothResiduals~DecTime | Country, data = PlotData, type  = "l", col = 1,        strip = function(bg = 'white',...)strip.default(bg = 'white',...)) 

Smoothed residuals

Explanation: lapply is used to iterate over the countries, applies a function and combines the function values in a list. The anonymous function does a loess fit (note my use of na.action=na.exclude to deal with NA values), does another loess fit on the residuals and returns a data.frame with the smoothed residuals. Subsequently, I combine all data.frames and create the plot.

Obviously you might want to adjust the parameters of the loess fits.

Similar Posts:

Rate this post

Leave a Comment