patch_Bettina/FLXMRglm.R

   1 FLXMRglm <- function(formula=.~., family=gaussian, offset=NULL)
   2 {
   3     if (is.character(family))
   4         family <- get(family, mode = "function", envir = parent.frame())
   5     if (is.function(family))
   6         family <- family()
   7     if (is.null(family$family)) {
   8         print(family)
   9         stop("'family' not recognized")
  10     }
  11     glmrefit <- function(x, y, w) {
  12       fit <- c(glm.fit(x, y, weights=w, offset=offset,
  13                        family=family),
  14                list(call = sys.call(), offset = offset,
  15                     control = eval(formals(glm.fit)$control),
  16                     method = "weighted.glm.fit"))
  17       fit$df.null <- sum(w) + fit$df.null - fit$df.residual - fit$rank
  18       fit$df.residual <- sum(w) - fit$rank
  19       fit$x <- x
  20       fit
  21     }
  22
  23     z <- new("FLXMRglm", weighted=TRUE, formula=formula,
  24              name=paste("FLXMRglm", family$family, sep=":"), offset = offset,
  25              family=family$family, refit=glmrefit)
  26     z@preproc.y <- function(x){
  27       if (ncol(x) > 1)
  28         stop(paste("for the", family$family, "family y must be univariate"))
  29       x
  30     }
  31
  32     if(family$family=="gaussian"){
  33       z@defineComponent <- function(para) {
  34         predict <- function(x, ...) {
  35           dotarg = list(...)
  36           if("offset" %in% names(dotarg)) offset <- dotarg$offset
  37           p <- x %*% para$coef
  38           if (!is.null(offset)) p <-  p + offset
  39           family$linkinv(p)
  40         }
  41
  42         logLik <- function(x, y, ...)
  43           dnorm(y, mean=predict(x, ...), sd=para$sigma, log=TRUE)
  44
  45         new("FLXcomponent",
  46             parameters=list(coef=para$coef, sigma=para$sigma),
  47             logLik=logLik, predict=predict,
  48             df=para$df)
  49       }
  50
  51       z@fit <- function(x, y, w, component){
  52           fit <- glm.fit(x, y, w=w, offset=offset, family = family)
  53           z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1,
  54                                         sigma =  sqrt(sum(fit$weights * fit$residuals^2 /
  55                                                           mean(fit$weights))/ (nrow(x)-fit$rank))))
  56       }
  57     }
  58     else if(family$family=="binomial"){
  59       z@preproc.y <- function(x){
  60         if (ncol(x) != 2)
  61           stop("for the binomial family, y must be a 2 column matrix\n",
  62                "where col 1 is no. successes and col 2 is no. failures")
  63         if (any(x < 0))
  64           stop("negative values are not allowed for the binomial family")
  65         x
  66       }
  67       z@defineComponent <- function(para) {
  68         predict <- function(x, ...) {
  69           dotarg = list(...)
  70           if("offset" %in% names(dotarg)) offset <- dotarg$offset
  71           p <- x %*% para$coef
  72           if (!is.null(offset)) p <- p + offset
  73           family$linkinv(p)
  74         }
  75         logLik <- function(x, y, ...)
  76           dbinom(y[,1], size=rowSums(y), prob=predict(x, ...), log=TRUE)
  77
  78         new("FLXcomponent",
  79             parameters=list(coef=para$coef),
  80             logLik=logLik, predict=predict,
  81             df=para$df)
  82       }
  83
  84       z@fit <- function(x, y, w, component){
  85         fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
  86         z@defineComponent(para = list(coef = coef(fit), df = ncol(x)))
  87       }
  88     }
  89     else if(family$family=="poisson"){
  90       z@defineComponent <- function(para) {
  91         predict <- function(x, ...) {
  92           dotarg = list(...)
  93           if("offset" %in% names(dotarg)) offset <- dotarg$offset
  94           p <- x %*% para$coef
  95           if (!is.null(offset)) p <- p + offset
  96           family$linkinv(p)
  97         }
  98         logLik <- function(x, y, ...)
  99           dpois(y, lambda=predict(x, ...), log=TRUE)
 100
 101         new("FLXcomponent",
 102             parameters=list(coef=para$coef),
 103             logLik=logLik, predict=predict,
 104             df=para$df)
 105       }
 106
 107       z@fit <- function(x, y, w, component){
 108         fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
 109         z@defineComponent(para = list(coef = coef(fit), df = ncol(x)))
 110       }
 111     }
 112     else if(family$family=="Gamma"){
 113       z@defineComponent <- function(para) {
 114         predict <- function(x, ...) {
 115           dotarg = list(...)
 116           if("offset" %in% names(dotarg)) offset <- dotarg$offset
 117           p <- x %*% para$coef
 118           if (!is.null(offset)) p <- p + offset
 119           family$linkinv(p)
 120         }
 121         logLik <- function(x, y, ...)
 122           dgamma(y, shape = para$shape, scale=predict(x, ...)/para$shape, log=TRUE)
 123
 124         new("FLXcomponent",
 125             parameters = list(coef = para$coef, shape = para$shape),
 126             predict = predict, logLik = logLik,
 127             df = para$df)
 128       }
 129
 130       z@fit <- function(x, y, w, component){
 131         fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
 132         z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1,
 133                               shape = sum(fit$prior.weights)/fit$deviance))
 134
 135       }
 136     }
 137     else stop(paste("Unknown family", family))
 138     z
 139 }