-
Notifications
You must be signed in to change notification settings - Fork 48
/
predab.resample.s
416 lines (347 loc) · 12.4 KB
/
predab.resample.s
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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
predab.resample <-
function(fit.orig,
fit,
measure,
method=c("boot","crossvalidation",".632","randomization"),
bw=FALSE,
B=50,
pr=FALSE, prmodsel=TRUE,
rule="aic",
type="residual",
sls=.05,
aics=0,
tol=1e-12,
force=NULL,
estimates=TRUE,
non.slopes.in.x=TRUE,
kint=1,
cluster,
subset,
group=NULL,
allow.varying.intercepts=FALSE,
debug=FALSE,
...)
{
method <- match.arg(method)
oldopt <- options('digits')
options(digits=4)
on.exit(options(oldopt))
efit <- function(...) list(fail=TRUE)
## Following logic prevents having to load a copy of a large x object
if(any(match(c("x", "y"), names(fit.orig), 0) == 0))
stop("must have specified x=T and y=T on original fit")
fparms <- fit.orig[c("terms", "Design")]
oassign <- fit.orig$assign
non.slopes <- num.intercepts(fit.orig, 'coef')
# x.index <- if(non.slopes==0 || non.slopes.in.x) function(i,...) i
# else
# function(i, ns) {
# if(any(i > ns)) i[i > ns] - ns
# else NULL
# }
x.index <- function(i, ns)
if(ns == 0) i else setdiff(i, 1 : ns) - ns
Xb <- function(x, b, non.slopes, n, kint=1) {
if(length(x)) matxv(x, b, kint=kint)
else if(non.slopes == 0 || ! length(kint)) rep(0, n)
else rep(b[kint], n)
}
# if(length(x)) {
# if(non.slopes == 0 || non.slopes.in.x) x %*% b
# else b[kint] + x %*% b[-(1 : non.slopes)]
# }
# else {
# if(non.slopes==0) rep(0, n)
# else
# rep(b[kint], n)
# }
# }
nac <- fit.orig$na.action
x <- as.matrix(fit.orig[['x']])
n <- nrow(x)
## Remove model.matrix class for subset operations later
attr(x,'class') <- NULL
y <- fit.orig[['y']]
if(is.factor(y)) y <- unclass(y)
if(! is.Surv(y)) y <- as.matrix(y)
## some subjects have multiple records now
multi <- ! missing(cluster)
if(length(group)) {
if(multi || method != 'boot')
stop('group is currently allowed only when method="boot" and cluster is not given')
if(length(group) > n) {
## Missing observations were deleted during fit
if(length(nac))
j <- ! is.na(naresid(nac, y) %*% rep(1, ncol(y)))
group <- group[j]
}
if(length(group) != n)
stop('length of group does not match # rows used in fit')
group.inds <- split(1 : n, group) # see bootstrap()
ngroup <- length(group.inds)
}
else ngroup <- 0
if(multi) {
if(method != 'boot')
stop('cluster only implemented for method="boot"')
if(length(cluster) > n) {
## Missing observations were deleted during fit
if(length(nac)) {
j <- ! is.na(naresid(nac, y) %*% rep(1, ncol(y)))
cluster <- cluster[j]
}
}
if(length(cluster) != n)
stop('length of cluster does not match # rows used in fit')
if(any(is.na(cluster)))
stop('cluster has NAs')
n.orig <- length(unique(cluster))
cl.samp <- split(1 : n, cluster)
}
else n.orig <- n
if(! missing(subset)) {
if(length(subset) > n && length(nac)) {
j <- ! is.na(naresid(nac, y) %*% rep(1, ncol(y)))
subset <- subset[j]
}
if(length(subset) != n && all(subset >= 0))
stop('length of subset does not match # rows used in fit')
if(any(is.na(subset))) stop('subset has NAs')
if(! is.logical(subset)) {
subset2 <- rep(FALSE, n)
subset2[subset] <- TRUE
subset <- subset2
subset2 <- NULL
}
}
stra <- fit.orig$strata
if(bw) {
if(fit.orig$fail) return()
if(prmodsel) cat("\n Backwards Step-down - Original Model\n")
fbw <- fastbw(fit.orig, rule=rule, type=type, sls=sls, aics=aics,
eps=tol, force=force)
if(prmodsel) print(fbw, estimates=estimates)
orig.col.kept <- fbw$parms.kept
if(! length(orig.col.kept))
stop("no variables kept in original model")
## Check that x.index works if allow.varying.intercepts
xcol <- x.index(orig.col.kept, non.slopes)
## Refit subset of predictors on whole sample
fit.orig <- fit(x[, xcol, drop=FALSE], y, strata=stra,
iter=0, tol=tol, xcol=xcol, ...)
}
else orig.col.kept <- seq(along=fit.orig$coef)
b <- fit.orig$coef
xcol <- x.index(orig.col.kept, non.slopes)
xb <- Xb(x[, xcol, drop=FALSE], b, non.slopes, n,
kint=kint)
index.orig <- if(missing(subset))
measure(xb, y, strata=stra, fit=fit.orig, iter=0, evalfit=TRUE,
fit.orig=fit.orig,
kint=kint, ...)
else
measure(xb[subset], y[subset,,drop=FALSE], strata=stra[subset],
fit=fit.orig,
iter=0, evalfit=FALSE, fit.orig=fit.orig, kint=kint, ...)
keepinfo <- attr(index.orig, 'keepinfo')
test.stat <- double(length(index.orig))
train.stat <- test.stat
name <- fparms$Design$name
if(bw) varin <- matrix(FALSE, nrow=B, ncol=length(name))
j <- 0
num <- 0
if(method == "crossvalidation") {
per.group <- n / B
if(per.group < 2) {
stop("B > n/2")
}
sb <- sample(n, replace=FALSE)
}
##Cross-val keeps using same random set of indexes, without replacement
ntest <- 0 #Used in getting weighted average for .632 estimator
if(method == ".632") {
## Must do assignments ahead of time so can weight estimates
## according to representation in bootstrap samples
S <- matrix(integer(1), nrow=n, ncol=B)
W <- matrix(TRUE, nrow=n, ncol=B)
for(i in 1 : B) {
S[, i] <- s <- sample(n, replace=TRUE)
W[s, i] <- FALSE #now these obs are NOT omitted
}
nomit <- drop(W %*% rep(1,ncol(W))) #no. boot samples omitting each obs
if(min(nomit) == 0)
stop("not every observation omitted at least once ",
"in bootstrap samples.\nRe--run with larger B")
W <- apply(W / nomit, 2, sum) / n
if(pr) {
cat("\n\nWeights for .632 method (ordinary bootstrap weights ",
format(1 / B), ")\n", sep="")
print(summary(W))
}
}
pb <- setPb(B, type=if(method == 'crossvalidation') 'Cross' else 'Boot',
onlytk=! pr,
every=1*(B < 20) + 5*(B >= 20 & B < 50) +
10*(B >= 50 & B < 100) + 20*(B >= 100 & B < 1000) +
50*(B >= 1000))
for(i in 1 : B) {
pb(i)
switch(method,
crossvalidation = {
is <- 1 + round((i - 1) * per.group)
ie <- min(n, round(is + per.group - 1))
test <- sb[is : ie]
train <- -test
}, #cross-val
boot = {
if(ngroup) {
train <- integer(n.orig)
for(si in 1 : ngroup) {
gi <- group.inds[[si]]
lgi <- length(gi)
train[gi] <- if(lgi == 1) gi
else {
## sample behaves differently when first arg is
## a single integer
sample(gi, lgi, replace=TRUE)
}
}
}
else {
train <- sample(n.orig, replace=TRUE)
if(multi) train <- unlist(cl.samp[train])
}
test <- 1 : n
}, #boot
".632" = {
train <- S[, i]
test <- -train
}, #boot .632
randomization =
{
train <- sample(n, replace=FALSE)
test <- 1 : n
}
) #randomization
xtrain <- if(method == "randomization") 1 : n
else train
if(debug) {
cat('\nSubscripts of training sample:\n')
print(train)
cat('\nSubscripts of test sample:\n')
print(test)
}
f <- tryCatch(fit(x[xtrain, , drop=FALSE], y[train, , drop=FALSE],
strata=stra[train], iter=i, tol=tol, ...), error=efit)
if(! length(f$fail)) f$fail <- FALSE
f$assign <- NULL #Some programs put a NULL assign (e.g. ols.val fit)
ni <- num.intercepts(f)
fail <- f$fail
if(! fail) {
## Following if..stop was before f$assign above
if(! allow.varying.intercepts && ni != non.slopes) {
stop('A training sample has a different number of intercepts (', ni ,')\n',
'than the original model fit (', non.slopes, ').\n',
'You probably fit an ordinal model with sparse cells and a re-sample\n',
'did not select at least one observation for each value of Y.\n',
'Add the argument group=y where y is the response variable.\n',
'This will force balanced sampling on levels of y.')
}
clf <- attr(f, "class") # class is removed by c() below
f[names(fparms)] <- fparms
assign <- oassign
## Slopes are shifted to the left when fewer unique values of Y
## occur (especially for orm models) resulting in fewer intercepts
if(non.slopes != ni) for(z in 1 : length(assign))
assign[[z]] <- assign[[z]] - (non.slopes - ni)
f$assign <- assign
attr(f, "class") <- clf
if(! bw) {
coef <- f$coef
col.kept <- seq(along=coef)
}
else {
f <- fastbw(f, rule=rule, type=type, sls=sls, aics=aics,
eps=tol, force=force)
if(pr && prmodsel) print(f, estimates=estimates)
varin[j + 1, f$factors.kept] <- TRUE
col.kept <- f$parms.kept
if(! length(col.kept))
f <- tryCatch(fit(NULL, y[train,, drop=FALSE], stra=stra[xtrain],
iter=i, tol=tol,...), error=efit)
else {
xcol <- x.index(col.kept, ni)
f <- tryCatch(fit(x[xtrain, xcol, drop=FALSE], strata=stra[xtrain],
y[train,, drop=FALSE],
iter=i, tol=tol, xcol=xcol, ...), error=efit)
}
if(length(f$fail) && f$fail) fail <- TRUE
else coef <- f$coef
}
}
if(! fail) {
j <- j + 1
xcol <- x.index(col.kept, ni)
xb <- Xb(x[,xcol,drop=FALSE], coef, ni, n, kint=kint)
if(missing(subset)) {
train.statj <-
measure(xb[xtrain], y[train,,drop=FALSE], strata=stra[xtrain],
fit=f, iter=i, fit.orig=fit.orig, evalfit=TRUE,
kint=kint, ...)
test.statj <- measure(xb[test], y[test,,drop=FALSE],
strata=stra[test],
fit=f, iter=i, fit.orig=fit.orig,
evalfit=FALSE,
kint=kint, ...)
}
else {
ii <- xtrain
if(any(ii < 0)) ii <- (1 : n)[ii]
ii <- ii[subset[ii]]
train.statj <- measure(xb[ii], y[ii,,drop=FALSE],
strata=stra[ii],
fit=f, iter=i, fit.orig=fit.orig,
evalfit=FALSE,
kint=kint, ...)
ii <- test
if(any(ii < 0)) ii <- (1 : n)[ii]
ii <- ii[subset[ii]]
test.statj <- measure(xb[ii], y[ii,,drop=FALSE], fit=f,
iter=i, strata=stra[ii],
fit.orig=fit.orig, evalfit=FALSE,
kint=kint, ...)
}
na <- is.na(train.statj + test.statj)
num <- num + ! na
if(pr)
print(cbind(training=train.statj, test=test.statj))
train.statj[na] <- 0
test.statj[na] <- 0
if(method == ".632") {
wt <- W[i]
if(any(na))
warning('method=".632" does not properly handle missing summary indexes')
}
else wt <- 1
train.stat <- train.stat + train.statj
test.stat <- test.stat + test.statj * wt
ntest <- ntest + 1
}
}
if(pr) cat("\n\n")
if(j != B)
cat("\nDivergence or singularity in", B - j, "samples\n")
train.stat <- train.stat / num
if(method != ".632") {
test.stat <- test.stat / num
optimism <- train.stat - test.stat
}
else optimism <- .632 * (index.orig - test.stat)
res <- cbind(index.orig=index.orig, training=train.stat, test=test.stat,
optimism=optimism, index.corrected=index.orig-optimism, n=num)
if(bw) {
varin <- varin[1 : j, , drop=FALSE]
dimnames(varin) <- list(rep("", j), name)
}
structure(res, class='validate', kept=if(bw) varin, keepinfo=keepinfo)
}