Chap4_Lab3.Rnw

\documentclass[b5paper,11pt]{book}
\input{preamble}

\ifthenelse{\equal{\jobname}{\detokenize{Chap4_Lab3}}}
{\standalonetrue}
{\standalonefalse}


\begin{document}

<<echo=FALSE>>=
opts_chunk$set(background = rgb(1,1,0.85),
               size='small')
@

\addcontentsline{toc}{section}{Lab 4.3: Temporal Basis Functions with {\bf SpatioTemporal}}
\section*{Lab 4.3: Temporal Basis Functions with {\bf SpatioTemporal}}
\markright{Lab 4.3: Temporal Basis Functions with {\bf SpatioTemporal}}

In this Lab we model the maximum temperature in the NOAA data set (\cc{Tmax}) using temporal basis functions and spatial random fields. Specifically, we use the model
\begin{equation}
Y(\bs;t) = \bx(\bs;t)' \bfbeta + \sum_{i=1}^{n_\alpha} \phi_{i}(t) \alpha_i(\bs) + \nu(\bs;t),
\label{eq:Yprocess_tbasis2}
\end{equation}
\noindent where $\bx(\bs;t)$ are the covariates; $\bfbeta$ are the regression coefficients; $\{\phi_i(t)\}$ are the tem\-poral basis functions; $\{\alpha_i(\bs)\}$ are coefficients of the temporal basis functions, modeled as multivariate (spatial) random fields; and $\nu(\bs;t)$ is a spatially correlated, but temporally independent, random process.

Spatio-temporal modeling using temporal basis functions can be carried out using the package {\bf SpatioTemporal}. For this Lab we require the following packages.

<<echo=TRUE,message=FALSE,results='hide',warning=FALSE>>=
library("dplyr")
library("ggplot2")
library("gstat")
library("RColorBrewer")
library("sp")
library("spacetime")
library("SpatioTemporal")
library("STRbook")
library("tidyr")
@

<<echo=FALSE,message=FALSE,results='hide'>>=
library("grid")
library("gridExtra")
@


The space-time object used by {\bf SpatioTemporal} is of class \cc{STdata} and is created using the function \fn{createSTdata}. This function takes the data either as a space-wide matrix with the row names containing the date and the column names the station ID, or as a data frame in long form. Here we use the latter. This data frame needs to have the station ID as characters in the field \cc{ID}, the data in the field \cc{obs}, and the date in the field \cc{date}. A new data frame of this form can be easily created using the function \fn{transmute} from the package {\bf dplyr}.

<<warning=FALSE>>=
data("NOAA_df_1990", package = "STRbook")   # load NOAA data
NOAA_sub <- filter(NOAA_df_1990,       # filter data to only
                   year == 1993 &      # contain July 1993
                   month == 7 &
                   proc == "Tmax")     # and just max. temp.

NOAA_sub_for_STdata <- NOAA_sub %>%
                       transmute(ID = as.character(id),
                                 obs = z,
                                 date = date)
@

The covariates that will be used to model the spatially varying effects also need to be supplied as a data frame. In our case we only consider the station coordinates as covariates. The station coordinates are extracted from the maximum temperature data as follows.

<<>>=
covars <-  dplyr::select(NOAA_sub, id, lat, lon) %>%
           unique() %>%
           dplyr::rename(ID = id)      # createSTdata expects "ID"
@

\noindent Now we can construct the \cc{STdata} object by calling the function \fn{createSTdata}.

<<message=FALSE,results='hide', warning = FALSE>>=
STdata <- createSTdata(NOAA_sub_for_STdata, covars = covars)
@

The model used in {\bf SpatioTemporal} assumes that $\nu(\bs;t)$ is temporally uncorrelated. Consequently, all temporal variability needs to be captured through the covariates or the basis functions. To check whether the data exhibit temporal autocorrelation (before adding any temporal basis functions), one can use the \fn{plot} function. For example, we plot the estimated autocorrelation function for station 3812 in the left panel of Figure \ref{fig:ACFs} (after the mean is removed from the data). The plot suggests that the data are correlated (the estimated lag-1 autocorrelation coefficient is larger than would be expected by chance at the 5\% level of significance).

