[R] Creating its own family

Stephane DRAY stephane.dray at umontreal.ca
Tue Mar 9 18:33:34 CET 2004


Hello,
I would like to create my own family for glm modelling

If we consider a matrix Y, I would like to model Yij/var(Yj) with an 
inverse variance link.
I have create my own family inspired by the negative.binomial of MASS:

#########################################################
myfamily=function (varcol)
{
     env <- new.env(parent = .GlobalEnv)

     assign(".varcol",varcol,envir=env)
     famname="myfamily"
     link="inv col var"
     linkfun=function(.varcol,mu){mu/(.varcol)}
     linkinv=function(.varcol,eta){eta*(.varcol)}
     variance <- function (mu) rep.int(1, length(mu))
     validmu <- function (mu) TRUE

     dev.resids <- function (y, mu, wt) wt * ((y - mu)^2)
     aic <- function (y, n, mu, wt, dev) sum(wt) * (log(dev/sum(wt) * 2 * 
pi) + 1) + 2
     mu.eta=function (eta) rep.int(1, length(eta))
     initialize <- expression({
     n <- rep.int(1, nobs)
     mustart <- y
     })

     environment(variance) <- environment(validmu) <- 
environment(dev.resids) <- environment(aic) <- env

     structure(list(family = famname, link = link, linkfun = linkfun,
         linkinv = linkinv, variance = variance, dev.resids = dev.resids,
         aic = aic, mu.eta = mu.eta, initialize = initialize,
         validmu = validmu), class = "family")
}
#####################################################

But it does not work when I try it

 > Y=matrix(runif(50),10,5)
 > glm(as.vector(t(Y))~x,family=myfamily(rep(apply(Y,2,var),10)))
Error in family$linkfun(mustart) : Argument "mu" is missing, with no default
 >

Hope that someone could help me,

Thanks in advances,
Sincerely.

Stéphane DRAY
-------------------------------------------------------------------------------------------------- 

Département des Sciences Biologiques
Université de Montréal, C.P. 6128, succursale centre-ville
Montréal, Québec H3C 3J7, Canada

Tel : 514 343 6111 poste 1233
E-mail : stephane.dray at umontreal.ca
-------------------------------------------------------------------------------------------------- 

Web                                          http://www.steph280.freesurf.fr/




More information about the R-help mailing list