[R] glmmPQL and variance structure
stephenb
stenka1 at go.com
Wed Aug 26 17:07:16 CEST 2009
this is very late, but I saw this now as I am dealing with it now:
I think varPower should not be needed here. The family should be one of the
quasi families eg quasibinomial and that will automatically allow
variance/dispersion to become a function of the fit. This is a feature of
glm inherently, which does not exist in lme, so it is called with varPower
(in lme).
If you see that, pls post whether it works as expected on the original
dataset.
regards
Stephen
Spencer Graves wrote:
>
> 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
>
> ______________________________________________
> 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
>
>
--
View this message in context: http://www.nabble.com/glmmPQL-and-variance-structure-tp2107860p25151412.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list