<<eval = FALSE>>=
plot(STdata, "acf", ID = "3812")
@

The role of the temporal basis functions is to adequately capture temporal modes of variation. When modeling data over a time interval that spans years, one of these is typically a seasonal component. As another example, when modeling trace-gas emissions, one basis function to use would be one that captures weekday/weekend cycles typically found in gaseous pollutants (e.g., due to vehicular traffic). The package {\bf SpatioTemporal} allows for user-defined basis functions (see the example at the end of this Lab) or data-driven basis functions (which we consider now). In both cases, the first temporal basis function, $\phi_1(t)$, is a constant; that is, $\phi_1(t) = 1$.

The basis functions extracted from the data are \emph{smoothed}, \emph{left singular vectors} (i.e., smoothed temporal EOFs) of the matrix $\widetilde{\bZ}$\ifstandalone. \else, described in Technical Note \ref{technote:EOFs}. \fi These make up the remaining $n_\alpha - 1$ basis functions, upon which smoothing is carried out using splines. In {\bf SpatioTemporal}, these basis functions are found (or set) using the function \fn{updateTrend}.

<<SpatioTemporal_setup,results='hide',message=FALSE>>=
STdata <- updateTrend(STdata, n.basis = 2)
@

<<echo=FALSE, results = 'hide', message=FALSE>>=
## See how many stations are no longer lag-1 correlated
STdata1 <- createSTdata(NOAA_sub_for_STdata, covars = covars)
STdata2 <- updateTrend(STdata, n.basis = 2)
DFcheck <- left_join(NOAA_sub_for_STdata, STdata2$trend)
Insig1 <- Insig2 <- NULL
for(i in unique(STdata1$obs$ID)) {
  ACF1 <- acf(residuals(lm(obs ~ 1, data = filter(DFcheck, ID == i))), plot = 0)
  ACF2 <- acf(residuals(lm(obs ~ V1 + V2, data = filter(DFcheck, ID == i))), plot = 0)
  Insig1[[i]] <- ACF1$acf[2] < qnorm((1 + 0.95)/2)/sqrt(ACF1$n.used) # Taken from plot.act
  Insig2[[i]] <- ACF2$acf[2] < qnorm((1 + 0.95)/2)/sqrt(ACF2$n.used)
}
Sig1 <- 100 - round(mean(unlist(Insig1)),2)*100
Sig2 <- 100 - round(mean(unlist(Insig2)),2)*100
@


\noindent We can see that the lag-1 autocorrelation coefficient is no longer significant (at the 5\% level) after adding in these basis functions; see the right panel of Figure \ref{fig:ACFs}. In practice, one should add basis functions until temporal autocorrelation in the data (at most stations) is considerably reduced. In this case study, it can be shown that \Sexpr{Sig1}\% of stations record maximum temperature data that have lag-1 autocorrelation coefficients that are significant at the 5\% level. On the other hand, with \args{n.basis}\cc{ = }\num{2} (i.e., with two temporal basis functions for capturing temporal variation),  the proportion of stations with residuals exhibiting a significant lag-1 autocorrelation coefficient is \Sexpr{Sig2}\%.

<<eval = FALSE>>=
plot(STdata, "acf", ID = "3812")
@

<<echo=FALSE, results='hide', message=FALSE>>=
png(filename = "img/Chapter_4/ACF1.png",width = 1500,height = 1500, res = 300)
STdata <- createSTdata(NOAA_sub_for_STdata, covars = covars)
plot(STdata, "acf", ID = "3812")
dev.off()

png(filename = "img/Chapter_4/ACF2.png",width = 1500,height = 1500, res = 300)
STdata <- updateTrend(STdata, n.basis = 2)
plot(STdata, "acf", ID = "3812")
dev.off()
@


\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.49\textwidth,angle=0]{img/Chapter_4/ACF1.png}
\includegraphics[width=0.49\textwidth,angle=0]{img/Chapter_4/ACF2.png}
\end{center}
\caption{Left: Estimated autocorrelation function for the time series of maximum temperature \cc{Tmax} at Station 3812. Right: Same as left panel, but with the data first detrended using an intercept and the two temporal basis functions shown in the top panel of Figure~\ref{fig:temp_basis_fn}.}
\label{fig:ACFs}
\end{figure}




