Commit | Line | Data |
---|---|---|
cbd88fe5 BA |
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 | } |