[R] power() function

halvorsen kjetilh at umsanet.edu.bo
Tue Jul 25 00:04:03 CEST 2000


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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list