[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