Fix computeW (was random - TODO)
[morpheus.git] / patch_Bettina / FLXMRglm.R
CommitLineData
cbd88fe5
BA
1FLXMRglm <- function(formula=.~., family=gaussian, offset=NULL)
2{
d294ece1
BA
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)
cbd88fe5 50 }
d294ece1
BA
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))))
cbd88fe5 57 }
d294ece1
BA
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")
cbd88fe5
BA
69 x
70 }
71
d294ece1
BA
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)
cbd88fe5 81 }
d294ece1
BA
82 logLik <- function(x, y, ...)
83 dbinom(y[,1], size=rowSums(y), prob=predict(x, ...), log=TRUE)
cbd88fe5 84
d294ece1
BA
85 new("FLXcomponent",
86 parameters=list(coef=para$coef),
87 logLik=logLik, predict=predict,
88 df=para$df)
cbd88fe5 89 }
cbd88fe5 90
d294ece1
BA
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)))
cbd88fe5 94 }
d294ece1
BA
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)
cbd88fe5 105 }
d294ece1
BA
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)
cbd88fe5 113 }
cbd88fe5 114
d294ece1
BA
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)
cbd88fe5 129 }
d294ece1
BA
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)
cbd88fe5 137 }
d294ece1
BA
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
cbd88fe5 148}