[R-sig-ME] Standard errors or t values from glmm.admb
Roberto Patuelli
roberto.patuelli at usi.ch
Wed Jun 23 12:44:47 CEST 2010
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
More information about the R-sig-mixed-models
mailing list