-
Notifications
You must be signed in to change notification settings - Fork 17
/
Copy pathfit_and_reporting.R
5014 lines (4701 loc) · 201 KB
/
fit_and_reporting.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
# Copyright 2007-2022 Timothy C. Bates
#
# 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
#
# https://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.
#' Quickly plot y ~ x with a regression line and R^2, and nice labels.
#'
#' @description
#' Want a figure for your paper or presentation but not the work of combining [ggplot2::ggplot()],
#' [ggplot2::geom_smooth()] and `method` options, plus [ggplot2::geom_point()]. Organizing [ggplot2::labs()] and its
#' `x`, `y`, and `title` components. Adding your preferred theme like [ggplot2::theme_gray()], plus recalling
#' for [cowplot::draw_label()], and/or [ggplot2::annotate()] to draw math-enabled labels on the plot, as well
#' as the required [bquote()], and extracting the relevant fit statistics from [lm()] and the subsidiary tasks
#' of [reformulate()] programmatic variables?
#'
#' `umxPlot` just takes `y ~ x` (or "x" and "y" as strings), and gives you a nicely labelled plot, with a fitted line, the
#' \ifelse{html}{\out{<em>R</em><sup>2</sup>}}{\eqn{R^2}} so readers can see how well this fitted. It knows how to put Greek symbols like \eqn{beta} into axes.
#' @param x formula or (alternatively) x as string
#' @param y variable as string.
#' @param data The data for the graph.
#' @param xlab X-axis label (default y).
#' @param ylab Y-axis label (default y).
#' @param title Graph title. Default = paste0(y, " as a function of ", x)
#' @param r2x x location for the fit summary (default 1).
#' @param r2y y location for the fit summary (default 2).
#' @param geom_point show points? (TRUE)
#' @param method Method for fitting curve (default = lm)
#' @param family for glm default = "gaussian"
#' @return - plot you can edit.
#' @export
#' @family Plotting functions
#' @seealso - [ggplot2::qplot()]
#' @md
#' @examples
#' data(mtcars)
#' umxPlot(mpg ~ wt, data = mtcars, r2x = 2, r2y = 10)
#' umxPlot(x = "wt", y = "mpg", mtcars, r2x = 2, r2y = 10)
umxPlot <- function(x, y= NULL, data, xlab= x, ylab = y, title = paste0(y, " as a function of ", x), r2x=NA, r2y=NA, geom_point = TRUE, method = c("lm", "auto", "loess", "glm", "gam"), family = c("gaussian","binomial", "Gamma", "inverse", "poisson", "quasi", "quasibinomial", "quasipoisson")) {
method = match.arg(method)
family = match.arg(family)
if(inherits(x, "formula")){
.formula = x
tmp = attr(terms(x), "factors")
tmp = dimnames(tmp)[[1]]
y = tmp[1]
x = tmp[2]
umx_check_names(c(x, y), data = data, die = TRUE)
} else if(is.null(y)){
# let qplot make a histogram
return(qplot(x, data = data, xlab =x, main = title))
} else {
# make formula from x and y strings
umx_check_names(c(x, y), data = data, die = TRUE)
.formula = reformulate(paste0(y, "~ ", x))
}
if(geom_point){
p = ggplot(data = data, aes_string(x, y)) + geom_smooth(method = method) + geom_point()
} else {
p = ggplot(data = data, aes_string(x, y)) + geom_smooth(method = method)
}
# data[,x]
if(is.na(r2x)){
r2x = min(data[,x])
r2y = max(data[,y])
}
if(method == "lm"){
m1 = lm(.formula, data = data)
r2 = round(summary(m1)$r.squared, 3)
lab = bquote(R^2 == .(r2))
p = p + labs(x= xlab, y= ylab, title= title)
p = p + cowplot::draw_label(lab, x = r2x, y = r2y, fontfamily = "Times", size = 12)
}else if (method == "glm"){
# m1 = glm(.formula, data = data, family=family)
message("polite note: Currently, I only know how to do method = lm")
}else{
message("polite note: Currently, I only know how to do method = lm or glm")
}
p = p + theme_minimal() # gray, bw, linedraw, light, dark, minimal, classic
p
# p + annotate("text", 3, 30, label = expression(R^2 == beta + 1 ~ hello), family="Optima")
}
#' `umxPlotPredict` takes a model and plots the y value against predicted(y)
#' @param model lm or other model that understands predict()
#' @param xlab X-axis label (default x).
#' @param ylab Y-axis label (default y).
#' @param r2x x location for the fit summary.
#' @param r2y y location for the fit summary.
#' @param font_size Default 13
#' @param font Default "Times"
#' @return - plot you can edit.
#' @export
#' @family Plotting functions
#' @seealso - [ggplot2::qplot()]
#' @md
#' @examples
#' data(mtcars)
#' tmp = lm(mpg ~ wt, data = mtcars)
#' umxPlotPredict(tmp, r2x = 2, r2y = 10)
umxPlotPredict <- function(model, xlab= "Predicted Y", ylab= "Observed Y", r2x= 1.5, r2y= 4.5, font_size = 13, font= "Times") {
r2 = paste0("r = ", round(summary(model)$adj.r.squared^.5, 2))
y_var = model.frame(model)[, 1]
p = ggplot() + geom_point(aes(x = predict(model), y = y_var))
p = p + geom_smooth(aes(x = predict(model), y = y_var), method = "lm", se = TRUE, color = "blue")
p = p + labs(x= xlab, y= ylab)
p = p + cowplot::draw_label(r2, x = r2x, y = r2y, fontfamily = font, size = (font_size+1))
p + theme_minimal(base_size = font_size, base_family= font)
}
# =====================
# = Model Diagnostics =
# =====================
#' Diagnose problems in a model - not working!
#'
#' The goal of this function **WILL BE** (not currently functional) to diagnose problems in
#' a model and return suggestions to the user.
#' It is a work in progress, and of no use as yet.
#'
#' Best diagnostics are:
#'
#' 1. Observed data variances and means
#' 2. Expected variances and means
#' 3. Difference of these?
#'
#' Try
#' * diagonalizeExpCov diagonal
#' * [umx_is_ordered()]
#'
#' Tricky, but reporting variances and standardized thresholds is ideal.
#' Guidance is to start with unit variances and thresholds within
#' +/- 2 SD of the mean. Like %p option in Classic Mx.
#' @param model an [OpenMx::mxModel()] to diagnose
#' @param tryHard whether I should try and fix it? (defaults to FALSE)
#' @param diagonalizeExpCov Whether to diagonalize the ExpCov
#' @return - helpful messages and perhaps a modified model
#' @export
#' @family Teaching and Testing functions
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#' manifests = names(demoOneFactor)
#'
#' m1 = umxRAM("OneFactor", data = demoOneFactor, type= "cov",
#' umxPath("G", to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = "G", fixedAt = 1)
#' )
#' m1 = mxRun(m1)
#' umxSummary(m1, std = TRUE)
#' umxDiagnose(m1)
#' }
umxDiagnose <- function(model, tryHard = FALSE, diagonalizeExpCov = FALSE){
# 1. First thing to check is whether the covariance matrix is positive definite.
minEigen = min(eigen(umxExpCov(model))$values)
if(minEigen<0){
message("The expected covariance matrix is not positive definite")
# now what?
}
}
# =============================
# = Fit and Reporting Helpers =
# =============================
#' AIC weight-based conditional probabilities.
#'
#' @description
#' Returns the best model by AIC, and computes the probabilities
#' according to AIC weight-based conditional probabilities (Wagenmakers & Farrell, 2004).
#'
#' @param models a list of models to compare.
#' @param digits (default 2)
#' @return - Best model
#' @export
#' @family Miscellaneous Stats Functions
#' @seealso - [AIC()]
#' @references - Wagenmakers E.J., Farrell S. (2004), 192-196. AIC model selection using Akaike weights. *Psychonomic Bulletin and Review*. **11**, 192-196. \url{https://pubmed.ncbi.nlm.nih.gov/15117008/}
#' @md
#' @examples
#' l1 = lm(mpg~ wt + disp, data=mtcars)
#' l2 = lm(mpg~ wt, data=mtcars)
#' umxWeightedAIC(models = list(l1, l2))
umxWeightedAIC <- function(models, digits= 2) {
if(inherits(models[[1]], "numeric")){
stop("Please input the list of models to compare as a list, i.e. models = list(model1, model2)")
}
AIClist = c()
for (i in models) {
AIClist = c(AIClist, AIC(i))
}
whichBest = which.min(AIClist)
bestModel = models[[whichBest]]
aic.weights = round(MuMIn::Weights(AIClist), 2)
if(isS4(models[[1]]) & is(models[[1]], "MxModel")){
message("The ", omxQuotes(bestModel$name), " model is the best fitting model according to AIC.")
# Probabilities according to AIC Weights (Wagenmakers et al https://pubmed.ncbi.nlm.nih.gov/15117008/ )
message("AIC weight-based conditional probabilities {Wagenmakers, 2004, 192-196} of being the best model for ",
omxQuotes(namez(models)), " respectively are: ",
omxQuotes(aic.weights), " Using MuMIn::Weights(AIC()).")
}else{
if("call" %in% names(bestModel)){
# ID = paste0("Model ", omxQuotes(bestModel$call))
ID = paste0("Model ", whichBest)
} else {
ID = paste0("Model ", whichBest)
}
message(ID, " is the best fitting model according to AIC.")
message("AIC weight-based conditional probabilities {Wagenmakers, 2004, 192-196} of being the best model are (for each model you gave me): ",
omxQuotes(aic.weights), " Using MuMIn::Weights(AIC()).")
}
invisible(bestModel)
}
#' Reduce models, and report the results.
#'
#' @description
#' Given a `umx` model (currently `umxACE` and `umxGxE` are supported - ask for more!)
#' `umxReduce` will conduct a formalised reduction process. It will also report
#' Akaike weights are also reported showing relative support across models.
#'
#' Specialized functions are called for different type of input:
#' 1. **GxE model reduction** For [umxGxE()] models [umxReduceGxE()] is called.
#' 2. **ACE model reduction** For [umxACE()] models,[umxReduceACE()] is called.
#'
#' `umxReduce` reports the results in a table. Set the format of the table with
#' [umx_set_table_format()], or set `report= "html"` to open a
#' table for pasting into a word processor.
#'
#' `umxReduce` can be extended to new cases as demand emerges.
#'
#' @param model The [OpenMx::mxModel()] which will be reduced.
#' @param report How to report the results. "html" = open in browser
#' @param intervals Recompute CIs (if any included) on the best model (default = TRUE)
#' @param testD Whether to test ADE and DE models (TRUE)
#' @param baseFileName (optional) custom filename for html output (defaults to "tmp")
#' @param tryHard Default = "yes"
#' @param silent Default = FALSE
#' @param ... Other parameters to control model summary
#' @family Model Summary and Comparison
#' @family Twin Modeling Functions
#' @seealso [umxReduceGxE()], [umxReduceACE()]
#' @references - Wagenmakers, E.J., & Farrell, S. (2004). AIC model selection using Akaike weights.
#' *Psychonomic Bulletin and Review*, **11**, 192-196. \doi{10.3758/BF03206482}
#' @export
#' @md
umxReduce <- function(model, report = c("markdown", "inline", "html"), intervals = TRUE, testD = TRUE, baseFileName = "tmp", tryHard = "yes", silent=FALSE, ...){
UseMethod("umxReduce", model)
}
#' @export
umxReduce.default <- function(model, report = c("markdown", "inline", "html"), intervals = FALSE, testD = TRUE, baseFileName = "tmp", tryHard = "yes", silent=FALSE, ...){
stop("umxReduce is not defined for objects of class:", class(model))
}
#' Reduce a GxE model.
#'
#' @description
#' This function can perform model reduction for [umxGxE()] models,
#' testing dropping a`,c` & e`, as well as c & c`, a & a` etc.
#'
#' It reports the results in a table. Set the format of the table with
#' [umx_set_table_format()]. Or set `report = "html"` to open a
#' table for pasting into a word processor.
#'
#' In addition to printing a table, the function returns the preferred model.
#'
#' @param model A [umxGxE()] to reduce.
#' @param report How to report the results. default = "markdown". "html" = open in browser.
#' @param intervals Recompute CIs (if any included) on the best model (default = TRUE)
#' @param testD Whether to test ADE and DE models (TRUE)
#' @param baseFileName (optional) custom filename for html output (default = "tmp").
#' @param tryHard Default ('no') uses normal mxRun. "yes" uses mxTryHard. Other options: "ordinal", "search"
#' @param silent Default (FALSE)
#' @param ... Other parameters to control model summary.
#' @return best model
#' @export
#' @family Twin Modeling Functions
#' @seealso [umxReduce()], [umxReduceACE()]
#' @references - Wagenmakers, E.J., & Farrell, S. (2004). AIC model selection using Akaike weights.
#' *Psychonomic Bulletin and Review*, **11**, 192-196. \doi{10.3758/BF03206482}.
#' @md
#' @examples
#' \dontrun{
#' model = umxReduce(model)
#' }
umxReduceGxE <- function(model, report = c("markdown", "inline", "html", "report"), intervals = TRUE, testD = TRUE, baseFileName = "tmp_gxe", tryHard = c("yes", "no", "ordinal", "search"), silent = FALSE, ...) {
report = match.arg(report)
umx_is_MxModel(model)
if(inherits(model, "MxModelGxE")){
# Reduce GxE Model
noAmod = umxModify(model, update = "am_r1c1", name = "No_mod_on_A", tryHard= tryHard)
noCmod = umxModify(model, update = "cm_r1c1", name = "No_mod_on_C", tryHard= tryHard)
noEmod = umxModify(model, update = "em_r1c1", name = "No_mod_on_E", tryHard= tryHard)
noACEmod = umxModify(model, regex = "[ace]m_r1c1" , name = "No_moderation", tryHard= tryHard)
no_a_no_am = umxModify(noAmod , update = "a_r1c1" , name = "No_A_no_mod_on_A" , tryHard= tryHard)
no_c_no_cm = umxModify(noCmod , update = "c_r1c1" , name = "No_C_no_mod_on_C" , tryHard= tryHard)
no_c_no_cem = umxModify(no_c_no_cm , update = "em_r1c1", name = "No_c_no_ce_mod" , tryHard= tryHard)
no_c_no_mod = umxModify(no_c_no_cem, update = "am_r1c1", name = "No_c_no_moderation", tryHard= tryHard)
# if("testFor quad11 in parameters"){
# no_sq_mean = umxModify(noACEmod, update = "quad11" , name = "No_mod_no_quad_mean", tryHard= tryHard)
# no_lin_mean = umxModify(noACEmod, update = "lin11" , name = "No_mod_no_lin_mean" , tryHard= tryHard)
# nomeans = umxModify(noACEmod, regex = "lin|quad", name = "No_mod_no_means_mod", tryHard= tryHard)
# } else {
# # likely has non-equal moderators don't bother trying to drop
# # betaSelf and betaCoTwin they are almost certainly necessary
# }
comparisons = c(
noAmod, noCmod, noEmod, noACEmod,
no_a_no_am, no_c_no_cm, no_c_no_cem,
no_c_no_mod
)
# ====================
# = everything table =
# ====================
umxCompare(model, comparisons, all = TRUE, report = report, file = paste0(baseFileName, "1.html"))
# umxCompare(no_c_no_cem, no_c_no_moderation, all = TRUE, report = report, file = paste0(baseFileName, "2.html"))
modelList = c(model, comparisons)
# get list of AICs
AIClist = c()
for (i in modelList) {
AIClist = c(AIClist, AIC(i))
}
whichBest = which.min(AIClist)
bestModel = modelList[[whichBest]]
message("The ", omxQuotes(bestModel$name), " model is the best fitting model according to AIC.")
# Probabilities according to AIC MuMIn::Weights (Wagenmakers et al https://pubmed.ncbi.nlm.nih.gov/15117008/ )
aic.weights = round(Weights(AIClist), 2)
message("AIC weight-based conditional probabilities {Wagenmakers, 2004, 192-196} of being the best model for ",
omxQuotes(namez(modelList)), " respectively are: ",
omxQuotes(aic.weights), " Using MuMIn::Weights(AIC())."
)
if(intervals){
bestModel = mxRun(bestModel, intervals = intervals)
}
invisible(bestModel)
} else {
stop("This function is for GxE. Feel free to let me know what you want...")
}
}
#' @export
umxReduce.MxModelGxE <- umxReduceGxE
#' Reduce an ACE model.
#'
#' This function can perform model reduction on [umxACE()] models,
#' testing dropping A and C, as well as an ADE or ACE model, displaying the results
#' in a table, and returning the best model.
#'
#' It is designed for testing univariate models. You can offer up either the ACE or ADE base model.
#'
#' Suggestions for more sophisticated automation welcomed!
#'
#' @param model an ACE or ADE [OpenMx::mxModel()] to reduce
#' @param report How to report the results. "html" = open in browser
#' @param intervals Recompute CIs (if any included) on the best model (default = TRUE)
#' @param testD Whether to test ADE and DE models (TRUE)
#' @param baseFileName (optional) custom filename for html output (defaults to "tmp")
#' @param tryHard (default = "yes")
#' @param silent Don't print the ACE models (default = FALSE)
#' @param digits rounding in printout (default = 2)
#' @param ... Other parameters to control model summary
#' @return Best fitting model
#' @export
#' @family Twin Modeling Functions
#' @seealso [umxReduceGxE()], [umxReduce()]
#' @references - Wagenmakers, E.J., & Farrell, S. (2004). AIC model selection using Akaike weights. *Psychonomic Bulletin and Review*, **11**, 192-196. \doi{10.3758/BF03206482}
#' @md
#' @examples
#' \dontrun{
#' data(twinData)
#' mzData = subset(twinData, zygosity == "MZFF")
#' dzData = subset(twinData, zygosity == "DZFF")
#' m1 = umxACE(selDVs = "bmi", dzData = dzData, mzData = mzData, sep = "")
#'
#' # ===========================================================================
#' # = Table of parameters + fit comparisons, ready too copy to word processor =
#' # ===========================================================================
#' umxReduce(m1, silent=TRUE, digits=2, repo="h")
#'
#' # ==========================================
#' # = Function captures the preferred model =
#' # ==========================================
#' m2 = umxReduce(m1)
#' umxSummary(m2)
#'
#' # works for ADE input also
#' m1 = umxACE(selDVs = "bmi", dzData = dzData, mzData = mzData, sep = "", dzCr = .25)
#'
#' }
umxReduceACE <- function(model, report = c("markdown", "inline", "html", "report"), intervals = TRUE, testD = TRUE, baseFileName = "tmp", tryHard = c("yes", "no", "ordinal", "search"), silent=FALSE, digits = 2, ...) {
# override umx_set_auto_run?
oldAutoRun = umx_set_auto_run(autoRun = TRUE)
report = match.arg(report)
tryHard = match.arg(tryHard)
if(silent){
oldSilent = umx_set_silent(TRUE)
}else{
oldSilent = FALSE
}
oldAutoPlot = umx_set_auto_plot(FALSE, silent = TRUE)
if(model$top$dzCr$values == 1){
ACE = model
ADE = umxModify(model, 'dzCr_r1c1', value = .25, name = "ADE", tryHard = tryHard)
if(-2*logLik(ACE) > -2*logLik(ADE) && testD){
CE = umxModify(ADE, regex = "a_r[0-9]+c[0-9]+" , name = "DE", tryHard = tryHard)
AE = umxModify(ADE, regex = "c_r[0-9]+c[0-9]+" , name = "AE", tryHard = tryHard)
E = umxModify( AE, regex = "a_r[0-9]+c[0-9]+" , name = "E", tryHard = tryHard)
message("A dominance model is preferred, set dzCr = 0.25")
}else{
CE = umxModify(ACE, regex = "a_r[0-9]+c[0-9]+" , name = "CE", tryHard = tryHard)
AE = umxModify(ACE, regex = "c_r[0-9]+c[0-9]+" , name = "AE", tryHard = tryHard)
E = umxModify( AE, regex = "a_r[0-9]+c[0-9]+" , name = "E", tryHard = tryHard)
}
}else if(model$top$dzCr$values == .25){
if(model$name=="ACE"){
message("You gave me an ADE model, but it was called 'ACE'. I have renamed it ADE for the purposes of clarity in model reduction.")
model = mxRename(model, newname = "ADE", oldname = "ACE")
} else {
message("You gave me an ADE model.")
}
ADE = model
ACE = umxModify(ADE, 'dzCr_r1c1', value = 1, name = "ACE", tryHard = tryHard)
AE = umxModify(ADE, regex = "c_r[0-9]+c[0-9]+" , name = "AE", tryHard = tryHard)
E = umxModify( AE, regex = "a_r[0-9]+c[0-9]+" , name = "E", tryHard = tryHard)
if(-2*logLik(ADE) > -2*logLik(ACE)){
CE = umxModify(ACE, regex = "a_r[0-9]+c[0-9]+" , name = "CE", tryHard=tryHard)
message("An ACE model is preferred, set dzCr = 1.0")
}else{
CE = umxModify(ADE, regex = "a_r[0-9]+c[0-9]+" , name = "DE", tryHard=tryHard)
}
}else{
stop(model$top$dzCr$values, " is an odd number for dzCr, isn't it? I was expecting 1 (C) or .25 (D)",
"\nPerhaps you're a John Loehlin (requiescat in pace :-( ) fan, and are doing an assortative mating test? e-mail me to get this added here.")
# TODO umxReduceACE handle odd values of dzCr as assortative mating etc.?
bestModel = model
}
# = Show fit table =
tmp = data.frame(matrix(nrow=5,ncol=4))
names(tmp) = c("a" , "c" , "e" , "d")
tmp[1,] = c(ACE$top$a_std$result, ACE$top$c_std$result, ACE$top$e_std$result, NA)
tmp[2,] = c(ADE$top$a_std$result, NA , ADE$top$e_std$result, ADE$top$c_std$result)
tmp[3,] = c( NA , CE$top$c_std$result, CE$top$e_std$result, NA)
tmp[4,] = c( AE$top$a_std$result, NA , AE$top$e_std$result, NA)
tmp[5,] = c( NA , NA , E$top$e_std$result, NA)
biggles = umxCompare(ACE, c(ADE, CE, AE, E), all = TRUE, report = report, silent=TRUE)
tmp2 = cbind(biggles[, 1, drop = FALSE], tmp, biggles[, 2:dim(biggles)[2] ] )
umx_print(tmp2, digits = digits, report = report)
whichBest = which.min(AIC(ACE, ADE, CE, AE)[,"AIC"])[1]
bestModel = list(ACE, ADE, CE, AE)[[whichBest]]
message("Among ACE, ADE, CE, and AE models ", omxQuotes(bestModel$name), " fit best according to AIC.\n")
# Probabilities according to AIC MuMIn::Weights (Wagenmakers et al https://pubmed.ncbi.nlm.nih.gov/15117008/ )
if(testD){
aic.weights = round(Weights(AIC(ACE, ADE, CE, AE)[,"AIC"]), 2)
aic.names = namez(c(ACE, ADE, CE, AE))
} else {
aic.weights = round(Weights(AIC(ACE, CE, AE)[,"AIC"]), 2)
aic.names = namez(c(ACE, CE, AE))
}
message("Conditional AIC probability {Wagenmakers, 2004, 192-196} using MuMIn::Weights(AIC()) indicated relative model support for ",
omxQuotes(aic.names), " respectively as: ", omxQuotes(aic.weights))
message(paste0(aic.names," (", aic.weights, "%)"))
if(intervals){
bestModel = mxRun(bestModel, intervals = intervals)
}
umx_set_auto_plot(oldAutoPlot, silent = TRUE)
umx_set_silent(oldSilent)
umx_set_auto_run(autoRun = oldAutoRun)
invisible(bestModel)
}
#' @export
umxReduce.MxModelACE <- umxReduceACE
#' Get residuals from an MxModel
#'
#' Return the [residuals()] from an OpenMx RAM model. You can format these (with digits), and suppress small values.
#'
#' @rdname residuals.MxModel
#' @param object An fitted [OpenMx::mxModel()] from which to get residuals
#' @param digits round to how many digits (default = 2)
#' @param suppress smallest deviation to print out (default = NULL = show all)
#' @param reorder optionally reorder the variables in the residuals matrix to show patterns
#' @param ... Optional parameters
#' @return - matrix of residuals
#' @export
#' @family Reporting functions
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#' manifests = names(demoOneFactor)
#'
#' m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
#' umxPath("G", to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = "G", fixedAt = 1.0)
#' )
#'
#'# ===================================
#'# = Show the residuals of the model =
#'# ===================================
#' residuals(m1)
#' # | |x1 |x2 |x3 |x4 |x5 |
#' # |:--|:----|:-----|:----|:-----|:--|
#' # |x1 |. |. |0.01 |. |. |
#' # |x2 |. |. |0.01 |-0.01 |. |
#' # |x3 |0.01 |0.01 |. |. |. |
#' # |x4 |. |-0.01 |. |. |. |
#' # |x5 |. |. |. |. |. |
#' # [1] "nb: You can zoom in on bad values with, e.g. suppress = .01, which
#' # will hide values smaller than this. Use digits = to round"
#'
#' residuals(m1, digits = 3)
#' residuals(m1, digits = 3, suppress = .005)
#' # residuals are returned as an invisible object you can capture in a variable
#' a = residuals(m1); a
#' }
#'
residuals.MxModel <- function(object, digits = 2, suppress = NULL, reorder=NULL, ...){
umx_check_model(object, type = NULL, hasData = TRUE)
expCov = umxExpCov(object, latents = FALSE)
if(object$data$type == "raw"){
obsCov = umxHetCor(object$data$observed)
} else {
obsCov = object$data$observed
}
resid = cov2cor(obsCov) - cov2cor(expCov)
if(!is.null(reorder)){
resid = umx_reorder(resid, newOrder = reorder, force = TRUE)
}
umx_print(data.frame(resid), digits = digits, zero.print = ".", suppress = suppress)
if(is.null(suppress)){
print("nb: Zoom in on bad values with, e.g. suppress = .1, which will hide smaller values. Use digits = to round")
}
invisible(resid)
}
# define generic loadings...
#' loadings
#' Generic loadings function to extract factor loadings from exploratory or confirmatory
#' factor analyses.
#'
#' See \code{\link[umx]{loadings.MxModel}} to access the loadings of OpenMx EFA models.
#'
#' Base \code{\link[stats]{loadings}} handles [factanal()] objects.
#'
#' @param x an object from which to get loadings
#' @param ... additional parameters
#' @return - matrix of loadings
#' @export
#' @family Reporting functions
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
loadings <- function(x, ...) UseMethod("loadings")
#' @export
loadings.default <- function(x, ...) stats::loadings(x, ...)
# TODO: alternative approach would be to use setGeneric("loadings")
#' Extract factor loadings from an EFA (factor analysis).
#'
#' loadings extracts the factor loadings from an EFA (factor analysis) model.
#' It behaves equivalently to stats::loadings, returning the loadings from an
#' EFA (factor analysis). However it does not store the rotation matrix.
#'
#' @param x A RAM model from which to get loadings.
#' @param ... Other parameters (currently unused)
#' @return - loadings matrix
#' @export
#' @family Reporting functions
#' @seealso - [factanal()], [loadings()]
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @md
#' @examples
#' \dontrun{
#' myVars = c("mpg", "disp", "hp", "wt", "qsec")
#' m1 = umxEFA(name = "test", factors = 2, data = mtcars[, myVars])
#' loadings(m1)
#' }
#'
loadings.MxModel <- function(x, ...) {
x$A$values[x@manifestVars, x@latentVars, drop = FALSE]
}
#' Get confidence intervals from a umx model
#'
#' Implements confidence interval function for umx models.
#'
#' Note: By default, requesting new CIs wipes the existing ones.
#' To keep these, set wipeExistingRequests = FALSE.
#'
#' Because CIs can take time to run, by default only already-computed CIs will be reported. To run new CIs, set run = TRUE .
#'
#' @details *Note*: OpenMx defines a `confint` function which will return SE-based CIs.
#'
#' If `parm` is empty, and `run = FALSE`, a message will alert you to set `run = TRUE`.
#'
#' @param object An [OpenMx::mxModel()], possibly already containing [OpenMx::mxCI()]s that have been [OpenMx::mxRun()] with intervals = TRUE))
#' @param parm Which parameters to get confidence intervals for. Can be "existing", "all", or one or more parameter names.
#' @param level The confidence level required (default = .95)
#' @param run Whether to run the model (defaults to FALSE)
#' @param wipeExistingRequests Whether to remove existing CIs when adding new ones (ignored if parm = 'existing').
#' @param optimizer For difficult CIs, trying other optimizers can help!
#' @param showErrorCodes (default = FALSE)
#' @export
#' @return - [OpenMx::mxModel()]
#' @family Reporting functions
#' @seealso - [stats::confint()], [OpenMx::mxSE()], [umxCI()], [OpenMx::mxCI()]
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#'
#' manifests = names(demoOneFactor)
#' m1 = umxRAM("OneFactor", data = demoOneFactor, type = "cov",
#' umxPath(from = "G", to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = "G", fixedAt = 1)
#' )
#'
#' m1 = umxConfint(m1, run = TRUE) # There are no existing CI requests...
#'
#' # Add a CI request for "G_to_x1", run, and report. Save with this CI computed
#' m2 = umxConfint(m1, parm = "G_to_x1", run = TRUE)
#'
#' # Just print out any existing CIs
#' umxConfint(m2)
#'
#' # CI requests added for free matrix parameters. User prompted to set run = TRUE
#' m3 = umxConfint(m1, "all")
#'
#' # Run the requested CIs
#' m3 = umxConfint(m3, run = TRUE)
#'
#' # Run CIs for free one-headed (asymmetric) paths in RAM model.
#' # note: Deletes other existing requests,
#' tmp = umxConfint(m1, parm = "A", run = TRUE)
#'
#' # Wipe existing CIs, add G_to_x1
#' tmp = umxConfint(m1, parm = "G_to_x1", run = TRUE, wipeExistingRequests = TRUE)
#'
#' # For some twin models, a "smart" mode is implemented
#' # note: only implemented for umxCP so far
#' m2 = umxConfint(m1, "smart")
#' }
#'
umxConfint <- function(object, parm = c("existing", "all", "or one or more labels", "smart"), wipeExistingRequests = TRUE, level = 0.95, run = FALSE, showErrorCodes = FALSE, optimizer= c("SLSQP", "NPSOL", "CSOLNP", "current")) {
optimizer = match.arg(optimizer)
if(optimizer == "current"){
# optimizer is sent to omxRunCI so set "current" to the actual name
optimizer = umx_set_optimizer(silent = TRUE)
}
parm = xmu_match.arg(parm, c("existing", "all", "or one or more labels", "smart"), check = FALSE)
# Upgrade "all" to "smart" for CP
if(class(object)[[1]] == "MxModelCP" && parm == "all"){
parm = "smart"
}
# 1. remove existing CIs if requested to
if(wipeExistingRequests && (parm != "existing")){
if(length(object$intervals)){
object = mxModel(object, remove = TRUE, object$intervals)
message("Removed existing CIs")
}
}
# 2. Add CIs if requested
if (length(parm) >1){
# Add requested CIs to model
# TODO umxConfint: Check that these are valid and not duplicates
object = mxModel(object, mxCI(parm, interval = level))
} else if (parm == "all") {
CIs_to_set = names(omxGetParameters(object, free = TRUE))
object = mxModel(object, mxCI(CIs_to_set, interval = level))
} else if (parm == "smart"){
if(class(object)[[1]] == "MxModelCP"){
# Add individual smart (only free cell) mxCI requests
# For CP model, these are the free cells in
# top.as_std, top.cs_std, top.es_std
# object = m1
this = object$top$as$free
template = umxMatrix("A", "Full", dim(this)[1], dim(this)[2])$labels
patt = "^.*_r([0-9]+)c([0-9]+)$"
as_free = gsub(pattern = patt, replacement= "top.as_std[\\1,\\2]", template)[which(object$top$as$free)]
cs_free = gsub(pattern = patt, replacement= "top.cs_std[\\1,\\2]", template)[which(object$top$cs$free)]
es_free = gsub(pattern = patt, replacement= "top.es_std[\\1,\\2]", template)[which(object$top$es$free)]
# Get labels for free cells in top.cp_loadings_std
this = object$top$cp_loadings$free
template = umxMatrix("A", "Full", dim(this)[1], dim(this)[2])$labels
cp_loadings_free = gsub(pattern = patt, replacement= "top.cp_loadings_std[\\1,\\2]", template)[which(this)]
# top.a_cp, top.c_cp, top.e_cp
this = object$top$a_cp$free
template = umxMatrix("A", "Full", dim(this)[1], dim(this)[2])$labels
a_cp_free = gsub(pattern = patt, replacement= "top.a_cp[\\1,\\2]", template)[which(object$top$a_cp$free)]
c_cp_free = gsub(pattern = patt, replacement= "top.c_cp[\\1,\\2]", template)[which(object$top$c_cp$free)]
e_cp_free = gsub(pattern = patt, replacement= "top.e_cp[\\1,\\2]", template)[which(object$top$e_cp$free)]
CIs2Add = c(a_cp_free, c_cp_free, e_cp_free, cp_loadings_free, as_free, cs_free, es_free)
umx_msg(CIs2Add)
object = mxModel(object, mxCI(CIs2Add, interval = level))
message("added ", length(CIs2Add), " CIs")
} else {
stop("I only know how to add smart CIs for CP models so far. Sorry")
}
} else if (parm == "existing"){
# nothing to do
} else {
# User requesting 1 new CI
# TODO umxConfint: Check that these are valid and not duplicates
object = mxModel(object, mxCI(parm, interval = level))
}
# 3. Run CIs if requested
if(run) {
# Check there are some in existence
if(!umx_has_CIs(object, "intervals")) {
message("Polite note: This model has no CIs yet. Perhaps you wanted to use parm = 'all' to request CIs on all free parameters? Or list some path labels?")
}else{
# object = mxRun(object, intervals = TRUE)
object = omxRunCI(object, optimizer = optimizer)
}
}
# 4. Report CIs
if(!umx_has_CIs(object, "both")) {
if(run == FALSE){
message("Polite note: Some CIs have been requested but not run. Add ", omxQuotes("run = TRUE"), " to your umxConfint() call to run them.\n",
"To store the model capture it from umxConfint like this:\n",
"m1 = umxConfint(m1, run = TRUE)")
} else if(length(object$intervals)==0){
message("No CIs requested...")
} else{
message("Polite note: hmmm... You wanted it run, but I don't see any computed CIs despite there being ", length(object$intervals), " requested...",
"\nThat's a bug. Please report it to timothy.c.bates@gmail.com")
}
} else {
# model has CIs and they have been run
# 1. Summarize model
model_summary = summary(object, verbose = TRUE)
# 2. Extract CIs and details, and arrange for merging
CIdetail = model_summary$CIdetail
CIdetail = CIdetail[, c("parameter", "value", "side", "diagnostic", "statusCode")]
CIdetail$diagnostic = as.character(CIdetail$diagnostic)
CIdetail$statusCode = as.character(CIdetail$statusCode)
CIdetail$diagnostic = namez(CIdetail$diagnostic, pattern = "alpha level not reached" , replacement = "alpha hi")
CIdetail$statusCode = namez(CIdetail$statusCode, pattern = "infeasible non-linear constraint", replacement = "constrained")
CIdetail$statusCode = namez(CIdetail$statusCode, pattern = "iteration limit/blue" , replacement = "blue")
CIs = model_summary$CI
CIs$parameter = row.names(CIs)
row.names(CIs) = NULL
CIs = CIs[, c("parameter", "estimate", "lbound", "ubound", "note")]
intersect(names(CIdetail), names(CIs))
tmp = merge(CIs, CIdetail[CIdetail$side == "lower", ], by = "parameter", all.x = TRUE)
tmp = merge(tmp, CIdetail[CIdetail$side == "upper", ], by = "parameter", all.x = TRUE, suffixes = c(".lower",".upper"))
tmp$side.lower = NULL
tmp$side.upper = NULL
# 3. Format CIs
model_CIs = round(CIs[,c("lbound", "estimate", "ubound")], 3)
model_CI_OK = object$output$confidenceIntervalCodes
colnames(model_CI_OK) = c("lbound Code", "ubound Code")
model_CIs = cbind(round(model_CIs, 3), model_CI_OK)
print(model_CIs)
npsolMessages = list(
'1' = 'The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).',
'2' = 'The linear constraints and bounds could not be satisfied. The problem has no feasible solution.',
'3' = 'The nonlinear constraints and bounds could not be satisfied. The problem may have no feasible solution.',
'4' = 'The major iteration limit was reached (Mx status BLUE).',
'5' = 'not used',
'6' = 'The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)',
'7' = 'The function derivatives returned by funcon or funobj appear to be incorrect.',
'8' = 'not used',
'9' = 'An input parameter was invalid')
if(!all(is.null(model_CI_OK)) && (any(model_CI_OK !=0) |any(is.na(model_CI_OK))) && showErrorCodes){
codeList = c(model_CI_OK[,"lbound Code"], model_CI_OK[,"ubound Code"])
relevantCodes = unique(codeList); relevantCodes = relevantCodes[relevantCodes !=0]
for(i in relevantCodes) {
print(paste0(i, ": ", npsolMessages[i][[1]]))
}
}
}
invisible(object)
}
# 1776(liberty), 1789(liberty+terror), 1815 (liberty+inequality)
#' Add (and, optionally, run) confidence intervals to a structural model.
#'
#' `umxCI` adds [OpenMx::mxCI()] calls for requested (default all) parameters in a model,
#' runs these CIs if necessary, and reports them in a neat summary.
#'
#' @details
#' `umxCI` also reports if any problems were encountered. The codes are standard OpenMx errors and warnings
#' \itemize{
#' \item 1: The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN)
#' \item 2: The linear constraints and bounds could not be satisfied. The problem has no feasible solution.
#' \item 3: The nonlinear constraints and bounds could not be satisfied. The problem may have no feasible solution.
#' \item 4: The major iteration limit was reached (Mx status BLUE).
#' \item 6: The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
#' \item 7: The function derivatives returned by funcon or funobj appear to be incorrect.
#' \item 9: An input parameter was invalid.
#' }
#' If `run = "no"`, the function simply adds the CI requests, but returns the model without running them.
#'
#' @param model The [OpenMx::mxModel()] you wish to report [OpenMx::mxCI()]s on
#' @param which What CIs to add: c("ALL", NA, "list of your making")
#' @param remove = FALSE (if set, removes existing specified CIs from the model)
#' @param run Whether or not to compute the CIs. Valid values = "no" (default), "yes", "if necessary".
#' 'show' means print the intervals if computed, or list their names if not.
#' @param interval The interval for newly added CIs (defaults to 0.95)
#' @param type The type of CI (defaults to "both", options are "lower" and "upper")
#' @param regex Add CIs for labels matching this regular expression (over-rides which)
#' @param showErrorCodes Whether to show errors (default == TRUE)
#' @return - [OpenMx::mxModel()]
#' @family Reporting functions
#' @seealso - [stats::confint()], [umxConfint()], [umxCI()], [umxModify()]
#' @references - <https://github.com/tbates/umx>
#' @export
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#' manifests = names(demoOneFactor)
#'
#' m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
#' umxPath("G", to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = "G", fixedAt = 1)
#' )
#' m1$intervals # none yet - empty list()
#' m1 = umxCI(m1)
#' m1$intervals # $G_to_x1...
#' m1 = umxCI(m1, remove = TRUE) # remove CIs from the model and return it
#' m1$intervals # none again
#'
#' # Add CIs by name
#' parameters(m1, patt="_with_")
#' m1 = umxCI(m1, which = "x1_with_x1")
#' m1 = umxCI(m1, which = c("x1_with_x1", "x2_with_x2"))
#' m1 = umxCI(m1, regex = "x1_with_", run= "yes")
#' # lbound estimate ubound lbound Code ubound Code
#' # x1_with_x1 0.036 0.041 0.047 0 0
#'
#' # ========================
#' # = A twin model example =
#' # ========================
#' data(twinData)
#' mzData = subset(twinData, zygosity == "MZFF")
#' dzData = subset(twinData, zygosity == "DZFF")
#' m1 = umxACE(selDVs = c("bmi1","bmi2"), dzData = dzData, mzData = mzData)
#' umxCI(m1, run = "show") # show what will be requested
#' umxCI(m1, run = "yes") # actually compute the CIs
#' # Don't force update of CIs, but if they were just added, then calculate them
#' umxCI(m1, run = "if necessary")
#' m1 = umxCI(m1, remove = TRUE) # remove them all
#' m1$intervals # none!
#' # Show what parameters are available to get CIs on
#' umxParameters(m1)
#' # Request a CI by label:
#' m1 = umxCI(m1, which = "a_r1c1", run = "yes")
#' }
umxCI <- function(model = NULL, which = c("ALL", NA, "list of your making"), remove = FALSE, run = c("no", "yes", "if necessary", "show"), interval = 0.95, type = c("both", "lower", "upper"), regex = NULL, showErrorCodes = TRUE) {
# Note: OpenMx now overloads confint, returning SE-based intervals.
run = match.arg(run)
which = xmu_match.arg(which, c("ALL", NA, "list of your making"), check = FALSE)
if(remove){
if(!is.null(regex)){
CIs = namez(model$intervals, pattern = regex)
} else if(any(which == "ALL")){
CIs = names(model$intervals)
} else {
CIs = which
}
if(length(names(model$intervals)) > 0){
model = mxModel(model, mxCI(CIs), remove = TRUE)
} else {
message("model has no intervals to remove")
}
invisible(model)
} else {
# Adding CIs
# TODO Avoid duplicating existing CIs
# TODO Add each CI individually
# TODO Break them out into separate models and reassemble if on cluster?
if(any(is.na(which)) && is.null(regex)){
# nothing to add
} else {
if(!is.null(regex)){
CIs = umxGetParameters(model, regex = regex, free = TRUE)
} else if(any(which == "ALL")){
CIs = names(omxGetParameters(model, free = TRUE))
} else {
CIs = which
}
model = mxModel(model, mxCI(CIs, interval = interval, type = type))
}
}
if(run == "yes" | (!umx_has_CIs(model) & run == "if necessary")) {
model = mxRun(model, intervals = TRUE)
} else {
message("Not running CIs, run == ", run)
}
if(run == "show") {
print("CI requests in the model:")
print(names(model$intervals))
}
if(umx_has_CIs(model)){
umxConfint(model, showErrorCodes = showErrorCodes)
message("Table: Computed CIs in model ", model$name)
}
invisible(model)
}
#' Shows a compact, publication-style, summary of umx models
#'
#' @description
#' Report the fit of a OpenMx model or specialized model class (such as ACE, CP etc.)
#' in a compact form suitable for reporting in a journal.
#'
#' See documentation for RAM models summary here: [umxSummary.MxModel()].
#'
#' View documentation on the ACE model subclass here: [umxSummaryACE()].
#'
#' View documentation on the ACEv model subclass here: [umxSummaryACEv()].
#'
#' View documentation on the IP model subclass here: [umxSummaryIP()].
#'
#' View documentation on the CP model subclass here: [umxSummaryCP()].
#'
#' View documentation on the GxE model subclass here: [umxSummaryGxE()].
#'
#' @param model The [OpenMx::mxModel()] whose fit will be reported
#' @param ... Other parameters to control model summary
#' @family Model Summary and Comparison
#' @md
#' @export
umxSummary <- function(model, ...){
UseMethod("umxSummary", model)
}
#' @export
umxSummary.default <- function(model, ...){
stop("umxSummary is not defined for objects of class:", class(model))
}
#' Shows a compact, publication-style, summary of a RAM model
#'
#' Report the fit of a model in a compact form suitable for a journal.
#' It reports parameters in a markdown or html table (optionally standardized), and fit indices
#' RMSEA (an absolute fit index, comparing the model to a perfect model) and CFI and TLI (incremental fit indices comparing model a model with the worst fit).
#'
#' `umxSummary` alerts you when model fit is worse than accepted criterion (TLI >= .95 and RMSEA <= .06; (Hu & Bentler, 1999; Yu, 2002).
#'
#' Note: For some (multi-group) models, you will need to fall back on [summary()]
#'
#' **CIs and Identification**
#' This function uses the standard errors reported by OpenMx to produce the CIs you see in umxSummary
#' These are used to derive confidence intervals based on the formula 95%CI = estimate +/- 1.96*SE)
#'
#' Sometimes SEs appear NA. This may reflect a model which is not identified (see <http://davidakenny.net/cm/identify.htm>).
#' This can include empirical under-identification - for instance two factors
#' that are essentially identical in structure. use [OpenMx::mxCheckIdentification()] to check identification.