\noindent The basis functions, available in \cc{STdata\$trend}, are shown in the top panel of Figure \ref{fig:temp_basis_fn}.


<<echo=FALSE>>=
data_to_plot <- gather(STdata$trend,trend, z, -date)
data_to_plot <- rbind(data=data.frame(date = as.Date(c("1993-07-01", "1993-07-31")),
                                      trend = "V0",
                                      z = c(1, 1)),
                      data_to_plot)
LineTrends <- ggplot(data_to_plot) +
    geom_line(aes(x=date,y=z,group=trend,linetype=trend)) +
    scale_linetype_manual(labels=c(expression(varphi[1](t)),
                                   expression(varphi[2](t)),
                                   expression(varphi[3](t))),
                          values=c(1,3,4),name="basis") +
    ylab("") + theme_bw()
ggsave(LineTrends, file="./img/Chapter_4/temp_basis_fns.png",height=2,width=6)
@

In {\bf SpatioTemporal}, the spatial quantities $\{\alpha_i(\bs)\}$ are themselves modeled as spatial fields. Once the $\{\phi_i(t)\}$ are declared, empirical estimates of $\{\alpha_i(\bs)\}$ can be found using the function \fn{estimateBetaFields}. Note that we use the Greek letter ``alpha'' to denote these fields, which differs from the name ``Beta'' inside the command. The following and all subsequent references to ``Beta'' and ``beta'' should be interpreted as representing spatial fields $\{\alpha_i(\bs)\}$.

<<SpatioTemporal_beta_empirical,echo=TRUE>>=
beta.lm <- estimateBetaFields(STdata)
@

\noindent The resulting object, \cc{beta.lm}, contains two fields; \cc{beta} (estimated coefficients) and \cc{beta.sd} (standard error of the estimates) with row names equal to the station ID, and three columns corresponding to estimates of $\alpha_1(\bs)$, $\alpha_2(\bs)$, and $\alpha_3(\bs)$, respectively. We are interested in seeing whether the empirical estimates are correlated with our covariate, latitude. To this end, the authors of {\bf SpatioTemporal} suggest using the package {\bf plotrix}, and the function \fn{plotCI}, to plot the estimates and covariance intervals against a covariate of choice. When plotting using \fn{plotCI}, care should be taken that the ordering of the stations in \cc{beta} and \cc{beta.sd} is the same as that if the covariate data frame. For example, consider

<<>>=
head(row.names(beta.lm$beta))
head(covars$ID)
@

\noindent This illustrates a discrepancy, since the ordering of strings is not necessarily that of the ordered integers. For this reason we recommend employing best practice and always merging (e.g., using \fn{left\_join}) on a column variable; in this case, we choose the integer version of the field \cc{ID}. In the following commands, we first convert the \cc{beta} and \cc{beta.sd} objects into data frames, add the column \cc{ID}, join into a data frame \cc{BETA}, and then combine with \cc{covars} containing the latitude data.

<<>>=
beta.lm$beta <- data.frame(beta.lm$beta)
beta.lm$beta.sd <- data.frame(beta.lm$beta.sd)
beta.lm$beta$ID <- as.integer(row.names(beta.lm$beta))
BETA <- cbind(beta.lm$beta, beta.lm$beta.sd)
colnames(BETA) <- c("alpha1", "alpha2", "alpha3", "ID",
                    "alpha1_CI", "alpha2_CI", "alpha3_CI")
BETA <- left_join(BETA, covars, by = "ID")
@

\noindent Once \cc{BETA} is constructed, the empirical estimates can be plotted using \fn{ggplot}, with \fn{geom\_errorbar} to also plot error bars, as follows.

<<eval=FALSE>>=
ggplot(BETA) + geom_point(aes(x = lat, y = alpha1)) +
    geom_errorbar(aes(x = lat,
                      ymin = alpha1 - 1.96*alpha1_CI,
                      ymax = alpha1 + 1.96*alpha1_CI)) +
    ylab(expression(alpha[1](s))) +
    xlab("lat (deg)") + theme_bw()
