« back to all changes in this revision
Viewing changes to R/GLMM.R
- browse files at revision 5
- compare with another revision
- download diff
- download tarball
- view history from revision 5
- Committer: Bazaar Package Importer
- Author(s): Douglas Bates
- Date: 2005-05-20 09:30:11 UTC
- mfrom: (1.1.3 upstream)
- Revision ID: james.westby@ubuntu.com-20050520093011-vp8clxugp8hbjejp
Tags: 0.95.8-1
* New upstream release
* Upstream release no longer requires r-cran-latticeextra. Closes: Bug#307497
* Upstream release no longer requires r-cran-latticeextra. Closes: Bug#307497
- files added:
- Done
- files removed:
- .Rbuildindex.1534
- data
- files modified:
- ChangeLog
added
removed
R/GLMM.R
1
setMethod("GLMM", signature(formula = "missing"),
2
function(formula, family, data, random,
3
method = c("PQL", "Laplace"),
4
control = list(),
5
subset,
6
weights,
7
na.action,
8
offset,
9
model = TRUE, x = FALSE, y = FALSE, ...)
10
{
11
Call = match.call()
12
resp = getResponseFormula(data)[[2]]
13
cov = getCovariateFormula(data)[[2]]
14
nCall$formula = eval(substitute(resp ~ cov))
15
.Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
16
mCall, PACKAGE = "Matrix")
17
})
18
19
20
setMethod("GLMM", signature(random = "formula"),
21
function(formula, family, data, random,
22
method = c("PQL", "Laplace"),
23
control = list(),
24
subset,
25
weights,
26
na.action,
27
offset,
28
model = TRUE, x = FALSE, y = FALSE, ...)
29
{
30
nCall = mCall = match.call()
31
cov = getCovariateFormula(random)
32
nCall$random <- lapply(getGroupsFormula(random, asList = TRUE),
33
function(f) cov)
34
.Call("nlme_replaceSlot", eval(nCall, parent.frame()), "call",
35
mCall, PACKAGE = "Matrix")
36
})
37
38
39
make.glm.call <-
40
function (mf, frm)
41
{
42
m <- match(c("formula", "family", "data", "weights",
43
"subset", "na.action", "offset"), names(mf), 0)
44
mf <- mf[c(1, m)]
45
mf[[1]] <- as.name("glm")
46
## environment(frm) = environment(formula) ???
47
mf$formula = frm
48
mf$model = FALSE
49
mf$x = FALSE
50
mf$y = TRUE
51
mf
52
}
53
54
setMethod("GLMM",
55
signature(formula = "formula",
56
random = "list"),
57
function(formula, family, data, random,
58
method = c("PQL", "Laplace"),
59
control = list(),
60
subset,
61
weights,
62
na.action,
63
offset,
64
model = TRUE, x = FALSE, y = FALSE, ...)
65
{
66
gVerb <- getOption("verbose")
67
random <-
68
lapply(random,
69
get("formula", pos = parent.frame(), mode = "function"))
70
controlvals <- do.call("lmeControl", control)
71
controlvals$REML <- FALSE
72
method <- match.arg(method)
73
74
## BEGIN glm fit without random effects
75
76
## several arguments are handled at this point.
77
## what about
78
79
## subset, na.action (arguments to model.frame)
80
81
## Are these already handled ?
82
## The rest could be part of ..., why have them explicitly ?
83
## Do we want possibility of supplying control in glm() ?
84
85
glm.fit <- eval(make.glm.call(match.call(expand.dots = TRUE),
86
formula), parent.frame())
87
family <- glm.fit$family
88
offset <- if (is.null(glm.fit$offset)) 0 else glm.fit$offset
89
weights <- sqrt(abs(glm.fit$prior.weights))
90
## initial 'fitted' values on linear scale
91
eta <- glm.fit$linear.predictors
92
## FIXME: does this include offsets (make sure it does)
93
etaold <- eta
94
95
## END using glm fit results
96
## Note: offset is on the linear scale
97
98
## FIXME: Not clear how (user specified) offset works. It
99
## doesn't make sense for it to be an offset for the
100
## response on the mu scale, so I'm assuming it's on the
101
## linear predictor scale
102
103
104
ndat <- eval(make.mf.call(match.call(expand.dots = FALSE),
105
formula, random), parent.frame())
106
facs <- lapply(names(random),
107
function(x) as.factor(eval(as.name(x), envir = ndat)))
108
names(facs) <- names(random)
109
## order factor list by decreasing number of levels
110
ford <- rev(order(sapply(facs, function(fac) length(levels(fac)))))
111
if (any(ford != seq(a = ford))) { # re-order both facs and random
112
facs <- facs[ford]
113
random <- random[ford]
114
}
115
## creates model matrices
116
Xmat <- model.matrix(formula, data)
117
mmats.unadjusted <-
118
c(lapply(random,
119
function(x) model.matrix(formula(x), data = ndat)),
120
list(.Xy = cbind(Xmat, .response = glm.fit$y)))
121
responseIndex <- ncol(mmats.unadjusted$.Xy)
122
obj <- ## creates ssclme structure
123
.Call("ssclme_create", facs,
124
unlist(lapply(mmats.unadjusted, ncol)),
125
PACKAGE = "Matrix")
126
facs = facshuffle(obj, facs)
127
obj = obj[[1]]
128
mmats <- mmats.unadjusted
129
## the next line is to force a copy of mmats, because we are
130
## going to use both mmats and mmats.unadjusted as arguments
131
## in a .Call where one of them will be modified (don't want
132
## the other to be modified as well)
133
mmats[[1]][1,1] <- mmats[[1]][1,1]
134
conv <- FALSE
135
firstIter <- TRUE
136
msMaxIter.orig <- controlvals$msMaxIter
137
138
for (iter in seq(length = controlvals$PQLmaxIt))
139
{
140
mu <- family$linkinv(eta)
141
dmu.deta <- family$mu.eta(eta)
142
## weights (note: weights is already square-rooted)
143
w <- weights * dmu.deta / sqrt(family$variance(mu))
144
## adjusted response (should be comparable to X \beta, not including offset
145
z <- eta - offset + (mmats.unadjusted$.Xy[, responseIndex] - mu) / dmu.deta
146
.Call("nlme_weight_matrix_list",
147
mmats.unadjusted, w, z, mmats, PACKAGE="Matrix")
148
.Call("ssclme_update_mm", obj, facs, mmats, PACKAGE="Matrix")
149
if (firstIter) {
150
.Call("ssclme_initial", obj, PACKAGE="Matrix")
151
if (gVerb) cat(" PQL iterations convergence criterion\n")
152
}
153
.Call("ssclme_EMsteps", obj,
154
controlvals$niterEM,
155
FALSE, #controlvals$REML,
156
controlvals$EMverbose,
157
PACKAGE = "Matrix")
158
LMEoptimize(obj) = controlvals
159
eta[] <- offset + ## FIXME: should the offset be here ?
160
.Call("ssclme_fitted", obj, facs,
161
mmats.unadjusted, TRUE, PACKAGE = "Matrix")
162
crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
163
if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
164
## use this to determine convergence
165
if (crit < controlvals$tolerance) {
166
conv <- TRUE
167
break
168
}
169
etaold[] <- eta
170
171
## Changing number of iterations on second and
172
## subsequent iterations.
173
if (firstIter)
174
{
175
controlvals$niterEM <- 2
176
controlvals$msMaxIter <- 10
177
firstIter <- FALSE
178
}
179
}
180
if (!conv) warning("IRLS iterations for glmm did not converge")
181
controlvals$msMaxIter <- msMaxIter.orig
182
183
## Need to optimize L(theta, beta) using Laplace approximation
184
185
## Things needed for that:
186
##
187
## 1. reduced ssclme object, offset, weighted model matrices
188
## 2. facs, reduced model matrices
189
190
## Of these, those in 2 will be fixed given theta and beta,
191
## and can be thought of arguments to the L(theta, beta)
192
## function. However, the ones in 1 will have the same
193
## structure. So the plan is to pre-allocate them and pass
194
## them in too so they can be used without creating/copying
195
## them more than once
196
197
198
## reduced ssclme
199
200
reducedObj <- .Call("ssclme_collapse", obj, PACKAGE = "Matrix")
201
reducedMmats.unadjusted <- mmats.unadjusted
202
reducedMmats.unadjusted$.Xy <-
203
reducedMmats.unadjusted$.Xy[, responseIndex, drop = FALSE]
204
reducedMmats <- mmats
205
reducedMmats$.Xy <-
206
reducedMmats$.Xy[, responseIndex, drop = FALSE]
207
208
## define function that calculates bhats given theta and beta
209
210
bhat <-
211
function(pars = NULL) # 1:(responseIndex-1) - beta, rest - theta
212
{
213
if (is.null(pars))
214
{
215
off <- drop(mmats.unadjusted$.Xy %*%
216
c(fixef(obj), 0)) + offset
217
}
218
else
219
{
220
.Call("ssclme_coefGets",
221
reducedObj,
222
as.double(pars[responseIndex:length(pars)]),
223
TRUE,
224
PACKAGE = "Matrix")
225
off <- drop(mmats.unadjusted$.Xy %*%
226
c(pars[1:(responseIndex-1)], 0)) + offset
227
}
228
229
niter <- 20
230
conv <- FALSE
231
232
eta <- offset +
233
.Call("ssclme_fitted", obj, facs,
234
mmats.unadjusted, TRUE, PACKAGE = "Matrix")
235
etaold <- eta
236
237
for (iter in seq(length = niter))
238
{
239
mu <- family$linkinv(eta)
240
dmu.deta <- family$mu.eta(eta)
241
w <- weights * dmu.deta / sqrt(family$variance(mu))
242
z <- eta - off + (reducedMmats.unadjusted$.Xy[, 1]
243
- mu) / dmu.deta
244
.Call("nlme_weight_matrix_list",
245
reducedMmats.unadjusted, w, z, reducedMmats,
246
PACKAGE="Matrix")
247
.Call("ssclme_update_mm", reducedObj, facs, reducedMmats,
248
PACKAGE="Matrix")
249
eta[] <- off +
250
.Call("ssclme_fitted", reducedObj, facs,
251
reducedMmats.unadjusted, TRUE,
252
PACKAGE = "Matrix")
253
##cat(paste("bhat Criterion:", max(abs(eta - etaold)) /
254
## (0.1 + max(abs(eta))), "\n"))
255
## use this to determine convergence
256
if (max(abs(eta - etaold)) <
257
(0.1 + max(abs(eta))) * controlvals$tolerance)
258
{
259
conv <- TRUE
260
break
261
}
262
etaold[] <- eta
263
264
}
265
if (!conv) warning("iterations for bhat did not converge")
266
267
## bhat doesn't really need to return anything, we
268
## just want the side-effect of modifying reducedObj
269
## In particular, we are interested in
270
## ranef(reducedObj) and reducedObj@bVar (?). But
271
## the mu-scale response will be useful for log-lik
272
## calculations later, so return them anyway
273
274
invisible(family$linkinv(eta))
275
}
276
277
## function that calculates log likelihood (the thing that
278
## needs to get evaluated at each Gauss-Hermite location)
279
280
## log scale ? worry about details later, get the pieces in place
281
282
## this is for the Laplace approximation only. GH is more
283
## complicated
284
285
devLaplace <- function(pars = NULL)
286
{
287
## FIXME: This actually returns half the deviance.
288
289
## gets correct values of bhat and bvars. As a side
290
## effect, mu now has fitted values
291
mu <- bhat(pars = pars)
292
293
## GLM family log likelihood (upto constant ?)(log scale)
294
## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)
295
296
## Keep everything on (log) likelihood scale
297
298
## log lik from observations given fixed and random effects
299
## get deviance, then multiply by -1/2 (since deviance = -2 log lik)
300
ans <- -sum(family$dev.resids(y = mmats.unadjusted$.Xy[,
301
responseIndex],
302
mu = mu,
303
wt = weights^2))/2
304
305
ranefs <- ranef(reducedObj)
306
# ans <- ans + reducedObj@devComp[2]/2 # log-determinant of Omega
307
Omega <- reducedObj@Omega
308
for (i in seq(along = ranefs))
309
{
310
## contribution for random effects (get it working,
311
## optimize later)
312
## symmetrize RE variance
313
Omega[[i]] <- Omega[[i]] + t(Omega[[i]])
314
diag(Omega[[i]]) <- diag(Omega[[i]]) / 2
315
316
## want log of `const det(Omega) exp(-1/2 b' Omega b )`
317
## i.e., const + log det(Omega) - .5 * (b' Omega b)
318
## FIXME: need to adjust for sigma^2 for appropriate models (easy)
319
## these are all the b'Omega b, summed as they eventually need to be
320
## think of this as sum(rowSums((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]]))
321
322
ranef.loglik.det <- nrow(ranefs[[i]]) *
323
determinant(Omega[[i]], logarithm = TRUE)$modulus/2
324
ranef.loglik.re <- -sum((ranefs[[i]] %*% Omega[[i]]) *
325
ranefs[[i]])/2
326
ranef.loglik <- ranef.loglik.det + ranef.loglik.re
327
328
## Jacobian adjustment
329
log.jacobian <- sum(log(abs(apply(reducedObj@bVar[[i]],
330
3,
331
function(x) sum(diag(x)))
332
)))
333
334
## the constant terms from the r.e. and the final
335
## Laplacian integral cancel out both being:
336
## ranef.loglik.constant <- 0.5 * length(ranefs[[i]]) * log(2 * base::pi)
337
338
ans <- ans + ranef.loglik + log.jacobian
339
}
340
## ans is (up to some constant) log of the Laplacian
341
## approximation of the likelihood. Return it's negative
342
## to be minimized
343
344
# cat("Parameters: ")
345
# print(pars)
346
347
# cat("Value: ")
348
# print(as.double(-ans))
349
350
-ans
351
}
352
353
if (method == "Laplace")
354
{
355
## no analytic gradients or hessians
356
optimRes <-
357
optim(fn = devLaplace,
358
par = c(fixef(obj), coef(obj, unconst = TRUE)),
359
method = "BFGS", hessian = TRUE,
360
control = list(trace = getOption("verbose"),
361
reltol = controlvals$msTol,
362
maxit = controlvals$msMaxIter))
363
if (optimRes$convergence != 0)
364
warning("optim failed to converge")
365
optpars <- optimRes$par
366
Hessian <- optimRes$hessian
367
368
##fixef(obj) <- optimRes$par[seq(length = responseIndex - 1)]
369
if (getOption("verbose")) {
370
cat(paste("optim convergence code",
371
optimRes$convergence, "\n"))
372
cat("Fixed effects:\n")
373
print(fixef(obj))
374
print(optimRes$par[seq(length = responseIndex - 1)])
375
cat("(Unconstrained) variance coefficients:\n")
376
print(coef(obj, unconst = TRUE))
377
coef(obj, unconst = TRUE) <-
378
optimRes$par[responseIndex:length(optimRes$par)]
379
print(coef(obj, unconst = TRUE))
380
}
381
# else if (controlvals$optimizer == "nlm")
382
# {
383
# typsize <- rep(1.0, length(coef(obj, unconst = TRUE)) +
384
# responseIndex - 1)
385
# if (is.null(controlvals$nlmStepMax))
386
# controlvals$nlmStepMax <-
387
# max(100 * sqrt(sum((c(fixef(obj),
388
# coef(obj, unconst = TRUE))/
389
# typsize)^2)), 100)
390
# nlmRes =
391
# nlm(f = devLaplace,
392
# p = c(fixef(obj), coef(obj, unconst = TRUE)),
393
# hessian = TRUE,
394
# print.level = if (getOption("verbose")) 2 else 0,
395
# steptol = controlvals$msTol,
396
# gradtol = controlvals$msTol,
397
# stepmax = controlvals$nlmStepMax,
398
# typsize=typsize,
399
# ## fscale=xval,
400
# iterlim = controlvals$msMaxIter)
401
# if (nlmRes$code > 3)
402
# warning("nlm probably failed to converge")
403
# optpars <- nlmRes$estimate
404
# Hessian <- nlmRes$hessian
405
406
# if (getOption("verbose")) {
407
# cat(paste("nlm convergence code", nlmRes$code, "\n"))
408
# cat("Fixed effects:\n")
409
# print(fixef(obj))
410
# print(nlmRes$estimate[seq(length = responseIndex - 1)])
411
# cat("(Unconstrained) variance coefficients:\n")
412
# print(coef(obj, unconst = TRUE))
413
# coef(obj, unconst = TRUE) <-
414
# nlmRes$estimate[responseIndex:
415
# length(nlmRes$estimate)]
416
# print(coef(obj, unconst = TRUE))
417
# }
418
# }
419
420
## need to calculate likelihood also need to store new
421
## estimates of fixed effects somewhere (probably cannot
422
## update standard errors)
423
} else {
424
optpars <- c(fixef(obj), coef(obj, unconst = TRUE))
425
Hessian <- new("matrix")
426
}
427
428
429
## Before finishing, we need to call devLaplace with the
430
## optimum pars to get the final log likelihood (still need
431
## to make sure it's the actual likelihood and not a
432
## multiple). This would automatically call bhat() and hence
433
## have the 'correct' random effects in reducedObj.
434
435
loglik <- devLaplace(optpars)
436
ff <- optpars[1:(responseIndex-1)]
437
names(ff) <- names(fixef(obj))
438
439
if (!x) mmats <- list()
440
## Things to include in returned object: new ranef
441
## estimates, new parameter estimates (fixed effects and
442
## coefs) (with standard errors from hessian ?) and
443
## (Laplace) approximate log likelihood.
444
445
new("GLMM",
446
family = family,
447
logLik = -loglik,
448
fixef = ff,
449
Hessian = Hessian,
450
method = method,
451
call = match.call(),
452
facs = facs,
453
x = mmats,
454
model = if(model) ndat else data.frame(list()),
455
REML = FALSE,
456
rep = if(method == 'Laplace') reducedObj else obj,
457
fitted = eta,
458
terms = attr(model.frame(formula, data), "terms"),
459
assign = attr(Xmat, "assign")
460
)
461
})
462
463
setMethod("logLik", signature(object = "GLMM"),
464
function(object, ...) {
465
val = object@logLik
466
rr = object@rep
467
nc = rr@nc[-seq(a = rr@Omega)]
468
attr(val, "nall") = attr(val, "nobs") = nc[2]
469
attr(val, "df") = nc[1] + length(coef(rr))
470
attr(val, "REML") = FALSE
471
class(val) <- "logLik"
472
val
473
})
474
475
setMethod("deviance", signature(object = "GLMM"),
476
function(object, ...) -2 * object@logLik)
477
478
setMethod("fixef", signature(object = "GLMM"),
479
function(object) object@fixef)
480
481
setMethod("ranef", signature(object = "GLMM"),
482
function(object, ...)
483
{
484
object = object@rep
485
callGeneric()
486
})
487
488
setMethod("coef", signature(object = "GLMM"),
489
function(object, ...)
490
{
491
object = object@rep
492
callGeneric()
493
})
494
495
setMethod("show", signature(object = "summary.GLMM"),
496
function(object)
497
{
498
rdig <- 5
499
cat("Generalized Linear Mixed Model\n\n")
500
cat("Family:", object@family$family, "family with",
501
object@family$link, "link\n")
502
if (!is.null(object@call$formula)) {
503
cat("Fixed:", deparse(object@call$formula),"\n")
504
}
505
if (!is.null(object@call$data)) {
506
cat("Data:", deparse(object@call$data), "\n")
507
}
508
if (!is.null(object@call$subset)) {
509
cat(" Subset:",
510
deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
511
}
512
llik = object@logLik
513
print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
514
logLik = c(object@logLik), row.names = ""))
515
cat("\n")
516
object@re@useScale = FALSE
517
object@re@showCorrelation = TRUE
518
show(object@re)
519
invisible(object)
520
})
521
522
setMethod("show", signature(object = "GLMM"),
523
function(object)
524
{
525
sumry = summary(object)
526
rdig <- 5
527
cat("Generalized Linear Mixed Model\n\n")
528
cat("Family:", object@family$family, "family with",
529
object@family$link, "link\n")
530
if (!is.null(object@call$formula)) {
531
cat("Fixed:", deparse(object@call$formula),"\n")
532
}
533
if (!is.null(object@call$data)) {
534
cat("Data:", deparse(object@call$data ), "\n")
535
}
536
if (!is.null(object@call$subset)) {
537
cat(" Subset:",
538
deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
539
}
540
cat(paste(" log-likelihood: ", sep = ''), logLik(object), "\n")
541
sumry@re@useScale = FALSE
542
sumry@re@showCorrelation = FALSE
543
saveopt = options(show.signif.stars=FALSE)
544
on.exit(options(saveopt))
545
show(sumry@re)
546
if (object@method != "PQL") { # fix up the fixed effects
547
print(t(object@fixef))
548
}
549
invisible(object)
550
})
551
552
setMethod("summary", signature(object="GLMM"),
553
function(object, ...) {
554
llik <- logLik(object)
555
resd <- residuals(object, type="pearson")
556
if (length(resd) > 5) {
557
resd <- quantile(resd)
558
names(resd) <- c("Min","Q1","Med","Q3","Max")
559
}
560
re = summary(object@rep, REML = FALSE, useScale = FALSE)
561
if (object@method != 'PQL') {
562
hess <- object@Hessian
563
corFixed <- as(as(solve(hess[-nrow(hess),-ncol(hess)]),
564
"pdmatrix"), "corrmatrix")
565
## FIXME: change the name of Hessian to info, to make it
566
## explicit that this is the information matrix
567
fcoef <- object@fixef
568
re@coefficients <- array(c(fcoef, corFixed@stdDev),
569
c(length(fcoef), 2),
570
list(names(fcoef),
571
c("Estimate", "Std. Error")))
572
dimnames(corFixed) <- list(names(fcoef), names(fcoef))
573
re@corFixed <- corFixed
574
}
575
new("summary.GLMM",
576
family = object@family,
577
call = object@call,
578
logLik = llik,
579
re = re,
580
residuals = resd)
581
})
Loggerhead is a web-based interface for Breezy
Version: 2.0.1