Predict & liability scale issue

[szymekdr]

@itchyshin summary of different ways to predict in MCMCglmm

Model: E(y) = exp(Xb+e)
Linear predictor: eta = Xb
Liability: l = eta+e = Xb+e (Liab = linear predictor + residual variance)
Prediction with rand marginalised: E_e(y) = exp(Xb+0.5*var_e)

Three ways of getting l-scale and predicted values, first (native prediction) and last match:

x <- rnorm(200, 5, 1.5)
ll <- 3 * x + rnorm(200, 0, 3)
toy <- data.frame(y = rpois(200,
    lambda = ifelse(ll<0, 0, ll)
), x = x)

model <- MCMCglmm(y ~ x,
    family = "poisson",
    pl = T, pr = T, data = toy,
    nitt = 500000, thin = 200, burnin = 50000,
    verbose = F,
    saveX = T, saveXL = T, saveZ = T
)
summary(model)

y1 <- predict(model, type = "response")
l1 <- predict(model, type = "terms")

l2 <- data.frame(l2 = apply(model$Liab, 2, mean))
y2 <- data.frame(y2 = exp(l2 + 0.5 * mean(model$VCV[, "units"])))

# posterior of fixed predictions
A <- apply(as.matrix(model$Sol), 1,
    function(x) as.matrix(model$X) %*% x,
    simplify = T
)

# mean of posterior of fixed predictions
l3 <- data.frame(l3 = rowMeans(A))

# posterior of predictions with marginalised units
B <- t(A) + 0.5 * kronecker(t(rep(1, 200)), model$VCV[, "units"])
# doable without kronecker - but this way was easiest to assemble for me
y3 <- data.frame(y3 = colMeans(exp(B)))

head(cbind(l1, y1, l2, y2, l3, y3))