[R] glmmPQL and variance structure
Spencer Graves
spencer.graves at pdf.com
Sun Jan 15 02:59:06 CET 2006
I just identified an error in my recent post on this subject: There
is a very good reason that Venables & Ripley's "glmmPQL" did NOT include
an argument like the "weights.lme" in the "glmmPQL." included in my
recent post: Their function calls "glm" first and then provides weights
computed by "glm" to "lme". If you want other weights, you need to
consider appropriately the weights computed by "glm". It may be
reasonable to make your other weights proportional to the "glm" weights,
but it would not be smart to do as I suggested a few hours ago, which
ignored the "glm" weights.
I hope this error doesn't seriously inconvenience you or anyone else
who might read these comments.
spencer graves
####################################
Thanks for providing a partially reproducible example. I believe the
error message you cite came from "lme". I say this, because I modified
your call to "glmmPQL2" to call lme and got the following:
> library(nlme)
> fit.lme <- lme(y ~ trt + I(week > 2), random = ~ 1 | ID,
+ data = bacteria, weights=varPower(~1))
Error in unlist(x, recursive, use.names) :
argument not a list
I consulted Pinheiro and Bates (2000) Mixed-Effects Models in S and
S-Plus (Springer, sec. 5.2, p.211) to see that the syntax for "varPower"
appears to be correct. I removed "~" and it worked, mostly:
> fit.lme <- lme(y ~ trt + I(week > 2), random = ~ 1 | ID,
+ data = bacteria, weights=varPower(1))
Warning message:
- not meaningful for factors in: Ops.factor(y[revOrder], Fitted)
I got "an answer", though with a warning and not for the problem you
want to solve. However, I then made this modification to a call to my
own modification of Venables and Ripley's glmmPQL and it worked:
> fit. <- glmmPQL.(y ~ trt + I(week > 2), random = ~ 1 | ID,
+ family = binomial, data = bacteria,
+ weights.lme=varPower(1))
iteration 1
iteration 2
iteration 3
> fit.
Linear mixed-effects model fit by maximum likelihood
Data: bacteria
Log-likelihood: -541.0882
Fixed: y ~ trt + I(week > 2)
(Intercept) trtdrug trtdrug+ I(week > 2)TRUE
2.7742329 -1.0852566 -0.5896635 -1.2682626
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 4.940885e-05 2.519018
Variance function:
Structure: Power of variance covariate
Formula: ~fitted(.)
Parameter estimates:
power
0.3926788
Number of Observations: 220
Number of Groups: 50
>
This function "glmmPQL." adds an argument "weights.lme" to Venables
and Ripley's "glmmPQL" and uses that in place of
'quote(varFixed(~invwt))' when provided; see below.
hope this helps.
spencer graves
glmmPQL. <-
function (fixed, random, family, data, correlation, weights,
weights.lme, control, niter = 10, verbose = TRUE, ...)
{
if (!require("nlme"))
stop("package 'nlme' is essential")
if (is.character(family))
family <- get(family)
if (is.function(family))
family <- family()
if (is.null(family$family)) {
print(family)
stop("'family' not recognized")
}
m <- mcall <- Call <- match.call()
nm <- names(m)[-1]
wts.lme <- mcall$weights.lme
keep <- is.element(nm, c("weights", "data", "subset", "na.action"))
for (i in nm[!keep]) m[[i]] <- NULL
allvars <- if (is.list(random))
allvars <- c(all.vars(fixed), names(random), unlist(lapply(random,
function(x) all.vars(formula(x)))))
else c(all.vars(fixed), all.vars(random))
Terms <- if (missing(data))
terms(fixed)
else terms(fixed, data = data)
off <- attr(Terms, "offset")
if (length(off <- attr(Terms, "offset")))
allvars <- c(allvars, as.character(attr(Terms, "variables"))[off +
1])
m$formula <- as.formula(paste("~", paste(allvars, collapse = "+")))
environment(m$formula) <- environment(fixed)
m$drop.unused.levels <- TRUE
m[[1]] <- as.name("model.frame")
mf <- eval.parent(m)
off <- model.offset(mf)
if (is.null(off))
off <- 0
w <- model.weights(mf)
if (is.null(w))
w <- rep(1, nrow(mf))
mf$wts <- w
fit0 <- glm(formula = fixed, family = family, data = mf,
weights = wts, ...)
w <- fit0$prior.weights
eta <- fit0$linear.predictor
zz <- eta + fit0$residuals - off
wz <- fit0$weights
fam <- family
nm <- names(mcall)[-1]
keep <- is.element(nm, c("fixed", "random", "data", "subset",
"na.action", "control"))
for (i in nm[!keep]) mcall[[i]] <- NULL
fixed[[2]] <- quote(zz)
mcall[["fixed"]] <- fixed
mcall[[1]] <- as.name("lme")
mcall$random <- random
mcall$method <- "ML"
if (!missing(correlation))
mcall$correlation <- correlation
# weights.lme
{
if(is.null(wts.lme))
mcall$weights <- quote(varFixed(~invwt))
else
mcall$weights <- wts.lme
}
mf$zz <- zz
mf$invwt <- 1/wz
mcall$data <- mf
for (i in 1:niter) {
if (verbose)
cat("iteration", i, "\n")
fit <- eval(mcall)
etaold <- eta
eta <- fitted(fit) + off
if (sum((eta - etaold)^2) < 1e-06 * sum(eta^2))
break
mu <- fam$linkinv(eta)
mu.eta.val <- fam$mu.eta(eta)
mf$zz <- eta + (fit0$y - mu)/mu.eta.val - off
wz <- w * mu.eta.val^2/fam$variance(mu)
mf$invwt <- 1/wz
mcall$data <- mf
}
attributes(fit$logLik) <- NULL
fit$call <- Call
fit$family <- family
oldClass(fit) <- c("glmmPQL", oldClass(fit))
fit
}
Patrick Giraudoux wrote:
> Dear listers,
>
> On the line of a last (unanswered) question about glmmPQL() of the
> library MASS, I am still wondering if it is possible to pass a variance
> structure object to the call to lme() within the functions (e.g.
> weights=varPower(1), etc...). The current weights argument of glmmPQL is
> actually used for a call to glm -and not for lme). I have tried to go
> through the code, and gathered that the variance structure passed to the
> call to lme() was:
>
> mcall$weights <- quote(varFixed(~invwt))
>
> and this cannot be modified by and argument of glmmPQL().
>
> I have tried to modify the script a bit wildly and changed varFixed
> into VarPower(~1), in a glmmPQL2 function. I get the following error:
>
> > glmmPQL2(y ~ trt + I(week > 2), random = ~ 1 | ID,
> + family = binomial, data = bacteria)
> iteration 1
> Error in unlist(x, recursive, use.names) :
> argument not a list
>
> I get the same error whatever the change in variance structure on this line.
>
> Beyond this I wonder why variance structure cannot be passed to lme via
> glmmPQL...
>
> Any idea?
>
> Patrick Giraudoux
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list