First commit
[morpheus.git] / 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 }