-
Notifications
You must be signed in to change notification settings - Fork 55
/
glmmTMB.R
1973 lines (1797 loc) · 87.5 KB
/
glmmTMB.R
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
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
## internal flag for debugging OpenMP behaviour
openmp_debug <- function() {
getOption("glmmTMB_openmp_debug", FALSE)
}
## glmmTMB openmp controller copied from TMB (Windows needs it).
openmp <- function (n = NULL) {
if (openmp_debug() && !is.null(n)) {
cat("setting OpenMP threads to ", n, "\n")
}
## FIXME: redundant with integer-setting within omp_num_threads C++ def in utils.cpp
null_arg <- is.null(n)
if (!null_arg) n <- as.integer(n)
## only want to warn if attempt to set >1 threads in absence
## of OpenMP support ..
if (null_arg || n <= 1) {
w <- options(warn = -1)
on.exit(options(warn = w[["warn"]]))
}
TMB::openmp(n, DLL="glmmTMB")
}
##' Change starting parameters, either by residual method or by user input (start)
##' @inheritParams mkTMBStruc
##' @param formula current formula, containing both fixed & random effects
##' @param ziformula a \emph{one-sided} (i.e., no response variable) formula for zero-inflation combining fixed and random effects: the default \code{~0} specifies no zero-inflation. Specifying \code{~.} sets the zero-inflation formula identical to the right-hand side of \code{formula} (i.e., the conditional effects formula); terms can also be added or subtracted. \strong{When using \code{~.} as the zero-inflation formula in models where the conditional effects formula contains an offset term, the offset term will automatically be dropped}. The zero-inflation model uses a logit link.
##' @param dispformula a \emph{one-sided} formula for dispersion containing only fixed effects: the default \code{~1} specifies the standard dispersion given any family. The argument is ignored for families that do not have a dispersion parameter. For an explanation of the dispersion parameter for each family, see \code{\link{sigma}}. The dispersion model uses a log link. In Gaussian mixed models, \code{dispformula=~0} fixes the residual variance to be 0 (actually a small non-zero value), forcing variance into the random effects. The precise value can be controlled via \code{control=glmmTMBControl(zero_dispval=...)}; the default value is \code{sqrt(.Machine$double.eps)}.
##' @param fr model frame
##' @param yobs observed y
##' @param size number of trials in binomial and betabinomial families
##' @param family family object
##' @param start starting values, expressed as a list with possible components \code{beta}, \code{betazi}, \code{betad} (fixed-effect parameters for conditional, zero-inflation, dispersion models); \code{b}, \code{bzi} (conditional modes for conditional and zero-inflation models); \code{theta}, \code{thetazi} (random-effect parameters, on the standard deviation/Cholesky scale, for conditional and z-i models); \code{psi} (extra family parameters, e.g., shape for Tweedie models).
##' @param sparseX see \code{\link{glmmTMB}}
##' @param start_method Options to initialise the starting values for rr parameters; jitter.sd adds variation to the starting values of latent variables when start = "res".
##' @keywords internal
##' @importFrom stats ppois pbinom rnorm
startParams <- function(parameters,
formula, ziformula, dispformula,
fr,
yobs,
weights,
size = NULL,
Xd = NULL,
XdS = NULL,
family,
condReStruc,
start = NULL,
sparseX = NULL,
start_method = list(method = NULL, jitter.sd = 0)) {
start.met <- start_method$method
jitter.sd <- ifelse(!is.null(start_method$jitter.sd),
start_method$jitter.sd, 0)
## rrValues calculates residuals from the fixed model,
## fits a reduced rank model to obtain starting values for the latent variables and the factor loadings
rrValues <- function(yobs, weights, fr, mu,
family, formula, ziformula, dispformula, condReStruc,
phi = NULL, jitter.sd = 0){
nobs <- length(yobs)
resid <- rep(NA, nobs)
# get either dunn-smyth residuals or
fam <- family$family
res.families <- c("poisson", "nbinom2", "binomial", "gaussian")
if (fam %in% res.families) {
#### Get the dunn smyth residuals
if (fam == "poisson") {
a <- ppois(yobs - 1, mu)
b <- ppois(yobs, mu)
u <- runif(n = nobs, min = a, max = b)
resid <- qnorm(u)
}
if (fam == "nbinom2") {
phi <- phi + 1e-05
a <- pnbinom(yobs - 1, mu = mu, size = phi)
b <- pnbinom(yobs, mu = mu, size = phi)
u <- runif(n = nobs, min = a, max = b)
resid <- qnorm(u)
}
if (fam == "nbinom1") {
phi <- phi + 1e-05
a <- pnbinom(yobs - 1, mu = mu, size = mu/phi)
b <- pnbinom(yobs, mu = mu, size = mu/phi)
u <- runif(n = nobs, min = a, max = b)
resid <- qnorm(u)
}
if (fam == "binomial"){
a <- pbinom(yobs - 1, 1, mu)
b <- pbinom(yobs, 1, mu)
u <- runif(n = nobs, min = a, max = b)
resid <- qnorm(u)
}
if (fam == "gaussian"){
resid <- yobs - mu
}
} else {
resid <- family$dev.resids(y = yobs, mu = mu, wt = weights)
}
resid[is.infinite(resid)] <- 0; resid[is.nan(resid)] <- 0
resid <- as.data.frame(resid)
get_rank <- function(x) {
if (x[["blockCode"]] != .valid_covstruct[["rr"]]) return(0)
p <- x$blockSize
nt <- x$blockNumTheta
rank <- (2*p + 1 - sqrt((2*p+1)^2 - 8*nt))/2
return(rank)
}
rank <- vapply(condReStruc,
get_rank,
FUN.VALUE=numeric(1))
nlv <- rank[rank > 0]
namBlk <- names(nlv)
par.list <- vector("list", length = 3)
names(par.list) <- c("theta", "b", "fact_load")
# Use glmmTMB to get initial starting values for factor loadings and latent variables
fr.res <- cbind(fr, resid)
ranForm <- no_specials(findbars_x(RHSForm(formula)))
nrr <- length(namBlk)
rrTrm <- lapply(1:length(namBlk), function(x) as.character(ranForm[ranForm == namBlk][[x]]))
x <- sapply(1:nrr, function(x) paste(rrTrm[[x]][2], rrTrm[[x]][1], rrTrm[[x]][3]))
resForm <- formula(paste("resid ~ 0 "))
for(i in 1:nrr){
rrForm <- formula(paste("~ rr(", x[i], ",", nlv[i], ")"))
resForm <- addForm(resForm, rrForm)
}
# residual model; assuming gaussian and fixing sd to 1
fit.res <- glmmTMB(resForm, data = fr.res, family = gaussian, start = list(betad = c(log(1))), map = list(betad = factor(c(NA))))
par.list$theta <- fit.res$obj$env$parList(fit.res$fit$par, fit.res$fit$parfull)$theta
par.list$b <- fit.res$obj$env$parList(fit.res$fit$par, fit.res$fit$parfull)$b
# Add jitter to latent variables
par.list$b <- par.list$b + rnorm(length(par.list$b), 0, jitter.sd)
return(par.list)
}
# Fit a fixed model to get the starting values for the fixed parameters, call rrValues to get starting parameters for the rr cov struct (theta and b)
startVals <- function(yobs, weights, fr, Xd, XdS, sparseX,
family, formula, ziformula, dispformula, condReStruc,
parameters, jitter.sd){
start <- parameters #starting parameters
fam <- family$family
### fit a fixed model
fixedform <- formula
RHSForm(fixedform) <- nobars(RHSForm(fixedform))
# FIX ME: Need to add offset?
fit.fixed <- glmmTMB(fixedform, data = fr, family = fam,
ziformula = ziformula, dispformula = dispformula,
weights = weights, sparseX = sparseX)
fixed.pars <- fit.fixed$obj$env$parList(fit.fixed$fit$par, fit.fixed$fit$parfull)
nu <- predict(fit.fixed)
mu <- family$linkinv(nu)
sparseXd <- ifelse(dim(Xd)[1] == 0 && dim(Xd)[2] == 0, 1, 0)
if(length(fixed.pars$betad) != 0){
if(!sparseXd)
phi <- as.matrix(Xd) %*% exp(fixed.pars$betad)
else
phi <- as.vector(XdS %*% exp(fixed.pars$betad))
}
# obtain residuals and get starting values for rr
rrStart <- rrValues(yobs, weights, fr, mu,
family, formula, ziformula, dispformula, condReStruc,
phi, jitter.sd)
start.fixed <- fixed.pars
# Set starting values for fixed parameters from model fit.fixed
fix.names <- !(names(start) %in% c("b", "theta"))
for (i in names(start)[fix.names]) {
if (length(start[[i]]) > 0 & (length(start.fixed[[i]]) == length(start[[i]])))
start[[i]] <- start.fixed[[i]]
}
# Change starting parameters for b and theta for the rr structure
tp <- trrp <- 1 #theta position for full model, and for model with only rr
bp <- brrp <- 1 #b position for full model, and for model with only rr
for (j in seq_along(condReStruc)) {
nt <- condReStruc[[j]]$blockNumTheta
nb <- condReStruc[[j]]$blockReps * condReStruc[[j]]$blockSize
if (condReStruc[[j]]$blockCode == .valid_covstruct[["rr"]]) {
start$b[bp:(bp + nb - 1)] <- rrStart$b[brrp:(brrp + nb - 1)]
start$theta[tp:(tp + nt - 1)] <- rrStart$theta[trrp:(trrp + nt - 1)]
brrp <- brrp + nb; trrp <- trrp + nt
}
bp <- bp + nb
tp <- tp + nt
}
return(start)
}
if(!is.null(start.met)){
start <- startVals(yobs, weights, fr, Xd, XdS, sparseX,
family, formula, ziformula, dispformula, condReStruc,
parameters, jitter.sd)
}
for (p in names(start)) {
if (!(p %in% names(parameters))) {
stop(sprintf("unrecognized vector '%s' in %s",p,sQuote("start")),
call. = FALSE)
}
if ((Lp <- length(parameters[[p]])) != (Ls <- length(start[[p]]))) {
stop(sprintf("parameter vector length mismatch: in %s, length(%s)==%d, should be %d", sQuote("start"), p, Ls, Lp),
call. = FALSE)
}
parameters[[p]] <- start[[p]]
}
return(parameters)
}
##' Extract info from formulas, reTrms, etc., format for TMB
##' @inheritParams glmmTMB
##' @param combForm combined formula
##' @param mf call to model frame
##' @param fr model frame
##' @param yobs observed y
##' @param respCol response column
##' @param weights model weights (for binomial-type models, used as size/number of trials)
##' @param family family object
##' @param se (logical) compute standard error?
##' @param call original \code{glmmTMB} call
##' @param ziPredictCode zero-inflation code
##' @param doPredict flag to enable sds of predictions
##' @param whichPredict which observations in model frame represent predictions
##' @param sparseX see \code{\link{glmmTMB}}
##' @param old_smooths (optional) smooth components from a previous fit: used when constructing a new model structure for prediction
##' from an existing model. A list of smooths for each model component (only cond and zi at present); each smooth has sm and re elements
##' @keywords internal
mkTMBStruc <- function(formula, ziformula, dispformula,
combForm,
mf, fr,
yobs,
respCol,
## no conditional offset argument
## (should be stored in model frame)
weights = NULL,
contrasts,
family,
se=NULL,
call=NULL,
verbose=NULL,
ziPredictCode="corrected",
doPredict=0,
whichPredict=integer(0),
REML=FALSE,
start=NULL,
map=NULL,
sparseX=NULL,
control=glmmTMBControl(),
old_smooths = NULL) {
if (is.null(sparseX)) sparseX <- logical(0)
for (component in c("cond", "zi", "disp")) {
if (!component %in% names(sparseX)) {
sparseX[[component]] <- FALSE
}
}
## make sure the order is right!
sparseX <- sparseX[c("cond","zi","disp")]
## FIXME: (1) should use proper tryCatch below
if (!is(family,"family")) {
## if family specified as list
## if specified as character or function, should have been converted
## to a list with class "family" upstream ...
if (is.list(family)) {
warning("specifying ",sQuote("family")," as a plain list is deprecated")
}
fname <- family$family
args <- family["link"]
ff <- try(do.call(fname,args),silent=TRUE)
if (!inherits(ff,"try-error")) {
family <- ff
} else {
## fallback: add link information to family
## FIXME: is this ever used?
if (is.null(family$linkfun)) {
family <- c(family,make.link(family$link))
}
}
}
## Handle ~0 dispersion for gaussian family.
mapArg <- map
dispformula.orig <- dispformula ## Restore when done
# family$family contains the *name* of the family
if ( usesDispersion(family$family) && (dispformula == ~0) ) {
if (family$family != "gaussian")
stop("~0 dispersion not implemented for ",
sQuote(family$family),
" family")
## FIXME: Depending on the final estimates, we should somehow
## check that this fixed dispersion is small enough.
betad_init <- control$zerodisp_val
dispformula[] <- ~1
mapArg <- c(mapArg,list(betad = factor(NA))) ## Fix betad
} else {
betad_init <- 0
}
## Ignore 'dispformula' argument for non-dispersion families.
if ( ! usesDispersion(family$family) ) {
## if ( dispformula != ~0 &&
## dispformula != ~1 )
## warning(sprintf("dispersion parameter ignored for family %s",
## sQuote(family)))
dispformula[] <- ~0
}
## fixme: may need to modify here, or modify getXReTrms, for smooth-term prediction
condList <- getXReTrms(formula, mf, fr, type="conditional", contrasts=contrasts, sparse=sparseX[["cond"]],
old_smooths = old_smooths$cond)
ziList <- getXReTrms(ziformula, mf, fr, type="zero-inflation", contrasts=contrasts, sparse=sparseX[["zi"]],
old_smooths = old_smooths$zi)
dispList <- getXReTrms(dispformula, mf, fr,
ranOK=FALSE, type="dispersion",
contrasts=contrasts, sparse=sparseX[["disp"]])
condReStruc <- with(condList, getReStruc(reTrms, ss, aa, reXterms, fr))
ziReStruc <- with(ziList, getReStruc(reTrms, ss, aa, reXterms, fr))
grpVar <- with(condList, getGrpVar(reTrms$flist))
nobs <- nrow(fr)
if (is.null(weights)) weights <- rep(1, nobs)
size <- numeric(0)
## binomial family:
## binomial()$initialize was only executed locally
## yobs could be a factor -> treat as binary following glm
## yobs could be cbind(success, failure)
## yobs could be binary
## (yobs, weights) could be (proportions, size)
## On the C++ side 'yobs' must be the number of successes.
if ( binomialType(family$family) ) {
if (is.factor(yobs)) {
## following glm, ‘success’ is interpreted as the factor not
## having the first level (and hence usually of having the
## second level).
yobs <- pmin(as.numeric(yobs)-1,1)
size <- rep(1, nobs)
} else {
if(is.matrix(yobs)) { # yobs=cbind(success, failure)
size <- yobs[,1] + yobs[,2]
yobs <- yobs[,1] #successes
} else {
## previously tested for binary data
## shouldn't need to do this, as weights is a vector of ones
## by default (see NEWS for 1.1.8-9000/1.1.9)
yobs <- weights*yobs
size <- weights
weights <- rep(1.0, nobs)
}
}
}
denseXval <- function(component,lst) if (sparseX[[component]]) matrix(nrow=0,ncol=0) else lst$X
## need a 'dgTMatrix' (double, general, Triplet representation)
sparseXval <- function(component,lst) {
if (sparseX[[component]]) lst$X else nullSparseMatrix()
}
data.tmb <- namedList(
X = denseXval("cond",condList),
XS = sparseXval("cond",condList),
Z = condList$Z,
Xzi = denseXval("zi",ziList),
XziS = sparseXval("zi",ziList),
Zzi = ziList$Z,
Xd = denseXval("disp",dispList),
XdS = sparseXval("disp",dispList),
## Zdisp=dispList$Z,
## use c() on yobs, size to strip attributes such as 'AsIs'
## (which confuse MakeADFun)
yobs = c(yobs),
respCol,
## strip attributes from offset terms
offset = c(condList$offset),
zioffset = c(ziList$offset),
doffset = c(dispList$offset),
weights,
size = c(size),
## information about random effects structure
terms = condReStruc,
termszi = ziReStruc,
family = .valid_family[family$family],
link = .valid_link[family$link],
ziPredictCode = .valid_zipredictcode[ziPredictCode],
doPredict = doPredict,
whichPredict = whichPredict
)
# function to set value for dorr
rrVal <- function(lst) if(any(lst$ss == "rr")) 1 else 0
dorr = rrVal(condList)
getVal <- function(obj, component)
vapply(obj, function(x) x[[component]], numeric(1))
## safer initialization for link functions that might give
## illegal predictions for certain families
## (sqrt() behaves weirdly for beta=0
## [because inverse-link is symmetric around 0?]
beta_init <- if (family$link %in% c("identity","inverse","sqrt")) 1 else 0
rr0 <- function(n) {
if (is.null(n)) numeric(0) else rep(0, n)
}
## Extra family specific parameters
## FIXME: switch/rewrite to be less ugly?
psiLength <- if (family$family %in% c("t", "tweedie"))
{ 1 } else if (family$family == "ordbeta") { 2 } else { 0 }
psi_init <- if (family$family == "ordbeta") c(-1, 1) else rr0(psiLength)
# theta is 0, except if dorr, theta is 1
t01 <- function(dorr, condReStruc){
theta <- rr0(sum(getVal(condReStruc,"blockNumTheta")))
if(dorr){
nt <- 1
blockNumTheta <- getVal(condReStruc,"blockNumTheta")
blockCode <- getVal(condReStruc, "blockCode")
for (i in 1:length(blockCode)) {
if(names(.valid_covstruct)[match(blockCode[i], .valid_covstruct)]=="rr") {
theta[nt:(nt + blockNumTheta[i] - 1)] <- rep(1, blockNumTheta[i])
}
nt <- nt + blockNumTheta[i]
}
}
theta
}
parameters <- with(data.tmb,
list(
beta = rep(beta_init, max(ncol(X), ncol(XS))),
betazi = rr0(max(ncol(Xzi),ncol(XziS))),
b = rep(beta_init, ncol(Z)),
bzi = rr0(ncol(Zzi)),
betad = rep(betad_init, max(ncol(Xd),ncol(XdS))),
theta = t01(dorr, condReStruc),
thetazi = rr0(sum(getVal(ziReStruc, "blockNumTheta"))),
psi = psi_init
))
if(!is.null(start) || !is.null(control$start_method$method)){
parameters <- startParams(parameters,
formula, ziformula, dispformula,
fr,
yobs = data.tmb$yobs,
weights = data.tmb$weights,
size = data.tmb$size,
Xd = data.tmb$Xd,
XdS = data.tmb$XdS,
family,
condReStruc,
start = start,
sparseX = sparseX,
start_method = control$start_method)
}
randomArg <- c(if(ncol(data.tmb$Z) > 0) "b",
if(ncol(data.tmb$Zzi) > 0) "bzi")
## REML
if (REML) randomArg <- c(randomArg, "beta")
dispformula <- dispformula.orig ## May have changed - restore
return(namedList(data.tmb, parameters, mapArg, randomArg, grpVar,
condList, ziList, dispList, condReStruc, ziReStruc,
family, contrasts, respCol,
allForm=namedList(combForm,formula,ziformula,dispformula),
fr, se, call, verbose, REML, map, sparseX))
}
##' Create X and random effect terms from formula
##' @param formula current formula, containing both fixed & random effects
##' @param mf matched call
##' @param fr full model frame
##' @param ranOK random effects allowed here?
##' @param type label for model type
##' @param contrasts a list of contrasts (see ?glmmTMB)
##' @param sparse (logical) return sparse model matrix?
##' @param old_smooths smooth information from a prior model fit (for prediction)
##' @return a list composed of
##' \item{X}{design matrix for fixed effects}
##' \item{Z}{design matrix for random effects}
##' \item{reTrms}{output from \code{\link{mkReTrms}} from \pkg{lme4}, possibly augmented with information about \code{mgcv}-style smooth terms}
##' \item{ss}{splitform of the formula}
##' \item{aa}{additional arguments, used to obtain rank}
##' \item{terms}{terms for the fixed effects}
##' \item{offset}{offset vector, or vector of zeros if offset not specified}
##' \item{reXterms}{terms for the model matrix in each RE term}
##'
##' @importFrom stats model.matrix contrasts
##' @importFrom methods new
##' @importFrom lme4 findbars nobars
##' @importFrom mgcv smoothCon smooth2random s PredictMat
getXReTrms <- function(formula, mf, fr, ranOK=TRUE, type="",
contrasts, sparse=FALSE, old_smooths = NULL) {
has_re <- !is.null(findbars_x(formula))
has_smooths <- anySpecial(formula, specials = "s")
## fixed-effects model matrix X -
## remove random effect parts from formula:
fixedform <- formula
RHSForm(fixedform) <- nobars(RHSForm(fixedform))
terms <- NULL ## make sure it's empty in case we don't set it
nobs <- nrow(fr)
## check for empty fixed form
## need to ignore environments when checking!
## ignore.environment= arg only works with closures
idfun <- function(x,y) {
environment(x) <- emptyenv()
environment(y) <- emptyenv()
return(identical(x,y))
}
if (idfun(RHSForm(fixedform, as.form=TRUE), ~ 0) ||
idfun(RHSForm(fixedform, as.form=TRUE), ~ -1)) {
X <- matrix(ncol=0, nrow=nobs)
offset <- rep(0,nobs)
} else {
## check for mgcv-style smooth terms, adjust accordingly ...
if (has_smooths) {
## extract s() terms
smooth_terms <- findbars_x(fixedform, default.special = NULL, target = "s")
## *remove* s() terms from fixed formula
fixedform <- noSpecials(fixedform, specials = "s")
## FIXME: could be fragile about eval environments (how
## far up do we have to go with eval.parent? Or do we
## use the environment of the formula?
if (!is.null(old_smooths)) {
smooth_terms2 <- lapply(old_smooths[lengths(old_smooths)>0],
function(s) {
if (is.null(s)) return(NULL)
X <- PredictMat(s$sm, fr) ## get prediction matrix for new data
## transform to r.e. parameterization
if (!is.null(s$re$trans.U)) X <- X%*%s$re$trans.U
X <- t(t(X)*s$re$trans.D)
## re-order columns according to random effect re-ordering...
X[,s$re$rind] <- X[,s$re$pen.ind!=0]
## re-order penalization index in same way
pen.ind <- s$re$pen.ind; s$pen.ind[s$re$rind] <- pen.ind[pen.ind>0]
## start return object...
s_new <- list(re = list(rand=list(), Xf=X[,which(s$re$pen.ind==0),drop=FALSE]))
for (i in 1:length(s$re$rand)) { ## loop over random effect matrices
s_new$re$rand[[i]] <- X[, which(pen.ind==i), drop=FALSE]
attr(s_new$re$rand[[i]], "s.label") <- attr(s$re$rand[[i]], "s.label")
}
names(s_new$re$rand) <- names(s$re$rand)
return(s_new)
})
} else {
smooth_terms2 <- lapply(smooth_terms,
function(tt) {
## ‘smoothCon’ returns a list of smooths because factor ‘by’
## variables result in multiple copies of a smooth, each multiplied
## by the dummy variable associated with one factor level.
sm <- eval(bquote(mgcv::smoothCon(.(tt),
## don't want intercept etc ...
absorb.cons = TRUE,
data=mf$data)))
if (length(sm)>1) stop("can't handle 'by' arguments in smooths yet")
sm <- sm[[1]]
list(sm = sm,
## FIXME: will vnames ("a vector of names
## to avoid as dummy variable names in the
## random effects form") ever be non-empty?
re = mgcv::smooth2random(sm, vnames = "", type = 2))
})
} ## create (new) smooth terms
} ## has_smooths
tt <- terms(fixedform)
pv <- attr(mf$formula,"predvars")
attr(tt, "predvars") <- fix_predvars(pv,tt)
mf$formula <- tt
terms_fixed <- terms(eval(mf,envir=environment(fixedform)))
terms <- list(fixed=terms(terms_fixed))
if (!sparse) {
X <- model.matrix(drop.special(fixedform), fr, contrasts)
} else {
X <- Matrix::sparse.model.matrix(drop.special(fixedform), fr, contrasts)
## FIXME? ?sparse.model.matrix recommends MatrixModels::model.Matrix(*,sparse=TRUE)
## (but we may not need it, and would add another dependency etc.)
}
if (has_smooths) {
if (sparse) warning("smooth terms may not be compatible with sparse X matrices")
for (s in smooth_terms2) {
cnm <- colnames(X)
snm <- attr(s$re$rand$Xr, "s.label")
X <- cbind(X, s$re$Xf)
colnames(X) <- c(cnm, paste0(snm, seq.int(ncol(s$re$Xf))))
}
}
## will be 0-column matrix if fixed formula is empty
offset <- rep(0,nobs)
if (inForm(fixedform,quote(offset))) {
## hate to match offset terms with model frame names
## via deparse, but since that what was presumably done
## internally to get the model frame names in the first place ...
for (o in extractForm(fixedform,quote(offset))) {
offset_nm <- deparse1(o)
## don't think this will happen, but ...
if (length(offset_nm)>1) {
stop("trouble reconstructing offset name")
}
offset <- offset + fr[[offset_nm]]
}
}
}
## ran-effects model frame (for predvars)
## important to COPY formula (and its environment)?
ranform <- formula
if (!has_re && !has_smooths) {
reTrms <- reXterms <- NULL
Z <- new("dgCMatrix",Dim=c(as.integer(nobs),0L)) ## matrix(0, ncol=0, nrow=nobs)
aa <- integer(0) #added for rr to get rank
ss <- integer(0)
} else {
## FIXME: check whether predvars are carried along correctly in terms
if (!ranOK) stop("no random effects allowed in ", type, " term")
RHSForm(ranform) <- subbars(RHSForm(reOnly(formula)))
if (has_re) {
mf$formula <- ranform
reTrms <- mkReTrms(no_specials(findbars_x(formula)),
fr, reorder.terms=FALSE)
} else {
## dummy elements
reTrms <- list(Ztlist = list(), flist = list(), cnms = list(),
theta = list())
}
## formula <- Reaction ~ s(Days) + (1|Subject)
ss <- splitForm(formula)
## contains: c("Zt", "theta", "Lind", "Gp", "lower", "Lambdat", "flist", "cnms", "Ztlist", "nl")
## we only need "Zt", "flist", "Gp", "cnms", "Ztlist" (I think)
## "theta", "Lind", "lower", "Lambdat", "nl" (?) are lme4-specific
## need to fill in smooth terms in **correct locations**
## do we need to reconstitute
## post-process mkReTrms to add smooths (incorporate in mkReTrms?)
if (has_smooths) {
ns <- length(ss$reTrmClasses)
augReTrms <- list(Ztlist = vector("list", ns),
flist = vector("list", ns),
cnms = vector("list", ns),
smooth_info = vector("list", ns))
barpos <- which(ss$reTrmClasses != "s")
nonbarpos <- which(ss$reTrmClasses == "s")
for (p in c("Ztlist", "flist", "cnms")) {
augReTrms[[p]][barpos] <- reTrms[[p]]
names(augReTrms[[p]])[barpos] <- names(reTrms[[p]])
}
## one theta value for smooths
## FIXME: do we actually use theta values?
## FIXME: can we set
augReTrms$theta <- rep(0, ## default value
sum(lengths(reTrms$theta)) +
length(nonbarpos))
augReTrms$theta[barpos] <- reTrms$theta
## only need one 'dummy' factor for all the smooth terms
ff <- factor(rep(1, nobs))
augReTrms$flist <- c(reTrms$flist, list(dummy = ff))
avec <- rep(NA_integer_, ns)
avec[barpos] <- attr(reTrms$flist, "assign")
## mkReTrms returns more than we need (some is for lme4)
## ... which bits are actually used hereafter?
avec[nonbarpos] <- length(augReTrms$flist)
attr(augReTrms$flist, "assign") <- avec
for (i in seq_along(smooth_terms2)) {
s <- smooth_terms2[[i]]
pos <- nonbarpos[i]
Zt <- as(t(s$re$rand$Xr), "dgCMatrix")
npar <- nrow(Zt)
augReTrms$Ztlist[[pos]] <- Zt
nm <- attr(s$re$rand$Xr, "s.label")
names(augReTrms$Ztlist)[pos] <- nm
## cnms
augReTrms$cnms[[pos]] <- paste0("dummy", seq(npar))
names(augReTrms$cnms)[pos] <- "dummy"
}
## store smooth info in relevant spots
for (i in seq_along(nonbarpos)) {
augReTrms$smooth_info[[i]] <- smooth_terms2[[nonbarpos[i]]]
}
## reconstitute other pieces
augReTrms$Zt <- do.call(rbind, augReTrms$Zt)
augReTrms$Gp <- cumsum(c(0, vapply(augReTrms$Ztlist, nrow, 0L)))
##
reTrms <- augReTrms
}
ss$reTrmClasses[ss$reTrmClasses == "s"] <- "homdiag"
# FIX ME: migrate this (or something like it) down to reTrms,
## allow for more different covstruct types that have additional arguments
## e.g. phylo(.,tree); fixed(.,Sigma)
# FIX ME: use NA rather than 0 as a placeholder in aa?
## FIXME: make sure that eval() happens in the right environment/
## document potential issues
get_num <- function(v) {
if (length(v) == 1) return(NA_real_)
payload <- v[[2]]
res <- tryCatch(eval(payload, envir = environment(formula)),
error = function(e)
stop("can't evaluate reduced-rank dimension ",
sQuote(deparse(payload)),
.call = FALSE))
if (is.na(suppressWarnings(as.numeric(res)))) {
stop("non-numeric value for reduced-rank dimension",
call. = FALSE)
}
return(res)
}
aa <- ifelse(ss$reTrmClass=="rr",
vapply(ss$reTrmAddArgs,
get_num,
FUN.VALUE=numeric(1)),
0)
## terms for the model matrix in each RE term
## this is imperfect: it should really be done in mkReTrms/mkBlist,
## where we are generating these terms anyway on the way
## to constructing Z, but that's in lme4 so we can't change it
## unless absolutely necessary
termsfun <- function(x) {
## this is a little magic: copying lme4:::mkBlist approach
ff <- eval(substitute( ~ foo, list(foo = x[[2]]))) ## make formula from LHS
tt <- try(terms(ff, data=fr), silent=TRUE) ## construct terms
if (inherits(tt,"try-error")) {
stop(
sprintf("can't evaluate RE term %s: simplify?",
sQuote(deparse(ff)))
)
}
tt
}
## HACK: should duplicate 'homdiag' definition, keep it as 's' (or call it 'mgcv_smooth")
## so we can recognize it.
## Here, we're using the fact that the ...AddArgs stuff is still in an unevaluated form
reXterms <- Map(function(f, a) {
if (identical(head(a), as.symbol('s'))) NA else termsfun(f)
}, ss$reTrmFormulas, ss$reTrmAddArgs)
ss <- unlist(ss$reTrmClasses)
Z <- t(reTrms$Zt) ## still sparse ...
}
## if(is.null(rankX.chk <- control[["check.rankX"]]))
## rankX.chk <- eval(formals(lmerControl)[["check.rankX"]])[[1]]
## X <- chkRank.drop.cols(X, kind=rankX.chk, tol = 1e-7)
## if(is.null(scaleX.chk <- control[["check.scaleX"]]))
## scaleX.chk <- eval(formals(lmerControl)[["check.scaleX"]])[[1]]
## X <- checkScaleX(X, kind=scaleX.chk)
## list(fr = fr, X = X, reTrms = reTrms, family = family, formula = formula,
## wmsgs = c(Nlev = wmsgNlev, Zdims = wmsgZdims, Zrank = wmsgZrank))
namedList(X, Z, reTrms, ss, aa, terms, offset, reXterms)
}
##' Extract grouping variables for random effect terms from a factor list
##' @title Get Grouping Variable
##' @param x "flist" object; a data frame of factors including an \code{assign} attribute
##' matching columns to random effect terms
##' @return character vector of grouping variables
##' @keywords internal
##' @examples
##' data(cbpp,package="lme4")
##' cbpp$obs <- factor(seq(nrow(cbpp)))
##' rt <- lme4::glFormula(cbind(size,incidence-size)~(1|herd)+(1|obs),
##' data=cbpp,family=binomial)$reTrms
##' getGrpVar(rt$flist)
##' @export
getGrpVar <- function(x)
{
assign <- attr(x,"assign")
names(x)[assign]
}
##' Calculate random effect structure
##' Calculates number of random effects, number of parameters,
##' block size and number of blocks. Mostly for internal use.
##' @param reTrms random-effects terms list
##' @param ss a character string indicating a valid covariance structure.
##' Must be one of \code{names(glmmTMB:::.valid_covstruct)};
##' default is to use an unstructured variance-covariance
##' matrix (\code{"us"}) for all blocks).
##' @param reXterms terms objects corresponding to each RE term
##' @param fr model frame
##' @param aa additional arguments (i.e. rank)
##' @return a list
##' \item{blockNumTheta}{number of variance covariance parameters per term}
##' \item{blockSize}{size (dimension) of one block}
##' \item{blockReps}{number of times the blocks are repeated (levels)}
##' \item{covCode}{structure code}
##' @examples
##' data(sleepstudy, package="lme4")
##' rt <- lme4::lFormula(Reaction~Days+(1|Subject)+(0+Days|Subject),
##' sleepstudy)$reTrms
##' rt2 <- lme4::lFormula(Reaction~Days+(Days|Subject),
##' sleepstudy)$reTrms
##' getReStruc(rt)
##' @importFrom stats setNames dist .getXlevels
##' @export
getReStruc <- function(reTrms, ss=NULL, aa=NULL, reXterms=NULL, fr=NULL) {
## information from ReTrms is contained in cnms, flist elements
## cnms: list of column-name vectors per term
## flist: data frame of grouping variables (factors)
## 'assign' attribute gives match between RE terms and factors
if (is.null(reTrms)) {
list()
} else {
## Get info on sizes of RE components
assign <- attr(reTrms$flist,"assign")
nreps <- vapply(assign,
function(i) length(levels(reTrms$flist[[i]])),
0)
blksize <- diff(reTrms$Gp) / nreps
## figure out number of parameters from block size + structure type
if (is.null(ss)) {
ss <- rep("us",length(blksize))
}
if ( any(is.na(aa[ss=="rr"]))) {
aa0 <- which(is.na(aa) & ss=="rr")
aa[aa0] <- 2 #set default rank to 2 if it's not specified
}
if ( is.null(aa)) {
aa <- rep(0,length(blksize)) #set rank to 0
}
blkrank <- aa
covCode <- .valid_covstruct[ss]
parFun <- function(struc, blksize, blkrank) {
switch(as.character(struc),
"diag" = blksize, # (heterogenous) diag
"us" = blksize * (blksize+1) / 2,
"cs" = blksize + 1,
"ar1" = 2,
"ou" = 2,
"exp" = 2,
"gau" = 2,
"mat" = 3,
"toep" = 2 * blksize - 1,
"rr" = blksize * blkrank - (blkrank - 1) * blkrank / 2, #rr
"homdiag" = 1 ## (homogeneous) diag
)
}
blockNumTheta <- mapply(parFun, ss, blksize, blkrank, SIMPLIFY=FALSE)
ans <- list()
for (i in seq_along(ss)) {
tmp <- list(blockReps = nreps[i],
blockSize = blksize[i],
blockNumTheta = blockNumTheta[[i]],
blockCode = covCode[i]
)
if(ss[i] == "ar1") {
## FIXME: Keep this warning ?
if (any(reTrms$cnms[[i]][1] == "(Intercept)") )
warning("AR1 not meaningful with intercept")
if (length(.getXlevels(reXterms[[i]],fr))!=1) {
stop("ar1() expects a single, factor variable as the time component")
}
} else if(ss[i] == "ou"){
times <- parseNumLevels(reTrms$cnms[[i]])
if (ncol(times) != 1)
stop("'ou' structure is for 1D coordinates only.")
if (is.unsorted(times, strictly=TRUE))
stop("'ou' is for strictly sorted times only.")
tmp$times <- drop(times)
} else if(ss[i] %in% c("exp", "gau", "mat")){
coords <- parseNumLevels(reTrms$cnms[[i]])
tmp$dist <- as.matrix( dist(coords) )
}
ans[[i]] <- tmp
}
setNames(ans, names(reTrms$Ztlist))
}
}
.noDispersionFamilies <- c("binomial", "poisson", "truncated_poisson")
## BMB: why not just sigma(x)!=1.0 ... ? (redundant with sigma.glmmTMB)
usesDispersion <- function(x) {
is.na(match(x, .noDispersionFamilies))
## !x %in% .noDispersionFamilies
}
.classicDispersionFamilies <- c("gaussian","Gamma","t")
## select only desired pieces from results of getXReTrms
stripReTrms <- function(xrt, whichReTrms = c("cnms","flist","smooth_info"), which="terms") {
whichReTrms <- intersect(whichReTrms, names(xrt$reTrms)) ## not all models will have smooth_info
c(xrt$reTrms[whichReTrms],setNames(xrt[which],which))
}
#.okWeightFamilies <- c("binomial", "betabinomial")
okWeights <- function(x) {
TRUE
#!is.na(match(x, .okWeightFamilies))
## x %in% .okWeightFamilies
}
## Families for which binomial()$initialize is used
.binomialFamilies <- c("binomial", "betabinomial")
binomialType <- function(x) {
!is.na(match(x, .binomialFamilies))
}
##' Fit Models with TMB
##'
##' Fit a generalized linear mixed model (GLMM) using Template Model Builder (TMB).
##' @param formula combined fixed and random effects formula, following lme4 syntax.
##' @param data data frame (tibbles are OK) containing model variables. Not required, but strongly recommended; if \code{data} is not specified, downstream methods such as prediction with new data (\code{predict(fitted_model, newdata = ...)}) will fail. If it is necessary to call \code{glmmTMB} with model variables taken from the environment rather than from a data frame, specifying \code{data=NULL} will suppress the warning message.
##' @param family a family function, a character string naming a family function, or the result of a call to a family function (variance/link function) information. See \code{\link{family}} for a generic discussion of families or \code{\link{family_glmmTMB}} for details of \code{glmmTMB}-specific families.
##' @param ziformula a \emph{one-sided} (i.e., no response variable) formula for zero-inflation combining fixed and random effects: the default \code{~0} specifies no zero-inflation. Specifying \code{~.} sets the zero-inflation formula identical to the right-hand side of \code{formula} (i.e., the conditional effects formula); terms can also be added or subtracted. \strong{When using \code{~.} as the zero-inflation formula in models where the conditional effects formula contains an offset term, the offset term will automatically be dropped}. The zero-inflation model uses a logit link.
##' @param dispformula a \emph{one-sided} formula for dispersion containing only fixed effects: the default \code{~1} specifies the standard dispersion given any family. The argument is ignored for families that do not have a dispersion parameter. For an explanation of the dispersion parameter for each family, see \code{\link{sigma}}. The dispersion model uses a log link. In Gaussian mixed models, \code{dispformula=~0} fixes the residual variance to be 0 (actually a small non-zero value), forcing variance into the random effects. The precise value can be controlled via \code{control=glmmTMBControl(zero_dispval=...)}; the default value is \code{sqrt(.Machine$double.eps)}.
##' @param weights weights, as in \code{glm}. Not automatically scaled to have sum 1.
##' @param offset offset for conditional model (only).
##' @param contrasts an optional list, e.g., \code{list(fac1="contr.sum")}. See the \code{contrasts.arg} of \code{\link{model.matrix.default}}.
##' @param na.action a function that specifies how to handle observations
##' containing \code{NA}s. The default action (\code{na.omit},
##' inherited from the 'factory fresh' value of
##' \code{getOption("na.action")}) strips any observations with any
##' missing values in any variables. Using \code{na.action = na.exclude}
##' will similarly drop observations with missing values while fitting the model,
##' but will fill in \code{NA} values for the predicted and residual
##' values for cases that were excluded during the fitting process
##' because of missingness.
##' @param se whether to return standard errors.
##' @param verbose whether progress indication should be printed to the console.
##' @param doFit whether to fit the full model, or (if FALSE) return the preprocessed data and parameter objects, without fitting the model.
##' @param control control parameters, see \code{\link{glmmTMBControl}}.
##' @param REML whether to use REML estimation rather than maximum likelihood.
##' @param start starting values, expressed as a list with possible components \code{beta}, \code{betazi}, \code{betad} (fixed-effect parameters for conditional, zero-inflation, dispersion models); \code{b}, \code{bzi} (conditional modes for conditional and zero-inflation models); \code{theta}, \code{thetazi} (random-effect parameters, on the standard deviation/Cholesky scale, for conditional and z-i models); \code{psi} (extra family parameters, e.g., shape for Tweedie models).
##' @param map a list specifying which parameter values should be fixed to a constant value rather than estimated. \code{map} should be a named list containing factors corresponding to a subset of the internal parameter names (see \code{start} parameter). Distinct factor values are fitted as separate parameter values, \code{NA} values are held fixed: e.g., \code{map=list(beta=factor(c(1,2,3,NA)))} would fit the first three fixed-effect parameters of the conditional model and fix the fourth parameter to its starting value. In general, users will probably want to use \code{start} to specify non-default starting values for fixed parameters. See \code{\link[TMB]{MakeADFun}} for more details.
##' @param sparseX a named logical vector containing (possibly) elements named "cond", "zi", "disp" to indicate whether fixed-effect model matrices for particular model components should be generated as sparse matrices, e.g. \code{c(cond=TRUE)}. Default is all \code{FALSE}
##' @importFrom stats gaussian binomial poisson nlminb as.formula terms model.weights
##' @importFrom lme4 subbars findbars mkReTrms nobars
##' @importFrom Matrix t
##' @importFrom TMB MakeADFun sdreport
##' @details
##' \itemize{
##' \item Binomial models with more than one trial (i.e., not binary/Bernoulli) can either be specified in the form \code{prob ~ ..., weights = N}, or in the more typical two-column matrix \code{cbind(successes,failures)~...} form.
##' \item Behavior of \code{REML=TRUE} for Gaussian responses matches \code{lme4::lmer}. It may also be useful in some cases with non-Gaussian responses (Millar 2011). Simulations should be done first to verify.
##' \item Because the \code{\link{df.residual}} method for \code{glmmTMB} currently counts the dispersion parameter, users should multiply this value by \code{sqrt(nobs(fit) / (1+df.residual(fit)))} when comparing with \code{lm}.
##' \item Although models can be fitted without specifying a \code{data} argument, its use is strongly recommended; drawing model components from the global environment, or using \code{df$var} notation within model formulae, can lead to confusing (and sometimes hard-to-detect) errors.
##' \item By default, vector-valued random effects are fitted with unstructured (general symmetric positive definite) variance-covariance matrices. Structured variance-covariance matrices can be specified in the form \code{struc(terms|group)}, where \code{struc} is one of
##' \itemize{
##' \item \code{diag} (diagonal, heterogeneous variance)
##' \item \code{ar1} (autoregressive order-1, homogeneous variance)
##' \item \code{cs} (compound symmetric, heterogeneous variance)
##' \item \code{ou} (* Ornstein-Uhlenbeck, homogeneous variance)
##' \item \code{exp} (* exponential autocorrelation)
##' \item \code{gau} (* Gaussian autocorrelation)
##' \item \code{mat} (* Matérn process correlation)
##' \item \code{toep} (* Toeplitz)
##' \item \code{rr} (reduced rank/factor-analytic model)
##' \item \code{homdiag} (diagonal, homogeneous variance)
##' }
##' Structures marked with * are experimental/untested. See \code{vignette("covstruct", package = "glmmTMB")} for more information.
##' \item For backward compatibility, the \code{family} argument can also be specified as a list comprising the name of the distribution and the link function (e.g. \code{list(family="binomial", link="logit")}). However, \strong{this alternative is now deprecated}; it produces a warning and will be removed at some point in the future. Furthermore, certain capabilities such as Pearson residuals or predictions on the data scale will only be possible if components such as \code{variance} and \code{linkfun} are present, see \code{\link{family}}.
##' \item Smooths taken from the \code{mgcv} package can be included in \code{glmmTMB} formulas using \code{s}; these terms will appear as additional components in both the fixed and the random-effects terms. This functionality is \emph{experimental} for now. We recommend using \code{REML=TRUE}. See \code{\link[mgcv]{s}} for details of specifying smooths (and \code{\link[mgcv]{smooth2random}} and the appendix of Wood (2004) for technical details).