www.phidot.org

I realised that the reason my previous code did not work in the post titled "Time-varying environmental covariates in marked" was because my ch is not made up of only ones and zeros. Therefore, I have switched model types and am using "Mscjs" rather than "probitcjs".

This is how I am adding the covariate (total number of open river days per state over 48 time periods) to the design data:

Code: Select all

# Create vectors of Open River values for each state

OR.a <- OR$OpenRiver[c(1:48)]
OR.b <- OR$OpenRiver[c(49:96)]
OR.c <- OR$OpenRiver[c(97:144)]
OR.d <- OR$OpenRiver[c(145:192)]
OR.e <- OR$OpenRiver[c(193:240)]
OR.f <- OR$OpenRiver[c(241:288)]
OR.g <- OR$OpenRiver[c(289:336)]
OR.h <- OR$OpenRiver[c(337:384)]

# Add or (Open River) to design data in pdfh1.ddl and assign each vector to respective strata

pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "A"] = OR.a
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "B"] = OR.b
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "C"] = OR.c
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "D"] = OR.d
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "E"] = OR.e
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "F"] = OR.f
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "G"] = OR.g
pdfh1.ddl$Psi$or[pdfh1.ddl$Psi$stratum == "H"] = OR.h

This is how I have coded my model:

Code: Select all

#### mod.ss.pt.psior ####
mod.ss.pt.psior <- crm(pdfh1.proc, ddl = pdfh1.ddl, model = "Mscjs",
                      model.parameters = list(p = list(formula = ~ time),
                                              S = list(formula = ~ stratum),
                                              Psi = list(formula = ~ or)),
                      iter = 1000, use.admb = TRUE,
                      hessian = TRUE,
                      se = TRUE)
mod.ss.pt.psior.real <- predict(mod.ss.pt.psior, ddl = pdfh1.ddl, se = TRUE)
mod.ss.pt.psior.real

And this is a snippet of the output:

Code: Select all

mod.ss.pt.psior.real$Psi
    stratum tostratum or occ    estimate           se          lcl         ucl
1         A         A  0  45 1.000000000 0.0000000000 1.0000000000 1.000000000
2         A         B  1  45 0.003108990 0.0007087873 0.0019881155 0.004858723
3         A         C  2  45 0.003108188 0.0006700570 0.0020365631 0.004741014
4         A         D  0  45 0.003109795 0.0007536991 0.0019332898 0.004998676
5         A         E  1  45 0.003108990 0.0007087873 0.0019881155 0.004858723
6         A         F  0  45 0.003109795 0.0007536991 0.0019332898 0.004998676
7         A         G  0  45 0.003109795 0.0007536991 0.0019332898 0.004998676
8         A         H  0  45 0.003109795 0.0007536991 0.0019332898 0.004998676
9         B         A 12  45 0.003100190 0.0007634269 0.0019126692 0.005021297
10        B         B 12  45 1.000000000 0.0000000000 1.0000000000 1.000000000
11        B         C  0  45 0.003109815 0.0007558671 0.0019306655 0.005005514
12        B         D  0  45 0.003109815 0.0007558671 0.0019306655 0.005005514

What doesn't make sense to me are the values in the "or" column in the output. They seem to be pulled randomly from the values in my open river input file.

Does the way that I have gone about this seem logical? If so, can someone help me in getting an open river coefficient?

If this is not the correct way to go about this, can someone offer advice on how to properly add the covariate to my data and use it in the model parameterization.

Thank you!

Not sure but here are a couple of things to consider. First the coefficient is in the betas and not reals. Second, you don't need to use predict to get reals. They are in your model$results. Finally, did you look at ddl$Psi to see if your or values were assigned correctly?

Hi Jeff,

Thank you for the pointers about the coefficients and the reals. Upon investigation, the open river values were not assigned correctly in the ddl$Psi. I am working to correct this issue and will hopefully have more success in running my models.

Thank you,
Katie

Hi Jeff,

Thank you again for addressing my issues with the environmental covariate. I was able to identify the problem and am working to address it by properly formatting my data frame for input in the ddl$Psi.

I have a question pertaining to the MCMC chain outputs. Where in the model results can I find the raw values for the chain estimates?

Best,
Katie

Hi Katie,
Jeff passed this on to me because I wrote the code for probitCJS(). I assume you mean the raw MCMC chain, correct? not the beta estimates for the parameters. If so, you can find them at mymodel$results$beta.mcmc for the beta values. If you fitted a random effects model, you can find the random effects chains in mymodel$alpha.mcmc for the effects themselves and mymodel$results$var.comp.mcmc for the associated variance components. These are mcmc class objects from the package coda So you can look at the docs for that package to see what you can do with them, e.g., trace plots. convergence stats, etc.
-Devin