atsa-es / safs-timeseries Goto Github PK

Hi all,

I'm having no trouble getting the simplified model in problem 1 to work, but when I try to add a covariate, my process error for the "slope" parameter converges to zero, so the parameter is not dynamic through time. Here's the code/output. I feel like I'm missing something obvious:

Recruits = matrix(log(SRdata$Rec), nrow=1)
Prop_Rec = matrix(log(SRdata$Rec/SRdata$Sp), nrow=1)
Spawners = matrix(log(SRdata$Sp), nrow=1)
TT = length(SRdata$brood.yr)
m = dim(Spawners)[1] + 1

#Process error#
B = diag(m)
U = matrix(0, nrow=m, ncol=1)		
Q = matrix(list(0),m,m)			
diag(Q) = c("q.alpha","q.beta")

#Observation error#
R = matrix("r")				
Z = array(NA, c(1,m,TT)) 		
Z[1,1,] = rep(1,TT) 
Z[1,2,] = Spawners
A = matrix(0)			
inits.list = list(x0=matrix(c(0,0), nrow=m))
mod.list = list(B=B, U=U, Q=Q, Z=Z, A=A, R=R)
DLM.salmon3 = MARSS(Prop_Rec, inits=inits.list, model=mod.list, control=list(maxit=2000))

Log-likelihood: -48.51628
AIC: 107.0326 AICc: 108.9076

      Estimate

R.r 0.269
Q.q.alpha 0.316
Q.q.beta 0.000
x0.X1 1.344
x0.X2 -0.051

Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.

The AICc value for this model is also greater than for the simplified model, so I get the feeling something's off. Same thing happens when I include PDO as a covariate. Am I setting up my matrices completely wrong?

You'd enter your question here...

When fitting a DLM, if our model output estimates the value of one value of the Q matrix as zero despite you allowing it to be non-zero, should this be interpreted as an issue with our model fitting or that our model is fitting correctly but determining that the best fit is without the time-varying effect of a covariate? I am happy to post a reprex if it's appropriate to do so for a homework assignment.

Here's what my Q matrix looks like and what the model output looks like - in this case it's the "Q.q.delta" value that's estimated as zero:

> Q
     [,1]      [,2]      [,3]
[1,] "q.alpha" 0         0   
[2,] 0         "q.delta" 0   
[3,] 0         0         0   

> fit_3

MARSS fit is
Estimation method: kem 
Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
Estimation converged in 679 iterations. 
Log-likelihood: -51.78324 
AIC: 115.5665   AICc: 117.8998  

Estimate
R.r       2.00e-01
Q.q.alpha 3.16e-01
Q.q.delta 0.00e+00
x0.X1     1.01e+00
x0.X2     2.85e-02
x0.X3     5.49e-06
Initial states (x0) defined at t=0

Standard errors have not been calculated. 
Use MARSSparamCIs to compute CIs and bias estimates.

Thanks!

I am struggling with creating the Z matrix for problem 5 and solving for the parameters. I presently have the following for the y and Z matrices-

y=matrix(homeworkdat$Ozone, nrow=1)

fit=lm(Ozone ~ -1 + region, data=homeworkdat)
Z=rbind(c(1,0,1,0,1), c(0,1,0,1,0))

I can't get the solve to run (solve(t(Z)%%Z)%%t(Z)%*%y) and I figure it's because I'm not setting the Z matrix up properly. Can anyone help?

Couple of questions on HW1.

First, should the adf.test() and ur.df() functions lead to the same conclusions as long as we specify the same assumptions/lags within each function? For example, I would expect that the following two lines of code should lead to the same conclusion (WRT rejecting or not rejecting the null). However, I am finding that that is not the case.
adf.test(dat, k = 0)
ur.df(dat, type = "trend", lags = 0)

For question 4b, when using the auto.arima() function with trace=TRUE, I am finding that none of the candidate models are within Delta AICc = 2 of the top model. Is this user error on my part, or am I missing something?

for question 1.3 in HW #2. It was confusing that you said m_t = 0. When we do decompose(pDat) the trend panel is not zero. In't that 2nd panel m_t so definitely not zero?

I'm having issues with Problems 6 and 7 on the homework, namely, neither the arima() call, nor the cross-validation technique are working with multivariate datasets. I can get both sets of code to work if I limit my analysis to one taxonomic group...but with the full matrix it's errors galore. I've done a pretty thorough search for answers on stack overflow, but nothing is coming up. Has anyone found a solution? Am I missing something easy here?

What have I gotten into?

x <- 1

I am stuck on how to simulate the sample trajectories, not even sure where to start. Please advise!

I'm confused about the lecture slides for lecture 10 (Hidden Markov Models) - not sure if this is a typo or if I'm missing something.

On slide 19, it says this:

There are two general flavours of transition matrices:
􏰀- Homogenous (or stationary), where transition probabilities don’t depend on t
􏰀- Non-homogeneous, where transition probabilities are time-varying
􏰀 In this class, we’ll only consider homogeneous cases

Then, on slide 21, there is an example where there is "a temporal trend in the transition probabilities" and there is a question in the homework that asks you to include "time varying parameters (e.g. year) in the means of each state, or state transition probabilities"

Is slide 19 a typo then, since the examples on slide 21 and in the homework are non-homogenous? Or am I misunderstanding the definition of homogenous vs. non-homogenous transition matrices.