@

The three empirical estimates, plotted as a function of latitude, are shown in Figure \ref{fig:alpha_emp_ests}. The function $\alpha_1(\bs)$ exhibits a strong latitudinal trend, as expected; $\alpha_2(\bs)$ shows a weak latitudinal trend; and $\alpha_3(\bs)$ exhibits no trend. For this reason we model the expectations of these fields as \ifstandalone
\begin{align}
E(\alpha_1(\bs)) &= \alpha_{11} + \alpha_{12}s_2,\label{eq:alpha1}\\
E(\alpha_2(\bs)) &= \alpha_{21} + \alpha_{22}s_2,\label{eq:alpha2} \\
E(\alpha_3(\bs)) &= \alpha_{31},\label{eq:alpha3} \end{align}
where $s_2$ denotes the latitude coordinate at $\bs = (s_1, s_2)'$.
\else in \eqref{eq:alpha1}--\eqref{eq:alpha3}. \fi  Note that in this model we do not consider any spatio-temporal covariates, and hence the term $\bx(\bs;t)'\bfbeta = 0$ in \eqref{eq:Yprocess_tbasis2}. This does not mean that we do not have an intercept in our model: although it is random, the spatial field $\alpha_1(\bs)$ acts as a temporally invariant spatial covariate and includes a global space-time mean ($\alpha_{11}$ in \eqref{eq:alpha1}), which is estimated.

<<echo=FALSE,results='hide',fig.keep='none'>>=
g1 <- ggplot(BETA) + geom_point(aes(x=lat,y=alpha1)) +
    geom_errorbar(aes(x = lat,
                      ymin = alpha1 - 1.96*alpha1_CI,
                      ymax = alpha1 + 1.96*alpha1_CI)) +
    ylab(expression(alpha[1](s))) +
    xlab("lat (deg)") +
    theme_bw()

g2 <- ggplot(BETA) + geom_point(aes(x=lat,y=alpha2)) +
    geom_errorbar(aes(x = lat,
                      ymin = alpha2 - 1.96*alpha2_CI,
                      ymax = alpha2 + 1.96*alpha2_CI)) +
    ylab(expression(alpha[2](s))) +
    xlab("lat (deg)") +
    theme_bw()


g3 <- ggplot(BETA) + geom_point(aes(x=lat,y=alpha3)) +
    geom_errorbar(aes(x = lat,
                      ymin = alpha3 - 1.96*alpha3_CI,
                      ymax = alpha3 + 1.96*alpha3_CI)) +
    ylab(expression(alpha[3](s))) +
    xlab("lat (deg)") +
    theme_bw()

g <- grid.arrange(g1,g2,g3,nrow=1)
ggsave(g, file="./img/Chapter_4/alpha_emp_ests.png",height=2.8,width=8)
@

\ifstandalone
\begin{figure}[t!]
\begin{center}
\includegraphics[width=\textwidth,angle=0]{img/Chapter_4/temp_basis_fns.png}
\includegraphics[width=0.3\textwidth,angle=0]{img/Chapter_4/alpha_field1.png}
\includegraphics[width=0.3\textwidth,angle=0]{img/Chapter_4/alpha_field2.png}
\includegraphics[width=0.3\textwidth,angle=0]{img/Chapter_4/alpha_field3.png}
\end{center}
\caption{Top: Basis functions $\phi_1(t), \phi_2(t)$, and $\phi_3(t)$, where the latter two were obtained from the left-singular vectors following a singular value decomposition of the data matrix (i.e., temporal EOFs). Bottom: $E(\alpha_1(\bs) \mid \bz), E(\alpha_2(\bs) \mid \bz)$, and $E(\alpha_3(\bs) \mid \bz)$.}
\label{fig:temp_basis_fn}
\end{figure}
\fi

\begin{figure}[t!]
\begin{center}
\includegraphics[width=\textwidth,angle=0]{img/Chapter_4/alpha_emp_ests.png}
\end{center}
\caption{Empirical estimates of $\alpha_1(\bs)$, $\alpha_2(\bs)$, and $\alpha_3(\bs)$  at each station, with 95\% confidence intervals, plotted as a function of latitude.}
\label{fig:alpha_emp_ests}
\end{figure}



