[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