Negative values for Q and R

[abby-caplan]

I'm fitting a DLM with TMB to pass to tmbstan, but am getting negative values for Q and R matrices.

Here's the data:
weekly_ne_anom.txt
weekly_sst_ne.txt

And here's the code:

B <- diag(2)
U <- matrix(0, 2, 1)
Q <- matrix("0", 2, 2)
diag(Q) <- c("q.alpha", "q.beta")
Z <- array(NA, c(1, 2, 1356))
Z[1, 1, ] <- rep(1, 1356)
Z[1, 2, ] <- na.approx(weekly_sst_ne$anom)
x0 <- matrix(c(0,0), nrow = 2)

model1 <- list(B = B, U = U, Q = Q,
               Z = Z, A = matrix(0), R = matrix("r"))
inits_list <- list(x0 = x0)

tmb_ne <- estimate_marss(MARSS(ts(weekly_ne_anom$anom), model = model1, inits = inits_list, fit = F))

This gives me initial values of -1.5 for Q and R diagonals and 0 for everything else. When optimized, the entire Q matrix is negative. If I use estimate_marss2, I get very different answers and a false convergence (8) warning. I'm also not sure why it's estimating the Q off diagonal parameter, as I've set it to zero.

[mdscheuerell]

Hi @abby-caplan,

I haven't tried your code yet, but I can see a problem with your Q matrix. Recall that MARSS() will treat any character class as an unknown parameter to be estimated. Thus, these lines

Q <- matrix("0", 2, 2)
diag(Q) <- c("q.alpha", "q.beta")

will create a matrix with 3 unknown elements to be estimated:

1. q.alpha
2. q.beta
3. 0

Q
     [,1]      [,2]    
[1,] "q.alpha" "0"     
[2,] "0"       "q.beta"

To fix that, you'll need to first create a matrix where each element is a list and then do the substitution. For example,

Q <- matrix(list(0), 2, 2)
diag(Q) <- c("q.alpha", "q.beta")

which yields

Q
     [,1]      [,2]    
[1,] "q.alpha" 0       
[2,] 0         "q.beta"

[ericward-noaa]

I think Mark's change to the structure is right. The other thing we should include in the documentation is that the parameters estimates for Q and R that are being returned are log(standard deviations) -- so can be negative -- example below:

library(MARSS)
library(marssTMB)
library(zoo)

weekly_ne_anom <- read.table("weekly_ne_anom.txt", sep = " ")
weekly_sst_ne <- read.table("weekly_sst_ne.txt", sep = " ")
  
B <- diag(2)
U <- matrix(0, 2, 1)
Q <- matrix(list(0), 2, 2)
diag(Q) <- c("q.alpha", "q.beta")
Z <- array(NA, c(1, 2, 1356))
Z[1, 1, ] <- rep(1, 1356)
Z[1, 2, ] <- na.approx(weekly_sst_ne$anom)
x0 <- matrix(c(0,0), nrow = 2)

model1 <- list(B = B, U = U, Q = Q,
               Z = Z, A = matrix(0), R = matrix("r"))
inits_list <- list(x0 = x0)

# construct MARSS model -- don't fit
marss_fit <- MARSS(ts(weekly_ne_anom$anom), 
                               model = model1, 
                               inits = inits_list,
                   fit=FALSE, method="TMB")

# Fit the model with marssTMB -- note, optimization 
# with TMB is being done in log space, log(standard devs)
# so the parameters need to be exponentiated and squared
fit_tmb <- marssTMB::estimate_marss(marss_fit)
exp(fit_tmb$opt$par[3:4])^2
#> Qdiag.q.alpha  Qdiag.q.beta 
#>  0.4025178739  0.0007256394

^{Created on 2024-09-23 with reprex v2.1.1}

Here's the place in the TMB code where the log sd is converted to sd,

marssTMB/src/TMB/marss.hpp

Line 66 in 3230896

vector<Type> sdQ=QdiagMat.col(0); /* log sd of Q (diag) */