# [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:
{

if (lambda >= 0){
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 {
eta^(1/lambda)
mu.eta <- function(eta)
(1/lambda) * eta^(1/lambda - 1)
valideta <- function(eta) all(eta>0)
}
}
else
"logit" = { ...

The simple change to power:
> power1
function(lambda = 1) {
if(lambda == 0)
else if(lambda == 1)
else
}

> testcar
function(pow){
ob <-
eval(substitute(glm(Pound~CG+Age+Vage,data=car,weights=No,
subset=No>0,
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```