calcSmoothTrends: Smooth Basis Functions for a STdata Object

Description

A front end function for calling SVDsmooth (and SVDsmoothCV), with either a STdata object or vectors containing observations, dates and locations.

Usage

calcSmoothTrends(STdata = NULL, obs = STdata$obs$obs,
  date = STdata$obs$date, ID = STdata$obs$ID, subset = NULL,
  extra.dates = NULL, n.basis = 2, cv = FALSE, ...)

Arguments

STdata

A STdata/STmodel data structure containing observations, see mesa.data.raw. Use either this or the obs, date, and ID inputs.

obs

A vector of observations.

date

A vector of observation times.

A vector of observation locations.

subset

A subset of locations to extract the data matrix for. A warning is given for each name not found in ID.

extra.dates

Additional dates for which smooth trends should be computed (any duplicates will be removed).

n.basis

Number of basis functions to compute, see SVDsmooth.

Also compute smooth functions using leave one out cross-validation, see SVDsmoothCV.

...

Additional parameters passed to SVDsmooth and SVDsmoothCV; except fnc, which is always TRUE.

Value

Returns a list with

trend

A data.frame containing the smooth trends and the dates. This can be used as the trend in STdata$trend.

trend.cv

If cv==TRUE a list of data.frames; each one containing the smooth trend obtained when leaving one site out. Similar to SVDsmoothCV(data)$smoothSVD[[1]]).

trend.fnc,trend.fnc.cv

Functions that produce the content of the above data.frames, see SVDsmooth.

Details

The function uses createDataMatrix to create a data matrix which is passed to SVDsmooth (and SVDsmoothCV). The output can be used as STdata$trend = calcSmoothTrends(...)$trend, or STdata$trend = calcSmoothTrends(...)$trend.cv[[i]]. However, it is recommended to use updateTrend.STdata.

Examples

Run this code

# NOT RUN {
##let's load some data
data(mesa.model)

##let's compute two smooth trend functions
trend <- calcSmoothTrends(mesa.model, n.basis=2)

##or with some other parameters for the splines
trend.alt <- calcSmoothTrends(mesa.model, n.basis=2, df=100)

##and study the trends
par(mfrow=c(2,1), mar=c(2.5,2.5,.5,.5))
plot(trend$trend$date, trend$trend$V1, type="l", ylab="", xlab="",
     ylim=range(c(trend$trend$V1, trend$trend$V2)))
lines(trend$trend$date, trend$trend$V2, col=2)
plot(trend.alt$trend$date, trend.alt$trend$V1,
     type="l", ylab="", xlab="",
     ylim=range(c(trend.alt$trend$V1, trend.alt$trend$V2)))
lines(trend.alt$trend$date, trend.alt$trend$V2, col=2)

##Let's exclude locations with fewer than 100 observations
IND <- names(which(table(mesa.model$obs$ID) >= 100))
##now we also compute the CV trends.
trend2 <- calcSmoothTrends(mesa.model, n.basis=2, subset=IND, cv=TRUE)

##Let's compare to the previous result
lines(trend2$trend$date, trend2$trend$V1, lty=2)
lines(trend2$trend$date, trend2$trend$V2, lty=2, col=2)

##we can also study the cross validated results to examine the
##possible variation in the estimated trends.
plot(trend$trend$date, trend2$trend$V1, type="n", ylab="", xlab="",
     ylim=range(c(trend2$trend$V1, trend2$trend$V2)))
for(i in 1:length(trend2$trend.cv)){
  lines(trend2$trend.cv[[i]]$date, trend2$trend.cv[[i]]$V1, col=1)
  lines(trend2$trend.cv[[i]]$date, trend2$trend.cv[[i]]$V2, col=2)
}

##trend functions for every week (mesa.model$obs$date is every 2-weekss)
trend.more <- calcSmoothTrends(mesa.model, n.basis=2,
                               extra.dates=seq(min(mesa.model$obs$date),
                                 max(mesa.model$obs$date), by=7))
##This results in a message detailing how many times that
##have hade interpolated trends (i.e. no data)

##compare to the earlier
plot(trend$trend$date, trend$trend$V1, pch=19,
     ylim=range(c(trend$trend$V1, trend.more$trend$V1)))
points(trend.more$trend$date, trend.more$trend$V1, col=2, pch=3)
# }

Run the code above in your browser using DataCamp Workspace