Use variable in GLM quasi specification

Question

I'm fitting a GLM to some data, using a quasi-likelihood approach (family=quasi(...)).

I'd like to use a variable, p in the variance specification, like so:

family = quasi(link=log, variance=mu^p)

This however doesn't work (it no longer recongises mu).

Is there any way to get R to just insert the value of p in the expression before it is evaluated, so I can use pinstead of a number?

Here's an example that doesn't work:

set.seed(1)
x <- runif(100)
y <- x^2+2*x+sin(2*pi*x) + rnorm(100)

fitModel <- function(x,y, p) {
  model <- glm(y~x, family=quasi(link=log, variance=mu^p))
  return(model)
}
fitModel(x,y,2)

Thanks!

Otherwise you need to provide "a list containing components `varfun`, `validmu`, `dev.resids`, `initialize` and `name`" (see `?family`). — Roland, Sep 23 '13 at 11:01
Ideally I'd like to use `variance = paste0("mu^", p)`this in a function that calls `glm`. This gives me an error: 'variance' "NA" is invalid: possible values are "mu(1-mu)", "mu", "mu^2", "mu^3" and "constant" - probably because p is NA at this point. Giving it a default value doesn't seem to help either. Any ideas? — user2249626, Sep 23 '13 at 11:13
Add all relevant information (including code and toy data for reproducibility) to your question. — Roland, Sep 23 '13 at 11:18

Ben Bolker · Accepted Answer · 2013-09-27T14:45:41.963

The family function does fancy parsing which means the paste0 solution suggested in the comments won't work without jumping through considerable hoops. Also, the following function fails if any of the y values are <= 0, so I changed the example a little bit (if you do have negative response values you'll have to think about what you want to do about this ...)

set.seed(1)
x <- seq(2,10,length=100)
y <- x^2+2*x+sin(2*pi*x) + rnorm(100,)

What I did was to create a quasi family object, then modify its variance function on the fly.

pfamily <- quasi(link="log",variance="mu")
fitModel <- function(x,y, p) {
    pfamily[["variance"]] <- function(mu) mu^p
    model <- glm(y~x, family=pfamily)
    model
}

fitModel(x,y,2)
fitModel(x,y,1)

For what it's worth, this variant should be able to do arbitrary values of p, so e.g. you can draw a curve over the variance power:

dfun <- function(p) {
   deviance(fitModel(x,y,p))
}
pvec <- seq(0.1,3,by=0.1)
dvec <- sapply(pvec,dfun)
par(las=1,bty="l")
plot(pvec,dvec,type="b",xlab="variance power",ylab="deviance")

enter image description here

Use variable in GLM quasi specification

1 Answers1