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!

Answer 1

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")