/
ModelFit.R
1137 lines (1014 loc) · 44.8 KB
/
ModelFit.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
# @file ModelFit.R
#
# Copyright 2016 Observational Health Data Sciences and Informatics
#
# This file is part of cyclops
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#' @title Fit a Cyclops model
#'
#' @description
#' \code{fitCyclopsModel} fits a Cyclops model data object
#'
#' @details
#' This function performs numerical optimization to fit a Cyclops model data object.
#'
#' @param cyclopsData A Cyclops data object
#' @template prior
#' @param control A \code{"cyclopsControl"} object constructed by \code{\link{createControl}}
#' @param weights Vector of 0/1 weights for each data row
#' @param forceNewObject Logical, forces the construction of a new Cyclops model fit object
#' @param returnEstimates Logical, return regression coefficient estimates in Cyclops model fit object
#' @param startingCoefficients Vector of starting values for optimization
#' @param fixedCoefficients Vector of booleans indicating if coefficient should be fix
#' @param warnings Logical, report regularization warnings
#' @param computeDevice String: Name of compute device to employ; defaults to \code{"native"} C++ on CPU
#'
#' @return
#' A list that contains a Cyclops model fit object pointer and an operation duration
#'
#' @references
#' Suchard MA, Simpson SE, Zorych I, Ryan P, Madigan D.
#' Massive parallelization of serial inference algorithms for complex generalized linear models.
#' ACM Transactions on Modeling and Computer Simulation, 23, 10, 2013.
#'
#' Simpson SE, Madigan D, Zorych I, Schuemie M, Ryan PB, Suchard MA.
#' Multiple self-controlled case series for large-scale longitudinal observational databases.
#' Biometrics, 69, 893-902, 2013.
#'
#' Mittal S, Madigan D, Burd RS, Suchard MA.
#' High-dimensional, massive sample-size Cox proportional hazards regression for survival analysis.
#' Biostatistics, 15, 207-221, 2014.
#'
#' @examples
#' ## Dobson (1990) Page 93: Randomized Controlled Trial :
#' counts <- c(18,17,15,20,10,20,25,13,12)
#' outcome <- gl(3,1,9)
#' treatment <- gl(3,3)
#' cyclopsData <- createCyclopsData(counts ~ outcome + treatment, modelType = "pr")
#' cyclopsFit <- fitCyclopsModel(cyclopsData, prior = createPrior("none"))
#' coef(cyclopsFit)
#' confint(cyclopsFit, c("outcome2","treatment3"))
#' predict(cyclopsFit)
#'
#' @export
fitCyclopsModel <- function(cyclopsData,
prior = createPrior("none"),
control = createControl(),
weights = NULL,
forceNewObject = FALSE,
returnEstimates = TRUE,
startingCoefficients = NULL,
fixedCoefficients = NULL,
warnings = TRUE,
computeDevice = "native") {
# Delegate to control$setHook if exists
if (!is.null(control$setHook)) {
return(control$setHook(cyclopsData, prior, control,
weights, forceNewObject, returnEstimates,
startingCoefficients, fixedCoefficients))
}
# Delegate to prior$fitHook if exists
if (!is.null(prior$fitHook)) {
return(prior$fitHook(cyclopsData, prior, control,
weights, forceNewObject, returnEstimates,
startingCoefficients, fixedCoefficients))
}
cl <- match.call()
# Check conditions
.checkData(cyclopsData)
if (getNumberOfRows(cyclopsData) < 1 ||
getNumberOfStrata(cyclopsData) < 1 ||
getNumberOfCovariates(cyclopsData) < 1) {
stop("Data are incompletely loaded")
}
.checkInterface(cyclopsData, computeDevice = computeDevice, forceNewObject = forceNewObject)
# Set up prior
stopifnot(inherits(prior, "cyclopsPrior"))
if (!is.null(prior$setHook)) {
prior$setHook(cyclopsData, warnings) # Call-back
} else {
prior$exclude <- .checkCovariates(cyclopsData, prior$exclude)
if (!is.null(prior$neighborhood)) {
prior$neighborhood <- lapply(prior$neighborhood,
function(element) {
list(.checkCovariates(cyclopsData, element[[1]]),
.checkCovariates(cyclopsData, element[[2]]))
})
}
if (prior$priorType[1] != "none" &&
is.null(prior$graph) && # TODO Ignore hierarchical models for now
.cyclopsGetHasIntercept(cyclopsData) &&
!prior$forceIntercept) {
interceptId <- bit64::as.integer64(.cyclopsGetInterceptLabel(cyclopsData))
warn <- FALSE
if (is.null(prior$exclude)) {
prior$exclude <- c(interceptId)
warn <- TRUE
} else {
if (!(interceptId %in% prior$exclude)) {
prior$exclude <- c(interceptId, prior$exclude)
warn <- TRUE
}
}
if (warn && warnings) {
warning("Excluding intercept from regularization")
}
}
if (is.null(prior$graph)) {
graph <- NULL
} else {
graph <- .makeHierarchyGraph(cyclopsData, prior$graph)
if (length(prior$priorType) != length(prior$variance)) {
stop("Prior types and variances have a dimensionality mismatch")
}
if (any(prior$priorType != "normal")) {
stop("Only normal-normal hierarchies are currently supported")
}
}
if (is.null(prior$neighborhood)) {
neighborhood <- NULL
} else {
neighborhood <- prior$neighborhood
if (length(prior$priorType) != length(prior$variance)) {
stop("Prior types and variances have a dimensionality mismatch");
}
if (any(prior$priorType != "laplace")) {
stop("Only Laplace-Laplace fused neighborhoods are currently supported")
}
}
if (is.null(graph) && is.null(neighborhood) && length(prior$priorType) > 1) {
if (length(prior$priorType) != getNumberOfCovariates(cyclopsData)) {
stop("Length of priors must equal the number of covariates")
}
}
if (any(prior$priorType == "jeffreys")) {
if (Cyclops::getNumberOfCovariates(cyclopsData) > 1) {
stop("Jeffreys prior is currently only implemented for 1 covariate")
}
covariate <- Cyclops::getCovariateIds(cyclopsData)
if (Cyclops::getCovariateTypes(cyclopsData, covariate) != "indicator") {
count <- reduce(cyclopsData, covariate, power = 0)
sum <- reduce(cyclopsData, covariate, power = 1)
mean <- sum / count
if (!(mean == 0.0 || mean == 1.0)) {
stop("Jeffreys prior is currently only implemented for indicator covariates")
}
}
}
.cyclopsSetPrior(cyclopsData$cyclopsInterfacePtr, prior$priorType, prior$variance,
prior$exclude, graph, neighborhood)
}
if (control$selectorType == "auto") {
if (cyclopsData$modelType %in% c("pr", "lr")) {
control$selectorType <- "byRow"
} else {
rowsPerStratum <- (getNumberOfRows(cyclopsData) / getNumberOfStrata(cyclopsData))
if (rowsPerStratum < getNumberOfStrata(cyclopsData)) {
control$selectorType <- "byPid"
} else {
control$selectorType <- "byRow"
}
}
if (prior$useCrossValidation && control$noiseLevel != "silent") {
writeLines(paste("Using cross-validation selector type", control$selectorType))
}
}
if (control$cvRepetitions == "auto") {
control$cvRepetitions <- .getNumberOfRepetitions(getNumberOfRows(cyclopsData))
}
control <- .setControl(cyclopsData$cyclopsInterfacePtr, control)
threads <- control$threads
if (!is.null(startingCoefficients)) {
if (length(startingCoefficients) != getNumberOfCovariates(cyclopsData)) {
stop("Must provide a value for each coefficient")
}
if (.cyclopsGetHasOffset(cyclopsData)) {
startingCoefficients <- c(1.0, startingCoefficients)
}
.cyclopsSetBeta(cyclopsData$cyclopsInterfacePtr, startingCoefficients)
.cyclopsSetStartingBeta(cyclopsData$cyclopsInterfacePtr, startingCoefficients)
}
if (!is.null(fixedCoefficients)) {
if (length(fixedCoefficients) != getNumberOfCovariates(cyclopsData)) {
stop("Must provide a boolean for each coefficient")
}
offset <- ifelse(.cyclopsGetHasOffset(cyclopsData), 1, 0)
for (i in 1:length(fixedCoefficients)) {
.cyclopsSetFixedBeta(cyclopsData$cyclopsInterfacePtr, offset + i, fixedCoefficients[i] == TRUE)
}
}
# Handle weights
weightsUnsorted <- TRUE
if (!is.null(cyclopsData$weights)) {
if (!is.null(weights)) {
warning("Using weights passed to fitCyclopsModel()")
} else {
weights <- cyclopsData$weights
weightsUnsorted <- FALSE
}
}
if (!is.null(weights)) {
if (length(weights) != getNumberOfRows(cyclopsData)) {
stop("Must provide a weight for each data row")
}
if (any(weights < 0)) {
stop("Only non-negative weights are allowed")
}
if (weightsUnsorted) {
if (!is.null(cyclopsData$sortOrder)) {
weights <- weights[cyclopsData$sortOrder]
}
}
.cyclopsSetWeights(cyclopsData$cyclopsInterfacePtr, weights)
}
# censorWeight check for the Fine-Gray model
if (cyclopsData$modelType == "fgr" & is.null(cyclopsData$censorWeights)) {
stop("Subject-specific censoring weights must be specified for modelType = 'fgr'.")
}
if (!is.null(cyclopsData$censorWeights)) {
if (cyclopsData$modelType != 'fgr' && warnings) {
warning(paste0("modelType = '", cyclopsData$modelType, "' does not use censorWeights. These weights will not be passed further."))
}
if (length(cyclopsData$censorWeights) != getNumberOfRows(cyclopsData)) {
stop("Must provide a censorWeight for each data row")
}
if (any(cyclopsData$censorWeights < 0) || any(cyclopsData$censorWeights > 1)) {
stop("Only weights between 0 and 1 are allowed for censorWeights")
}
.cyclopsSetCensorWeights(cyclopsData$cyclopsInterfacePtr, cyclopsData$censorWeights)
}
if (prior$useCrossValidation) {
minCVData <- control$minCVData
if (control$selectorType == "byRow" && minCVData > getNumberOfRows(cyclopsData)) {
stop("Insufficient data count for cross validation")
}
if (control$selectorType == "byPid" && minCVData > getNumberOfStrata(cyclopsData)) {
stop("Insufficient data count for cross validation")
}
fit <- .cyclopsRunCrossValidation(cyclopsData$cyclopsInterfacePtr)
} else {
fit <- .cyclopsFitModel(cyclopsData$cyclopsInterfacePtr)
}
if (fit$return_flag == "POOR_BLR_STEP" && control$convergenceType == "gradient") {
if (warnings) {
warning("BLR convergence criterion failed; coefficient may be infinite")
}
control$convergenceType <- "lange"
return(fitCyclopsModel(cyclopsData = cyclopsData,
prior = prior,
control = control,
weights = weights,
forceNewObject = forceNewObject,
returnEstimates = returnEstimates,
startingCoefficients = startingCoefficients,
fixedCoefficients = fixedCoefficients,
computeDevice = computeDevice))
}
if (returnEstimates) {
estimates <- .cyclopsLogModel(cyclopsData$cyclopsInterfacePtr)
fit <- c(fit, estimates)
fit$estimation <- as.data.frame(fit$estimation)
}
fit$call <- cl
fit$cyclopsData <- cyclopsData
fit$coefficientNames <- cyclopsData$coefficientNames
if (!is.null(fixedCoefficients)) {
fit$fixedCoefficients <- fixedCoefficients
}
fit$rowNames <- cyclopsData$rowNames
fit$scale <- cyclopsData$scale
fit$threads <- threads
fit$seed <- control$seed
if (prior$useCrossValidation) {
fit$cvRepetitions <- control$cvRepetitions
}
class(fit) <- "cyclopsFit"
return(fit)
}
.checkCovariates <- function(cyclopsData, covariates) {
if (!is.null(covariates)) {
saved <- covariates
indices <- NULL
if (inherits(covariates, "character")) {
# Try to match names
indices <- match(covariates, cyclopsData$coefficientNames)
covariates <- getCovariateIds(cyclopsData)[indices]
}
if (!bit64::is.integer64(covariates)) {
covariates <- bit64::as.integer64(covariates)
}
if (any(is.na(covariates))) {
stop("Unable to match all covariates: ", paste(saved, collapse = ", "))
}
attr(covariates, "indices") <- indices
}
covariates
}
.checkData <- function(x) {
# Check conditions
if (missing(x) || is.null(x$cyclopsDataPtr) || !inherits(x$cyclopsDataPtr, "externalptr")) {
stop("Improperly constructed cyclopsData object")
}
if (.isRcppPtrNull(x$cyclopsDataPtr)) {
stop("Data object is no longer initialized")
}
}
.checkInterface <- function(x, computeDevice = "native", forceNewObject = FALSE, testOnly = FALSE) {
if (forceNewObject
|| is.null(x$cyclopsInterfacePtr)
|| !inherits(x$cyclopsInterfacePtr, "externalptr")
|| .isRcppPtrNull(x$cyclopsInterfacePtr)
|| .cyclopsGetComputeDevice(x$cyclopsInterfacePtr) != computeDevice
) {
if (testOnly == TRUE) {
stop("Interface object is not initialized")
}
# if (computeDevice != "native") {
# stopifnot(computeDevice %in% listOpenCLDevices())
# }
# Build interface
interface <- .cyclopsInitializeModel(x$cyclopsDataPtr, modelType = x$modelType, computeDevice, computeMLE = TRUE)
# TODO Check for errors
assign("cyclopsInterfacePtr", interface$interface, x)
}
}
#' @title Extract model coefficients
#'
#' @description
#' \code{coef.cyclopsFit} extracts model coefficients from an Cyclops model fit object
#'
#' @param object Cyclops model fit object
#' @param rescale Boolean: rescale coefficients for unnormalized covariate values
#' @param ignoreConvergence Boolean: return coefficients even if fit object did not converge
#' @param ... Other arguments
#'
#' @return Named numeric vector of model coefficients.
#'
#' @export
coef.cyclopsFit <- function(object, rescale = FALSE, ignoreConvergence = FALSE, ...) {
if (object$return_flag != "SUCCESS" && !ignoreConvergence) {
stop("Cyclops estimation is null; suspect that estimation did not converge.")
}
result <- object$estimation$estimate
if (is.null(object$coefficientNames)) {
names(result) <- object$estimation$column_label
if ("0" %in% names(result)) {
names(result)[which(names(result) == "0")] <- "(Intercept)"
}
} else {
names(result) <- object$coefficientNames
}
if (!is.null(object$scale) && rescale) {
result <- result * object$scale
}
result
}
#' @title Get hyperparameter
#'
#' @description
#' \code{getHyperParameter} returns the current hyper parameter in a Cyclops model fit object
#'
#' @param object A Cyclops model fit object
#'
#' @template elaborateExample
#'
#' @export
getHyperParameter <- function(object) {
if (inherits(object, "cyclopsFit")) {
object$variance
} else {
NULL
}
}
#' @title Extract log-likelihood
#'
#' @description
#' \code{logLik} returns the current log-likelihood of the fit in a Cyclops model fit object
#'
#' @param object A Cyclops model fit object
#' @param ... Additional arguments
#'
#' @template elaborateExample
#'
#' @export
logLik.cyclopsFit <- function(object, ...) {
out <- object$log_likelihood
attr(out, 'df') <- sum(!is.na(coefficients(object)))
attr(out, 'nobs') <- getNumberOfRows(object$cyclopsData)
class(out) <- 'logLik'
out
}
#' @method print cyclopsFit
#' @title Print a Cyclops model fit object
#'
#' @description
#' \code{print.cyclopsFit} displays information about a Cyclops model fit object
#'
#' @param x A Cyclops model fit object
#' @param show.call Logical: display last call to update the Cyclops model fit object
#' @param ... Additional arguments
#'
#' @export
print.cyclopsFit <- function(x, show.call=TRUE ,...) {
cat("Cyclops model fit object\n\n")
if (show.call && !is.null(x$call)) {
cat("Call: ",paste(deparse(x$call),sep="\n",collapse="\n"),"\n\n",sep="")
}
cat(" Model: ", x$cyclopsData$modelType, "\n", sep="")
cat(" Prior: ", x$prior_info, "\n", sep="")
cat(" Hyperparameter: ", paste(x$variance, collapse=" "), "\n", sep="")
cat(" Return flag: ", x$return_flag, "\n", sep="")
if (x$return_flag == "SUCCESS") {
cat("Log likelikehood: ", x$log_likelihood, "\n", sep="")
cat(" Log prior: ", x$log_prior, "\n", sep="")
}
invisible(x)
}
#' @title Create a Cyclops control object
#'
#' @description
#' \code{createControl} creates a Cyclops control object for use with \code{\link{fitCyclopsModel}}.
#'
#' @param maxIterations Integer: maximum iterations of Cyclops to attempt before returning a failed-to-converge error
#' @param tolerance Numeric: maximum relative change in convergence criterion from successive iterations to achieve convergence
#' @param convergenceType String: name of convergence criterion to employ (described in more detail below)
#' @param cvType String: name of cross validation search.
#' Option \code{"auto"} selects an auto-search following BBR.
#' Option \code{"grid"} selects a grid-search cross validation
#' @param fold Numeric: Number of random folds to employ in cross validation
#' @param lowerLimit Numeric: Lower prior variance limit for grid-search
#' @param upperLimit Numeric: Upper prior variance limit for grid-search
#' @param gridSteps Numeric: Number of steps in grid-search
#' @param cvRepetitions Numeric: Number of repetitions of X-fold cross validation
#' @param minCVData Numeric: Minimum number of data for cross validation
#' @param noiseLevel String: level of Cyclops screen output (\code{"silent"}, \code{"quiet"}, \code{"noisy"})
#' @param threads Numeric: Specify number of CPU threads to employ in cross-validation; default = 1 (auto = -1)
#' @param seed Numeric: Specify random number generator seed. A null value sets seed via \code{\link{Sys.time}}.
#' @param resetCoefficients Logical: Reset all coefficients to 0 between model fits under cross-validation
#' @param startingVariance Numeric: Starting variance for auto-search cross-validation; default = -1 (use estimate based on data)
#' @param useKKTSwindle Logical: Use the Karush-Kuhn-Tucker conditions to limit search
#' @param tuneSwindle Numeric: Size multiplier for active set
#' @param selectorType String: name of exchangeable sampling unit.
#' Option \code{"byPid"} selects entire strata.
#' Option \code{"byRow"} selects single rows.
#' If set to \code{"auto"}, \code{"byRow"} will be used for all models except conditional models where
#' the average number of rows per stratum is smaller than the number of strata.
#' @param initialBound Numeric: Starting trust-region size
#' @param maxBoundCount Numeric: Maximum number of tries to decrease initial trust-region size
#' @param algorithm String: name of fitting algorithm to employ; default is `ccd`
#'
#' Todo: Describe convegence types
#'
#' @return
#' A Cyclops control object of class inheriting from \code{"cyclopsControl"} for use with \code{\link{fitCyclopsModel}}.
#'
#' @template elaborateExample
#'
#' @export
createControl <- function(maxIterations = 1000,
tolerance = 1E-6,
convergenceType = "gradient",
cvType = "auto",
fold = 10,
lowerLimit = 0.01,
upperLimit = 20.0,
gridSteps = 10,
cvRepetitions = 1,
minCVData = 100,
noiseLevel = "silent",
threads = 1,
seed = NULL,
resetCoefficients = FALSE,
startingVariance = -1,
useKKTSwindle = FALSE,
tuneSwindle = 10,
selectorType = "auto",
initialBound = 2.0,
maxBoundCount = 5,
algorithm = "ccd") {
validCVNames = c("grid", "auto")
stopifnot(cvType %in% validCVNames)
validNLNames = c("silent", "quiet", "noisy")
stopifnot(noiseLevel %in% validNLNames)
stopifnot(threads == -1 || threads >= 1)
stopifnot(startingVariance == -1 || startingVariance > 0)
stopifnot(selectorType %in% c("auto","byPid", "byRow"))
validAlgorithmNames = c("ccd", "mm")
stopifnot(algorithm %in% validAlgorithmNames)
structure(list(maxIterations = maxIterations,
tolerance = tolerance,
convergenceType = convergenceType,
autoSearch = (cvType == "auto"),
fold = fold,
lowerLimit = lowerLimit,
upperLimit = upperLimit,
gridSteps = gridSteps,
minCVData = minCVData,
cvRepetitions = cvRepetitions,
noiseLevel = noiseLevel,
threads = threads,
seed = seed,
resetCoefficients = resetCoefficients,
startingVariance = startingVariance,
useKKTSwindle = useKKTSwindle,
tuneSwindle = tuneSwindle,
selectorType = selectorType,
initialBound = initialBound,
maxBoundCount = maxBoundCount,
algorithm = algorithm),
class = "cyclopsControl")
}
#' @title Create a Cyclops prior object
#'
#' @description
#' \code{createPrior} creates a Cyclops prior object for use with \code{\link{fitCyclopsModel}}.
#'
#' @param priorType Character: specifies prior distribution. See below for options
#' @param variance Numeric: prior distribution variance
#' @param exclude A vector of numbers or covariateId names to exclude from prior
#' @param graph Child-to-parent mapping for a hierarchical prior
#' @param neighborhood A list of first-order neighborhoods for a partially fused prior
#' @param useCrossValidation Logical: Perform cross-validation to determine prior \code{variance}.
#' @param forceIntercept Logical: Force intercept coefficient into prior
#'
#' @section Prior types:
#'
#' We specify all priors in terms of their variance parameters.
#' Similar fitting tools for regularized regression often parameterize the Laplace distribution
#' in terms of a rate \code{"lambda"} per observation.
#' See \code{"glmnet"}, for example.
#'
#' variance = 2 * / (nobs * lambda)^2 or lambda = sqrt(2 / variance) / nobs
#'
#' @template elaborateExample
#'
#' @return
#' A Cyclops prior object of class inheriting from \code{"cyclopsPrior"} for use with \code{fitCyclopsModel}.
#'
#' @export
createPrior <- function(priorType,
variance = 1,
exclude = c(),
graph = NULL,
neighborhood = NULL,
useCrossValidation = FALSE,
forceIntercept = FALSE) {
validNames = c("none", "laplace","normal", "barupdate", "hierarchical", "jeffreys")
stopifnot(priorType %in% validNames)
if (!is.null(exclude)) {
if (!inherits(exclude, "character") &&
!inherits(exclude, "numeric") &&
!inherits(exclude, "integer")
) {
stop(cat("Unable to parse excluded covariates:"), exclude)
}
}
if (length(priorType) != length(variance)) {
stop("Prior types and variances have a dimensionality mismatch")
}
if (all(priorType == "none") && useCrossValidation) {
stop("Cannot perform cross validation with a flat prior")
}
if (any(priorType == "barupdate") && useCrossValidation) {
stop("Cannot perform cross valudation with BAR updates")
}
if (any(priorType == "hierarchical") && missing(graph)) {
stop("Must provide a graph for a hierarchical prior")
}
if (!is.null(neighborhood)) {
allNames <- unlist(neighborhood)
if (!inherits(allNames, "character") &&
!inherits(allNames, "numeric") &&
!inherits(allNames, "integer")) {
stop(cat("Unable to parse neighborhood covariates:"), allNames)
}
}
structure(list(priorType = priorType, variance = variance, exclude = exclude,
graph = graph,
neighborhood = neighborhood,
useCrossValidation = useCrossValidation, forceIntercept = forceIntercept),
class = "cyclopsPrior")
}
# .cyclopsSetCoefficients <- function(object, coefficients) {
# .checkInterface(object$cyclopsData, testOnly = TRUE)
#
# if (length(coefficients) != getNumberOfCovariates(object$cyclopsData)) {
# stop("Must provide a value for each coefficient")
# }
#
# if (.cyclopsGetHasOffset(object$cyclopsData)) {
# coefficients <- c(1.0, coefficients)
# }
#
# .cyclopsSetBeta(object$cyclopsData$cyclopsInterfacePtr, coefficients)
# }
#' @title Compute predictive log-likelihood from a Cyclops model fit
#'
#' @description
#' \code{getCyclopsPredictiveLogLikelihood} returns the log-likelihood of a subset of the data in a Cyclops model fit object.
#'
#' @param object A Cyclops model fit object
#' @param weights Numeric vector: vector of 0/1 identifying subset (=1) of rows from \code{object} to use in computing the log-likelihood
#' @return The predictive log-likelihood
#'
#' @keywords internal
getCyclopsPredictiveLogLikelihood <- function(object, weights) {
.checkInterface(object$cyclopsData, testOnly = TRUE)
if (length(weights) != getNumberOfRows(object$cyclopsData)) {
stop("Must provide a weight for each data row")
}
if (any(weights < 0)) {
stop("Only non-negative weights are allowed")
}
if(!is.null(object$cyclopsData$sortOrder)) {
weights <- weights[object$cyclopsData$sortOrder]
}
# TODO Remove code duplication with weights section of fitCyclopsModel
.cyclopsGetNewPredictiveLogLikelihood(object$cyclopsData$cyclopsInterfacePtr, weights)
}
#' @title Get cross-validation information from a Cyclops model fit
#'
#' @description {getCrossValidationInfo} returns the predicted optimal cross-validation point and ordinate
#'
#' @param object A Cyclops model fit object
#'
#' @keywords internal
getCrossValidationInfo <- function(object) {
info <- object$cross_validation
if (is.na(info) || info == "") {
stop("No cross-validation information is available")
}
values <- as.numeric(unlist(strsplit(info, " ")))
list(ordinate = values[1],
point = values[-1])
}
.setControl <- function(cyclopsInterfacePtr, control) {
if (!missing(control)) {
stopifnot(inherits(control, "cyclopsControl"))
if (is.null(control$seed)) {
control$seed <- as.integer(Sys.time())
}
if (is.null(control$algorithm) || is.na(control$algorithm)) { # Provide backwards compatibility
control$algorithm <- "ccd"
}
.cyclopsSetControl(cyclopsInterfacePtr, control$maxIterations, control$tolerance,
control$convergenceType, control$autoSearch, control$fold,
(control$fold * control$cvRepetitions),
control$lowerLimit, control$upperLimit, control$gridSteps,
control$noiseLevel, control$threads, control$seed, control$resetCoefficients,
control$startingVariance, control$useKKTSwindle, control$tuneSwindle,
control$selectorType, control$initialBound, control$maxBoundCount,
control$algorithm
)
return(control)
}
return(NULL)
}
#' @title Extract standard errors
#'
#' @description
#' \code{getSEs} extracts asymptotic standard errors for specific covariates from a Cyclops model fit object.
#'
#' @details This function first computes the (partial) Fisher information matrix for
#' just the requested covariates and then returns the square root of the diagonal elements of
#' the inverse of the Fisher information matrix. These are the asymptotic standard errors
#' when all possible covariates are included.
#' When the requested covariates do not equate to all coefficients in the model,
#' then interpretation is more challenging.
#'
#' @param object A Cyclops model fit object
#' @param covariates Integer or string vector: list of covariates for which asymptotic standard errors are wanted
#'
#' @return Vector of standard error estimates
#'
#' @keywords internal
getSEs <- function(object, covariates) {
.checkInterface(object$cyclopsData, testOnly = TRUE)
covariates <- .checkCovariates(object$cyclopsData, covariates)
if (getNumberOfCovariates(object$cyclopsData) != length(covariates)) {
warning("Asymptotic standard errors are only valid if computed for all covariates simultaneously")
}
fisherInformation <- .cyclopsGetFisherInformation(object$cyclopsData$cyclopsInterfacePtr, covariates)
ses <- sqrt(diag(solve(fisherInformation)))
names(ses) <- object$coefficientNames[as.integer(covariates)]
ses
}
#' @title Confidence intervals for Cyclops model parameters
#'
#' @description
#' \code{confinit.cyclopsFit} profiles the data likelihood to construct confidence intervals of
#' arbitrary level. Usually it only makes sense to do this for variables that have not been regularized.
#'
#' @param object A fitted Cyclops model object
#' @param parm A specification of which parameters require confidence intervals,
#' either a vector of numbers of covariateId names
#' @param level Numeric: confidence level required
## @param control A Cyclops \code{\link{control}} object
#' @param overrideNoRegularization Logical: Enable confidence interval estimation for regularized parameters
#' @param includePenalty Logical: Include regularized covariate penalty in profile
#' @param rescale Boolean: rescale coefficients for unnormalized covariate values
#' @param ... Additional argument(s) for methods
#'
#' @return
#' A matrix with columns reporting lower and upper confidence limits for each parameter.
#' These columns are labelled as (1-level) / 2 and 1 - (1 - level) / 2 in percent
#' (by default 2.5 percent and 97.5 percent)
#'
#' @template elaborateExample
#'
#' @export
confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
overrideNoRegularization = FALSE,
includePenalty = TRUE,
rescale = FALSE, ...) {
.checkInterface(object$cyclopsData, testOnly = TRUE)
#.setControl(object$cyclopsData$cyclopsInterfacePtr, control)
parm <- .checkCovariates(object$cyclopsData, parm)
if (level < 0.01 || level > 0.99) {
stop("level must be between 0 and 1")
}
threshold <- qchisq(level, df = 1) / 2
threads <- object$threads
if (!is.null(object$fixedCoefficients)) {
if (any(object$fixedCoefficients[as.integer(parm)])) {
stop("Cannot estimate confidence interval for a fixed coefficient")
}
}
prof <- .cyclopsProfileModel(object$cyclopsData$cyclopsInterfacePtr, parm,
threads, threshold,
overrideNoRegularization,
includePenalty)
indices <- match(parm, getCovariateIds(object$cyclopsData))
if (!is.null(object$scale) && rescale) {
prof$lower <- prof$lower * object$scale[indices]
prof$upper <- prof$upper * object$scale[indices]
}
prof <- as.matrix(as.data.frame(prof))
rownames(prof) <- object$coefficientNames[indices]
qs <- c((1 - level) / 2, 1 - (1 - level) / 2) * 100
colnames(prof)[2:3] <- paste(sprintf("%.1f", qs), "%")
# Change NaN to NA
prof[which(is.nan(prof[, 2])), 2] <- NA
prof[which(is.nan(prof[, 3])), 3] <- NA
prof
}
.initAdaptiveProfile <- function(object, parm, bounds, includePenalty) {
# If an MLE was found, let's not throw that bit of important information away:
if (object$return_flag == "SUCCESS" &&
coef(object)[as.character(parm)] > bounds[1] &&
coef(object)[as.character(parm)] < bounds[2]) {
profile <- tibble(point = coef(object)[as.character(parm)],
value = fixedGridProfileLogLikelihood(object, parm, coef(object)[as.character(parm)], includePenalty)$value)
} else {
profile <- tibble()
}
}
#' @title Profile likelihood for Cyclops model parameters
#'
#' @description
#' \code{getCyclopsProfileLogLikelihood} evaluates the profile likelihood at a grid of parameter values.
#'
#' @param object Fitted Cyclops model object
#' @param parm Specification of which parameter requires profiling,
#' either a vector of numbers of covariateId names
#' @param x Vector of values of the parameter
#' @param bounds Pair of values to bound adaptive profiling
#' @param tolerance Absolute tolerance allowed for adaptive profiling
#' @param initialGridSize Initial grid size for adaptive profiling
#' @param includePenalty Logical: Include regularized covariate penalty in profile
#'
#' @return
#' A data frame containing the profile log likelihood. Returns NULL when the adaptive profiling fails
#' to converge.
#'
#' @export
getCyclopsProfileLogLikelihood <- function(object,
parm,
x = NULL,
bounds = NULL,
tolerance = 1E-3,
initialGridSize = 10,
includePenalty = TRUE) {
maxResets <- 10
if (!xor(is.null(x), is.null(bounds))) {
stop("Must provide either `x` or `bounds`, but not both.")
}
if (!is.null(bounds)) { # Adaptive profiling using recursive calls
if (length(bounds) != 2 || bounds[1] >= bounds[2]) {
stop("Must provide bounds[1] < bounds[2]")
}
profile <- .initAdaptiveProfile(object, parm, bounds, includePenalty)
# Start with sparse grid:
grid <- seq(bounds[1], bounds[2], length.out = initialGridSize)
# Iterate until stopping criteria met:
priorMaxMaxError <- Inf
resetsPerformed <- 0
while (length(grid) != 0) {
ll <- fixedGridProfileLogLikelihood(object, parm, grid, includePenalty)
profile <- bind_rows(profile, ll) %>% arrange(.data$point)
if (any(is.nan(profile$value))) {
if (all(is.nan(profile$value))) {
warning("Failing to compute likelihood at entire initial grid.")
return(NULL)
}
start <- min(which(!is.nan(profile$value)))
end <- max(which(!is.nan(profile$value)))
if (start == end) {
warning("Failing to compute likelihood at entire grid except one. Giving up")
return(NULL)
}
profile <- profile[start:end, ]
if (any(is.nan(profile$value))) {
warning("Failing to compute likelihood in non-extreme regions. Giving up.")
return(NULL)
}
warning("Failing to compute likelihood at extremes. Truncating bounds.")
}
deltaX <- profile$point[2:nrow(profile)] - profile$point[1:(nrow(profile) - 1)]
deltaY <- profile$value[2:nrow(profile)] - profile$value[1:(nrow(profile) - 1)]
slopes <- deltaY / deltaX
if (resetsPerformed < maxResets && !all(slopes[2:length(slopes)] < slopes[1:(length(slopes)-1)])) {
warning("Coefficient drift detected. Resetting Cyclops object and recomputing all likelihood values computed so far.")
grid <- profile$point
profile <- tibble()
interface <- .cyclopsInitializeModel(object$cyclopsData$cyclopsDataPtr, modelType = object$cyclopsData$modelType, "native", computeMLE = TRUE)
assign("cyclopsInterfacePtr", interface$interface, object)
resetsPerformed <- resetsPerformed + 1
next
}
# Compute where prior and posterior slopes intersect
slopes <- c(slopes[1] + (slopes[2] - slopes[3]),
slopes,
slopes[length(slopes)] - (slopes[length(slopes) - 1] - slopes[length(slopes)]))
interceptX <- (profile$value[2:nrow(profile)] -
profile$point[2:nrow(profile)] * slopes[3:length(slopes)] -
profile$value[1:(nrow(profile) - 1)] +
profile$point[1:(nrow(profile) - 1)] * slopes[1:(length(slopes) - 2)]) /
(slopes[1:(length(slopes) - 2)] - slopes[3:length(slopes)])
# Compute absolute difference between linear interpolation and worst case scenario (which is at the intercept):
maxError <- abs((profile$value[1:(nrow(profile) - 1)] + (interceptX - profile$point[1:(nrow(profile) - 1)]) * slopes[1:(length(slopes) - 2)]) -
(profile$value[1:(nrow(profile) - 1)] + (interceptX - profile$point[1:(nrow(profile) - 1)]) * slopes[2:(length(slopes) - 1)]))
maxMaxError <- max(maxError, na.rm = TRUE)
if (is.na(maxMaxError) || maxMaxError > priorMaxMaxError) {
warning("Failing to converge when using adaptive profiling.")
return(NULL)
}
priorMaxMaxError <- maxMaxError
exceed <- which(maxError > tolerance)
grid <- (profile$point[exceed] + profile$point[exceed + 1]) / 2
}
} else { # Use x
profile <- fixedGridProfileLogLikelihood(object, parm, x, includePenalty)
}
return(profile)
}
fixedGridProfileLogLikelihood <- function(object, parm, x, includePenalty) {
.checkInterface(object$cyclopsData, testOnly = TRUE)
parm <- .checkCovariates(object$cyclopsData, parm)
threads <- object$threads
if (getNumberOfCovariates(object$cyclopsData) == 1 || length(x) == 1) {
grid <- .cyclopsGetProfileLikelihood(object$cyclopsData$cyclopsInterfacePtr, parm, x,
threads, includePenalty)
} else {
# Partition sequence
y <- sort(x)
midPt <- floor(length(x) / 2)
lower <- y[midPt:1]