Hola!
To my surprise, the power() function cannot be used to make links with
nehgative
powers (it returns the log link if one tries). I can see no good reason
for this, and have tried with simple changes to make it work usefully
for negative arguments. Are there any good reasons I have overlooked
not to allow negative arguments?
The first part of make.link must be augmented as
follows:> make.link
function (link)
{
if (is.character(link) && length(grep("^power", link) >
0))
return(eval(parse(text = link)))
else if(!is.character(link) && !is.na(lambda <-
as.numeric(link))) {
if (lambda >= 0){
linkfun <- function(mu) mu^lambda
linkinv <- function(eta)
pmax(.Machine$double.eps, eta^(1/lambda))
mu.eta <- function(eta)
pmax(.Machine$double.eps, (1/lambda) * eta^(1/lambda - 1))
valideta <- function(eta) all(eta>0)
} else {
linkfun <- function(mu) mu^lambda
linkinv <- function(eta)
eta^(1/lambda)
mu.eta <- function(eta)
(1/lambda) * eta^(1/lambda - 1)
valideta <- function(eta) all(eta>0)
}
}
else
switch(link,
"logit" = { ...
The simple change to power:> power1
function(lambda = 1) {
if(lambda == 0)
make.link("log")
else if(lambda == 1)
make.link("identity")
else
make.link(lambda)
}
> testcar
function(pow){
ob <-
eval(substitute(glm(Pound~CG+Age+Vage,data=car,weights=No,
subset=No>0,
family=quasi(link=power1(pow),
var=mu^2))),list(pow=pow))
deviance(ob)
}
The history:
devs <- numeric(31)
for (i in seq(along=devs)) {devs[i] <- testcar(-2.1+i/10)}
plot(x=-2.1+(1:31)/10, y=devs,xlab="power",ylab="dev")
The data to remake the graphics:
car <- cbind(expand.grid(CG=c("A",
"B","C","D"),
Age=c("17-20","21-24","25-29","30-34",
"35-39","40-49","50-59","60+"),
Vage=c("0-3","4-7","8-9","10+")),
Pound=c(289,372,189,763,302,420,268,407,268,275,334,
383,236,259,340,400,207,208,251,233,254,218,
239,387,251,196,268,391,264,224,269,385,282,
249,288,850,194,243,343,320,285,243,274,305,
270,226,260,349,129,214,232,325,213,209,250,
299,227,229,250,228,198,193,258,324,133,288,179,
0, 135,196,293,205,181,179,208,116,160,161,
189,147,157,149,204,207,149,172,174,325,172,
164,175,346,167,178,227,192,160,11, 0 ,0 ,
166,135,104,0 ,110,264,150,636,110,107,104,65,
113,137,141,0 ,98 ,110,129,137,98 ,132,152,
167,114,101,119,123),
No=c(8,10,9,3,18,59,44,24,56,125,163,72,43,179,197,
104,43,191,210,119,90,380,401,199,69,366,310,
105,64,228,183,62,8,28,13,2,31,96,39,18,55,172,
129,50,53,211,125,55,73,219,131,43,98,434,253,
88,120,353,148,46,100,233,103,22,4,1,1,0,10,13,7,2,17,36,18,
6,15,39,30,8,21,46,32,4,35,97,50,8,42,95,33,10,43,
73,20,6,1,1,0,0,4,3,2,0,12,10,8,1,12,19,9,2,14,23,8,0,22,59,15,9,35,45,
13,1,53,44,6,6))
(This is the car insurance data from page 298 of McCullagh&Nelder). I
have tried other examples also, and
it seems to work well.)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._