For homework 4 on dynamic linear models, I'm getting stuck on question 2. I want the observation and process equations to be the following:

$$ log(R_t/S_t) = \alpha_t - bS_t + v_t \\ \alpha_t = \alpha_{t-1} + w_t $$

My code is below. I'm getting an error that Q is not properly constrained, but I only want alpha to vary with time, not b. Please let me know if I set up the initial equations incorrectly or if there's a way to fix the MARSS error. Thanks!

Code:

library(MARSS)
data(KvichakSockeye, package = "atsalibrary")
SR_data <- KvichakSockeye
# response variable
data <- ts(data = log(SR_data$Rec/SR_data$Sp), start = 1952, end = 1989, frequency = 1)
# z-score the predictor variable
Sp_z <- matrix((SR_data$Sp-mean(SR_data$Sp))/sqrt(var(SR_data$Sp)), nrow = 1)
m <- 2
TT <- length(data)

B <-  diag(m)
U <-  matrix(0, nrow = m, ncol = 1)
Q <- matrix(c("q.alpha", 0, 0, 0), nrow = 2)
Z <- array(NA, c(1, m, TT))  ## NxMxT; empty for now
Z[1,1,] <- rep(1, TT)        ## Nx1; 1's for intercept
Z[1,2,] <- Sp_z              ## Nx1; predictor variable 
A = matrix(0)
R <-  matrix("r")

inits_list <- list(x0 = matrix(c(0, 0), nrow = m))
mod_list <- list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)

## fit the model with MARSS
dlm_fit <- MARSS(data, mod_list, inits = inits_list)

Error message:

Errors were caught in is.marssMODEL()
Errors were caught in is.marssMODEL_marxss()
The variance-covariance matrix Q is not properly constrained.
Covariances and variances cannot be shared at t=1

Let's see who gets this.

Hi, I was absent from lab this week (conflicting meeting). The syllabus lists our HW as TBD, and I'm wondering what the assigned problems are? Thanks!

Hi there,

I wanted to jump in and mention that in this case I do NOT recommend transforming (eg, z-score) your covariates in the DLMs you fit for this week's homework. I can talk about the reasons in class on Tues.

Mark

Hi,

I was working on homework 4 and I stumbled upon some issues on the scope of some questions.

• When we are asked to plot the residuals and acf for each model, does it mean both state and model residuals? Should we show them for all the models (or only the best one for each question)? Should we comment what we see in the residuals and acfs? I am thinking yes to every question but I just wanted to check before I start coding and writing too much.
• I don't know how to tackle questions 0.5 and 0.6. How should we compare the results? Are we looking more at the coefficients? To what should we compare the lm and arima models? The best model out of questions 0.2 to 0.4? I am not sure how to be able to compare the lm models (or the arima ones) to anything I did in questions 0.3 and 0.4. 0.3 looked at the difference between observation and process when here we have either one or the other (lm or arima). 0.4 asked for different or the same coefficients for the functional groups. How can we achieve this in terms of regression or arima? Also, should we include a lagged variable in the regression?

I am probably over thinking everything but a little guidance would be much appreciated.

Thanks!

Hey Eric,

R is having issues installing the STATS package to my machine from github. It gives the following error message:

• installing source package 'STATS' ...
  ** R
  ** exec
  ** preparing package for lazy loading
  ** help
  *** installing help indices
  ** building package indices
  ** testing if installed package can be loaded
  *** arch - i386
  During startup - Warning message:
  package 'stats' in options("defaultPackages") was not found
  Error : package 'STATS' is not installed for 'arch = i386'
  Error: loading failed
  Execution halted
  *** arch - x64
  ERROR: loading failed for 'i386'
• removing 'C:/Users/melaniedavis/Documents/R/R-3.3.0/library/STATS'
• restoring previous 'C:/Users/melaniedavis/Documents/R/R-3.3.0/library/STATS'
  Warning in file.copy(lp, dirname(pkgdir), recursive = TRUE, copy.date = TRUE) :
  problem copying C:\Users\melaniedavis\Documents\R\R-3.3.0\library\00LOCK-STATS\STATS\libs\x64\stats.dll to C:\Users\melaniedavis\Documents\R\R-3.3.0\library\STATS\libs\x64\stats.dll: Permission denied
  Error: Command failed (1)

Any idea what's going on here?

Thanks,
Melanie

I've tried various ways of installing the atsar package (restarting R studio, updating other packages first, updating R, selecting different options when prompted during the install), but I keep getting this error message:

devtools::install_github("nwfsc-timeseries/atsar")

ERROR: failed to lock directory ‘/Library/Frameworks/R.framework/Versions/4.0/Resources/library’ for modifying Try removing ‘/Library/Frameworks/R.framework/Versions/4.0/Resources/library/00LOCK-atsar’ Warning message: In i.p(...) : installation of package ‘/var/folders/0v/gmf2p3cj6jl70lg78pwyqhzc0000gn/T//RtmpNvBT1o/file1ae920c806f2/atsar_0.1.5.tar.gz’ had non-zero exit status

Is there a simple way to remedy this?
Thank you!

sessionInfo()
R version 4.0.4 (2021-02-15)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
other attached packages:
[1] rstantools_2.1.1 rstan_2.21.1