mdscheuerell / gretadfa

Using the greta package for R to fit Dynamic Factor Models

---

I'm the Assistant Unit Leader in the U.S. Geological Survey Washington Cooperative Fish and Wildlife Research Unit. I'm also an associate professor in the School of Aquatic and Fishery Sciences at the University of Washington.

I am an applied ecologist who integrates different data sources and analytical methods to study a variety of problems related to the conservation and management of aquatic resources, particularly along the west coast of North America. Much of my research is focused on the development and application of statistical methods for analyzing temporal and spatial data. Examples of recent projects include integrated population models for Pacific salmon, evaluation of the risks and rewards of ecological portfolios, and assessing the effects of large-scale disturbances from natural and anthropogenic causes. You can learn more about me here.

gretadfa's Issues

update constraints for Z

Update the code so that the diagonal of Z (z₁₁, z₂₂, ..., z_MM) is constrained to be positive.

initial values for diagonal of factor loadings matrix (Z)

@goldingn greta only allows one to set an initials() object for variable greta_array's, but the R6 defn we wrote for the special prior distribution of diag(Z) has not yet been incorporated into greta as one of the distribution() classes.

diagonal constraint for Z

The diagonal element of Z must be positive to insure convergence when sampling via Bayesian methods. In particular, Leung & Drton (2014) demonstrate a method that is invariant to any ordering of the observations.

problem with `model()` and `mcmc()`

Hey @goldingn, I finally found some time to return to the DFA model in greta, and I ran into a bunch of errors associated with the code changes in the new ver. I think I've got most of them ironed out now, but I can't seem to figure out the following problem with model()

Specifically, this line causes the following error (with traceback):

Error in parse(text = self$operation) : cannot coerce type 'closure' to vector of type 'character'
18. parse(text = self$operation)
17. eval(parse(text = self$operation), envir = self$tf_function_env)
16. self$tf(dag)
15. x$define_tf(dag)
14. FUN(X[[i]], ...)
13. lapply(self$children[which(!children_defined)], function(x) x$define_tf(dag))
12. x$define_tf(self)
11. FUN(X[[i]], ...)
10. lapply(target_nodes, function(x) x$define_tf(self))
9. force(expr)
8. tryCatchList(expr, classes, parentenv, handlers)
7. tryCatch(force(expr), finally = { data$`__exit__`(NULL, NULL, NULL) })
6. with.python.builtin.object(self$tf_graph$as_default(), expr)
5. with(self$tf_graph$as_default(), expr)
4. self$on_graph(lapply(target_nodes, function(x) x$define_tf(self)))
3. self$define_tf_body()
2. dag$define_tf()
1. model(xx_est, ZZ_est, RR_est, sigma_est)

Simplifying the call to

mod_fit <- model(sigma_est)

seems to work (i.e., it doesn't return an error), but then the following call to mcmc() causes the following error (with traceback):

Error in py_get_attr_impl(x, name, silent) : AttributeError: 'module' object has no attribute 'distributions'
29. stop(structure(list(message = "AttributeError: 'module' object has no attribute 'distributions'", call = py_get_attr_impl(x, name, silent), cppstack = structure(list( file = "", line = -1L, stack = c("1 reticulate.so 0x0000000108b0af9b _ZN4Rcpp9exceptionC2EPKcb + 219", "2 reticulate.so 0x0000000108b11a35 _ZN4Rcpp4stopERKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEE + 53", ...
28. py_get_attr_impl(x, name, silent)
27. py_get_attr(x, name)
26. py_get_attr_or_item(x, name, TRUE)
25. `$.python.builtin.object`(x, name)
24. `$.python.builtin.module`(tfp, distributions)
23. tfp$distributions
22. self$tf_distrib(parameters, dag)
21. self$tf_log_density_function(tf_target, tf_parameters, dag)
20. (function (tf_target) { tf_parameters <- self$tf_fetch_parameters(dag) target_params <- match_batches(c(list(tf_target), tf_parameters)) ...
19. (function (what, args, quote = FALSE, envir = parent.frame()) { if (!is.list(args)) stop("second argument must be a list") ...
18. mapply(do.call, density_functions, target_lists, MoreArgs = list(envir = tfe), SIMPLIFY = FALSE)
17. dag$define_joint_density()
16. force(expr)
15. tryCatchList(expr, classes, parentenv, handlers)
14. tryCatch(force(expr), finally = { data$`__exit__`(NULL, NULL, NULL) })
13. with.python.builtin.object(self$tf_graph$as_default(), expr)
12. with(self$tf_graph$as_default(), expr)
11. dag$on_graph(dag$define_joint_density())
10. self$valid_parameters(inits)
9. FUN(X[[i]], ...)
8. lapply(init_list, self$check_initial_values)
7. self$set_initial_values(initial_values)
6. super$initialize(initial_values = initial_values, model = model, parameters = parameters, seed = seed)
5. .subset2(public_bind_env, "initialize")(...)
4. sampler$class$new(initial_values, model, sampler$parameters, seed = seed)
3. FUN(X[[i]], ...)
2. lapply(initial_values_split, build_sampler, sampler, model)
1. mcmc(mod_fit, sampler = hmc(Lmin = 5, Lmax = 10, epsilon = 0.1, diag_sd = 1), warmup = 2000, n_samples = 5000, thin = 10, chains = 1, verbose = FALSE)

Ideas?

inverse rotation for factor

Need to apply an inverse factor rotation to the factors.

version problem?

Hi Mark,

I was trying to replicate the DFA example using greta: https://mdscheuerell.github.io/gretaDFA/#parameter_estimation

All codes are the same, except my greta and tensorflow versions are different.

tensorflow::tf_version()
[1] ‘1.10’
installed.packages()[c("greta", "tensorflow"), "Version"]
greta tensorflow
"0.2.5" "1.8"

I kept getting this warning message:
In utils::compareVersion("1.8", tf_version) : NAs introduced by coercion

Additionally, when I ran the line: model(xx_est,ZZ_est,RR_est_raw)
I got an error message:
Error in parse(text = self$operation) :
cannot coerce type 'closure' to vector of type 'character'

Do you think this is a version problem, and I should go back to the older versions as shown in the example?

Thanks,
Fan

priors for Z and R

If I worked out the math properly, scaling each of the variates (y) to have zero mean and unit variance should allow for constraining the elements of Z and R to be less than 1.