We let the covariance functions $\cov(\alpha_i(\bs),\alpha_i(\bs + \bh))$, $i = 1,2,3$, be exponential covariance functions without a nugget-effect term. In {\bf SpatioTemporal} these are declared as follows.

<<SpatioTemporal_modelling>>=
cov.beta <- list(covf = "exp", nugget = FALSE)
@

All that remains for constructing the spatio-temporal model is to define the spatial covariance function of the zero-mean, temporally independent, residual process $\nu(\bs;t)$\ifstandalone. \else ; see  \eqref{eq:Yprocess_tbasis2}. \fi We choose this to be an exponential covariance function with a nugget effect to account for measurement error. The argument \args{random.effect} \cc{= } \num{FALSE} is used to indicate that there is no random mean offset for the field at each time point.

<<>>=
cov.nu <- list(covf = "exp",
               nugget = ~1,
               random.effect = FALSE) # No random mean
                                      # for each nu
@

The function to create the spatio-temporal model is \fn{createSTmodel}. This takes as data the object \cc{STdata}, the covariates for the $\alpha$-fields (an intercept and latitude for $\alpha_1(\bs)$ and $\alpha_2(\bs)$, and just an intercept for $\alpha_3(\bs)$; see \eqref{eq:alpha1}--\eqref{eq:alpha3}, the covariance functions of the $\alpha$-fields and the $\nu$-field, and a list containing the names of station coordinate fields (\cc{lon} and \cc{lat}).

<<>>=
locations <- list(coords = c("lon", "lat"))
LUR <- list(~lat, ~lat, ~1)  # lat trend for phi1 and phi2 only
STmodel <- createSTmodel(STdata,              # data
                         LUR = LUR,           # spatial covariates
                         cov.beta = cov.beta, # cov. of alphas
                         cov.nu = cov.nu,     # cov. of nu
                         locations = locations) # coord. names
@

In order to fit the spatio-temporal model to the data, we need to provide initial values of the parameter estimates. The required parameter names can be extracted using the function \fn{loglikeSTnames} and, for our model, are as follows.

<<SpatioTemporal_estimation,>>=
parnames <- loglikeSTnames(STmodel, all = FALSE)
print(parnames)
@

\noindent Noting that all parameters are log-transforms of the quantities of interest, we let all of the initial values be equal to 3 (so that all initial ranges and sills are $e^3 \approx 20$). This seems reasonable when the temperature is varying on the order of several degrees Fahrenheit, and where the domain also spans several degrees (in latitude and longitude).

We use the function \fn{estimate} below to fit the spatio-temporal model to the data. This may take several minutes on a standard desktop computer. In this instance, the resulting object \cc{SpatioTemporalfit1} has been pre-computed and can be loaded directly from {\bf STRbook} by typing \fn{data}\cc{(}\strn{"SpatioTemporalfit1"}\cc{,} \args{package}\cc{ = }\strn{"STRbook"}\cc{)}.



<<eval=FALSE>>=
x.init <- matrix(3, 9, 1)
rownames(x.init) <- loglikeSTnames(STmodel, all = FALSE)
SpatioTemporalfit1 <- estimate(STmodel, x.init)
@

<<eval=FALSE, echo=FALSE>>=
save(SpatioTemporalfit1, file="data/SpatioTemporalfit1.rda")
@

<<echo=FALSE>>=
data("SpatioTemporalfit1", package = "STRbook")
@


\noindent The fitted coefficients for the parameters described by \cc{parnames} above can be extracted from the fitted object using the function \fn{coef}.

<<>>=
x.final <- coef(SpatioTemporalfit1, pars = "cov")$par
@


Having fitted the model, we now predict at unobserved locations. First, we establish the spatial and temporal grid upon which to predict; this proceeds by first initializing an \cc{STdata} object on a grid. We construct the grid following a very similar approach to what was done in Lab 4.1.

<<>>=
## Define space-time grid
spat_pred_grid <- expand.grid(lon = seq(-100, -80, length = 20),
                      lat = seq(32, 46, length = 20))
spat_pred_grid$id <- 1:nrow(spat_pred_grid)
temp_pred_grid <- as.Date("1993-07-01") + seq(3, 28, length = 6)

## Initialize data matrix
obs_pred_wide <- matrix(0, nrow = 6, ncol = 400)

## Set row names and column names
rownames(obs_pred_wide) <- as.character(temp_pred_grid)
colnames(obs_pred_wide) <- spat_pred_grid$id

covars_pred <- spat_pred_grid                     # covariates
STdata_pred <- createSTdata(obs = obs_pred_wide,  # ST object
                            covars = covars_pred)
@

Now prediction proceeds using the function \fn{predict}, which requires as arguments the model, the fitted model parameters, and the data matrix \cc{STdata\_pred}.

<<>>=
E <- predict(STmodel, x.final, STdata = STdata_pred)
@

\noindent The returned object \cc{E} contains both the $\alpha$-fields predictions as well as the $Y$-field prediction at the unobserved locations. For example, \cc{E\$beta\$EX} contains the conditional expectations of $\alpha_1(\bs),\alpha_2(\bs)$, and $\alpha_3(\bs)$ given the data. For conciseness, we do not illustrate the plotting commands here. In the bottom panels of Figure \ref{fig:temp_basis_fn}, we show the conditional expectations, while in Figures \ref{fig:SpatioTemporal_pred} and \ref{fig:SpatioTemporal_se} we show the predictions and prediction standard errors of maximum temperature over six days of interest in July 1993.


<<echo=FALSE,results='hide',fig.keep = 'none', message=FALSE>>=
spat_pred_grid <- expand.grid(lon = seq(-100,-80,length=20),
                      lat = seq(32,46,length=20)) %>%
            SpatialPoints(proj4string=CRS("+proj=longlat"))
gridded(spat_pred_grid) <- TRUE
DE_pred <- STFDF(sp=spat_pred_grid,
               time=temp_pred_grid,
               data = data.frame(pred=c(t(E$EX)),
                                 se = sqrt(c(t(E$VX)))))
color_pal <- rev(colorRampPalette(brewer.pal(11, "Spectral"))(16))
png("./img/Chapter_4/SpatioTemporal_pred.png",width = 1600,height=1000,res=300)
stplot(DE_pred[,,1],main="Predictions (degrees Fahrenheit)",
       layout=c(3,2),
       col.regions=color_pal)
dev.off()
png("./img/Chapter_4/SpatioTemporal_se.png",width = 1600,height=1000,res=300)
stplot(DE_pred[,,2],main="Prediction errors (degrees Fahrenheit)",
       layout=c(3,2),
       col.regions=color_pal)
dev.off()
@

<<echo=FALSE,fig.keep='none',message=FALSE,results='hide'>>=
XX <- E$beta$EX
XX <- cbind(XX,covars_pred)
g0 <- ggplot(XX) + geom_raster(aes(lon,lat,fill=const)) + scale_fill_distiller(palette="Spectral",name=expression(alpha[1](s)))  + theme_bw() + xlab("lon (deg)") + ylab("lat (deg)") + coord_fixed(ratio=1.5)
g1 <- ggplot(XX) + geom_raster(aes(lon,lat,fill=V1)) + scale_fill_distiller(palette="Spectral",name=expression(alpha[2](s)))  + theme_bw() + xlab("lon (deg)") + ylab("lat (deg)") + coord_fixed(ratio=1.5)
g2 <- ggplot(XX) + geom_raster(aes(lon,lat,fill=V2)) + scale_fill_distiller(palette="Spectral",name=expression(alpha[3](s)))  + theme_bw() + xlab("lon (deg)") + ylab("lat (deg)") + coord_fixed(ratio=1.5)
ggsave(g0, file="./img/Chapter_4/alpha_field1.png",height=2.8,width=3)
ggsave(g1, file="./img/Chapter_4/alpha_field2.png",height=2.8,width=3)
ggsave(g2, file="./img/Chapter_4/alpha_field3.png",height=2.8,width=3)
@

\begin{figure}[t!]
\begin{center}
\includegraphics[width=\textwidth,angle=0]{img/Chapter_4/SpatioTemporal_pred.png}
\end{center}
\caption{Predictions of \cc{Tmax} in degrees Fahrenheit within a square lat-lon box defining the spatial domain of interest, for six days in July 1993, using temporal basis functions. Data for 14 July 1993 were deliberately omitted from the original data set.}
\label{fig:SpatioTemporal_pred}
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[width=\textwidth,angle=0]{img/Chapter_4/SpatioTemporal_se.png}
\end{center}
\caption{Prediction standard errors of \cc{Tmax} in degrees Fahrenheit within a square lat-lon box enclosing the spatial domain of interest, for six days in July 1993, using temporal basis functions. Data for 14 July 1993 were deliberately omitted from the original data set.}
\label{fig:SpatioTemporal_se}
\end{figure}

\subsection*{Using {\bf SpatioTemporal} for Modeling Spatial Effects of Temporal Covariates}

In the first part of this Lab, we extracted the temporal basis functions from the data. However, {\bf SpatioTemporal} can also be used to model the spatially varying effect of exo\-gen\-ous temporal covariates. This can be done by manually setting the \cc{STdata\$trend} data frame. When modeling temperature, interesting covariates may include a periodic signal with period equal to one year, or an index such as the El Ni{\~n}o Southern Oscillation (ENSO) Index.

To use a pre-existing covariate, we need to use the \cc{fnc} argument in \fn{updateTrend} to define a function that takes a \cc{Date} object as an input and returns the covariate at these dates. The easiest way to do this in this example is to specify a look-up table in the function containing the covariate for each date, but an interpolant can also be used when the covariate has missing information for one or more dates.

As an exercise, repeat the Lab above, but this time use a single linear temporal trend as a temporal covariate. The look-up table we need is just a two-column data frame containing the \cc{date} in the first column, and \cc{V1} (first covariate) in the second column. This can be set up as follows.

<<>>=
all_dates <- NOAA_sub$date %>% unique()     # dates
lookup <- data.frame(date = all_dates,      # covariate (linear)
                     V1 = scale(as.numeric(all_dates)))
@

Type \fn{plot}\cc{(lookup)} to see the temporal covariate that we have just created. Now we need to create the function that takes a \cc{Date} object as input and returns the required covariate values. This can be done using \fn{left\_join}.

<<>>=
## Function that returns the covariates in a data frame
## at the required dates
fnc <- function(dates) {
  left_join(data.frame(date = dates),
            lookup, by = "date") %>%
  select(-date)
}
@

\noindent Now we can call \fn{updateTrend} with our covariate function as argument.

<<message=FALSE,results='hide'>>=
STdata <- updateTrend(STdata, fnc = fnc)
@


The rest of the code remains largely similar, except that now we are considering only two temporal basis functions and not three (the first basis function is constant in time, and the second one is linear in time). Changing the required parts of the code is left as an exercise.


<<echo=FALSE,eval=FALSE>>=
data_to_plot <- gather(STdata$trend,trend, z, -date)
data_to_plot <- rbind(data=data.frame(date = as.Date(c("1993-07-01", "1993-07-31")),
                                      trend = "V0",
                                      z = c(1, 1)),
                      data_to_plot)
LineTrends <- ggplot(data_to_plot) +
    geom_line(aes(x=date,y=z,group=trend,linetype=trend)) +
    scale_linetype_manual(labels=c(expression(phi[1](t)),
                                   expression(phi[2](t))),
                          values=c(1,3,4),name="basis") +
    ylab("") + theme_bw()
beta.lm <- estimateBetaFields(STdata)
beta.lm$beta <- data.frame(beta.lm$beta)
beta.lm$beta.sd <- data.frame(beta.lm$beta.sd)
beta.lm$beta$ID <-as.integer(row.names(beta.lm$beta))
BETA <- cbind(beta.lm$beta,beta.lm$beta.sd)
colnames(BETA) <- c("alpha1","alpha2","ID",
                    "alpha1_CI","alpha2_CI")
BETA <- left_join(BETA, covars, by = "ID")
g1 <- ggplot(BETA) + geom_point(aes(x=lat,y=alpha1)) +
    geom_errorbar(aes(x = lat,
                      ymin = alpha1 - 1.96*alpha1_CI,
                      ymax = alpha1 + 1.96*alpha1_CI)) +
    ylab(expression(alpha[1](s))) +
    xlab("lat (deg)") +
    theme_bw()

g2 <- ggplot(BETA) + geom_point(aes(x=lat,y=alpha2)) +
    geom_errorbar(aes(x = lat,
                      ymin = alpha2 - 1.96*alpha2_CI,
                      ymax = alpha2 + 1.96*alpha2_CI)) +
    ylab(expression(alpha[2](s))) +
    xlab("lat (deg)") +
    theme_bw()
g <- grid.arrange(g1,g2,nrow=1)
cov.beta <- list(covf = "exp", nugget = FALSE)
cov.nu <- list(covf = "exp",
               nugget = ~1,
               random.effect = FALSE) # No random mean
                                      # for each nu
locations <- list(coords=c("lon", "lat"))
LUR <- list(~lat, ~1)   # latitude trend for phi1 and phi2 only
STmodel <- createSTmodel(STdata,              # data
                         LUR = LUR,           # spatial covariates
                         cov.beta = cov.beta, # cov. of alpha fields
                         cov.nu = cov.nu,     # cov. of nu
                         locations = locations) # coord. names
parnames <- loglikeSTnames(STmodel, all = FALSE)
print(parnames)
x.init <- matrix(3,7,1)
rownames(x.init) <- loglikeSTnames(STmodel, all = FALSE)
if(0) {
  SpatioTemporalfit2 <- estimate(STmodel, x.init)
  save(SpatioTemporalfit2, file="data/SpatioTemporalfit2.rda")
} else {
  data(SpatioTemporalfit2)
}
x.final <- coef(SpatioTemporalfit2, pars="cov")$par

## Define space-time grid
spat_pred_grid <- expand.grid(lon = seq(-100,-80,length=20),
                      lat = seq(32,46,length=20))
spat_pred_grid$id <- 1:nrow(spat_pred_grid)
temp_pred_grid <- as.Date("1993-07-01") + seq(3,28,length=6)

## Construct empty data matrix
obs_pred_wide <- matrix(0,nrow=6,ncol=400)    # empty data matrix
rownames(obs_pred_wide) <- as.character(temp_pred_grid)     # time as rownames
colnames(obs_pred_wide) <- spat_pred_grid$id  # spatial id as names

covars_pred <- spat_pred_grid              # covariates
STdata_pred <- createSTdata(obs = obs_pred_wide,
                            covars = covars_pred)
E <- predict(STmodel, x.final,STdata = STdata_pred)

spat_pred_grid <- expand.grid(lon = seq(-100,-80,length=20),
                      lat = seq(32,46,length=20)) %>%
            SpatialPoints(proj4string=CRS("+proj=longlat"))
gridded(spat_pred_grid) <- TRUE
DE_pred <- STFDF(sp=spat_pred_grid,
               time=temp_pred_grid,
               data = data.frame(pred=c(t(E$EX)),
                                 se = sqrt(c(t(E$VX)))))
color_pal <- rev(colorRampPalette(brewer.pal(11, "Spectral"))(16))
stplot(DE_pred[,,1],main="Predictions (degrees Fahrenheit)",
       layout=c(3,2),
       col.regions=color_pal)
stplot(DE_pred[,,2],main="Prediction errors (degrees Fahrenheit)",
       layout=c(3,2),
       col.regions=color_pal)
XX <- E$beta$EX
XX <- cbind(XX,covars_pred)
g0 <- ggplot(XX) + geom_raster(aes(lon,lat,fill=const)) + scale_fill_distiller(palette="Spectral",name=expression(alpha[1](s)))  + theme_bw() + xlab("lon (deg)") + ylab("lat (deg)") + coord_fixed(ratio=1.5)
g1 <- ggplot(XX) + geom_raster(aes(lon,lat,fill=V1)) + scale_fill_distiller(palette="Spectral",name=expression(alpha[2](s)))  + theme_bw() + xlab("lon (deg)") + ylab("lat (deg)") + coord_fixed(ratio=1.5)
@


\end{document}