- 1.1 Creating matrices in R
- 1.2 Matrix multiplication, addition and transpose @@ -435,13 +461,13 @@
- Authors
- Citation -
- Preface
- 1 Basic matrix math in R
- 1.1 Creating matrices in R
- 1.2 Matrix multiplication, addition and transpose @@ -435,7 +461,7 @@
Applied Time Series Analysis for Fisheries and Environmental Sciences
-2019-03-28
+2019-11-26
@@ -467,29 +494,29 @@Preface
@@ -443,8 +469,8 @@Preface
Book package
The book uses a number of R packages and a variety of fisheries data sets. The packages and data sets can be installed by installing our atsalibrary package which is hosted on GitHub:
-+library(devtools) -devtools::install_github("nwfsc-timeseries/atsalibrary")Citation
-Holmes, E. E., M. D. Scheuerell, and E. J. Ward. Applied time series analysis for fisheries and environmental data. NOAA Fisheries, Northwest Fisheries Science Center, 2725 Montlake Blvd E., Seattle, WA 98112. Contacts eli.holmes@noaa.gov, eric.ward@noaa.gov, and mark.scheuerell@noaa.gov
+Holmes, E. E., M. D. Scheuerell, and E. J. Ward. Applied time series analysis for fisheries and environmental data. NOAA Fisheries, Northwest Fisheries Science Center, 2725 Montlake Blvd E., Seattle, WA 98112. Contacts eli.holmes@noaa.gov, eric.ward@noaa.gov, and mark.scheuerell@noaa.gov
+Citation
- Authors
- Citation
11.2 Example: Snotel Data
Let’s see an example using the Washington SNOTEL data. The data we will use is the snow water equivalent percent of normal. This represents the snow water equivalent compared to the average value for that site on the same day. We will look at a subset of sites in the Central Cascades in our snotel dataset (Figure 11.1).
-load("snotel.RData")
-y <- snotelmeta
-# Just use a subset
-y = y[which(y$Longitude < -121.4), ]
-y = y[which(y$Longitude > -122.5), ]
-y = y[which(y$Latitude < 47.5), ]
-y = y[which(y$Latitude > 46.5), ]
+
+y <- snotelmeta
+# Just use a subset
+y = y[which(y$Longitude < -121.4),]
+y = y[which(y$Longitude > -122.5),]
+y = y[which(y$Latitude < 47.5),]
+y = y[which(y$Latitude > 46.5),]
@@ -451,11 +477,10 @@ 11.2 Example: Snotel Data
For the first analysis, we are just going to look at February Snow Water Equivalent (SWE). Our subset of stations is y$Station.Id. There are many missing years among some of our stations (Figure 11.2).
-swe.feb <- snotel
-swe.feb <- swe.feb[swe.feb$Station.Id %in% y$Station.Id & swe.feb$Month ==
- "Feb", ]
-p <- ggplot(swe.feb, aes(x = Date, y = SWE)) + geom_line()
-p + facet_wrap(~Station)
+swe.feb <- snotel
+swe.feb <- swe.feb[swe.feb$Station.Id %in% y$Station.Id & swe.feb$Month=="Feb",]
+p <- ggplot(swe.feb, aes(x=Date, y=SWE)) + geom_line()
+p + facet_wrap(~Station)
@@ -465,7 +490,7 @@
11.2 Example: Snotel Data
11.2.1 Estimate Feb SWE using AR(1) with spatial correlation
Imagine that for our study we need an estimate of SWE for all sites. We will use the information from the sites with full data to estimate the missing SWE for other sites. We will use a MARSS model to use all the available data.
-\[\begin{equation}
+\[\begin{equation}
\begin{bmatrix}
x_1 \\ x_2 \\ \dots \\ x_{15}
\end{bmatrix}_t =
@@ -494,36 +519,35 @@ 11.2.1 Estimate Feb SWE using AR(
v_1 \\ v_2 \\ \dots \\ v_{15}
\end{bmatrix}_t
\tag{11.5}
-\end{equation}\]
+\end{equation}\]
We will use an unconstrained variance-covariance structure for \(\mathbf{w}\) and assume that \(\mathbf{v}\) is identical and independent and very low (SNOTEL instrument variability). The \(a_i\) determine the level of the \(x_i\).
We need our data to be in rows. We will use reshape2::acast().
-dat.feb <- reshape2::acast(swe.feb, Station ~ Year, value.var = "SWE")
+
We set up the model for MARSS so that it is the same as (11.5). We will fix the measurement error to be small; we could use 0 but the fitting is more stable if we use a small variance instead. When estimating \(\mathbf{B}\), setting the initial value to be at \(t=1\) instead of \(t=0\) works better.
-ns <- length(unique(swe.feb$Station))
-B <- "diagonal and equal"
-Q <- "unconstrained"
-R <- diag(0.01, ns)
-U <- "zero"
-A <- "unequal"
-x0 <- "unequal"
-mod.list.ar1 = list(B = B, Q = Q, R = R, U = U, x0 = x0, A = A,
- tinitx = 1)
+ns <- length(unique(swe.feb$Station))
+B <- "diagonal and equal"
+Q <- "unconstrained"
+R <- diag(0.01,ns)
+U <- "zero"
+A <- "unequal"
+x0 <- "unequal"
+mod.list.ar1 = list(B=B, Q=Q, R=R, U=U, x0=x0, A=A, tinitx=1)
Now we can fit a MARSS model and get estimates of the missing SWEs. Convergence is slow. We set \(\mathbf{a}\) equal to the mean of the time series to speed convergence.
-library(MARSS)
-m <- apply(dat.feb, 1, mean, na.rm = TRUE)
-fit.ar1 <- MARSS(dat.feb, model = mod.list.ar1, control = list(maxit = 5000),
- inits = list(A = matrix(m, ns, 1)))
-The \(b\) estimate is ````.
+library(MARSS)
+m <- apply(dat.feb, 1, mean, na.rm=TRUE)
+fit.ar1 <- MARSS(dat.feb, model=mod.list.ar1, control=list(maxit=5000), inits=list(A=matrix(m,ns,1)))
+The \(b\) estimate is 0.4494841.
Let’s plot the estimated SWEs for the missing years (Figure 11.3). These estimates use all the information about the correlation with other sites and uses information about correlation with the prior and subsequent years. We will use the tidy() function from the broom package to get the estimated 95% confidence intervals for the estimated states. Notice that for some sites, CIs are low in early years as these sites are highly correlated with site for which there are data. In other sites, the uncertainty is high in early years because the sites with data in those years are not highly correlated.
-fit <- fit.ar1
-d <- augment(fit, interval = "confidence")
-d$Year <- d$t + 1980
-d$Station <- d$.rownames
-p <- ggplot(data = d) + geom_line(aes(Year, .fitted)) + geom_ribbon(aes(x = Year,
- ymin = .conf.low, ymax = .conf.up), linetype = 2, alpha = 0.5)
-p <- p + geom_point(data = swe.feb, mapping = aes(x = Year, y = SWE))
-p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")
+fit <- fit.ar1
+d <- augment(fit, interval="confidence")
+d$Year <- d$t + 1980
+d$Station <- d$.rownames
+p <- ggplot(data = d) +
+ geom_line(aes(Year, .fitted)) +
+ geom_ribbon(aes(x=Year, ymin=.conf.low, ymax=.conf.up), linetype=2, alpha=0.5)
+p <- p + geom_point(data=swe.feb, mapping = aes(x=Year, y=SWE))
+p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")
@@ -536,9 +560,9 @@
11.2.1 Estimate Feb SWE using AR(
11.2.1.1 Diagnostics
The state residuals have a tendency for negative autocorrelation at lag-1 (Figure 11.4).
-fit <- fit.ar1
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30], 1, acf)
+fit <- fit.ar1
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30], 1, acf)
@@ -550,7 +574,7 @@
11.2.1.1 Diagnostics
11.2.2 Estimate Feb SWE using only correlation
Another approach is to treat the February data as temporally uncorrelated. The two longest time series (Paradise and Olallie Meadows) show minimal autocorrelation so we might decide to just use the correlation across stations for our estimates. In this case, the state of the missing SWE values at time \(t\) is the expected value conditioned on all the stations with data at time \(t\) given the estimated variance-covariance matrix \(\mathbf{Q}\).
-We could set this model up as
+We could set this model up as
\[\begin{equation}
\begin{bmatrix}
y_1 \\ y_2 \\ \dots \\ y_{15}
@@ -568,9 +592,9 @@ 11.2.2 Estimate Feb SWE using onl
\zeta_{15,1}&\zeta_{15,2}&\dots&\sigma_{15}
\end{bmatrix}
\tag{11.6}
-\end{equation}\]
+\end{equation}\]
However the EM algorithm used by MARSS() runs into numerical issues. Instead we will set the model up as follows. Allowing a hidden state observed with small error makes the estimation more stable.
-\[\begin{equation}
+\[\begin{equation}
\begin{bmatrix}
x_1 \\ x_2 \\ \dots \\ x_{15}
\end{bmatrix}_t =
@@ -605,30 +629,29 @@ 11.2.2 Estimate Feb SWE using onl
\end{bmatrix}
\tag{11.7}
\end{equation}\]
-
Again \(\mathbf{a}\) is the mean level in the time series. Note that the expected value of \(\mathbf{x}\) is zero if there are no data, so \(E(\mathbf{x}_0)=0\).
-ns <- length(unique(swe.feb$Station))
-B <- "zero"
-Q <- "unconstrained"
-R <- diag(0.01, ns)
-U <- "zero"
-A <- "unequal"
-x0 <- "zero"
-mod.list.corr = list(B = B, Q = Q, R = R, U = U, x0 = x0, A = A,
- tinitx = 0)
+Again \(\mathbf{a}\) is the mean level in the time series. Note that the expected value of \(\mathbf{x}\) is zero if there are no data, so \(E(\mathbf{x}_0)=0\).
+ns <- length(unique(swe.feb$Station))
+B <- "zero"
+Q <- "unconstrained"
+R <- diag(0.01,ns)
+U <- "zero"
+A <- "unequal"
+x0 <- "zero"
+mod.list.corr = list(B=B, Q=Q, R=R, U=U, x0=x0, A=A, tinitx=0)
Now we can fit a MARSS model and get estimates of the missing SWEs. Convergence is slow. We set \(\mathbf{a}\) equal to the mean of the time series to speed convergence.
-m <- apply(dat.feb, 1, mean, na.rm = TRUE)
-fit.corr <- MARSS(dat.feb, model = mod.list.corr, control = list(maxit = 5000),
- inits = list(A = matrix(m, ns, 1)))
+m <- apply(dat.feb, 1, mean, na.rm=TRUE)
+fit.corr <- MARSS(dat.feb, model=mod.list.corr, control=list(maxit=5000), inits=list(A=matrix(m,ns,1)))
The estimated SWEs for the missing years uses the information about the correlation with other sites only.
-fit <- fit.corr
-d <- broom::augment(fit, interval = "confidence")
-d$Year <- d$t + 1980
-d$Station <- d$.rownames
-p <- ggplot(data = d) + geom_line(aes(Year, .fitted)) + geom_ribbon(aes(x = Year,
- ymin = .conf.low, ymax = .conf.up), linetype = 2, alpha = 0.5)
-p <- p + geom_point(data = swe.feb, mapping = aes(x = Year, y = SWE))
-p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")
+fit <- fit.corr
+d <- broom::augment(fit, interval="confidence")
+d$Year <- d$t + 1980
+d$Station <- d$.rownames
+p <- ggplot(data = d) +
+ geom_line(aes(Year, .fitted)) +
+ geom_ribbon(aes(x=Year, ymin=.conf.low, ymax=.conf.up), linetype=2, alpha=0.5)
+p <- p + geom_point(data=swe.feb, mapping = aes(x=Year, y=SWE))
+p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")
@@ -638,15 +661,15 @@
11.2.2 Estimate Feb SWE using onl
11.2.2.1 Diagnostics
The state and model residuals have no tendency towards negative autocorrelation now that we removed the autoregressive component from the process (\(x\)) model.
-fit <- fit.corr
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30], 1, acf)
-mtext("State Residuals ACF", outer = TRUE, side = 3)
+fit <- fit.corr
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30], 1, acf)
+mtext("State Residuals ACF", outer=TRUE, side=3)
-fit <- fit.corr
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$model.residuals[, 1:30], 1, acf)
-mtext("Model Residuals ACF", outer = TRUE, side = 3)
+fit <- fit.corr
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$model.residuals[,1:30], 1, acf)
+mtext("Model Residuals ACF", outer=TRUE, side=3)
@@ -681,40 +704,34 @@ 11.2.3 Estimate Feb SWE using DFA
\end{bmatrix}_t
\]
The model is set up as follows:
-ns <- dim(dat.feb)[1]
-B <- matrix(list(0), 2, 2)
-B[1, 1] <- "b1"
-B[2, 2] <- "b2"
-Q <- diag(1, 2)
-R <- "diagonal and unequal"
-U <- "zero"
-x0 <- "zero"
-Z <- matrix(list(0), ns, 2)
-Z[1:(ns * 2)] <- c(paste0("z1", 1:ns), paste0("z2", 1:ns))
-Z[1, 2] <- 0
-A <- "unequal"
-mod.list.dfa = list(B = B, Z = Z, Q = Q, R = R, U = U, A = A,
- x0 = x0)
+ns <- dim(dat.feb)[1]
+B <- matrix(list(0),2,2)
+B[1,1] <- "b1"; B[2,2] <- "b2"
+Q <- diag(1,2)
+R <- "diagonal and unequal"
+U <- "zero"
+x0 <- "zero"
+Z <- matrix(list(0),ns,2)
+Z[1:(ns*2)] <- c(paste0("z1",1:ns),paste0("z2",1:ns))
+Z[1,2] <- 0
+A <- "unequal"
+mod.list.dfa = list(B=B, Z=Z, Q=Q, R=R, U=U, A=A, x0=x0)
Now we can fit a MARSS model and get estimates of the missing SWEs. We pass in the initial value for \(\mathbf{a}\) as the mean level so it fits easier.
-library(MARSS)
-m <- apply(dat.feb, 1, mean, na.rm = TRUE)
-fit.dfa <- MARSS(dat.feb, model = mod.list.dfa, control = list(maxit = 1000),
- inits = list(A = matrix(m, ns, 1)))
+library(MARSS)
+m <- apply(dat.feb, 1, mean, na.rm=TRUE)
+fit.dfa <- MARSS(dat.feb, model=mod.list.dfa, control=list(maxit=1000), inits=list(A=matrix(m,ns,1)))
11.2.4 Diagnostics
The state residuals are uncorrelated.
-fit <- fit.dfa
-par(mfrow = c(1, 2), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30, drop = FALSE], 1,
- acf)
+fit <- fit.dfa
+par(mfrow=c(1,2),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30,drop=FALSE], 1, acf)
As are the model residuals:
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$model.residual, 1, function(x) {
- acf(na.omit(x))
-})
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$model.residual, 1, function(x){acf(na.omit(x))})
@@ -729,24 +746,24 @@ 11.2.4 Diagnostics
+
+
@@ -41,6 +40,9 @@
+
+
+
@@ -51,41 +53,66 @@
@@ -109,7 +136,6 @@
-
+
+
+
@@ -41,6 +40,9 @@
+
+
+
@@ -51,41 +53,66 @@
@@ -109,7 +136,6 @@
11.2 Example: Snotel Data
Let’s see an example using the Washington SNOTEL data. The data we will use is the snow water equivalent percent of normal. This represents the snow water equivalent compared to the average value for that site on the same day. We will look at a subset of sites in the Central Cascades in our snotel dataset (Figure 11.1).
load("snotel.RData")y <- snotelmeta
-# Just use a subset
-y = y[which(y$Longitude < -121.4), ]
-y = y[which(y$Longitude > -122.5), ]
-y = y[which(y$Latitude < 47.5), ]
-y = y[which(y$Latitude > 46.5), ]y <- snotelmeta
+# Just use a subset
+y = y[which(y$Longitude < -121.4),]
+y = y[which(y$Longitude > -122.5),]
+y = y[which(y$Latitude < 47.5),]
+y = y[which(y$Latitude > 46.5),]11.2 Example: Snotel Data
For the first analysis, we are just going to look at February Snow Water Equivalent (SWE). Our subset of stations is y$Station.Id. There are many missing years among some of our stations (Figure 11.2).
swe.feb <- snotel
-swe.feb <- swe.feb[swe.feb$Station.Id %in% y$Station.Id & swe.feb$Month ==
- "Feb", ]
-p <- ggplot(swe.feb, aes(x = Date, y = SWE)) + geom_line()
-p + facet_wrap(~Station)swe.feb <- snotel
+swe.feb <- swe.feb[swe.feb$Station.Id %in% y$Station.Id & swe.feb$Month=="Feb",]
+p <- ggplot(swe.feb, aes(x=Date, y=SWE)) + geom_line()
+p + facet_wrap(~Station)@@ -465,7 +490,7 @@
11.2 Example: Snotel Data
11.2.1 Estimate Feb SWE using AR(1) with spatial correlation
Imagine that for our study we need an estimate of SWE for all sites. We will use the information from the sites with full data to estimate the missing SWE for other sites. We will use a MARSS model to use all the available data.
-\[\begin{equation} +\[\begin{equation}
\begin{bmatrix}
x_1 \\ x_2 \\ \dots \\ x_{15}
\end{bmatrix}_t =
@@ -494,36 +519,35 @@ 11.2.1 Estimate Feb SWE using AR(
v_1 \\ v_2 \\ \dots \\ v_{15}
\end{bmatrix}_t
\tag{11.5}
-\end{equation}\]
+\end{equation}\]
We will use an unconstrained variance-covariance structure for \(\mathbf{w}\) and assume that \(\mathbf{v}\) is identical and independent and very low (SNOTEL instrument variability). The \(a_i\) determine the level of the \(x_i\).
We need our data to be in rows. We will use reshape2::acast().
dat.feb <- reshape2::acast(swe.feb, Station ~ Year, value.var = "SWE")We set up the model for MARSS so that it is the same as (11.5). We will fix the measurement error to be small; we could use 0 but the fitting is more stable if we use a small variance instead. When estimating \(\mathbf{B}\), setting the initial value to be at \(t=1\) instead of \(t=0\) works better.
-ns <- length(unique(swe.feb$Station))
-B <- "diagonal and equal"
-Q <- "unconstrained"
-R <- diag(0.01, ns)
-U <- "zero"
-A <- "unequal"
-x0 <- "unequal"
-mod.list.ar1 = list(B = B, Q = Q, R = R, U = U, x0 = x0, A = A,
- tinitx = 1)ns <- length(unique(swe.feb$Station))
+B <- "diagonal and equal"
+Q <- "unconstrained"
+R <- diag(0.01,ns)
+U <- "zero"
+A <- "unequal"
+x0 <- "unequal"
+mod.list.ar1 = list(B=B, Q=Q, R=R, U=U, x0=x0, A=A, tinitx=1)Now we can fit a MARSS model and get estimates of the missing SWEs. Convergence is slow. We set \(\mathbf{a}\) equal to the mean of the time series to speed convergence.
-library(MARSS)
-m <- apply(dat.feb, 1, mean, na.rm = TRUE)
-fit.ar1 <- MARSS(dat.feb, model = mod.list.ar1, control = list(maxit = 5000),
- inits = list(A = matrix(m, ns, 1)))The \(b\) estimate is ````.
+library(MARSS)
+m <- apply(dat.feb, 1, mean, na.rm=TRUE)
+fit.ar1 <- MARSS(dat.feb, model=mod.list.ar1, control=list(maxit=5000), inits=list(A=matrix(m,ns,1)))The \(b\) estimate is 0.4494841.
Let’s plot the estimated SWEs for the missing years (Figure 11.3). These estimates use all the information about the correlation with other sites and uses information about correlation with the prior and subsequent years. We will use the tidy() function from the broom package to get the estimated 95% confidence intervals for the estimated states. Notice that for some sites, CIs are low in early years as these sites are highly correlated with site for which there are data. In other sites, the uncertainty is high in early years because the sites with data in those years are not highly correlated.
fit <- fit.ar1
-d <- augment(fit, interval = "confidence")
-d$Year <- d$t + 1980
-d$Station <- d$.rownames
-p <- ggplot(data = d) + geom_line(aes(Year, .fitted)) + geom_ribbon(aes(x = Year,
- ymin = .conf.low, ymax = .conf.up), linetype = 2, alpha = 0.5)
-p <- p + geom_point(data = swe.feb, mapping = aes(x = Year, y = SWE))
-p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")fit <- fit.ar1
+d <- augment(fit, interval="confidence")
+d$Year <- d$t + 1980
+d$Station <- d$.rownames
+p <- ggplot(data = d) +
+ geom_line(aes(Year, .fitted)) +
+ geom_ribbon(aes(x=Year, ymin=.conf.low, ymax=.conf.up), linetype=2, alpha=0.5)
+p <- p + geom_point(data=swe.feb, mapping = aes(x=Year, y=SWE))
+p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")@@ -536,9 +560,9 @@
11.2.1 Estimate Feb SWE using AR(
11.2.1.1 Diagnostics
The state residuals have a tendency for negative autocorrelation at lag-1 (Figure 11.4).
-fit <- fit.ar1
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30], 1, acf)fit <- fit.ar1
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30], 1, acf)@@ -550,7 +574,7 @@
11.2.1.1 Diagnostics
11.2.2 Estimate Feb SWE using only correlation
Another approach is to treat the February data as temporally uncorrelated. The two longest time series (Paradise and Olallie Meadows) show minimal autocorrelation so we might decide to just use the correlation across stations for our estimates. In this case, the state of the missing SWE values at time \(t\) is the expected value conditioned on all the stations with data at time \(t\) given the estimated variance-covariance matrix \(\mathbf{Q}\).
-We could set this model up as +We could set this model up as
\[\begin{equation}
\begin{bmatrix}
y_1 \\ y_2 \\ \dots \\ y_{15}
@@ -568,9 +592,9 @@ 11.2.2 Estimate Feb SWE using onl
\zeta_{15,1}&\zeta_{15,2}&\dots&\sigma_{15}
\end{bmatrix}
\tag{11.6}
-\end{equation}\]
+\end{equation}\]
However the EM algorithm used by MARSS() runs into numerical issues. Instead we will set the model up as follows. Allowing a hidden state observed with small error makes the estimation more stable.
\[\begin{equation}
\begin{bmatrix}
x_1 \\ x_2 \\ \dots \\ x_{15}
\end{bmatrix}_t =
@@ -605,30 +629,29 @@ 11.2.2 Estimate Feb SWE using onl
\end{bmatrix}
\tag{11.7}
\end{equation}\]
-
Again \(\mathbf{a}\) is the mean level in the time series. Note that the expected value of \(\mathbf{x}\) is zero if there are no data, so \(E(\mathbf{x}_0)=0\).
-ns <- length(unique(swe.feb$Station))
-B <- "zero"
-Q <- "unconstrained"
-R <- diag(0.01, ns)
-U <- "zero"
-A <- "unequal"
-x0 <- "zero"
-mod.list.corr = list(B = B, Q = Q, R = R, U = U, x0 = x0, A = A,
- tinitx = 0)ns <- length(unique(swe.feb$Station))
+B <- "zero"
+Q <- "unconstrained"
+R <- diag(0.01,ns)
+U <- "zero"
+A <- "unequal"
+x0 <- "zero"
+mod.list.corr = list(B=B, Q=Q, R=R, U=U, x0=x0, A=A, tinitx=0)Now we can fit a MARSS model and get estimates of the missing SWEs. Convergence is slow. We set \(\mathbf{a}\) equal to the mean of the time series to speed convergence.
-m <- apply(dat.feb, 1, mean, na.rm = TRUE)
-fit.corr <- MARSS(dat.feb, model = mod.list.corr, control = list(maxit = 5000),
- inits = list(A = matrix(m, ns, 1)))m <- apply(dat.feb, 1, mean, na.rm=TRUE)
+fit.corr <- MARSS(dat.feb, model=mod.list.corr, control=list(maxit=5000), inits=list(A=matrix(m,ns,1)))The estimated SWEs for the missing years uses the information about the correlation with other sites only.
-fit <- fit.corr
-d <- broom::augment(fit, interval = "confidence")
-d$Year <- d$t + 1980
-d$Station <- d$.rownames
-p <- ggplot(data = d) + geom_line(aes(Year, .fitted)) + geom_ribbon(aes(x = Year,
- ymin = .conf.low, ymax = .conf.up), linetype = 2, alpha = 0.5)
-p <- p + geom_point(data = swe.feb, mapping = aes(x = Year, y = SWE))
-p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")fit <- fit.corr
+d <- broom::augment(fit, interval="confidence")
+d$Year <- d$t + 1980
+d$Station <- d$.rownames
+p <- ggplot(data = d) +
+ geom_line(aes(Year, .fitted)) +
+ geom_ribbon(aes(x=Year, ymin=.conf.low, ymax=.conf.up), linetype=2, alpha=0.5)
+p <- p + geom_point(data=swe.feb, mapping = aes(x=Year, y=SWE))
+p + facet_wrap(~Station) + xlab("") + ylab("SWE (demeaned)")@@ -638,15 +661,15 @@
11.2.2 Estimate Feb SWE using onl
11.2.2.1 Diagnostics
The state and model residuals have no tendency towards negative autocorrelation now that we removed the autoregressive component from the process (\(x\)) model.
-fit <- fit.corr
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30], 1, acf)
-mtext("State Residuals ACF", outer = TRUE, side = 3)
+fit <- fit.corr
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30], 1, acf)
+mtext("State Residuals ACF", outer=TRUE, side=3)
-fit <- fit.corr
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$model.residuals[, 1:30], 1, acf)
-mtext("Model Residuals ACF", outer = TRUE, side = 3)
+fit <- fit.corr
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$model.residuals[,1:30], 1, acf)
+mtext("Model Residuals ACF", outer=TRUE, side=3)
11.2.2.1 Diagnostics
The state and model residuals have no tendency towards negative autocorrelation now that we removed the autoregressive component from the process (\(x\)) model.
-fit <- fit.corr
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30], 1, acf)
-mtext("State Residuals ACF", outer = TRUE, side = 3)fit <- fit.corr
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30], 1, acf)
+mtext("State Residuals ACF", outer=TRUE, side=3)
fit <- fit.corr
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$model.residuals[, 1:30], 1, acf)
-mtext("Model Residuals ACF", outer = TRUE, side = 3)fit <- fit.corr
+par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$model.residuals[,1:30], 1, acf)
+mtext("Model Residuals ACF", outer=TRUE, side=3)
11.2.3 Estimate Feb SWE using DFA \end{bmatrix}_t \]
The model is set up as follows:
-ns <- dim(dat.feb)[1]
-B <- matrix(list(0), 2, 2)
-B[1, 1] <- "b1"
-B[2, 2] <- "b2"
-Q <- diag(1, 2)
-R <- "diagonal and unequal"
-U <- "zero"
-x0 <- "zero"
-Z <- matrix(list(0), ns, 2)
-Z[1:(ns * 2)] <- c(paste0("z1", 1:ns), paste0("z2", 1:ns))
-Z[1, 2] <- 0
-A <- "unequal"
-mod.list.dfa = list(B = B, Z = Z, Q = Q, R = R, U = U, A = A,
- x0 = x0)ns <- dim(dat.feb)[1]
+B <- matrix(list(0),2,2)
+B[1,1] <- "b1"; B[2,2] <- "b2"
+Q <- diag(1,2)
+R <- "diagonal and unequal"
+U <- "zero"
+x0 <- "zero"
+Z <- matrix(list(0),ns,2)
+Z[1:(ns*2)] <- c(paste0("z1",1:ns),paste0("z2",1:ns))
+Z[1,2] <- 0
+A <- "unequal"
+mod.list.dfa = list(B=B, Z=Z, Q=Q, R=R, U=U, A=A, x0=x0)Now we can fit a MARSS model and get estimates of the missing SWEs. We pass in the initial value for \(\mathbf{a}\) as the mean level so it fits easier.
-library(MARSS)
-m <- apply(dat.feb, 1, mean, na.rm = TRUE)
-fit.dfa <- MARSS(dat.feb, model = mod.list.dfa, control = list(maxit = 1000),
- inits = list(A = matrix(m, ns, 1)))library(MARSS)
+m <- apply(dat.feb, 1, mean, na.rm=TRUE)
+fit.dfa <- MARSS(dat.feb, model=mod.list.dfa, control=list(maxit=1000), inits=list(A=matrix(m,ns,1)))
11.2.4 Diagnostics
The state residuals are uncorrelated.
-fit <- fit.dfa
-par(mfrow = c(1, 2), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$state.residuals[, 1:30, drop = FALSE], 1,
- acf)fit <- fit.dfa
+par(mfrow=c(1,2),mar=c(2,2,1,1))
+apply(residuals(fit)$state.residuals[,1:30,drop=FALSE], 1, acf)
As are the model residuals:
-par(mfrow = c(4, 4), mar = c(2, 2, 1, 1))
-apply(residuals(fit)$model.residual, 1, function(x) {
- acf(na.omit(x))
-})par(mfrow=c(4,4),mar=c(2,2,1,1))
+apply(residuals(fit)$model.residual, 1, function(x){acf(na.omit(x))})
It is easy enough to get the expected value of 
@@ -521,24 +548,24 @@ 
@@ -536,7 +568,8 @@ 
@@ -558,24 +591,24 @@ 
@@ -512,16 +533,16 @@ 
@@ -562,24 +590,24 @@ 
@@ -494,24 +520,24 @@ 
@@ -460,24 +486,24 @@ 
+
@@ -464,24 +490,24 @@ 
@@ -462,24 +488,24 @@ 
@@ -481,26 +507,26 @@ 
@@ -510,32 +536,32 @@ 
@@ -557,24 +583,24 @@ 
@@ -472,22 +498,22 @@ 
@@ -539,16 +565,16 @@ 
@@ -562,30 +588,30 @@ 
@@ -595,41 +621,41 @@ 
@@ -638,8 +664,8 @@ 
@@ -661,24 +687,24 @@ 
@@ -492,24 +518,24 @@ 
@@ -548,20 +574,20 @@ 
@@ -570,7 +596,7 @@ 
@@ -475,11 +501,11 @@ 
@@ -501,24 +527,24 @@ 
@@ -492,26 +518,26 @@ 
@@ -533,24 +559,24 @@ 
@@ -470,7 +496,7 @@ 
@@ -523,24 +549,24 @@ 
@@ -459,22 +485,22 @@ 
@@ -495,24 +521,24 @@ 
@@ -475,12 +501,12 @@ 
@@ -502,24 +528,24 @@ 
@@ -452,14 +478,14 @@ 
@@ -480,24 +506,24 @@ 
@@ -446,14 +472,13 @@ 
@@ -481,24 +507,24 @@