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
  12   glmrefit <- function(x, y, w) {
  13     fit <- c(glm.fit(x, y, weights=w, offset=offset, 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
  27   z@preproc.y <- function(x) {
  28     if (ncol(x) > 1)
  29       stop(paste("for the", family$family, "family y must be univariate"))
  30     x
  31   }
  32
  33   if (family$family=="gaussian") {
  34     z@defineComponent <- function(para) {
  35       predict <- function(x, ...) {
  36         dotarg = list(...)
  37         if("offset" %in% names(dotarg)) offset <- dotarg$offset
  38         p <- x %*% para$coef
  39         if (!is.null(offset)) p <-  p + offset
  40         family$linkinv(p)
  41       }
  42
  43       logLik <- function(x, y, ...)
  44         dnorm(y, mean=predict(x, ...), sd=para$sigma, log=TRUE)
  45
  46       new("FLXcomponent",
  47         parameters=list(coef=para$coef, sigma=para$sigma),
  48         logLik=logLik, predict=predict,
  49         df=para$df)
  50     }
  51
  52     z@fit <- function(x, y, w, component){
  53       fit <- glm.fit(x, y, w=w, offset=offset, family = family)
  54       z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1,
  55         sigma =  sqrt(sum(fit$weights * fit$residuals^2 /
  56         mean(fit$weights))/ (nrow(x)-fit$rank))))
  57     }
  58   }
  59
  60   else if (family$family=="binomial") {
  61     z@preproc.y <- function(x) {
  62       if (ncol(x) != 2)
  63       {
  64         stop("for the binomial family, y must be a 2 column matrix\n",
  65           "where col 1 is no. successes and col 2 is no. failures")
  66       }
  67       if (any(x < 0))
  68         stop("negative values are not allowed for the binomial family")
  69       x
  70     }
  71
  72     z@defineComponent <- function(para) {
  73       predict <- function(x, ...) {
  74         dotarg = list(...)
  75         if("offset" %in% names(dotarg))
  76           offset <- dotarg$offset
  77         p <- x %*% para$coef
  78         if (!is.null(offset))
  79           p <- p + offset
  80         family$linkinv(p)
  81       }
  82       logLik <- function(x, y, ...)
  83       dbinom(y[,1], size=rowSums(y), prob=predict(x, ...), log=TRUE)
  84
  85       new("FLXcomponent",
  86         parameters=list(coef=para$coef),
  87         logLik=logLik, predict=predict,
  88         df=para$df)
  89     }
  90
  91     z@fit <- function(x, y, w, component) {
  92       fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
  93       z@defineComponent(para = list(coef = coef(fit), df = ncol(x)))
  94     }
  95   }
  96
  97   else if (family$family=="poisson") {
  98     z@defineComponent <- function(para) {
  99       predict <- function(x, ...) {
 100         dotarg = list(...)
 101         if("offset" %in% names(dotarg)) offset <- dotarg$offset
 102         p <- x %*% para$coef
 103         if (!is.null(offset)) p <- p + offset
 104         family$linkinv(p)
 105       }
 106       logLik <- function(x, y, ...)
 107         dpois(y, lambda=predict(x, ...), log=TRUE)
 108
 109       new("FLXcomponent",
 110         parameters=list(coef=para$coef),
 111         logLik=logLik, predict=predict,
 112         df=para$df)
 113     }
 114
 115     z@fit <- function(x, y, w, component) {
 116       fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
 117       z@defineComponent(para = list(coef = coef(fit), df = ncol(x)))
 118     }
 119   }
 120
 121   else if (family$family=="Gamma") {
 122     z@defineComponent <- function(para) {
 123       predict <- function(x, ...) {
 124         dotarg = list(...)
 125         if("offset" %in% names(dotarg)) offset <- dotarg$offset
 126         p <- x %*% para$coef
 127         if (!is.null(offset)) p <- p + offset
 128         family$linkinv(p)
 129       }
 130       logLik <- function(x, y, ...)
 131         dgamma(y, shape = para$shape, scale=predict(x, ...)/para$shape, log=TRUE)
 132
 133       new("FLXcomponent",
 134         parameters = list(coef = para$coef, shape = para$shape),
 135         predict = predict, logLik = logLik,
 136         df = para$df)
 137     }
 138
 139     z@fit <- function(x, y, w, component) {
 140       fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef)
 141       z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1,
 142         shape = sum(fit$prior.weights)/fit$deviance))
 143     }
 144   }
 145
 146   else stop(paste("Unknown family", family))
 147   z
 148 }