Playing with `mcglm` for a multivariate poisson model. Not sure how to extract the inter-DV covariance.

mcglm.r

# From the docs for mcglm::ahs

require(mcglm)
data(ahs, package="mcglm")
form1 <- Ndoc ~ income + age
form2 <- Nndoc ~ income + age

Z0 <- mc_id(ahs)
fit.ahs <- mcglm(linear_pred = c(form1, form2),
                 matrix_pred = list(Z0, Z0), link = c("log","log"),
                 variance = c("poisson_tweedie","poisson_tweedie"),
                 data = ahs)
summary(fit.ahs)

fit.ahs$Covariance

#####################################################################
# Example above extended to the full model, fitting 10 IVs to 5 DVs #
#####################################################################

require(mcglm)
data(ahs, package="mcglm")

iv = names(ahs)[1:10]
dv = names(ahs)[11:15]

rhs = paste(iv, collapse=" + ")
formulas_str = paste(dv, rhs, sep = " ~ ")

formulas = sapply(formulas_str, as.formula)

Z0 <- mc_id(ahs)

# `control_algorithm` settings lifted from the Poisson-Tweedie demo in the GLMExamples.R
# Not sure about the rest of it, but 4.1 of the article specifies that they use bias correction,
# which presumably refers to "correct=TRUE", and that the fitting algorithm should be the 
# "modified chaser" algorithm.
fit.ahs.full <- mcglm(linear_pred = formulas,
                 matrix_pred = list(Z0, Z0, Z0, Z0, Z0), 
                 link = rep("log", 5),
                 variance = rep("poisson_tweedie", 5),
                 data = ahs,
                 control_algorithm = list("correct" = TRUE, 
                                          #"verbose" = TRUE,
                                          "tol" = 1e-5,
                                          "max_iter" = 100, 
                                          "method" = "chaser", 
                                          "tunning" = 1)
                 )

# Presumably this is how we're supposed to use this...?
length(fit.ahs.full$Covariance) # 15
n = length(dv)
M = matrix(0, n, n)
sum(lower.tri(M, diag=TRUE)) # 15
M[lower.tri(M, diag=TRUE)] = fit.ahs.full$Covariance
M

# Diagonals aren't even close to 1. 
# off diagonals don't align with the paper, especially (6,4) which is 19.1
# This doesn't align with the \hat{\Sigma_b} given at the top of page 13 of the paper.

####

# This is another way we could go, but then we have a whole dimension on top of the 
# table in the paper, on top of the values being off...
M = matrix(0, n+1, n+1)
sum(lower.tri(M, diag=FALSE))
M[lower.tri(M, diag=FALSE)] = fit.ahs.full$Covariance
diag(M) = 1
M