SpatioTemporal (version 1.1.2)

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. Use either this or the obs, date, and ID inputs.

obs

A vector of observations.

date

A vector of observation times.

ID

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.

n.basis

Number of basis functions to compute, see SVDsmooth.

cv

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

...

Additional parameters passed to SVDsmooth and SVDsmoothCV.

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]]).

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 updateSTdataTrend.

See Also

Other SVD for missing data: plot.SVDcv, print.SVDcv, SVDmiss, SVDsmooth, SVDsmoothCV, updateSTdataTrend

Examples

##let's load some data
data(mesa.data)

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

##or with some other parameters for the splines
trend.alt <- calcSmoothTrends(mesa.data, 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.data$obs$ID) >= 100))
##now we also compute the CV trends.
trend2 <- calcSmoothTrends(mesa.data, 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.data$obs$date is every 2-weekss)
trend.more <- calcSmoothTrends(mesa.data, n.basis=2,
                               extra.dates=seq(min(mesa.data$obs$date),
                                 max(mesa.data$obs$date), by=7))
##This results in a message detailing 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)


##Let's use our new trend functions in the mesa.data
##NOTE: IT IS PREFERABLE TO USE updateSTdataTrend INSTEAD!!!
mesa.data$trend <- trend2$trend