1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
|
.onLoad <- function(lib, pkg) {
if(is.null(getOption("max.print")))
options(max.print = 10000)#-> show() of large matrices
}
## Model Matrix
setClass("modelMatrix",
representation(assign = "integer",
contrasts = "list", "VIRTUAL"),
contains = "Matrix",
validity = function(object) {
if(length(object@assign) != (p <- ncol(object)))
return(gettextf("'%s' slot must be integer of length %d",
"assign", p))
contr <- object@contrasts
c.cl <- vapply(contr, function(x) class(x)[[1L]], "<class>", USE.NAMES=FALSE)
if(length(nc <- names(contr)) != length(c.cl) || !all(nchar(nc) > 0))
return(gettextf("'%s' slot must be a correctly named list",
"contrasts"))
## TODO? length(contrasts) < maximal value in 'assign' <= p
contrCls <- c("character", "function", "matrix", "Matrix")
if(any(vapply(c.cl, function(cl) all(is.na(match(extends(cl), contrCls))), NA)))
return(gettextf(
"'%s' slot must be named list of contrast functions or their names, or matrices",
"contrasts"))
TRUE
})
setClass("sparseModelMatrix", representation("VIRTUAL"),
contains = c("CsparseMatrix", "modelMatrix"))
setClass("denseModelMatrix", representation("VIRTUAL"),
contains = c("unpackedMatrix", "modelMatrix"))
## The currently only *actual* denseModelMatrix class:
setClass("ddenseModelMatrix", contains = c("dgeMatrix", "denseModelMatrix"))
' The following gives --- with latest R-devel (2020-01-07 r77634) ----
Found the following significant warnings:
Warning: unable to find a consistent ordering of superclasses for class
"dsparseModelMatrix": order chosen is inconsistent with the superclasses of
"sparseModelMatrix"
'
## The currently only *actual* sparseModelMatrix class:
setClass("dsparseModelMatrix", contains = c("dgCMatrix", "sparseModelMatrix"))
###------ Modules related to modelling --- basically in two parts --------------
###------ 1) "prediction-Module" -- currently in a sparse and dense flavor
###------ 2) "response-Module"
## Linear predictor modules, which consist of the model matrix, the
## coefficient vector and a triangular factor of the weighted model matrix.
## the super class contains the slots already;
setClass("predModule",
representation(X = "modelMatrix", coef = "numeric", Vtr = "numeric",
fac = "CholeskyFactorization",
"VIRTUAL"))
## the sub classes specify more specific classes for the two non-trivial slots:
setClass("dPredModule", contains = "predModule",
representation(X = "ddenseModelMatrix", fac = "Cholesky"))
setClass("sPredModule", contains = "predModule",
representation(X = "dsparseModelMatrix", fac = "CHMfactor"))
## Response modules for models with a linear predictor, which can
## include linear models, generalized linear models, nonlinear models
## and generalized nonlinear models.
## y, offset and mu are as expected. Note that length(offset) can be a multiple of length(y)
## weights are the prior weights
## sqrtrwt and sqrtXwt are the square roots of residual and X weights
setClass("respModule",
representation(mu = "numeric", # of length n
offset = "numeric", # of length n * s
sqrtXwt = "matrix", # of dim(.) == (n, s)
sqrtrwt = "numeric", # sqrt(residual weights)
weights = "numeric", # prior weights
wtres = "numeric",
y = "numeric"),
validity = function(object) {
n <- length(object@y)
if (any(n != sapply(lapply(c("weights","sqrtrwt","mu","wtres"
), slot, object = object), length)))
return("lengths of weights, sqrtwt and mu must match length(y)")
lo <- length(object@offset)
if (!lo || lo %% n)
return("length(offset) must be a positive multiple of length(y)")
if (length(object@sqrtXwt) != lo)
return("length(sqrtXwt) must equal length(offset)")
if (nrow(object@sqrtXwt) != n)
return("nrow(sqrtXwt) != length(y)")
TRUE
})
setOldClass("family")
##' glm response module
setClass("glmRespMod",
representation(family = "family",
eta = "numeric",
n = "numeric"), # for evaluation of the aic
contains = "respModule",
validity = function(object) {
if (length(object@eta) != length(object@y))
return("lengths of eta and y must match")
})
##' nls response module
setClass("nlsRespMod",
representation(nlenv = "environment",
nlmod = "call",
pnames = "character"),
contains = "respModule",
validity = function(object) {
n <- length(object@y)
N <- length(object@offset)
s <- N %/% n
lpn <- length(object@pnames)
if (lpn != s) return(sprintf("length(pnames) = %d != s = %d", lpn, s))
dd <- dim(object@sqrtXwt)
if (!all(dd == c(n, s))) {
return(sprintf("dim(gradient) = (%d, %d), n = %d, s = %d",
dd[1], dd[2], n, s))
}
TRUE
})
##' nglm response module
setClass("nglmRespMod", contains = c("glmRespMod", "nlsRespMod"))
### FIXME: move this eventually to 'methods':
## -----
##' The mother class of all (S4 based) (statistical / physical / ...) models in R:
setClass("Model", representation(call = "call", fitProps = "list",
"VIRTUAL"))
##' Statistical models based on linear predictors
##' "glpModel" := General Linear Prediction Models
setClass("glpModel", representation(resp = "respModule", pred = "predModule"),
contains = "Model")
rMod <-
setRefClass("RespModule",
fields = list(
mu = "numeric", # of length n
n = "integer", # for evaluation of the aic
offset = "numeric", # of length n * s
sqrtXwt = "matrix", # of dim(.) == (n, s)
sqrtrwt = "numeric", # sqrt(residual weights)
weights = "numeric", # prior weights
wtres = "numeric",
y = "numeric"),
methods = list(
initialize = function(...) {
initFields(...)
if (length(n) == 0L) n <<- length(y)
s <- 0L
##currently fails at pkg INSTALL time: stopifnot(n > 0L)
if (length(weights) == 0L) weights <<- numeric(n) + 1
sqrtrwt <<- sqrt(weights)
if (any((dd <- dim(sqrtXwt)) < 1L))
sqrtXwt <<- matrix(1, ncol = 1L, nrow = n)
else {
stopifnot(nrow(sqrtXwt) == n)
s <- ncol(sqrtXwt)
}
swrk <- max(s, 1L)
if (length(offset) == 0) offset <<- numeric(n * swrk)
else {
so <- length(offset) %/% n
stopifnot(length(offset) %% n == 0, s == 0 || so == s)
}
wtres <<- mu <<- numeric(n) * NA_real_
.self
},
updateMu = function(gamma) {
gamma <- as.numeric(gamma)
stopifnot(length(gamma) == length(offset))
mu <<- gamma + offset
wtres <<- sqrtrwt * (y - mu)
},
updateWts = function() {}
))
rMod$lock("y","n","weights","offset")
glrMod <-
setRefClass("GLMrespMod",
fields = list(
family = "family",
eta = "numeric"),
contains = "RespModule",
methods = list(
initialize = function(...) {
callSuper(...)
args <- list(...)
stopifnot("family" %in% names(args), is(args$family, "family"))
family <<- args$family
.self
},
updateMu = function(gamma) {
gamma <- as.numeric(gamma)
stopifnot(length(gamma) == length(offset))
mu <<- family$linkinv(eta <<- offset + gamma)
wtres <<- sqrtrwt * (y - mu)
},
updateWts = function() {
sqrtrwt <<- rtrwt <- sqrt(weights/family$variance(mu))
sqrtXwt[] <<- rtrwt * family$mu.eta(eta)
wtres <<- rtrwt * (y - mu)
}
))
glrMod$lock("family")
|
