# Toy example to demonstrate that we cannot adjust mGamma to enforce offsetting
# using the setPriors() function

library(Hmsc)

# Predictor matrix: an intercept plus an offset column
X <- matrix(c(rep(1, 10), log(2^(-2:7))), ncol=2)
dimnames(X)[[2]] <- c("Intercept", "Offset")
# Species counts - Twenty species from an identical distribution
Y <- matrix(rpois(10*20, exp(0)*exp(X[,2])), ncol=20)

# Desired outcome:
# - estimate for gamma[1] = 0 with s.e. > 0
# - estimate for gamma[2] = 1 with s.e. = 0 (Offset term)

# Define two sampling groups
studyDesign <- data.frame(Group = factor(rep(LETTERS[1:5], each=2)))
rL = HmscRandomLevel(units = studyDesign$Group)


### Set up the model
m.setup = Hmsc(Y = Y,
               X = X,
               XScale = FALSE,
               studyDesign = studyDesign,
               ranLevels = list("Group" = rL),
               ranLevelsUsed = "Group",
               distr = "lognormal poisson"
)

# Declare an IW(V0, f0) prior for the V matrix
#   and a N(mGamma, UGamma) prior for the gamma parameters
# The following should declare for the offset parameter that
#   Beta[2] ~ N(Gamma[2], 1e-8) and
#   Gamma[2] ~ N(1, 1e-8)
K <- ncol(X)
m.setup <- setPriors(m.setup,
                     V0 = diag(c(rep(1, K-1), 1e-8)),
                     f0 = K+1,
                     mGamma = c(rep(0, K-1), 1),
                     UGamma = diag(c(rep(1, K-1), 1e-8))
)

# Perform MCMC sampling
samples <- 400
thin <- 1
hmsc.m = sampleMcmc(m.setup,
                    samples = samples,
                    thin = thin,
                    adaptNf = 0,
                    transient = 0,   # Look at burn-in samples
                    nChains = 4, nParallel = 2)

fit <- convertToCodaObject(hmsc.m)


# Visualize results
par(mfrow=c(2,1))
# Gamma[,2] traceplot initiates at one and rapidly centers to zero with tiny
# variance indicating that the prior for Gamma[,2] is N(0,tiny) and that
# sampling uses an initial value of 1, whereas what we desire is a prior of
# N(1,tiny) with a spread of initial values.
traceplot(fit$Gamma)
traceplot(fit$Beta[,c(1,2)])
# Beta[,2] traceplot is centered at zero with tiny variance showing that the
# prior on V[2,2] is correctly applied

### Other issues:

# (1) Labeling of Etas: The correct labels should be: "Eta1[A, factor 1]"
# "Eta1[B, factor 1]" "Eta1[C, factor 1]" "Eta1[D, factor 1]", etc., but Hmsc
# is labeling random effect levels (A-E) with rL$pi instead of levels(rL$pi)
colnames(as.matrix(fit$Eta[[1]]))


# (2) Lambda matrix
# Two separate ways to extract the n-th sample of the Lambda matrix for the
# first random effect Version 1:
n <- 7
V1 <- as.matrix(fit$Lambda[[1]])[n,]

# Version 2:
postList <- poolMcmcChains(hmsc.m$postList)
V2 <- postList[[n]]$Lambda[[1]]

# V1 is a named vector.  V2 is a matrix.  If we trust the names from V1, we
# would convert V1 into a matrix via: matrix(V1, nrow=2, byrow=TRUE), but if we
# look at the values from V1, we would convert it into a matrix via: matrix(V1,
# nrow=2, byrow=FALSE).

# I am pretty sure V2 is correct and that the issue is in the construction of V1.


# (3) Finally, it would be nice to have an informative error message in the case
# where the user accidentally declares adaptNf > transient:
hmsc.m = sampleMcmc(m.setup,
                    samples = samples,
                    thin = thin,
                    adaptNf = 25,
                    transient = 0,
                    nChains = 4, nParallel = 2)