Re: [R] glmmPQL and variance structure

[Spencer Graves]

  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) trtdrugtrtdrug+ 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(~invw

Re: [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) trtdrugtrtdrug+ 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 Girau

Re: [R] glmmPQL and variance structure

[Patrick Giraudoux]

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@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

Re: [R] glmmPQL and variance structure

[Spencer Graves]

  Have you received a reply to this post?  I haven't seen one.  I don't 
have an answer for you, but if you'd still like help from this list, I 
suggest you prepare the simplest possible toy example that you can 
conceive and send it to this list, restating your question in terms of 
that example.  You question indicates you've read the code for glmmPQL 
and seem to know enough to experiment with modifying the "glmmPQL" code 
or with extracting crucial pieces to make part of the simplified example 
illustrating your question (consistent with the posting guide, 
"www.R-project.org/posting-guide.html").

  If you are not sure what "glmmPQL" does, you can say "debug(glmmPQL)" 
before executing a command that invokes "glmmPQL".  That will open a 
browser that will allow you to look at and change anything in the 
environment of "glmmPQL" before and after any command;  if a command 
commits a fatal error, you will be evicted from "glmmPQL" and will have 
to start over.  This is the quickest way I know to understand and debug 
R code.

  hope this helps.
  spencer graves

Patrick Giraudoux wrote:

> Dear listers,
> 
> glmmPQL (package MASS) is given to work by repeated call to lme. In the 
> classical outputs glmmPQL  the Variance Structure is given  as " fixed 
> weights,  Formula: ~invwt".  The script shows that the function 
> varFixed() is used, though the place where 'invwt' is defined remains 
> unclear to me.  I wonder if there is an easy way to specify another 
> variance structure (eg varPower, etc..), preferably using an lme object 
> of the varFunc classes ? Some trials show that the 'weights' argument of 
> glmmPQL is just the same as in glm (which is clearly stated in the help) 
> and I wonder actually, if not a nonsense, how to pass eg a 'weights' 
> arguments as used in lme (eg weights=varPower()) to specify a variance 
> function (in the same way as a correlation structure can be passed easy).
> 
> Thanks in advance for any hint,
> 
> Patrick
> 
> __
> R-help@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

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

[EMAIL PROTECTED]
www.pdf.com 
Tel:  408-938-4420
Fax: 408-280-7915

__
R-help@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] glmmPQL and variance structure

[Patrick Giraudoux]

Dear listers,

glmmPQL (package MASS) is given to work by repeated call to lme. In the 
classical outputs glmmPQL  the Variance Structure is given  as " fixed 
weights,  Formula: ~invwt".  The script shows that the function 
varFixed() is used, though the place where 'invwt' is defined remains 
unclear to me.  I wonder if there is an easy way to specify another 
variance structure (eg varPower, etc..), preferably using an lme object 
of the varFunc classes ? Some trials show that the 'weights' argument of 
glmmPQL is just the same as in glm (which is clearly stated in the help) 
and I wonder actually, if not a nonsense, how to pass eg a 'weights' 
arguments as used in lme (eg weights=varPower()) to specify a variance 
function (in the same way as a correlation structure can be passed easy).

Thanks in advance for any hint,

Patrick

__
R-help@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