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