Actually AD Model Builder has been freely available and
open source for some time.
http://admb-project.org
But the easiest way to see everything
is to run glmmADMB with the option of keeping the temporary files
and then using the ADMB documentation to interpret these files.
There are a lot of interesting options that are not available
using the R interface. However first running the model in R is
useful to create the model matrix etc in a way familiar to r
users.
Standard errors or t values from glmm.admb
2 messages · dave fournier, Roberto Patuelli
Dear Paul, Thanks a lot for the "dedication" in investigating the function. Indeed, I should have looked more carefully at the attributes of the regression object. This is also confirmed to me off-list by Dave Fournier, who works for ADMB. So at least standard errors are not a problem, as well as probably t-values. P-values are a different story of course... I'll run the model again (it takes 12 hours) and see if everything hidden looks fine. Cheers Roberto ----- Original Message ----- From: "Paul Johnson" <pauljohn32 at gmail.com> To: "Patuelli Roberto" <roberto.patuelli at usi.ch>; "R-SIG-Mixed-Models at r-project.org" <r-sig-mixed-models at r-project.org> Sent: Tuesday, June 22, 2010 9:24 PM Subject: Re: [R-sig-ME] Standard errors or t values from glmm.admb On Tue, Jun 22, 2010 at 3:42 AM, Roberto Patuelli
<roberto.patuelli at usi.ch> wrote:
Dear All, I have experimented with estimating a negative binomial mixed model through the glmmADMB package (that is, using the glmm.admb function). After an entire night of computing :) I have results (though with a warning message saying that proper convergence could not be reached). But when I use summary.glmm.admb, I can only see, with regard to the fixed effects, the betas estimates, but no standard errors or other measures regarding the precision of the regression coefficients. See below:
I've not tried glmmADMB before, but I went and found it at http://otter-rsch.com/admbre/examples/glmmadmb/glmmADMB.html. I agree, this is lacking some creature comforts, summary is not exhaustive, vcov() is not implemented. Appears no way to get the Hessian out of it. So, well, I agree this would be a frustrating user experience. Here's where being old helps. :) I went old school on it! Suppose your output object is ad1, run attributes(ad1), I see:
attributes(ad1)
$names [1] "n" "q" "fixed" "call" [5] "group" "family" "corStruct" "impSamp" [9] "easyFlag" "zeroInflation" "b" "stdbeta" [13] "alpha" "sd_alpha" "S" "sd_S" [17] "random" "U" "sd_U" "fitted" [21] "residuals" "npar" "loglik" "gradloglik" $class [1] "glmm.admb" I wondered if ad1$stdbeta might give me what you want. From the name, it appears promising. You'd need to run some model in some program that we trust, like "lmer", and see if the values match up. Without documentation from the ADMB folks, its hard to guess. ## Create your own t values for the fixed effects.
ad1$b/ad1$stdbeta
(Intercept) Base trtprogabide Age
-1.129856 6.759541 -2.328344 1.374022
Visit Base:trtprogabide
-1.619971 1.665581
Now, I wonder if that is the correct number. You can read the code if
you do this (this used to work for almost all R commands, still works
for many). Read down a ways to find "stdbeta"
glmm.admb
function (fixed, random, group, data, family = "poisson", link,
corStruct = "diag", impSamp = 0, easyFlag = TRUE, zeroInflation = FALSE,
imaxfn = 10, save.dir = NULL, offset)
{
olddir <- getwd()
dirname <- ifelse(is.null(save.dir), "_glmm_ADMB_temp_dir_",
save.dir)
dir.create(dirname)
setwd(dirname)
if (is.null(save.dir))
on.exit(unlink(dirname, recursive = TRUE))
call <- match.call()
if (!(corStruct == "diag" | corStruct == "full"))
stop("Argument \"corStruct\" must be either \"diag\" or \"full\"")
if (!(family == "poisson" | family == "nbinom" | family ==
"binomial"))
stop("Argument \"family\" must be either \"poisson\",
\"binomial\" or \"nbinom\"")
if (missing(random) & ((!missing(impSamp)) | (!missing(corStruct))))
stop("When \"random\" is missing neither of \"impSamp\" or
\"corStruct\" make sense")
if (missing(group) || (!(is.character(group) & (length(group) ==
1))) || !member(group, names(data)))
stop("Argument \"group\" must be a character string spesifying
the name of the grouping variable (also when \"random\" is missing)")
if (missing(link) & family == "binomial")
stop("Argument \"link\" must be provided for the \"binomial\"
family)")
if ((!missing(offset)) && ((!is.character(offset)) || (length(offset) !=
1) || !member(offset, names(data))))
stop("An error occured in relatation to the argument
\"offset\". It must be a character string spesifying of the variable
holding the offset")
if (!missing(offset))
Offset = data[[offset]]
else Offset = rep(0, nrow(data))
tmpu <- table(data[[group]])
tmpu[] <- 1:length(tmpu)
tmpI <- tmpu[as.character(data[[group]])]
data <- data[order(tmpI), ]
n <- nrow(data)
y <- data[[as.character(fixed[2])]]
X <- model.matrix(fixed, data)
p <- ncol(X)
II <- as.numeric(tmpu[as.character(data[[group]])])
group_d <- c(1, (2:n)[diff(II) == 1], n + 1)
if (missing(random))
Z <- model.matrix(~1, data)
else Z <- model.matrix(random, data)
m <- ncol(Z)
q <- length(unique(II))
cmdoptions = paste("-maxfn 500", ifelse(impSamp == 0, "",
paste("-is", impSamp)))
dat_list = list(n = n, y = y, p = p, X = X, q = q, m = m,
Z = Z, group_d = group_d, II = II, cor_flag = ifelse(corStruct ==
"full", 1, 0), like_type_flag = as.numeric(family ==
"poisson" || ((family == "binomial") && (link ==
"logit"))), no_rand_flag = as.numeric(missing(random)),
easy_flag = as.numeric(easyFlag), zi_flag =
as.numeric(zeroInflation),
intermediate_maxfn = 10, offset = Offset)
pin_list = list(pz = ifelse(zeroInflation, 0.02, 1e-04),
b = numeric(p), tmpL = 0.25 + numeric(m), tmpL1 = 1e-04 +
numeric(m * (m - 1)/2), logalpha = 2, kkludge = 0,
u = matrix(numeric(q * m), q, m))
if (family == "binomial") {
if (!all(member(y, 0:1)))
stop("Only Bernoulli responce (0 or 1) allowed currently")
file_name = "bvprobit"
dat_list = dat_list[c("n", "y", "p", "X", "q", "m", "Z",
"group_d", "cor_flag", "like_type_flag", "no_rand_flag")]
pin_list = pin_list[c("b", "tmpL", "tmpL1", "kkludge",
"u")]
}
else {
file_name = "nbmm"
dat_list = dat_list[c("n", "y", "p", "X", "q", "m", "Z",
"II", "cor_flag", "easy_flag", "zi_flag", "like_type_flag",
"no_rand_flag", "intermediate_maxfn", "offset")]
pin_list = pin_list[c("pz", "b", "tmpL", "tmpL1", "logalpha",
"kkludge", "u")]
if (!missing(link))
warning("Argument \"link\" ignored for other familes than
\"binomial\" (exponential link used)")
}
dat_write(file_name, dat_list)
pin_write(file_name, pin_list)
std_file = paste(file_name, ".std", sep = "")
file.remove(std_file)
if (.Platform$OS.type == "windows") {
shell(paste(.path.package("glmmADMB"), "/admb/", file_name,
".exe", " ", cmdoptions, sep = ""), invisible = TRUE)
}
else {
system(paste("cp ", .path.package("glmmADMB"), "/admb/",
file_name, " .", sep = ""))
system(paste("./", file_name, " ", cmdoptions, sep = ""))
unlink(file_name)
}
if (!file.exists(std_file))
stop("The function maximizer failed")
tmp <- read.table(paste(file_name, ".std", sep = ""), skip = 1)
tmpindex <- as.character(tmp[, 2])
out <- list(n = n, q = q, fixed = fixed, call = call, group = group,
family = family, corStruct = corStruct, impSamp = impSamp,
easyFlag = easyFlag, zeroInflation = zeroInflation)
if (zeroInflation)
out$pz = as.numeric(tmp[tmpindex == "pz", 3])
out$b = as.numeric(tmp[tmpindex == "real_b", 3])
names(out$b) = colnames(X)
out$stdbeta = as.numeric(tmp[tmpindex == "real_b", 4])
names(out$stdbeta) = colnames(X)
if (family == "nbinom") {
out$alpha <- as.numeric(tmp[tmpindex == "alpha", 3])
out$sd_alpha <- as.numeric(tmp[tmpindex == "alpha", 3])
}
if (!missing(random)) {
out$S = matrix(as.numeric(tmp[tmpindex == "S", 3]), nrow = m,
dimnames = list(colnames(Z), colnames(Z)))
out$sd_S = matrix(as.numeric(tmp[tmpindex == "S", 4]),
nrow = m, dimnames = list(colnames(Z), colnames(Z)))
if (corStruct == "diag") {
out$S[below(m, TRUE) | t(below(m, TRUE))] = 0
out$sd_S[below(m, TRUE) | t(below(m, TRUE))] = 0
}
out$random = random
}
if (!missing(link))
out$link = link
if (!missing(random)) {
U = matrix(as.numeric(tmp[tmpindex == "u", 3]), ncol = m,
byrow = TRUE, dimnames = list(data[, group][group_d[-1] -
1], colnames(Z)))
out$U = U
out$sd_U = matrix(as.numeric(tmp[tmpindex == "u", 4]),
ncol = m, byrow = TRUE, dimnames = list(data[,
group][group_d[-1] -
1], colnames(Z)))
}
else U = matrix(rep(0, q), ncol = 1, byrow = TRUE)
mu <- as.numeric(X %*% out$b)
lambda <- 0
for (i in 1:n) lambda[i] <- exp(mu[i] + sum(Z[i, ] * U[II[i],
]))
if (family == "binomial")
out$fitted <- lambda/(1 + lambda)
else out$fitted <- lambda
tmpsd <- switch(family, poisson = sqrt(lambda), nbinom = sqrt(lambda *
(1 + lambda/out$alpha)), binomial = sqrt(out$fitted *
(1 - out$fitted)))
out$residuals <- as.numeric(y - lambda)/tmpsd
tmp = par_read(file_name)
out$npar <- as.numeric(scan(paste(file_name, ".par", sep = ""),
what = "", quiet = TRUE)[6])
out$loglik = tmp$loglik
out$gradloglik = tmp$gradient
if (abs(out$gradloglik) >= 0.001)
warning("Proper convergence could not be reached")
setwd(olddir)
class(out) <- "glmm.admb"
return(out)
}
It doesn't tell us exactly what stdbeta represents, but it is still
interesting.
I'm able to see that this routine does not actually do calculations,
rather it calls a separate program that ADMB installed in your R lib
folder, "admb.exe" or "admb" depending on your OS. This writes
commands to a file, runs that external program, writes the output into
a temp file. I suppose I would not go further to deconstruct the ADMB
software to see what they mean, since they have claimed it is secret
in the past.
Well, that's all I can do. But it does remind me of the dangers of
using packages that don't exist in CRAN and aren't exactly "open
source". You probably need to directly contact the owners....
Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas