[R] glmmPQL and variance structure

[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