# 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.
# devtools::document("~/bin/umx"); devtools::install("~/bin/umx");
# utility naming convention: "umx_" prefix, lowercase, and "_" (not camel case) e.g. xmu_data_swap_a_block()
# ===================================
# = Poems one should know by heart: =
# ===================================
# William Shakespeare
# [Macbeth Tomorrow and tomorrow soliloquy](https://www.poetryfoundation.org/poems/56964/speech-tomorrow-and-tomorrow-and-tomorrow)
# [Othello To be or not to be](https://www.poetryfoundation.org/poems/56965/speech-to-be-or-not-to-be-that-is-the-question)
# [The Merchant of Venice](https://www.goodreads.com/work/quotes/2682703-the-merchant-of-venice)
# * "How far that little candle throws his beams! So shines a good deed in a weary world."
# * The quality of mercy is not strained.
# * "One half of me is yours, the other half is yours,
# Mine own, I would say; but if mine, then yours,
# And so all yours."
# * If to do were as easy as to know what were good to do,
# chapels had been churches,
# and poor men's cottages princes' palaces.
# * This above all: to thine own self be true,
# # PERCY BYSSHE SHELLEY
# [Ozymandias](https://www.poetryfoundation.org/poems/46565/ozymandias)
#
# [Yeats](https://en.wikipedia.org/wiki/W._B._Yeats)
# * [The Second Coming](https://en.wikipedia.org/wiki/The_Second_Coming_(poem))
# * [The Wild Swans at Coole](https://en.wikipedia.org/wiki/The_Wild_Swans_at_Coole_(poem))
#
# * [Invictus](https://en.wikipedia.org/wiki/Invictus)
# Brevia
# * [Abou ben Adhem](https://www.poetryfoundation.org/poems/44433/abou-ben-adhem)
# * [Odi et amo](https://en.wikipedia.org/wiki/Catullus_85)
# I hate and I love.
# Why?, perhaps you ask.
# I know not, but I feel it so,
# and it tortures me.
#' Justified P/E Ratio
#'
#' @description
#' Compute the Justified P/E of a stock.
#' Justified P/E = ( (DPS / EPS) * (1 + g)) / (k – g)
#' DPS is the dividend per share, EPS is the earnings per share,
#' g is the sustainable growth rate, and k is the required rate of return.
#' @param Dividend The dividend.
#' @param EPS The Earnings per Share.
#' @param growthRate The growth rate.
#' @param discountRate Your chosen discount rate.
#' @param basePE The base PE.
#' @param yrs Years.
#' @return - A PE that is justified for this stock.
#' @export
#' @family Miscellaneous Functions
#' @seealso - [fin_interest()], [fin_percent()], [fin_NI()]
#' @md
#' @examples
#' # fin_JustifiedPE(Dividend= .8, EPS = 2, growthRate = .06, discountRate = .1)
#'
fin_JustifiedPE <- function(Dividend= .02, EPS = 1, growthRate = .08, discountRate = .12, basePE= 20, yrs=10) {
paste0("Based on growth (", growthRate*100, "% expected growth for ", yrs, " years and a base P/E of ",
basePE, "), the justified P/E would be: ", (growthRate * yrs) + basePE )
# ((0.4 * 2) * (1 + 0.06)) / (0.1 - 0.06)
# ((Dividend/EPS) * (1 + growthRate)) / (k-growthRate)
# Justified P/E Ratio = 16.8
}
#' Open a ticker in yahoo finance.
#'
#' @description
#' Open a stock ticker, currently in yahoo finance
#'
#' @param ticker A stock symbol to look up, e.g., "OXY"
#' @return - Open a ticker in a finance site online
#' @export
#' @family Miscellaneous Functions
#' @seealso - [fin_interest()], [fin_percent()], [fin_NI()]
#' @md
#' @examples
#' # Open $INTC in yahoo finance.
#' \dontrun{
#' fin_ticker("INTC")
#' }
fin_ticker <- function(ticker= "INTC") {
url =paste0("https://finance.yahoo.com/quote/", ticker)
# https://www.google.com/finance/quote/SWBI:NASDAQ
browseURL(url, browser = getOption("browser"))
}
#' Add a fit statistic to a ggplot
#'
#' @description
#' Add a fit statistic to a ggplot
#'
#' @param model a statistical model which contains a fit measure.
#' @param effect optional hard coded fit/effect.
#' @param xloc x location of R.
#' @param yloc y location of R.
#' @return - plot
#' @export
#' @family Plotting functions
#' @seealso - [umxPlot()], [umxPlotFun()]
#' @md
#' @examples
#' \dontrun{
#' m1 = lm(mpg ~ wt, data = mtcars)
#' p = ggplot2::ggplot(data = mtcars, aes(x = wt, y = mpg))+ geom_point() +geom_smooth()+
#' ggAddR(m1, effect = NA, xloc=2, yloc= 10); p
#' }
ggAddR <- function(model, effect = NA, xloc=8, yloc= 10) {
if(is.na(effect)){
r2 = round(summary(model)$r.squared, 3)
lab = bquote(R^2 == .(r2))
return(cowplot::draw_label(lab, x = xloc, y = yloc, fontfamily = "Times", size = 12))
} else {
return(cowplot::draw_label(effect, x = xloc, y = yloc, fontfamily = "Times", size = 12))
}
}
# #' Easily use the Box-Cox transform
# #'
# #' @description
# #' Applies Box-Cox to the input columns. Box-Cox finds (x ^ (lambda - 1)) / lambda
# #'
# #' | \eqn{\lambda} | Transformation |
# #' |:--------------|:----------------|
# #' | -2 | \eqn{1/x^2} |
# #' | -1 | 1/x |
# #' | -.5 | 1/sqrt(x) |
# #' | 0 | log(x) |
# #' | 0.5 | sqrt(x) |
# #'
# #' This is a crucial processing step. You should graph (`plotit =TRUE`) and examine.
# #'
# #' Improvements would allow taking more than one column and (optionally)#
# #' snapping to common values.
# #' @details
# #'
# #' @param x A variable to transform.
# #' @param verbose Plot the likelihood of lambda and show nearest "common" transform.
# #' @return - The Box-Cox transformed variable.
# #' @export
# #' @family Data Functions
# #' @seealso - [MASS::boxcox()]
# #' @references - Box, G. E. P. and Cox, D. R. (1964) An analysis of transformations (with discussion).
# #' *Journal of the Royal Statistical Society B*, **26**, 211–252. <https://www.jstor.org/stable/2984418>
# #' @md
# #' @examples
# #' tmp = log(1:100)
# #' hist(tmp)
# #' hist(umx_boxcox(tmp))
# umx_boxcox <- function(x, verbose = TRUE, shiftMin=FALSE) {
# if(min(x)<0){
# con.expr
# } else {
# alt.expr
# }
# xbox = MASS::boxcox(lm(x ~ 1), plotit = verbose)
# lambda = xbox$x[which.max(xbox$y)]
# if(lambda == 0){
# newX = log(x)
# } else {
# newX = ((x^lambda) - 1) / lambda
# }
# if(verbose){
# nearest = c("1/x^2", "1/x", "1/sqrt(x)", "log(x)", "sqrt(x)")[which.min(abs(lambda - c(-2, -1, -.5, 0, 0.5)))]
# cat("lambda = ", lambda, "\nTransform was: (x ^ ", lambda, " - 1) / ", lambda, "\nThe nearest common transform is:", omxQuotes(nearest))
# }
# return(newX)
# }
#' load libraries
#'
#' @description
#' `libs` allows loading multiple libraries in one call
#'
#' @param ... library names as strings, e.g. "pwr"
#' @param force.update install.package even if present (to get new version) FALSE
#' @return - nothing.
#' @export
#' @family Miscellaneous Utility Functions
#' @seealso - [library()], [install.packages()], [remove.packages()]
#' @md
#' @examples
#' \dontrun{
#' libs("umx", "OpenMx", "car")
#' libs("umx", c("OpenMx", "car"))
#' remove.packages()
#' }
libs <- function(... , force.update = FALSE) {
dot.items = list(...) # grab all the dot items
dot.items = unlist(dot.items) # In case any dot items are lists
for (pack in dot.items) {
result = tryCatch({
if(force.update){
install.packages(pack)
}
library(pack, character.only = TRUE)
}, warning = function(warn) {
umx_msg("Who's, Z?")
}, error = function(err) {
cat(paste0("I'll try and install.packages(", omxQuotes(pack), ") for you"))
install.packages(pack)
library(pack, character.only = TRUE)
}, finally={
})
}
}
#' Succinctly select complete rows from a dataframe
#'
#' @description
#' Succinctly select complete rows from a dataframe.
#'
#' @param df an [data.frame()] to select on
#' @param rows Rows to keep (optional, incomplete rows still discarded)
#' @param cols Cols to keep
#' @param drop Whether to return a vector when only 1 column is selected (default TRUE)
#' @return - Complete rows and (optionally) selected columns
#' @export
#' @family Data Functions
#' @md
#' @examples
#' tmp = mtcars
#' tmp[2,1] = NA
#' noNAs(tmp, cols="mpg")
#' noNAs(tmp, cols="mpg", drop = FALSE)
#' noNAs(tmp) # no Mazda RX4 Wag
#'
noNAs <- function(df, rows = NULL, cols = NULL, drop = TRUE) {
if(is.null(rows)){
if(is.null(cols)){
# every complete row, all cols
df[complete.cases(df), , drop = drop]
}else{
# every complete row, selected cols
df[complete.cases(df[, cols]), cols, drop = drop]
}
} else {
# remove discluded rows
df = df[rows, ]
if(is.null(cols)){
# every remaining complete row, all cols
df[complete.cases(df), , drop = drop]
}else{
# every remaining complete row, selected cols
df[complete.cases(df[, cols]), cols, drop = drop]
}
}
}
#' Return names of models found within a model
#'
#' @description
#' `umxModelNames` returns the names of each model contained in the model provided to it
#' (optionally excluding the out model itself).
#'
#' @param model an [OpenMx::mxModel()] to search for model names.
#' @param includeOuterModelName FALSE
#' @return - All models names
#' @export
#' @family Miscellaneous Utility Functions
#' @seealso - [OpenMx::mxRename()], [umxSuperModel()]
#' @md
#' @examples
#' \dontrun{
#' data(GFF)
#' mzData = subset(GFF, zyg_2grp == "MZ")
#' dzData = subset(GFF, zyg_2grp == "DZ")
#' selDVs = c("gff", "fc", "qol")
#' m1 = umxCP(selDVs= selDVs, nFac= 1, dzData= dzData, mzData= mzData, sep= "_T", autoRun= TRUE)
#' m2 = mxRename(m1, "model2")
#' umxModelNames(m1) # "top" "MZ" "DZ"
#' umxModelNames(m2) # "top" "MZ" "DZ"
#'
#' super = umxSuperModel("myModel", m1, m2, autoRun = TRUE)
#' umxModelNames(super)
#'
#' plot(super$CP1fac)
#' }
umxModelNames <- function(model, includeOuterModelName = FALSE) {
nameList = c()
if(includeOuterModelName){
nameList = c(nameList, model$name)
}
if(length(model$submodels) > 0){
# get the names of the submodels
subNames = names(model$submodels)
nameList = c(nameList, subNames)
# iterate overthem finding sub-sub-models
for (thisSubModel in subNames) {
newNames = names(eval(parse(text = paste0("model$", thisSubModel, "$submodels"))))
nameList = c(nameList, newNames)
}
}
return(nameList)
}
# =========================================================
# = Obscure enough to be in xmu internal not for end user =
# =========================================================
#' Rename a umxMatrix (even in a model)
#'
#' @description
#' Rename a [umxMatrix()], including updating its labels to match the new name.
#'
#' @param x A model or matrix
#' @param matrixName Name of the matrix
#' @param name The new name
#' @return - updated matrix or model with updated matrix in it.
#' @export
#' @family xmu internal not for end user
#' @md
#' @examples
#' \dontrun{
#' data(twinData) # ?twinData from Australian twins.
#' twinData[, c("ht1", "ht2")] = twinData[, c("ht1", "ht2")] * 10
#' mzData = twinData[twinData$zygosity %in% "MZFF", ]
#' dzData = twinData[twinData$zygosity %in% "DZFF", ]
#' m1 = umxACE(selDVs= "ht", sep= "", dzData= dzData, mzData= mzData, autoRun= FALSE)
#' tmp = umxRenameMatrix(m1$top, matrixName = "a", name="hello")
#' umx_check(tmp$hello$labels == "hello_r1c1") # new is there
#' umx_check(is.null(tmp$a)) # old is gone
#' }
#'
umxRenameMatrix <- function(x, matrixName, name) {
if(umx_is_MxModel(x)){
# 1. Grab a copy of the matrix
tmp = x[[matrixName]]
umx_check(!is.null(tmp), "stop", paste0("matrix ", matrixName, " not found (calling umxRenameMatrix in model ", omxQuotes(x$name), ")"))
# 2. Update the new copy
tmp$name = name
tmp$labels = namez(tmp$labels, pattern = paste0(matrixName), replacement = paste0(name))
# 2. Delete the old one from the model
x = mxModel(x, matrixName, remove= TRUE)
# 4. Add back to top
x = mxModel(x, tmp)
return(x)
} else {
stop("Haven't implemented umxRenameMatrix for matrices")
}
}
# ==============================
# = Get and set OpenMx options =
# ==============================
#' Display umx options
#'
#' Show the umx options. Useful for beginners to discover, or people like me to remember :-)
#'
#' @return - message
#' @export
#' @family Get and set
#' @examples
#' umx_get_options()
umx_get_options <- function() {
umx_set_auto_plot()
umx_set_plot_format()
umx_set_plot_file_suffix()
umx_set_plot_use_hrbrthemes()
umx_set_table_format()
umx_set_optimizer()
umx_set_auto_run()
umx_set_condensed_slots()
message(umx_set_cores(silent = TRUE), " cores will be used")
}
#' Set theme system to use for plots.
#'
#' Set output file suffix (default = "gv", alternative is "dot"). If you call this with no
#' value, it will return the current setting. If you call it with TRUE, it toggles the setting.
#'
#' @param umx.plot.use_hrbrthemes whether to them plots with hrbrthemes (if empty returns the current value)
#' @param silent If TRUE, no message will be printed.
#' @return - Current setting
#' @export
#' @family Get and set
#' @md
#' @examples
#' umx_set_plot_use_hrbrthemes() # print current state
#' old = umx_set_plot_use_hrbrthemes(silent = TRUE) # store current value
#' umx_set_plot_use_hrbrthemes(TRUE)
#' umx_set_plot_use_hrbrthemes(old) # reinstate
umx_set_plot_use_hrbrthemes <- function(umx.plot.use_hrbrthemes = NULL, silent = FALSE) {
if(is.null(umx.plot.use_hrbrthemes)) {
if(!silent){
message("Currently option to use hrbrthemes for plots is",
omxQuotes(getOption("umx.plot.use_hrbrthemes")),
". Valid options are TRUE or FALSE"
)
}
invisible(getOption("umx.plot.use_hrbrthemes"))
} else {
umx_check(umx.plot.use_hrbrthemes %in% c(TRUE, FALSE), "stop", "valid options are TRUE or FALSE)")
options("umx.plot.use_hrbrthemes" = umx.plot.use_hrbrthemes)
}
}
#' Set output suffix used in umx SEM diagram files saved to disk.
#'
#' `umx` SEM diagram files can have a suffix of "gv" (default) or "dot".
#' Interrogate the setting by calling with no value: it will return the current setting.
#' To change the setting call with "gv" or "dot". Or use TRUE to toggle the setting.
#'
#' @param umx.plot.suffix The suffix for plot files (if empty current value is returned). "TRUE", toggles setting.
#' @param silent If TRUE, no message will be printed.
#' @return - Current setting
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' umx_set_plot_file_suffix() # print current state
#' old = umx_set_plot_file_suffix(silent = TRUE) # store current value
#' umx_set_plot_file_suffix("dot")
#' umx_set_plot_file_suffix("gv")
#' umx_set_plot_file_suffix(old) # reinstate
umx_set_plot_file_suffix <- function(umx.plot.suffix = NULL, silent = FALSE) {
if(is.null(umx.plot.suffix)) {
if(!silent){
message("Current format is",
omxQuotes(getOption("umx.plot.suffix")),
". Valid options are 'gv' or 'dot'. Use TRUE to toggle"
)
}
invisible(getOption("umx.plot.suffix"))
} else {
if(umx.plot.suffix == TRUE){
# if T then toggle
if(getOption("umx.plot.suffix") == "gv"){
umx.plot.suffix = "dot"
} else {
umx.plot.suffix = "gv"
}
} else {
umx_check(umx.plot.suffix %in% c("gv", "dot"), "stop", "valid options are 'gv' or 'dot'. Use TRUE to toggle)")
}
options("umx.plot.suffix" = umx.plot.suffix)
}
}
#' Set output format of plots (structural diagrams) in umx
#'
#' Set output format of plots (default = "[DiagrammeR::DiagrammeR()]", alternatives are graphviz, svg, png, pdf). If you call this with no
#' value, it will return the current setting. If you call it with TRUE, it toggles the setting.
#'
#' @param umx.plot.format format for plots (if empty, returns the current value of umx.plot.format). If "TRUE", then toggles
#' @param silent If TRUE, no message will be printed.
#' @return - Current umx.plot.format setting
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' umx_set_plot_format() # print current state
#' old = umx_set_plot_format(silent = TRUE) # store current value
#' umx_set_plot_format("graphviz")
#' umx_set_plot_format("DiagrammeR")
#' umx_set_plot_format("png")
#' umx_set_plot_format("pdf")
#' umx_set_plot_format(old) # reinstate
umx_set_plot_format <- function(umx.plot.format = NULL, silent = FALSE) {
if(is.null(umx.plot.format)) {
if(!silent){
message("Current format is ",
omxQuotes(getOption("umx.plot.format")),
". Valid options are 'DiagrammeR', 'pdf', 'png', 'svg', or 'graphviz'. TRUE toggles between DiagrammeR and graphviz"
)
}
invisible(getOption("umx.plot.format"))
} else {
if(umx.plot.format == TRUE){
# if T then toggle
if(getOption("umx.plot.format") == "graphviz"){
umx.plot.format = "DiagrammeR"
} else {
umx.plot.format = "graphviz"
}
} else {
umx_check(umx.plot.format %in% c("DiagrammeR", "pdf", "png", 'svg', "graphviz"), "stop", "valid options are 'DiagrammeR', 'pdf', 'png', 'svg', 'graphviz'. TRUE toggles between DiagrammeR and graphviz)")
}
options("umx.plot.format" = umx.plot.format)
}
}
#' Set the symbol for money
#'
#' Set umx_set_dollar_symbol (used in e.g. [fin_interest()]
#'
#' @param umx.dollar.symbol symbol for money calculations.
#' @param silent If TRUE, no message will be printed.
#' @return - Current umx.dollar.symbol
#' @export
#' @family Get and set
#' @examples
#' library(umx)
#' umx_set_dollar_symbol() # show current state
#' old = umx_set_dollar_symbol(silent=TRUE) # store existing value
#' fin_interest(100)
#' umx_set_dollar_symbol(old) # reinstate
umx_set_dollar_symbol <- function(umx.dollar.symbol = NULL, silent = FALSE) {
if(is.null(umx.dollar.symbol)) {
if(!silent){ message("Current format is ", omxQuotes(getOption("umx.dollar.symbol")) ) }
invisible(getOption("umx.dollar.symbol"))
} else {
options("umx.dollar.symbol" = umx.dollar.symbol)
}
}
#' Set the separator
#'
#' Set umx_default_separator (used in CI\[low sep high\] ). Default = ","
#'
#' @param umx_default_separator separator for CIs etc. (if empty, returns the current value)
#' @param silent If TRUE, no message will be printed.
#' @return - Current umx_default_separator
#' @export
#' @family Get and set
#' @examples
#' library(umx)
#' umx_set_separator() # show current state
#' old = umx_set_separator(silent=TRUE) # store existing value
#' umx_set_separator("|")
#' umxAPA(.3, .2)
#' umx_set_separator(old) # reinstate
umx_set_separator <- function( umx_default_separator = NULL, silent = FALSE) {
if(is.null( umx_default_separator)) {
if(!silent){
message("Current separator is", omxQuotes(getOption(" umx_default_separator")) )
}
invisible(getOption(" umx_default_separator"))
} else {
options(" umx_default_separator" = umx_default_separator)
}
} # end umx_set_separator
#' umx_set_table_format
#'
#' Set knitr.table.format default (output style for tables). Legal values are
#' "latex", "html", "markdown", "pandoc", and "rst".
#'
#' @param knitr.table.format format for tables (if empty, returns the current value of knitr.table.format)
#' @param silent If TRUE, no message will be printed.
#' @return - Current knitr.table.format setting
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' umx_set_table_format() # show current state
#' old = umx_set_table_format() # store existing value
#' umx_set_table_format("latex")
#' umx_set_table_format("html")
#' umx_set_table_format("markdown")
#' umx_set_table_format("") # get available options
#' umx_set_table_format(old) # reinstate
umx_set_table_format <- function(knitr.table.format = NULL, silent = FALSE) {
legal = c('latex', 'html', 'pipe', 'simple', 'markdown', 'rst')
if(is.null(knitr.table.format)) {
if(!silent){
message("Current format is", omxQuotes(getOption("knitr.table.format")),
". Valid options are ", omxQuotes(legal)
)
}
invisible(getOption("knitr.table.format"))
} else {
if(!knitr.table.format %in% legal){
message("legal options are ", omxQuotes(legal))
} else {
options("knitr.table.format" = knitr.table.format)
}
}
} # end umx_set_table_format
#' umx_set_auto_plot
#'
#' Set autoPlot default for models like umxACE umxGxE etc.
#'
#' @param autoPlot If TRUE, sets the umx_auto_plot option. Else returns the current value of umx_auto_plot
#' @param silent If TRUE, no message will be printed.
#' @return - Current umx_auto_plot setting
#' @export
#' @family Get and set
#' @return - existing value
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' umx_set_auto_plot() # print current state
#' old = umx_set_auto_plot(silent = TRUE) # store existing value
#' old
#' umx_set_auto_plot(TRUE) # set to on (internally stored as "name")
#' umx_set_auto_plot(FALSE) # set to off (internally stored as NA)
#' umx_set_auto_plot(old) # reinstate
umx_set_auto_plot <- function(autoPlot = NULL, silent = FALSE) {
oldAutoPlot = getOption("umx_auto_plot")
if(is.null(oldAutoPlot)){
# not initialised yet
options("umx_auto_plot" = NA)
oldAutoPlot = NA
}
if(is.na(oldAutoPlot)){
oldAutoPlot = FALSE
} else {
oldAutoPlot = TRUE
}
if(is.null(autoPlot)){
# just asking...
if(!silent){
message("Current plot format is ", omxQuotes(oldAutoPlot), "\n Use TRUE to turn on, FALSE to turn off.")
}
}else{
if(is.na(autoPlot) || autoPlot %in% FALSE){
# Set to NA (which is interpreted as FALSE
options("umx_auto_plot" = NA)
} else if(autoPlot == 'name' || autoPlot){
options("umx_auto_plot" = "name")
}
}
invisible(oldAutoPlot)
}
#' umx_set_data_variance_check
#'
#' Set default for data checking in models like umxACE umxGxE etc.
#'
#' @param minVar Set the threshold at which to warn user about variables with too-small variance. Else returns the current value of umx_minVar
#' @param maxVarRatio Set the option for threshold at which to warn user variances differ too much. Else returns the current value of umx_maxVarRatio
#' @param silent If TRUE, no message will be printed.
#' @return - list of umx_minVar and umx_maxVarRatio settings
#' @export
#' @seealso xmu_check_variance which uses these to check sanity in the variances of a data frame.
#' @family Get and set
#' @examples
#' library(umx)
#' umx_set_data_variance_check() # print current state
#' old = umx_set_data_variance_check(silent = TRUE) # store existing value
#' umx_set_data_variance_check(minVar = .01)
#' umx_set_data_variance_check(maxVarRatio = 500)
#' umx_set_data_variance_check(minVar = old$minVar, maxVarRatio = old$maxVarRatio) # reinstate
umx_set_data_variance_check <- function(minVar = NULL, maxVarRatio = NULL, silent = FALSE) {
if(is.null(minVar)){
minVar = getOption("umx_minVar")
if(!silent){
message("Current threshold for small variance warning in umx functions is ", omxQuotes(getOption("umx_minVar")))
}
}else{
options("umx_minVar" = minVar)
}
if(is.null(maxVarRatio)){
maxVarRatio = getOption("umx_maxVarRatio")
if(!silent){
message("Current threshold for excess ratio among variances warning in umx functions is ", omxQuotes(getOption("umx_maxVarRatio")))
}
}else{
options("umx_maxVarRatio" = maxVarRatio)
}
invisible(list(minVar = minVar, maxVarRatio = maxVarRatio))
}
#' Turn off most console and summary output from umx
#'
#' Running multiple analyses or simulations, it can be handy to turn off the automatic summary,
#' graphing, and printing that umx does to help interactive sessions. `umx_set_silent` does this.
#' Summary and graph output, as well as progress and durable console output will be suppressed.
#'
#' Not every function knows about silent, but most, like [umx::umxRAM()] etc do.
#'
#' @details
#' Under the hood, `umx_set_silent` sets options("umx_silent"). This can be set to either `TRUE` or `FALSE`.
#' If TRUE, then the progress messages from model runs are suppressed. Useful for power simulations etc.
#'
#' @param value Boolean stating if umx Models should run silently (TRUE).
#' @param silent If TRUE, this function itself will just return the state of the option, with no user message.
#' @return - Current silent value
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' old = umx_set_silent() # print & store existing value
#' umx_set_silent(FALSE, silent = TRUE) # set to FALSE
#' umx_set_silent(old) # reinstate
#' umx_set_silent() # print existing value
umx_set_silent <- function(value = NA, silent = FALSE) {
# initialize if needed
oldValue = getOption("umx_silent")
if(is.null(oldValue)) {
oldValue = FALSE
options("umx_silent" = FALSE)
}
if(is.na(value)) {
# get value
if(!silent){
message("Current silent setting is ", omxQuotes(oldValue),".\nValid options are TRUE or FALSE.")
}
} else {
umx_check(value %in% c(TRUE, FALSE), "stop")
options("umx_silent" = value)
}
invisible(oldValue)
}
#' Automatically run models?
#'
#' Set `autoRun` default for models like [umxRAM()], [umxACE()] etc.
#'
#' @param autoRun If TRUE or FALSE, sets the umx_auto_run option. Else returns the current value of umx_auto_run
#' @param silent If TRUE, no message will be printed.
#' @return - Current umx_auto_run setting
#' @export
#' @family Get and set
#' @md
#' @examples
#' library(umx)
#' umx_set_auto_run() # print existing value
#' old = umx_set_auto_run(silent = TRUE) # store existing value
#' umx_set_auto_run(FALSE) # set to FALSE
#' umx_set_auto_run(old) # reinstate
umx_set_auto_run <- function(autoRun = NA, silent = FALSE) {
oldValue = getOption("umx_auto_run")
if(is.null(oldValue)) {
oldValue = FALSE
options("umx_auto_run" = TRUE)
}
if(is.na(autoRun)) {
if(!silent){
message("Current auto-run setting is ", omxQuotes(oldValue),".\nValid options are TRUE or FALSE.")
}
invisible(getOption("umx_auto_run"))
} else {
umx_check(autoRun %in% c(TRUE, FALSE), "stop")
options("umx_auto_run" = autoRun)
}
}
#' umx_set_condensed_slots
#'
#' Sets whether newly-created mxMatrices are to be condensed (set to NULL if not being used) or not.
#'
#' @param state what state (TRUE or FALSE) to set condensed slots (default NA returns current value).
#' @param silent If TRUE, no message will be printed.
#' @return - current value of condensed slots
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' umx_set_condensed_slots() # print
#' old = umx_set_condensed_slots(silent = TRUE) # store the existing state
#' umx_set_condensed_slots(TRUE) # update globally
#' umx_set_condensed_slots(old) # set back
umx_set_condensed_slots <- function(state = NA, silent = FALSE) {
if(is.na(state)){
if(!silent){
message("mxCondenseMatrixSlots is currently: ",
omxQuotes(getOption('mxCondenseMatrixSlots'))
)
}
invisible(getOption('mxCondenseMatrixSlots'))
} else {
if(!is.logical(state)){
stop("mxCondenseMatrixSlots must be TRUE or FALSE you tried ", omxQuotes(state))
}else{
options(mxCondenseMatrixSlots = state)
}
}
}
#' Set options that affect optimization in OpenMx
#'
#' `umx_set_optimization_options` provides access to get and set options affecting optimization.
#'
#' *note*: For `mvnRelEps`, values between .0001 to .01 are conventional. Smaller values slow optimization.
#'
#' @param opt default returns current values of the options listed. Currently
#' "mvnRelEps", "mvnMaxPointsA", and "Parallel diagnostics".
#' @param value If not NULL, the value to set the opt to (can be a list of length(opt))
#' @param silent If TRUE, no message will be printed.
#' @param model A model for which to set the optimizer. Default (NULL) sets the optimizer globally.
#' @return - current values if no value set.
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' # show current value for selected or all options
#' umx_set_optimization_options() # print the existing state(s)
#' umx_set_optimization_options("mvnRelEps")
#' \dontrun{
#' umx_set_optimization_options("mvnRelEps", .01) # update globally
#' umx_set_optimization_options("Parallel diagnostics", value = "Yes")
#' }
umx_set_optimization_options <- function(opt = c("mvnRelEps", "mvnMaxPointsA", "Parallel diagnostics"), value = NULL, model = NULL, silent = FALSE) {
if(is.null(value)){
# print current values for each item in opt
for (this in opt) {
if(is.null(model)){
o = mxOption(NULL, this)
} else {
o = mxOption(model, this)
}
message(paste0("Current ", this , " is: ", omxQuotes(o)))
}
invisible(o)
} else {
# Set options
if(length(opt)!=length(value)){
stop("For safe coding, please match opt and value lengths")
} else {
i = 1
for (this in opt) {
if(is.null(model)){
o = mxOption(NULL, this, value[i])
} else {
o = mxOption(model, this, value[i])
}
i = i + 1
}
}
}
}
#' Set the optimizer in OpenMx
#'
#' `umx_set_optimizer` provides an easy way to get and set the default optimizer.
#'
#' @param opt default (NA) returns current value. Current alternatives are
#' "NPSOL" "SLSQP" and "CSOLNP".
#' @param model A model for which to set the optimizer. Default (NULL) sets the optimizer globally.
#' @param silent If TRUE, no message will be printed.
#' @return - current optimizer if nothing requested to be set.
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' umx_set_optimizer() # print the existing state
#' old = umx_set_optimizer(silent = TRUE) # store the existing state
#' umx_set_optimizer("SLSQP") # update globally
#' umx_set_optimizer(old) # set back
umx_set_optimizer <- function(opt = NA, model = NULL, silent = FALSE) {
if(is.na(opt)){
if(is.null(model)){
o = mxOption(NULL, "Default optimizer")
} else {
o = mxOption(model, "Default optimizer")
}
if(!silent){
quoteOptions = omxQuotes(mxAvailableOptimizers())
message("Current Optimizer is: ", omxQuotes(o), ". Options are: ", quoteOptions)
}
invisible(o)
} else {
if(!opt %in% mxAvailableOptimizers()){
stop("The Optimizer ", omxQuotes(opt), " is not legal. Legal values (from mxAvailableOptimizers() ) are:",
omxQuotes(mxAvailableOptimizers()))
}
if(is.null(model)){
mxOption(NULL, "Default optimizer", opt)
} else {
stop(paste0("'Default optimizer' is a global option and cannot be set on models. just say:\n",
"umx_set_optimizer(", omxQuotes(opt), ")"))
}
}
}
#' umx_set_cores
#'
#' set the number of cores (threads) used by OpenMx
#'
#' @param cores number of cores to use. NA (the default) returns current value. "-1" will set to `omxDetectCores()`.
#' @param model an (optional) model to set. If left NULL, the global option is updated.
#' @param silent If TRUE, no message will be printed.
#' @return - number of cores
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' library(umx)
#' manifests = c("mpg", "disp", "gear")
#' m1 = mxModel("ind", type = "RAM",
#' manifestVars = manifests,
#' mxPath(from = manifests, arrows = 2),
#' mxPath(from = "one", to = manifests),
#' mxData(mtcars[, manifests], type = "raw")
#' )
#' umx_set_cores() # print current value
#' oldCores = umx_set_cores(silent = TRUE) # store existing value
#' umx_set_cores(omxDetectCores()) # set to max
#' umx_set_cores(-1); umx_set_cores() # set to max
#' m1 = umx_set_cores(1, m1) # set m1 usage to 1 core
#' umx_set_cores(model = m1) # show new value for m1
#' umx_set_cores(oldCores) # reinstate old global value
umx_set_cores <- function(cores = NA, model = NULL, silent = FALSE) {
if(is.na(cores)){
n = mxOption(model, "Number of Threads") # get the old value
if(!silent){
message(n, "/", omxDetectCores() )
}
return(n)
} else if(umx_is_MxModel(cores)) {
stop("Call this as umx_set_cores(cores, model), not the other way around")
}else{
if(!is.numeric(cores)){
stop("cores must be a number. You gave me ", cores)
}
umx_check(isTRUE(all.equal(cores, as.integer(cores))), message = paste0("cores must be an integer. You gave me: ", cores))
if(cores > omxDetectCores() ){
message("cores set to maximum available (request (", cores, ") exceeds number possible: ", omxDetectCores() )
cores = omxDetectCores()
} else if (cores < 1){
cores = omxDetectCores()
}
mxOption(model, "Number of Threads", cores)
}
}
#' umx_set_checkpoint
#'
#' Set the checkpoint status for a model or global options
#'
#' @aliases umx_set_checkpoint umx_checkpoint
#' @param interval How many units between checkpoints: Default = 1.
#' A value of zero sets always to 'No' (i.e., do not checkpoint all models during optimization)
#' @param units units to count in: Default unit is 'evaluations' ('minutes' is also legal)
#' @param prefix string prefix to add to all checkpoint filenames (default = "")
#' @param directory a directory, i.e "~/Desktop" (defaults to getwd())
#' @param model (optional) model to set options in (default = NULL)
#' @return - mxModel if provided
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#'
#' \dontrun{
#' umx_set_checkpoint(interval = 1, "evaluations", dir = "~/Desktop/")
#' # Turn off checkpointing with interval = 0
#' umx_set_checkpoint(interval = 0)
#' umx_set_checkpoint(2, "evaluations", prefix="SNP_1")
#' 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 = umx_set_checkpoint(model = m1)
#' m1 = mxRun(m1)
#' umx_checkpoint(0)
#' }
umx_set_checkpoint <- function(interval = 1, units = c("evaluations", "iterations", "minutes"), prefix = "", directory = getwd(), model = NULL) {
if(umx_is_MxModel(interval)){
stop("You passed in a model as the first parameter. You probably want:\n",
"umx_set_checkpoint(model=yourModel)")
}
units = match.arg(units)
if(interval == 0){
always = "No"
} else {
always = "Yes"
}
if(is.null(model)){
# Whether to checkpoint all models during optimization.
mxOption(NULL, "Always Checkpoint" , always)
# The number of units between checkpoint intervals
mxOption(NULL, "Checkpoint Count" , interval)
# The type of units for checkpointing: 'minutes', 'iterations', or 'evaluations'.
mxOption(NULL, "Checkpoint Units" , units)
# The string prefix to add to all checkpoint filenames.
mxOption(NULL, "Checkpoint Prefix" , prefix)
# the directory into which checkpoint files are written.
mxOption(NULL, "Checkpoint Directory", directory)
} else {
model = mxOption(model, "Always Checkpoint" , always)
model = mxOption(model, "Checkpoint Count" , interval)
model = mxOption(model, "Checkpoint Units" , units)
model = mxOption(model, "Checkpoint Prefix" , prefix)
model = mxOption(model, "Checkpoint Directory", directory)
return(model)
}
}
#' @export
umx_checkpoint <- umx_set_checkpoint
#' Get or set checkpointing for a model
#'
#' Get the checkpoint status for a model or global options
#'
#' @param model an optional model to get options from
#' @return None
#' @export
#' @family Get and set
#' @references - <https://tbates.github.io>
#' @md
#' @examples
#' \dontrun{
#' umx_get_checkpoint() # current global default
#' 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)
#' )
#' umx_get_checkpoint(model = m1)
#' }
#'
umx_get_checkpoint <- function(model = NULL) {
message("Always Checkpoint: " , mxOption(model, "Always Checkpoint") )
message("Checkpoint Count: " , mxOption(model, "Checkpoint Count" ) )
message("Checkpoint Units: " , mxOption(model, "Checkpoint Units" ) )
message("Checkpoint Prefix: " , mxOption(model, "Checkpoint Prefix" ) )
message("Checkpoint Directory: ", mxOption(model, "Checkpoint Directory" ) )
}
#' Check if OpenMx is using OpenMP, test cores, and get timings
#'
#' Shows how many cores you are using, and runs a test script so user can check CPU usage.
#'
#' @details
#' Some historical (starting 2017-09-06) speeds on my late 2015 iMac, 3.3 GHz
#' Quad-core i7 desktop and then a quad i7 2018 MacBook Pro
#'
#' \tabular{rllll}{
#' Date \tab Version \tab Cores \tab Time \tab Notes \cr
#' 2021-07-28 \tab 2.19.6.19 (git) \tab 8 \tab 00 min 42.98 sec \tab \eqn{\Delta}:-80 (SLSQP laptop (55 sec under NPSOL)) \cr
#' 2021-07-28 \tab 2.19.6.19 (git) \tab 1 \tab 02 min 03 sec \tab (SLSQP on laptop) \cr
#' 2020-08-09 \tab 2.17.3 (git) \tab 1 \tab 01 min 52 sec \tab (CSOLNP on laptop) \cr
#' 2020-08-09 \tab 2.17.3 (git) \tab 4 \tab 00 min 40.18 sec \tab (CSOLNP on laptop) \cr
#' 2019-06-13 \tab v2.13.2 (git) \tab 1 \tab 01 min, 11 sec \tab (NPSOL) \cr
#' 2019-06-13 \tab v2.13.2 (git) \tab 4 \tab 00 min, 22 sec \tab (NPSOL) \cr
#' 2019-06-13 \tab v2.13.2 (git) \tab 6 \tab 00 min, 21 sec \tab (NPSOL) \cr
#' 2018-10-14 \tab v2.11.5 (CRAN) \tab 4 \tab 00 min, 36 sec \tab \eqn{\Delta}:-39.598) \cr
#' 2018-09-17 \tab v2.11.3 \tab 1 \tab 01 min, 31 sec \tab \cr
#' 2018-09-17 \tab v2.11.3 \tab 4 \tab 00 min, 30.6 sec \tab \eqn{\Delta}: -61.49) \cr
#' 2017-10-16 \tab v2.7.18-9 \tab 1 \tab 01 min, 07.30 sec \tab \cr
#' 2017-10-16 \tab v2.7.18-9 \tab 4 \tab 00 min, 22.63 sec \tab \eqn{\Delta}: -44.68) \cr
#' 2017-10-16 \tab Clang OpenMP \tab 1 \tab 01 min, 08.38 sec \tab \cr
#' 2017-10-16 \tab Clang OpenMP \tab 4 \tab 00 min, 24.89 sec \tab \eqn{\Delta}: -43.49) \cr
#' 2017-09-07 \tab Clang OpenMP \tab 1 \tab 01 min, 12.90 sec \tab \cr
#' 2017-09-07 \tab Clang OpenMP \tab 4 \tab 00 min, 32.20 sec \tab \eqn{\Delta}: -40.70 \cr
#' 2017-09-07 \tab Clang notOpenMP \tab 1 \tab 01 min, 09.90 sec \tab \cr
#' 2017-09-07 \tab TRAVIS \tab 1 \tab 01 min, 06.20 sec \tab \cr
#' 2017-09-07 \tab TRAVIS \tab 4 \tab 00 min, 21.10 sec \tab \eqn{\Delta}: -45.00 \cr
#' }
#'
#' @param nCores How many cores to run (defaults to c(1, max). -1 = all available.
#' @param testScript A user-provided script to run (NULL)
#' @param rowwiseParallel Whether to parallel-ize rows (default) or gradient computation
#' @param nSubjects Number of rows to model (Default = 1000) Reduce for quicker runs.
#' @param optimizer Set optimizer, e.g., "NPSOL")
#' @return None
#' @export
#' @family Test
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' # In 2016 1core took 1 minute
#' umx_check_parallel()
#' }
umx_check_parallel <- function(nCores = c(1, omxDetectCores()), testScript = NULL, rowwiseParallel = TRUE, nSubjects = 1000, optimizer=NULL) {
if(!is.null(optimizer)){
oldOpt = umx_set_optimizer()
umx_set_optimizer(optimizer)
}
if(!is.null(testScript)){
stop("test script not implemented yet - beat on tim to do it!")
}
oldCores = umx_set_cores()
if( (length(nCores) == 1) && (nCores == -1)){
nCores = omxDetectCores()
}
message("You have been using ", oldCores, " of ", omxDetectCores(), " available cores (0 means max - 1)")
message("I will now set cores to ", omxQuotes(nCores), " (they will be reset after) and run a script that hits that many cores if possible.\n",
"Check CPU while it's running and see if R is pegging the processor.")
set.seed(10)
# nSubjects = 1000
numberIndicators = 12
numberFactors = 3
fixedBMatrixF = matrix(c(.4, .2), 2, 1, byrow = TRUE)
randomBMatrixF = matrix(c(.3, .5), 2, 1, byrow = TRUE)
XMatrixF = matrix(rnorm(nSubjects * 2, mean = 0, sd = 1), nSubjects, 2)
UMatrixF = matrix(rnorm(nSubjects * 1, mean = 0, sd = 1), nSubjects, 1)
Z = matrix(rnorm(nSubjects, mean = 0, sd = 1), nrow=nSubjects, ncol = 2)
XMatrix = cbind(XMatrixF, XMatrixF %*% fixedBMatrixF + (XMatrixF*Z) %*% randomBMatrixF + UMatrixF)
BMatrix = matrix(c( 1, .6, .7, .8, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, .5, .6, .7, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, .7, .6, .5), numberFactors, numberIndicators, byrow=TRUE)
UMatrix = matrix(rnorm(nSubjects*numberIndicators, mean=0, sd=1), nSubjects, numberIndicators)
YMatrix = XMatrix %*% BMatrix + UMatrix
dimnames(YMatrix) = list(NULL, paste("X", 1:numberIndicators, sep=""))
latentMultiRegModerated1 = cbind(YMatrix,Z=Z[,1])
latentMultiRegModerated1[,'Z'] = latentMultiRegModerated1[,'Z'] - mean(latentMultiRegModerated1[,'Z'])
numberFactors = 3
numberIndicators = 12
numberModerators = 1
indicators = paste("X", 1:numberIndicators, sep="")
moderators = c("Z")
totalVars = numberIndicators + numberFactors + numberModerators
# Build orthogonal simple structure factor model
latents = paste0("F", 1:numberFactors)
latents1 = latents[1]
indicators1 = indicators[1:4]
latents2 = latents[2]
indicators2 = indicators[5:8]
latents3 = latents[3]
indicators3 = indicators[9:12]
# Create model with both direct and moderated paths
test1 = mxModel("threeLatentWithModerator", type = "RAM",
manifestVars = c(indicators),
latentVars = c(latents, "dummy1"),
umxPath(latents1 , to = indicators1, connect = "all.pairs", values = .2),
umxPath(latents2 , to = indicators2, connect = "all.pairs", values = .2),
umxPath(latents3 , to = indicators3, connect = "all.pairs", values = .2),
umxPath(latents1, to = indicators1[1], fixedAt = 1),
umxPath(latents2, to = indicators2[1], fixedAt = 1),
umxPath(latents3, to = indicators3[1], fixedAt = 1),
umxPath(var = latents , values = .8),
umxPath(var = indicators, values = .8),
umxPath(c("F1", "F2"), to = "F3", values = .2, labels = c("b11", "b12")),
umxPath("F1",to = "F2", values = .1, labels = "cF1F2"),
umxPath(c("F1", "F2"),to = "dummy1", values = .2, labels = c("b21", "b22")),
umxPath("dummy1",to = "F3", free = FALSE, labels = "data.Z"),
umxPath(means = indicators, fixedAt = 0),
umxPath(means = latents, values = .1),
mxData(latentMultiRegModerated1, type = "raw")
)
# set rowwiseParallel
if(packageVersion("OpenMx") >= "2.6.1"){
test1$fitfunction$rowwiseParallel = rowwiseParallel
} else {
message("ignored rowwiseParallel: upgrade to OpenMx 2.6.1 or better to use this")
# ignore: this is not supported by versions before 2.6.1
}
# nCores = 4
n = 1
for (thisCores in nCores) {
if(n == 1){
models = list(test1) # initialize
} else {
models = append(models, test1)
}
n = n + 1
}
n = 1
# run each model
for (thisCores in nCores) {
umx_set_cores(thisCores)
thisModel = mxRename(models[[n]], paste0("nCores_equals_", thisCores))
thisModel = mxRun(thisModel)
# umx_time(thisModel, autoRun= F)
models[[n]] = thisModel
n = n + 1
}
umx_set_cores(oldCores)
if(!is.null(optimizer)){
umx_set_optimizer(oldOpt)
}
invisible(umx_time(models, formatStr = "%M %OS3", autoRun = FALSE))
}
# ======================================
# = Lower-level Model building helpers =
# ======================================
#' umxJiggle
#'
#' umxJiggle takes values in a matrix and jiggles them
#'
#' @param matrixIn an [OpenMx::mxMatrix()] to jiggle the values of
#' @param mean the mean value to add to each value
#' @param sd the sd of the jiggle noise
#' @param dontTouch A value, which, if found, will be left as-is (defaults to 0)
#' @return - [OpenMx::mxMatrix()]
#' @family Advanced Model Building Functions
#' @references - <https://github.com/tbates/umx>
#' @export
#' @md
#' @examples
#' \dontrun{
#' mat1 = umxJiggle(mat1)
#' }
umxJiggle <- function(matrixIn, mean = 0, sd = .1, dontTouch = 0) {
mask = (matrixIn != dontTouch);
newValues = mask;
matrixIn[mask == TRUE] = matrixIn[mask == TRUE] + rnorm(length(mask[mask == TRUE]), mean = mean, sd = sd);
return (matrixIn);
}
# ===============
# = RAM helpers =
# ===============
#' umx_is_exogenous
#'
#' Return a list of all the exogenous variables (variables with no incoming single-arrow path) in a model.
#'
#' @param model an [OpenMx::mxModel()] from which to get exogenous variables
#' @param manifests_only Whether to check only manifests (default = TRUE)
#' @return - list of exogenous variables
#' @export
#' @family Check or test
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#' m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
#' umxPath("g", to = names(demoOneFactor)),
#' umxPath(var = "g", fixedAt = 1),
#' umxPath(var = names(demoOneFactor))
#' )
#' umx_is_exogenous(m1, manifests_only = TRUE)
#' umx_is_exogenous(m1, manifests_only = FALSE)
#'
#' }
umx_is_exogenous <- function(model, manifests_only = TRUE) {
umx_check_model(model, type = "RAM")
checkThese = model@manifestVars
if(!manifests_only){
checkThese = c(checkThese, model@latentVars)
}
if(length(checkThese) < 1){
return(c())
}
exog = c()
n = 1
for (i in checkThese) {
if(!any(model$matrices$A$free[i, ])){
exog[n] = i
n = n + 1
}
}
return(exog)
}
#' List endogenous variables in a model
#'
#' Return a list of all the endogenous variables (variables with at least one incoming single-arrow path) in a model.
#'
#' @param model an [OpenMx::mxModel()] from which to get endogenous variables
#' @param manifests_only Whether to check only manifests (default = TRUE)
#' @return - list of endogenous variables
#' @export
#' @family Check or test
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' data(demoOneFactor)
#' m1 = umxRAM("umx_is_endogenous", data = demoOneFactor, type = "cov",
#' umxPath("g", to = names(demoOneFactor)),
#' umxPath(var = "g", fixedAt = 1),
#' umxPath(var = names(demoOneFactor))
#' )
#' umx_is_endogenous(m1, manifests_only = TRUE)
#' umx_is_endogenous(m1, manifests_only = FALSE)
#'
#' }
umx_is_endogenous <- function(model, manifests_only = TRUE) {
# has_no_incoming_single_arrow
umx_check_model(model, type = "RAM")
checkThese = model@manifestVars
if(!manifests_only){
checkThese = c(checkThese, model@latentVars)
}
if(length(checkThese) < 1){
return(c())
}
endog = c()
n = 1
for (i in checkThese) {
# any free paths in this row?
if(any(model$matrices$A$free[i, ])){
endog[n] = i
n = n + 1
}
}
return(endog)
}
# ====================
# = Parallel Helpers =
# ====================
eddie_AddCIbyNumber <- function(model, labelRegex = "") {
# eddie_AddCIbyNumber(model, labelRegex="[ace][1-9]")
args = commandArgs(trailingOnly=TRUE)
CInumber = as.numeric(args[1]); # get the 1st argument from the cmdline arguments (this is called from a script)
CIlist = umxGetParameters(model ,regex= "[ace][0-9]", verbose= FALSE)
thisCI = CIlist[CInumber]
model = mxModel(model, mxCI(thisCI) )
return (model)
}
# ===================================
# = Ordinal/Threshold Model Helpers =
# ===================================
#' umxFactor
#'
#' A convenient version of [OpenMx::mxFactor()] supporting the common
#' case in which the factor levels are those in the variable.
#'
#' @aliases umx_factor
#' @param x A variable to recode as an mxFactor (see [OpenMx::mxFactor()])
#' @param levels (default NULL). Like [factor()] but UNLIKE [OpenMx::mxFactor()],
#' unique values will be used if levels not specified.
#' @param labels = levels (see [OpenMx::mxFactor()])
#' @param exclude = NA (see [OpenMx::mxFactor()])
#' @param ordered = TRUE By default return an ordered mxFactor
#' @param collapse = FALSE (see [OpenMx::mxFactor()])
#' @param verbose Whether to tell user about such things as coercing to factor
#' @param sep If twin data are being used, the string that separates the base from twin index
#' will try and ensure factor levels same across all twins.
#' @return - [OpenMx::mxFactor()]
#' @export
#' @family Data Functions
#' @seealso - [umxFactanal()], [OpenMx::mxFactor()]
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @md
#' @examples
#' umxFactor(letters)
#' umxFactor(letters, verbose = TRUE) # report coercions
#' umxFactor(letters, ordered = FALSE) # non-ordered factor like factor(x)
#' # Dataframe example:
#' x = umx_factor(mtcars[,c("cyl", "am")], ordered = FALSE); str(x)
#' # =================
#' # = Twin example: =
#' # =================
#' data(twinData)
#' tmp = twinData[, c("bmi1", "bmi2")]
#' tmp$bmi1[tmp$bmi1 <= 22] = 22
#' tmp$bmi2[tmp$bmi2 <= 22] = 22
#' # remember to factor _before_ breaking into MZ and DZ groups
#' x = umxFactor(tmp, sep = ""); str(x)
#' xmu_check_levels_identical(x, "bmi", sep="")
#'
#' # Simple example to check behavior
#' x = round(10 * rnorm(1000, mean = -.2))
#' y = round(5 * rnorm(1000))
#' x[x < 0] = 0; y[y < 0] = 0
#' jnk = umxFactor(x); str(jnk)
#' df = data.frame(x = x, y = y)
#' jnk = umxFactor(df); str(jnk)
umxFactor <- function(x = character(), levels= NULL, labels = levels, exclude = NA, ordered = TRUE, collapse = FALSE, verbose = FALSE, sep = NA){
if(is.data.frame(x)){
# x = tmp; sep = NA; sep = ""; thisName = "bmi"; levels = NA
ncols = ncol(x)
if(!is.na(sep)){
if(!is.null(levels)){
stop("leave levels = NA: I don't handle setting levels within data.frames AND sep. You set them to ", omxQuotes(levels))
}
tmp = umx_explode_twin_names(x, sep = sep)
sep = tmp$sep
baseNames = tmp$baseNames
twinIndexes = tmp$twinIndexes
for (thisName in baseNames) {
theseNames = umx_paste_names(thisName, sep, twinIndexes)
a = x[, theseNames]
allLevels = unique(as.vector(as.matrix(a)))
allLevels = sort(allLevels)
allLevels = allLevels[!is.na(allLevels)] # drop NA if present
# z = umxFactor(x = x[,theseNames], levels = allLevels, ordered = TRUE, verbose = TRUE, collapse=FALSE)
# z = umxFactor(x = x[,theseNames], levels = allLevels, labels = allLevels, ordered = TRUE, verbose = TRUE)
x[, theseNames] = umxFactor(x = x[, theseNames, drop = FALSE], levels = allLevels, labels = allLevels, exclude = exclude, collapse = collapse, ordered = ordered, verbose = verbose)
}
} else {
for (c in 1:ncols) {
x[,c] = umxFactor(x = x[,c], levels = levels, labels = labels, exclude = exclude, collapse = collapse, ordered = ordered, verbose = verbose)
}
}
} else {
if(!is.factor(x)){
if(!is.null(levels)) {
x = factor(x, levels = levels, labels = labels, exclude = exclude, ordered = ordered)
} else {
x = factor(x, exclude = exclude, ordered = ordered)
}
levels = levels(x)
if(verbose){
if(length(levels) > 20){
feedback = paste0(length(levels), " levels:", paste(c(levels[1:10], "..."), collapse = "', '"))
} else {
feedback = paste0("levels:", omxQuotes(levels))
}
message("Your variable was not a factor: I made it into one, with ", feedback)
}
}else{
# Already a factor
if(is.null(levels)) {
levels = levels(x)
} else {
if(!levels(x) == levels){
message("the levels you provided are not those I see in the data")
}
}
}
if(ordered){
x = mxFactor(x = x, levels = levels, labels = levels, exclude = exclude, ordered = ordered, collapse = collapse)
}
}
return(x)
}
#' @export
umx_factor <- umxFactor
# ===========
# = Utility =
# ===========
#' A wrapper to map columns of strings to numeric.
#'
#' If you give one column name, this is changed to numeric, and returned as a **vector**.
#' If you give multiple column names, or don't set cols, each is changed to numeric, and the updated **data.frame** is returned.
#'
#' @param df The df
#' @param cols (optional) list of columns (default = use all)
#' @param mapStrings legal strings which will be mapped in order to numbers.
#' @return - df
#' @export
#' @family Data Functions
#' @md
#' @examples
#' tmp = data.frame(x=letters)
#' umx_strings2numeric(tmp, mapStrings = letters)
#' umx_strings2numeric(tmp, cols= "x", mapStrings = letters)
umx_strings2numeric <- function(df, cols= NA, mapStrings = NULL) {
if(!all(is.na(cols))){
umx_check_names(cols, data = df)
df = df[, cols, drop=FALSE]
}else{
cols = names(df)
df = df[, cols, drop = FALSE]
}
for (thisCol in cols){
# check values
unique_values = unique(df[, thisCol, drop = TRUE])
unique_values = unique_values[!is.na(unique_values)]
if(is.null(mapStrings)){
# use table to find valid strings in some order... (not good, tell the user!)
mapStrings = unique_values
message("You didn't set mapStrings. This is unwise! I found the following responses, and used them in this order:", omxQuotes(mapStrings))
}else{
if(any(!(unique_values %in% mapStrings))){
notFound = unique_values[which(!(unique_values %in% mapStrings))]
stop("Some values in column ", omxQuotes(thisCol), " not in mapStrings, e.g.. :", omxQuotes(notFound))
}
}
# string 2 numeric
tmp = factor(df[, thisCol, drop = TRUE], levels = mapStrings, labels = 1: length(mapStrings))
df[, thisCol] = as.numeric(as.character(tmp))
}
if(length(cols) == 1){
return(df[, cols, drop = TRUE])
} else {
return(df)
}
}
#' A wrapper to make paran easier to use.
#' Just automates applying [complete.cases()]
#'
#' @param df The df (just the relevant columns)
#' @param cols (optional) list of columns (default = use all)
#' @param graph Whether to graph.
#' @param mapStrings optional mapping if cols are strings
#' @return - nothing
#' @export
#' @family Miscellaneous Stats Functions
#' @md
#' @examples
#' library(psych)
#' library(psychTools)
#' data(bfi)
#' umxParan(bfi[, paste0("A", 1:5)])
#' umxParan(bfi, cols= paste0("A", 1:5))
#' # umxParan(bfi, paste0("AB", 1))
umxParan <- function(df, cols= NA, graph = TRUE, mapStrings = NULL) {
if(!all(is.na(cols))){
umx_check_names(cols, data = df)
df = df[, cols]
}
if(!is.null(mapStrings)){
for (thisCol in names(df)){
unique_values = unique(df[, thisCol, drop = TRUE])
unique_values = unique_values[!is.na(unique_values)]
if(any(!(unique_values %in% mapStrings))){
notFound = unique_values[which(!(unique_values %in% mapStrings))]
stop("Some values in column ", omxQuotes(thisCol), " not in mapStrings, e.g.. :", omxQuotes(notFound))
}
tmp = factor(df[, thisCol, drop = TRUE], levels = mapStrings, labels = 1: length(mapStrings))
df[, thisCol] = as.numeric(as.character(tmp))
}
}
paran::paran(df[complete.cases(df), ], graph = graph)
}
#' Score a psychometric scale by summing normal and reversed items.
#'
#' In the presence of NAs, `score= "mean"` and `score = "totals"` both return NA unless na.rm = TRUE.
#' `score = "max"`, ignores NAs no matter what.
#'
#' @description
#' Use this function to generate scores as the appropriate sum of responses to the normal and reversed items in a scale.
#'
#' Items must be named on the pattern `basename + N + suffix`, where `base` is the prefix common to all item (column) names, N is item number in the scale, and suffix an optional trail (like "_T1").
#'
#' `pos` and `rev` are vectors of the item numbers for the normal and reverse-scored item numbers.
#'
#' To reverse items, the function uses `max` and `min` as the lowest and highest possible response scores to compute how to reverse items.
#'
#' *note*: `min` defaults to 1.
#' **TIP**: If you have strings, `umx_score_scale` will work (use `mapStrings = `). *BUT* if you want to make a numeric copy, use `umx_strings2numeric`
#'
#' @param base String common to all item names.
#' @param pos The positive-scored item numbers.
#' @param rev The reverse-scored item numbers.
#' @param min Minimum legal response value (default = 1). Not implemented for values other than 1 so far...
#' @param max Maximum legal response value (also used to compute reversed item values).
#' @param data The data frame
#' @param score Score total (default), proportionCorrect, errors, mean, max, or factor scores
#' @param name The name of the scale to be returned. Defaults to "`base`_score"
#' @param na.rm Whether to delete NAs when computing scores (Default = TRUE) Note: Choice affects mean!
#' @param minManifests How many missing items to tolerate for an individual (when score = factor)
#' @param alpha print Reliability (omega and Cronbach's alpha) (TRUE)
#' @param mapStrings Recoding input like "No"/"Maybe"/"Yes" into numeric values (0,1,2)
#' @param correctAnswer Use when scoring items with one correct response (1/0).
#' @param omegaNfactors Number of factors for the omega reliability (default = 1)
#' @param verbose Whether to print the whole omega output (FALSE)
#' @param digits Rounding for omega etc. (default 2)
#' @param suffix (if dealing with, e.g. "_T1")
#' @return - scores
#' @export
#' @family Data Functions
#' @seealso umx_strings2numeric
#' @references - Revelle, W. (2022) psych: Procedures for Personality and Psychological Research, Northwestern University, Evanston, Illinois, USA, <https://CRAN.R-project.org/package=psych> Version = 2.2.9.
#' * McNeish, D. (2018). Thanks coefficient alpha, we’ll take it from here. *Psychological Methods*, **23**, 412-433. \doi{10.1037/met0000144}.
#' @md
#' @examples
#' library(psych)
#' library(psychTools)
#' data(bfi)
#'
#' # ==============================
#' # = Score Agreeableness totals =
#' # ==============================
#'
#' # Handscore subject 1
#' # A1(R)+A2+A3+A4+A5 = (6+1)-2 +4+3+4+4 = 20
#'
#' tmp = umx_score_scale(base = "A", pos = 2:5, rev = 1, max = 6, data= bfi, name = "A")
#' tmp[1, namez(tmp, "A",ignore.case = FALSE)]
#' # A1 A2 A3 A4 A5 A
#' # 2 4 3 4 4 20
#'
#' # ====================
#' # = Request the mean =
#' # ====================
#' tmp = umx_score_scale(name = "A", base = "A",
#' pos = 2:5, rev = 1, max = 6, data= bfi, score="mean")
#' tmp$A[1] # = 4
#'
#' # ========================
#' # = Request factor score =
#' # ========================
#' \dontrun{
#'tmp = umx_score_scale(name = "A", base = "A", pos = 2:5, rev = 1,
#' max = 6, score = "factor", minManifests = 4, data= bfi)
#' # g
#' # A2 0.6574826
#' # A3 0.7581274
#' # A4 0.4814788
#' # A5 0.6272332
#' # A1 0.3736021
#'
#' # ==================
#' # = Request alpha =
#' # ==================
#'
#' tmp=umx_score_scale(base="A", pos=2:5, rev=1, max=6, data=bfi, alpha=TRUE)
#' # omega t = 0.72
#' }
#'
#' # ==================
#' # = na.rm = TRUE ! =
#' # ==================
#' tmpDF = bfi
#' tmpDF[1, "A1"] = NA
#' tmp = umx_score_scale("A", pos = 2:5, rev = 1, max = 6, data= tmpDF, score="mean")
#' tmp$A_score[1] # 3.75
#'
#' tmp= umx_score_scale("A", pos= 2:5, rev= 1, max = 6, data = tmpDF,
#' score="mean", na.rm=FALSE)
#' tmp$A_score[1] # NA (reject cases with missing items)
#'
#' # ===============
#' # = Score = max =
#' # ===============
#' tmp = umx_score_scale("A", pos = 2:5, rev = 1, max = 6,
#' data = bfi, name = "A", score = "max")
#' tmp$A[1] # Subject 1 max = 5 (reversed) item 1
#'
#' # Default scale name
#' tmp = umx_score_scale("E", pos = 3:5, rev = 1:2, max = 6,
#' data= tmp, score = "mean", na.rm = FALSE)
#' tmp$E_score[1]
#'
#' # Using @BillRevelle's psych package: More diagnostics, including alpha
#' scores= psych::scoreItems(items = bfi, min = 1, max = 6, keys = list(
#' E = c("-E1","-E2", "E3", "E4", "E5"),
#' A = c("-A1", "A2", "A3", "A4", "A5")
#' ))
#' summary(scores)
#' scores$scores[1, ]
#' # E A
#' # 3.8 4.0
#'
#' # Compare output
#' # (note, by default psych::scoreItems replaces NAs with the sample median...)
#' RevelleE = as.numeric(scores$scores[,"E"])
#' RevelleE == tmp[,"E_score"]
#'
#' # =======================
#' # = MapStrings examples =
#' # =======================
#' mapStrings = c(
#' "Very Inaccurate", "Moderately Inaccurate",
#' "Slightly Inaccurate", "Slightly Accurate",
#' "Moderately Accurate", "Very Accurate")
#' bfi$As1 = factor(bfi$A1, levels = 1:6, labels = mapStrings)
#' bfi$As2 = factor(bfi$A2, levels = 1:6, labels = mapStrings)
#' bfi$As3 = factor(bfi$A3, levels = 1:6, labels = mapStrings)
#' bfi$As4 = factor(bfi$A4, levels = 1:6, labels = mapStrings)
#' bfi$As5 = factor(bfi$A5, levels = 1:6, labels = mapStrings)
#' bfi= umx_score_scale(name="A" , base="A", pos=2:5, rev=1, max=6, data=bfi)
#' bfi= umx_score_scale(name="As", base="As", pos=2:5, rev=1, mapStrings = mapStrings, data= bfi)
umx_score_scale <- function(base= NULL, pos = NULL, rev = NULL, min= 1, max = NULL, data= NULL, score = c("total", "proportionCorrect", "errors", "mean", "max", "factor"), name = NULL, na.rm=TRUE, minManifests = NA, alpha = FALSE, mapStrings= NULL, correctAnswer = NULL, omegaNfactors = 1, digits = 2, verbose = FALSE, suffix = "") {
score = match.arg(score)
if(is.null(name)){ name = paste0(base, "_score", suffix) }
oldData = data
umx_check_names(namesNeeded= paste0(base, c(pos, rev), suffix), data=data)
if(!is.null(mapStrings)){
if(!is.null(max)){
# Check min max matches mapStrings
if(!(length(mapStrings) == length(min:max))){
stop(paste0("You set the max and min, but ", min, " to ", max, " must equal the number of map strings: ", length(mapStrings)))
}
}else{
min = 1
max = length(mapStrings)
}
relevantColumns = paste0(base, c(pos, rev), suffix)
for (thisCol in relevantColumns){
unique_values = unique(data[, thisCol, drop = TRUE])
unique_values = unique_values[!is.na(unique_values)]
if(any(!(unique_values %in% mapStrings))){
notFound = unique_values[which(!(unique_values %in% mapStrings))]
stop("Some values in column ", omxQuotes(thisCol), " not in mapStrings, e.g.. :", omxQuotes(notFound))
}
tmp = factor(data[, thisCol, drop = TRUE], levels = mapStrings, labels = min:max)
data[, thisCol] = as.numeric(as.character(tmp))
}
}
mins = umx_apply("min", data[ , paste0(base, c(pos, rev), suffix), drop = FALSE], by = "columns", na.rm=TRUE)
maxs = umx_apply("max", data[ , paste0(base, c(pos, rev), suffix), drop = FALSE], by = "columns", na.rm=TRUE)
if(any(mins < min)){
msg = paste0("Polite warning: the following columns had responses less than the min response you set (", omxQuotes(min), "):",
omxQuotes(names(mins)[(mins<min)]), "e.g.", omxQuotes(unique(mins[which(mins < min)]))
)
umx_msg(msg)
}
if(any(maxs > max)){
msg = paste0("Polite warning: the following columns had responses greater than the max response you set (", omxQuotes(max), "):", omxQuotes(names(max)[(maxs>max)]))
umx_msg(msg)
}
# ==================================
# = Reverse any items needing this =
# ==================================
if(!is.null(rev)){
if(is.null(max)){
maxs = umx_apply("max", data[ , paste0(base, rev, suffix), drop = FALSE], by = "columns", na.rm= TRUE)
message("If there are reverse items, you must set 'max' (the highest possible score for an item) in umx_score_scale (note: min defaults to 1)")
print(table(data[ , paste0(base, rev[1], suffix)] ))
stop("FYI, the max appears to be ", max(maxs))
}
revItems = data[,paste0(base, rev, suffix), drop = FALSE]
revItems = (max + min) - revItems
data[ , paste0(base, rev)] = revItems
}
allColNames = paste0(base, c(pos, rev), suffix)
df = data[ , allColNames, drop = FALSE]
if(!is.null(correctAnswer)){
dfDummy = matrix(nrow = nrow(df), ncol= ncol(df))
for (i in 1:nrow(df)) {
dfDummy[i,] = (df[i, ] == correctAnswer)
}
df = dfDummy
if(verbose){
print(str(df))
}
}
if(alpha){
print(reliability(cov(df, use = "pairwise.complete.obs")))
suppressWarnings({omegaOut = psych::omega(df, nfactors = omegaNfactors)})
if(verbose){
print(omegaOut)
print("\n")
}else{
if(omegaNfactors == 1){
# Omega_h for 1 factor is not meaningful, just omega_t
cat(paste0("\u03C9 t = ", round(omegaOut$omega.tot, digits), "\n"))
} else {
cat(paste0("\u03C9 h = ", round(omegaOut$omega_h, digits), "; \u03C9 t = ", round(omegaOut$omega.tot, digits), "\n"))
}
}
}
if(!is.null(correctAnswer)){
message("\nPolite note: I returned the sum of correct Answers scored 1/0.")
scaleScore = rep(NA, nrow(df))
for (i in 1:nrow(df)) {
scaleScore[i] = sum(df[i, ], na.rm = TRUE)
}
} else if(score == "max"){
scaleScore = rep(NA, nrow(df))
for (i in 1:nrow(df)) {
scaleScore[i] = max(df[i,], na.rm = TRUE)
}
}else if(score == "total"){
if(any(is.na(df))){
message("\nPolite note: You asked for scale totals, but some subjects have missing data: Perhaps use means?...")
}
scaleScore = rowSums(df, na.rm = na.rm)
}else if(score == "errors"){
if(any(is.na(df))){
message("\nPolite note: You asked for errors. Just to let you know, some subjects have NA data: I ignored those.")
}
scaleScore = rowSums(df, na.rm = na.rm)
attempted = rowSums(!is.na(df))
scaleScore = attempted - scaleScore
}else if(score == "proportionCorrect"){
if(any(is.na(df))){
message("\nPolite note: You asked for proportions (scaleScore/attempted). Just to let you know, some subjects have missing data.")
}
attempted = rowSums(!is.na(df))
scaleScore = rowSums(df, na.rm = na.rm)
scaleScore = scaleScore/attempted
}else if(score == "mean"){
scaleScore = rowMeans(df, na.rm = na.rm)
# Replace NaN with NA (generated by mean() when all items are NA)
scaleScore[is.nan(scaleScore)] = NA
}else if(score == "factor"){
tmp = umxEFA(name = "score", factors = "g", data = df, scores= "Regression", minManifests= minManifests)
scaleScore = tmp$g
}
oldData[, name] = scaleScore
return(oldData)
}
#' Get or print the version of umx, along with detail from OpenMx and general system info.
#'
#' @description
#' umxVersion returns the version information for umx, and for OpenMx and R.
#' Essential for bug-reports! This function can also test for a minimum version.
#'
#' @param model Optional to show optimizer in this model
#' @param min Optional minimum version string to test for, e.g. '2.7.0' (Default = NULL).
#' @param verbose = TRUE
#' @param return Which package (umx or OpenMx) to 'return' version info for (Default = umx).
#' @return - [OpenMx::mxModel()]
#' @export
#' @family Miscellaneous Utility Functions
#' @seealso - [packageVersion()], [install.OpenMx()]
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @md
#' @examples
#' x = umxVersion(); x
umxVersion <- function (model = NULL, min = NULL, verbose = TRUE, return = c("umx_vers", "OpenMx_vers")) {
return = match.arg(return)
umx_vers = try(packageVersion("umx"))
if (verbose) {
msg = paste0("umx version: ", umx_vers)
message(msg)
}
if(!is.null(min)){
if(umx_vers >= min){
message("umx version is recent enough")
} else {
stop("umx version is not recent enough to run this script! (min is ", min, "). You have ", umx_vers,
"\n You can run umx_open_CRAN_page('umx') to see the most recent version of umx on CRAN")
}
}
if(!is.null(model) && !umx_is_MxModel(model)){
message("Polite message - you should call umxVersion() with no parameters, or the first parameter should be a model")
model = NULL
}
OpenMx_vers = mxVersion(model = model, verbose = verbose)
if (verbose) {
message('Open the CRAN page for any package with umx_open_CRAN_page()\nYou can update OpenMx with:\ninstall.OpenMx(c("NPSOL", "travis", "CRAN", "open travis build page")')
}
if(return == "umx"){
invisible(umx_vers)
} else {
invisible(OpenMx_vers)
}
}
#' Open the CRAN page for a package
#'
#' On MacOS, this function opens the CRAN page for a package.
#' Useful for looking up documentation, checking you have an
#' up-to-date version, showing the package to people etc.
#' @param package An \R package name.
#' @param inst Install and load if not already installed?
#' @return None
#' @export
#' @family Miscellaneous Utility Functions
#' @md
#' @examples
#' \dontrun{
#' umx_open_CRAN_page("umx")
#' }
umx_open_CRAN_page <- function(package = "umx", inst=FALSE) {
for (p in package) {
# asString = deparse(substitute(parameter))
# if(!exists(asString)){
# p = asString
# }
# umx_msg(p)
# 1. Open the web pages
system(paste0("open 'https://cran.r-project.org/package=", p, "'"))
# 2. install if requested
if(inst){
install.packages(p)
} else {
# print("cleanup-code")
}
# 3. print data on current version and load
result = tryCatch({
ver = packageVersion(p)
print(ver)
library(p, character.only = TRUE)
}, warning = function(x) {
cat(p, "might not be installed locally:\n")
print(x)
}, error = function(x) {
cat(p, "might not be installed locally:\n")
print(x)
}, finally = {
#
})
}
}
#' Pad an Object with NAs
#'
#' This function pads an R object (list, data.frame, matrix, atomic vector)
#' with \code{NA}s. For matrices, lists and data.frames, this occurs by extending
#' each (column) vector in the object.
#' @param x An \R object (list, data.frame, matrix, atomic vector).
#' @param n The final length of each object.
#' @return - padded object
#' @export
#' @family Miscellaneous Utility Functions
#' @references - \url{https://github.com/kevinushey/Kmisc/tree/master/man}
#' @examples
#' umx_pad(1:3, 4)
#' umx_pad(1:3, 3)
umx_pad <- function(x, n) {
if (is.data.frame(x)) {
nrow = nrow(x)
attr(x, "row.names") = 1:n
for( i in 1:ncol(x) ) {
x[[i]] = c( x[[i]], rep(NA, times = n - nrow) )
}
return(x)
} else if (is.list(x)) {
if (missing(n)) {
max_len = max( sapply( x, length ) )
return( lapply(x, function(xx) {
return( c(xx, rep(NA, times=max_len-length(xx))) )
}))
} else {
return( lapply(x, function(xx) {
if (n > length(xx)) {
return( c(xx, rep(NA, times=n-length(xx))) )
} else {
return(xx)
}
}))
}
} else if (is.matrix(x)) {
return( rbind( x, matrix(NA, nrow=n-nrow(x), ncol=ncol(x)) ) )
} else {
if (n > length(x)) {
return( c( x, rep(NA, n-length(x)) ) )
} else {
return(x)
}
}
}
#' umx_apply
#'
#' Tries to make apply more readable. so "mean of x by columns", instead of "of x, by 2, mean"
#' Other functions to think of include:
#' [cumsum()], [rowSums()], [colMeans()], etc.
#'
#' @param FUN The function to apply.
#' @param of The dataframe to work with.
#' @param by Apply the function to columns or to rows (default = "columns")
#' @param ... optional arguments to FUN, e.g., na.rm = TRUE.
#' @return - object
#' @export
#' @seealso - [umx_aggregate()]
#' @family Miscellaneous Stats Functions
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' umx_apply(mean, mtcars, by = "columns")
#' umx_apply("mean", of = mtcars, by = "columns")
#' tmp = mtcars[1:3,]; tmp[1,1] = NA
#' umx_apply("mean", by = "rows", of = tmp)
#' umx_apply("mean", by = "rows", of = tmp, na.rm = TRUE)
umx_apply <- function(FUN, of, by = c("columns", "rows"), ...) {
by = match.arg(by)
by = ifelse(by == "rows", 1, 2)
apply(of, by, FUN, ...)
}
#' umx_as_numeric
#'
#' Convert each column of a dataframe to numeric
#'
#' @param df A [data.frame()] to convert
#' @param which which columns to convert (default (null) selects all)
#' @param force Whether to force conversion to numeric for non-numeric columns (defaults to FALSE)
#' @return - data.frame
#' @family Data Functions
#' @export
#' @references - <https://github.com/tbates/umx>
#' @examples
#' # make mpg into string, and cyl into a factor
#' df = mtcars
#' df$mpg = as.character(df$mpg)
#' df$cyl = factor(df$cyl)
#' df$am = df$am==1
#' df = umx_as_numeric(df); str(df) # mpg not touched
#' df = umx_as_numeric(df, force=TRUE); str(df) # mpg coerced back to numeric
#' \dontrun{
#' # coercing a real string will cause NAs
#' df$mpg = c(letters[1:16]); str(df) # replace mpg with letters.
#' df = umx_as_numeric(df, force=TRUE); str(df)
#' }
umx_as_numeric <- function(df, which = NULL, force = FALSE) {
# TODO umx_as_numeric: Handle matrices, vectors...
if(is.null(which)){
colsToConvert = names(df)
} else {
colsToConvert = which
}
umx_check_names(colsToConvert, data = df, die = TRUE)
if(!force){
# just the numeric names
colsToConvert = colsToConvert[umx_is_numeric(df)]
}
for (i in colsToConvert) {
df[ ,i] = as.numeric(df[ ,i, drop = TRUE])
}
return(df)
}
#' umx_find_object
#'
#' Find objects of a given class, whose name matches a search string.
#' The string (pattern) is grep-enabled, so you can match wild-cards
#'
#' @param pattern the pattern that matching objects must contain
#' @param requiredClass the class of object that will be matched
#' @return - a list of objects matching the class and name
#' @export
#' @references -
#' @family Miscellaneous Utility Functions
#' @examples
#' \dontrun{
#' umx_find_object("^m[0-9]") # mxModels beginning "m1" etc.
#' umx_find_object("", "MxModel") # all MxModels
#' }
umx_find_object <- function(pattern = ".*", requiredClass = "MxModel") {
# Use case: umxFindObject("Chol*", "MxModel")
matchingNames = ls(envir = sys.frame(-1), pattern = pattern) # envir
matchingObjects = c()
for (obj in matchingNames) {
if(class(get(obj))[1] == requiredClass){
matchingObjects = c(matchingObjects, obj)
}
}
return(matchingObjects)
}
#' umx_rename
#'
#' Returns a dataframe with variables renamed as desired.
#'
#' Unlike similar functions in other packages, it checks that the variables exist, and that the new names do not.
#'
#' Importantly, it also supports [regular expressions][regex]. This allows you to find and replace
#' text based on patterns and replacements. so to change "replacement" to "in place",
#' `grep=re(place)ment`, `replace= in \\1`.
#'
#' *note*:To use replace list, you must say c(old = "new"), not c(old -> "new")
#'
#' @param data The dataframe in which to rename variables
#' @param from List of existing names that will be found and replaced by the contents of replace. (optional: Defaults to NULL).
#' @param to If used alone, a named collection of c(oldName = "newName") pairs.
#' OR, if "from" is a list of existing names, the list of new names)
#' OR, if "regex" is a regular expression, the replace string)
#' @param regex Regular expression with matches will be replaced using replace as the replace string. (Optional: Defaults to NULL).
#' @param test Whether to report a "dry run", not changing anything. (Default = FALSE).
#' @param old deprecated: use from
#' @param replace deprecated: use to
#' @return - dataframe with columns renamed.
#' @export
#' @seealso [namez] to filter (and replace) names, Also [umx_check_names] to check for existence of names in a dataframe.
#' @family Data Functions
#' @md
#' @examples
#' tmp = mtcars
#'
#' tmp = umx_rename(tmp, to = c(cyl = "cylinder"))
#' # let's check cyl has been changed to cylinder...
#' namez(tmp, "c")
#'
#' # Alternate style: from->to, first with a test-run
#' # Dry run
#' tmp = umx_rename(tmp, from = "disp", to = "displacement", test= TRUE)
#' # Actually do it
#' tmp = umx_rename(tmp, from = c("disp"), to = c("displacement"))
#' umx_check_names("displacement", data = tmp, die = TRUE)
#' namez(tmp, "disp")
#'
#' # This will warn that "disp" does not exist (anymore)
#' new = c("auto", "displacement", "rear_axle_ratio")
#' tmp = umx_rename(tmp, from = c("am", "disp", "drat"), to = new)
#' namez(tmp, "a") # still updated am to auto (and rear_axle_ratio)
#'
#' # Test using regex (in this case to revert "displacement" to "disp")
#' tmp = umx_rename(tmp, regex = "lacement", to = "", test= TRUE)
#' tmp = umx_rename(tmp, regex = "lacement", to = "") # revert to disp
#' umx_names(tmp, "^d") # all names beginning with a d
#'
#' # dev: check deprecated format handled...
#' tmp = umx_rename(tmp, old = c("am", "disp", "drat"), replace = new)
umx_rename <- function(data, from = NULL, to = NULL, regex = NULL, test = FALSE, old = "deprecated_from", replace= "deprecated_to") {
# See also gdata::rename.vars(data, from, to)
if(any(old != "deprecated_from")){from = old; message("Polite message: Please use 'from' instead of 'old' in umx_rename()") }
if(any(replace != "deprecated_to")){to = replace; message("Polite message: Please use 'to' instead of 'replace' in umx_rename()") }
if(!is.null(attributes(from)$names)){
stop("You gave a list to from in umx_rename(). Lists (old='new') only allowed in to")
}
if(!is.null(from) && !is.null(regex)){
stop("Only one of from and regex can be used")
}
if(!is.null(regex)){
if(is.null(to)){
stop("Please set to to a valid replacement string!")
}
nameVector = umx_names(data)
if (is.null(nameVector)) {
stop(paste0("umx_rename requires a dataframe or something else with names(), ",
umx_str_from_object(data), " is a ", typeof(data)))
}
new_names = gsub(regex, to, nameVector)
if(test){
message("The following changes would be made (set test =FALSE to actually make them)")
message(length(nameVector), " names found. ",
length(nameVector[!(nameVector == new_names)]), " would be changed. Old:")
print(nameVector[!(nameVector == new_names)])
message("New:")
print(new_names[!(nameVector == new_names)])
} else {
if(class(data)[[1]] == "character"){
data = new_names
} else {
names(data) = new_names
}
}
invisible(data)
} else {
# Not regex
if(!is.null(from)){
# message("replacing from with to")
if(length(from) != length(to)){
stop("You are trying to replace ", length(from), " old names with ", length(to), " new names: Lengths must match")
}
names_to_replace = from
new_names_to_try = to
} else {
# to is a key-value list of names and replacements
names_to_replace = unlist(names(to))
new_names_to_try = unlist(unname(to))
}
from_names = names(data)
if(!all(names_to_replace %in% from_names)) {
warning("The following names did not appear in the dataframe:",
paste(names_to_replace[!names_to_replace %in% from_names], collapse= ", "), "\nPerhaps you already updated them")
}
if(anyDuplicated(names_to_replace)) {
err = paste("You are trying to update the following names more than once:",
paste(names_to_replace[duplicated(names_to_replace)], collapse = ", "))
stop(err)
}
if(anyDuplicated(new_names_to_try)) {
err = paste("You have the following duplicates in your to list:",
paste(new_names_to_try[duplicated(new_names_to_try)], collapse = ", ")
)
stop(err)
}
new_names = new_names_to_try[match(from_names, names_to_replace)]
if(test){
message("The following changes would be made (set test =FALSE to actually make them")
message("Names to be replaced")
print(names_to_replace)
message("replacement names:")
print(new_names[!is.na(new_names)])
invisible(data)
} else {
names(data) = new_names
setNames(data, ifelse(is.na(new_names), from_names, new_names)) # Also returns the new object
}
}
}
#' Search for text
#'
#' Search names if given a data.frame, or strings if given a vector of strings.
#'
#' The `namez` function is more flexible. A handy feature of `umx_grep` is that it can
#' search the labels of data imported from SPSS.
#'
#' *nb:* To simply grep for a pattern in a string use R's built-in [grep()] functions, e.g.:
#' `grepl("^NA\\[0-9]", "NA.3")`
#' @param df The [data.frame()] or string to search.
#' @param grepString the search string.
#' @param output the column name, the label, or both (default).
#' @param ignore.case whether to be case sensitive or not (default TRUE = ignore case).
#' @param useNames whether to search the names as well as the labels (for SPSS files with label metadata).
#' @return - list of matched column names and/or labels.
#' @seealso - [namez()], [umx_aggregate()], [grep()]
#' @family String Functions
#' @export
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' umx_grep(mtcars, "hp", output="both", ignore.case= TRUE)
#' umx_grep(c("hp", "ph"), "hp")
#' umx_grep(mtcars, "^h.*", output="both", ignore.case= TRUE)
#' \dontrun{
#' umx_grep(spss_df, "labeltext", output = "label")
#' umx_grep(spss_df, "labeltext", output = "name")
#' }
umx_grep <- function(df, grepString, output = c("both", "label", "name"), ignore.case=TRUE, useNames= FALSE) {
output = match.arg(output)
# if(length(grepString > 1)){
# for (i in grepString) {
# umx_grep_labels(df, i, output=output, ignore.case=ignore.case, useNames=useNames)
# }
if(is.data.frame(df)){
vLabels = attr(df, "variable.labels") # list of descriptive labels?
a = names(df)
if(is.null(vLabels)){
# message("No labels found")
return(grep(grepString, names(df), value=TRUE, ignore.case= ignore.case))
}
if(useNames) {
findIndex = grep(grepString,a, value=FALSE, ignore.case=ignore.case)
return( as.matrix(vLabels[findIndex]))
} else {
# need to cope with finding nothing
findIndex = grep(grepString,vLabels, value=FALSE, ignore.case=ignore.case)
if(output=="both") {
theResult = as.matrix(vLabels[findIndex])
} else if(output=="label"){
vLabels= as.vector(vLabels[findIndex])
theResult = (vLabels)
} else if(output=="name"){
theResult = names(vLabels)[findIndex]
}else{
stop(paste("bad choice of output:", output))
}
if(dim(theResult)[1]==0 |is.null(theResult)){
cat("using names!!!!\n")
findIndex = grep(grepString,a, value=FALSE, ignore.case=ignore.case)
return(as.matrix(vLabels[findIndex]))
} else {
return(theResult)
}
}
} else {
return(grep(grepString, df, value = TRUE, ignore.case = ignore.case))
}
}
#' A recipe Easter-egg for umx
#'
#' @description
#' How to cook steak.
#' @details Equipment matters. You should buy a heavy cast-iron skillet, and a digital internal thermometer.
#' Preferably cook over a gas flame.
#'
#' *note*: Cheaper cuts like blade steak can come out fine.
#'
#' @export
#' @family Miscellaneous Functions
#' @seealso - [OpenMx::omxBrownie()]
#' @examples
#' umxBrownie()
#' @md
umxBrownie <- function() {
message("Rub steak in a table spoon of salt, put it back in the fridge for an hour (longer is fine).\n",
"Place steak on a hot cast-iron skillet, with a little olive oil.\n",
"Turn steak as often as you wish. Control heat to below smoke point.\n",
"Remove and eat when internal temp reaches 110 \u0080 F.\n"
)
}
# ===========================
# = File handling functions =
# ===========================
#' Rename files
#'
#' Rename files. On OS X, the function can access the current front-most Finder window.
#' The file renaming is fast and, because you can use regular expressions too change names.
#'
#' @param findStr The pattern to find, i.e., "cats"
#' @param replaceStr The replacement pattern "\1 are not dogs"
#' @param baseFolder Folder to search in. Default ("Finder") will use the current front-most Finder window (on MacOS).
#' Set to NA for a "choose folder" dialog.
#' @param test Boolean determining whether to change files on disk, or just report on what would have happened (Defaults to test = TRUE)
#' @param ignoreSuffix Whether to ignore (don't search in) the suffix (file-type like .mpg) TRUE.
#' @param listPattern A pre-filter for files
#' @param overwrite Boolean determining if an existing file will be overwritten (Defaults to the safe FALSE)
#' @family File Functions
#' @return None
#' @export
#' @md
#' @examples
#' \dontrun{
#' # "Season 01" --> "S01" in current folder in MacOS Finder
#' umx_rename_file("[Ss]eason +([0-9]+)", replaceStr="S\\1", test = TRUE)
#'
#' # move date to end of file name
#' umx_rename_file("^(.*) *([0-9]{2}\\.[0-9]{2}\\.[0-9]+) *(.*)", replaceStr="\\1 \\3 \\2")
#'
#' }
umx_rename_file <- function(findStr = "old", replaceStr = NA, baseFolder = "Finder", test = TRUE, ignoreSuffix = TRUE, listPattern = NULL, overwrite = FALSE) {
umx_check(!is.na(replaceStr), "stop", "Please set a replaceStr to the replacement string you desire.")
# ==============================
# = 1. Set folder to search in =
# ==============================
if(baseFolder == "Finder"){
baseFolder = system(intern = TRUE, "osascript -e 'tell application \"Finder\" to get the POSIX path of (target of front window as alias)'")
message("Using front-most Finder window:", baseFolder)
} else if(baseFolder == "") {
baseFolder = paste(dirname(file.choose(new = FALSE)), "/", sep = "") # choose a directory
message("Using selected folder:", baseFolder)
}
# =================================================
# = 2. Find files matching listPattern or findStr =
# =================================================
fileList = list.files(baseFolder, pattern = listPattern)
message("found ", length(fileList), " possible files")
# if(test){
# omxQuotes(fileList)
# }
# return(fileList)
changed = 0
for (fn in fileList) {
if(ignoreSuffix){
baseName = sub(pattern = "(.*)(\\..*)$", x = fn, replacement = "\\1")
}else{
baseName = fn
}
if(grepl(pattern = findStr, baseName, perl= TRUE)){
# found pattern
if(ignoreSuffix){
# pull suffix and baseName (without suffix)
baseName = sub(pattern = "(.*)(\\..*)$", x = fn, replacement = "\\1")
suffix = sub(pattern = "(.*)(\\..*)$", x = fn, replacement = "\\2")
fnew = gsub(findStr, replacement = replaceStr, x = baseName, perl= TRUE) # replace all instances
fnew = paste0(fnew, suffix)
} else {
fnew = gsub(findStr, replacement = replaceStr, x = fn, perl= TRUE) # replace all instances
}
if(test){
message(fn, " would be changed to: ", omxQuotes(fnew))
} else {
if((!overwrite) & file.exists(paste(baseFolder, fnew, sep = ""))){
message("renaming ", fn, "to", fnew, "failed as already exists. To overwrite set T")
} else {
file.rename(paste0(baseFolder, fn), paste0(baseFolder, fnew))
changed = changed + 1;
}
}
}
}
if(test & (changed == 0)){
message("no matches for change (PS: once you get some hits, set test = FALSE to actually change files.")
} else {
umx_msg(changed)
}
}
#' Move files
#'
#' On OS X, `umx_move_file` can access the current front-most Finder window.
#' The file moves are fast and, because you can use regular expressions, powerful.
#'
#' @param baseFolder The folder to search in. If set to "Finder" (and you are on OS X) it will use the current
#' front-most Finder window. If it is blank, a choose folder dialog will be thrown.
#' @param regex string to select files to process within the selected folder.
#' @param fileNameList List of files to move.
#' @param destFolder Folder to move files to.
#' @param test Boolean determining whether to change the names, or just report a dry run.
#' @param overwrite Boolean determining whether to overwrite files or not (default = FALSE (safe)).
#' @return None
#' @seealso [file.rename()], [regex()]
#' @family File Functions
#' @md
#' @export
#' @examples
#' \dontrun{
#' base = "~/Desktop/"
#' dest = "~/Music/iTunes/iTunes Music/Music/"
#' umx_move_file(baseFolder = base, fileNameList = toMove, destFolder = dest, test= TRUE)
#'
#' # ============================================================
#' # = Move all files in downloads ending in ".jpeg" to Desktop =
#' # ============================================================
#' umx_move_file(baseFolder = "~/Downloads/", regex=".jpeg",
#' destFolder = "~/Desktop/", test= TRUE)
#' }
#'
umx_move_file <- function(baseFolder = NA, regex = NULL, fileNameList = NA, destFolder = NA, test = TRUE, overwrite = FALSE) {
if(!is.null(regex)){
if(!is.na(fileNameList)){
stop("Can't use regex and a fileNameList")
} else {
fileNameList = list.files(baseFolder, pattern = regex)
}
}
if(is.na(destFolder)){
stop("destFolder can't be NA")
}
if(baseFolder == "Finder"){
baseFolder = system(intern = TRUE, "osascript -e 'tell application \"Finder\" to get the POSIX path of (target of front window as alias)'")
message("Using front-most Finder window:", baseFolder)
} else if(baseFolder == "") {
baseFolder = paste(dirname(file.choose(new = FALSE)), "/", sep="") # choose a directory
message("Using selected folder:", baseFolder)
}
moved = 0
for (fn in fileNameList) {
if(test){
message("would move ", fn, " to ", destFolder)
moved = moved + 1;
message("Would have moved ", moved)
} else {
if((!overwrite) & file.exists(paste0(destFolder, fn))){
message("moving ", fn, "to", destFolder, "failed as already exists. To overwrite set overwrite= TRUE")
} else {
file.rename(paste0(baseFolder, fn), paste0(destFolder, fn))
moved = moved + 1;
}
}
message("Moved ", moved)
}
}
#' Open a file or folder
#'
#' Open a file or folder. Works on OS X, mostly on windows, and hopefully on unix.
#'
#' NOTE: Your filepath is [shQuote()]'d by this function.
#' @param filepath The file to open
#' @return None
#' @export
#' @family File Functions
#' @md
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @examples
#' \dontrun{
#' umx_open() # Default is to open working directory getwd()
#' umx_open("~/bin/umx/R/misc_and_utility copy.r")
#' }
umx_open <- function(filepath = getwd()) {
filepath = normalizePath(filepath)
if (umx_check_OS("Windows")){
shell(shQuote(filepath, type='cmd'), 'cmd.exe')
} else {
if(umx_check_OS("OSX")){
opener = "open "
} else { # *nix?
opener = "xdg-open "
}
system(paste(opener, shQuote(filepath)))
# system2(opener, shQuote(filepath)) # possibly more robust.
# check when around multiple machine types
}
}
#' umx_check_OS
#'
#' Check what OS we are running on (current default is OS X). Returns a boolean.
#' Optionally warn or die on failure of the test
#'
#' @param target Which OS(s) you wish to check for (default = "OSX")
#' @param action What to do on failure of the test: nothing (default), warn or die
#' @return - TRUE if on the specified OS (else FALSE)
#' @export
#' @family Test
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @md
#' @examples
#' umx_check_OS()
umx_check_OS <- function(target=c("OSX", "SunOS", "Linux", "Windows"), action = c("ignore", "warn", "die")) {
action = match.arg(action)
target = match.arg(target)
# OSX == Darwin
# Solaris == SunOS
sysinfo = Sys.info()
if (!is.null(sysinfo)){
os = sysinfo['sysname']
if (os == 'Darwin'){
os = "OSX"
}
} else {
os = .Platform$OS.type
if (grepl("linux-gnu", R.version$os)){
os = "Linux"
}
}
isTarget = (target == os)
if(!isTarget){
if(action == "die"){
stop("Sorry: You must be running on ", target, " OS. You're on ", os)
} else if(action == "warn"){
message("i was expecting the OS to be ", target, " not ", os)
}
}
return(isTarget)
}
#' Convert an excel spreadsheet in a text file on sql statements.
#'
#' Unlikely to be of use to anyone but the package author :-)
#'
#' On OS X, by default, the file selected in the front-most Finder window will be chosen.
#' If it is blank, a choose file dialog will be thrown.
#'
#' Read an xlsx file and convert into SQL insert statements (placed on the clipboard)
#' On MacOS, the function can access the current front-most Finder window.
#'
#' The file name should be the name of the test.
#' Columns should be headed:
#' itemText direction scale type \[optional response options\]
#'
#' The SQL fields generated are:
#' itemID, test, native_item_number, item_text, direction, scale, format, author
#' @details
#' tabbedPlus: list scored from 0 to n-1
#'
#' tabbedVertPlus: tabbed, but vertical lay-out
#'
#' number 2+2\<itemBreak\>min='0' max='7' step='1'
#'
#' 5fm Scored 1-5, anchored: Strongly Disagree | Disagree | Neutral | Agree | Strongly Agree
#'
#' intro (not) scored, and sequenced as item 0
#'
#' @param theFile The xlsx file to read. Default = "Finder")
#' @family File Functions
#' @return None
#' @export
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' # An example Excel spreadsheet
#' # local uncompiled path
#' fp = system.file("inst/extdata", "GQ6.sql.xlsx", package = "umx")
#' # installed path
#' fp = system.file("extdata", "GQ6.sql.xlsx", package = "umx")
#' umx_open(fp)
#' umx_make_sql_from_excel() # Using file selected in front-most Finder window
#' umx_make_sql_from_excel("~/Desktop/test.xlsx") # provide a path
#' }
umx_make_sql_from_excel <- function(theFile = "Finder") {
if(theFile == "Finder"){
umx_check_OS("OSX")
theFile = system(intern = TRUE, "osascript -e 'tell application \"Finder\" to get the POSIX path of (selection as alias)'")
message("Using file selected in front-most Finder window:", theFile)
} else if(theFile == "") {
theFile = file.choose(new = FALSE) # choose a file
message("Using selected file:", theFile)
} else if(theFile == "make") {
theFile = system.file("extdata", "GQ6.sql.xlsx", package = "umx")
}
umx_check(file.exists(theFile), message= paste0("file:'", theFile, "' does not exist..."))
# remove suffix (i.e., .xlsx )
testName = umx_trim(basename(theFile), "\\..+$")
df = openxlsx::read.xlsx(theFile, sheet = 1, stringsAsFactors= FALSE)
expect8 = c("itemText", "direction", "scale", "type")
if(!all(expect8 %in% names(df))){
stop(paste("I expected the following required column names:\n", omxQuotes(expect8), "\nYou gave me:",
omxQuotes(names(df))), call. = FALSE)
}
nItems = dim(df)[1]
nCols = dim(df)[2]
for (i in 1:nCols) {
df[,i] = as.character(df[,i])
}
df[df == ""] = NA
pre = "INSERT INTO Items VALUES ('"
end = paste0("');")
o = data.frame(sql="junk", stringsAsFactors = FALSE) ;
itemNumber = 1
for (lineNumber in 1:nItems) {
direction = df[lineNumber, "direction"]
scale = df[lineNumber, "scale"]
type = df[lineNumber, "type"]
if (type=="info" & itemNumber == 1){
# this will fail if there are two info questions at the top
itemNumber = 0
}
itemText = df[lineNumber, "itemText"]
# Any more cells in <itemBreak>?
if(nCols > 5){
items = df[lineNumber, 5:nCols]
if(any(!is.na(items))){
itemText = paste0(itemText, "<itemBreak>", paste(items[!is.na(items)], collapse = "<itemBreak>"))
}
}
thisRow = paste(pre, testName, itemNumber, itemText, direction, scale, type, testName, end, sep = "', '")
thisRow = umx_names(thisRow, pattern = ", '');", replacement = ");")
o[itemNumber, ] = thisRow
itemNumber = itemNumber + 1
}
umx_write_to_clipboard(x = o)
message("sql is on clipboard")
}
#' umx_write_to_clipboard
#'
#' @description
#' umx_write_to_clipboard writes data to the clipboard
#'
#' @details
#' Works on Mac. Let me know if it fails on windows or Unix.
#' @param x something to paste to the clipboard
#' @return None
#' @export
#' @family File Functions
#' @examples
#' \dontrun{
#' umx_write_to_clipboard("hello")
#' }
umx_write_to_clipboard <- function(x) {
if(umx_check_OS("OSX")){
clipboard = pipe("pbcopy", "w")
write.table(x, file = clipboard, sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
close(clipboard)
} else if (umx_check_OS("Windows")){
write.table(x, file = "clipboard", sep = "\t", col.names = NA)
}else{
message("clipboard not implemented for *nix - awaiting a reliable solution.
file(description='clipboard') might help. See:
https://stackoverflow.com/questions/13438556/how-do-i-copy-and-paste-data-into-r-from-the-clipboard#13438558")
}
}
#' dl_from_dropbox
#'
#' Download a file from Dropbox, given either the url, or the name and key
#'
#' Improvements would include error handling...
#' @param x Either the file name, or full dropbox URL (see example below)
#' @param key the code after s/ and before the file name in the dropbox url
#' @return None
#' @export
#' @family File Functions
#' @references - \url{https://thebiobucket.blogspot.kr/2013/04/download-files-from-dropbox.html}
#' @examples
#' \dontrun{
#' dl_from_dropbox("https://dl.dropboxusercontent.com/s/7kauod48r9cfhwc/tinytwinData.rda")
#' dl_from_dropbox("tinytwinData.rda", key = "7kauod48r9cfhwc")
#' }
dl_from_dropbox <- function(x, key=NULL){
# depends on RCurl::getBinaryURL
if(is.null(key)){
bin = RCurl::getBinaryURL(x, ssl.verifypeer = FALSE)
x = sub("^.+/(.*)$", "\\1", x, ignore.case = FALSE, perl = FALSE, fixed = FALSE, useBytes = FALSE)
} else {
# user has provided key and file name, so concatenate with https...
bin = RCurl::getBinaryURL(paste0("https://dl.dropboxusercontent.com/s/", key, "/", x), ssl.verifypeer = FALSE)
}
con = file(x, open = "wb")
writeBin(bin, con)
close(con)
message(noquote(paste(x, "read into", getwd())))
}
# ===================
# = Stats functions =
# ===================
# =======================
# = Financial utilities =
# =======================
#' Work the valuation of a company
#'
#' @description
#' `fin_valuation` uses the revenue, operating margin, expenses and PE to compute a market capitalization.
#' Better to use a more powerful online site.
#'
#' @details
#' Revenue is multiplied by opmargin to get a gross profit. From this the proportion specified in `expenses` is subtracted
#' and the resulting earnings turned into a price via the `PE`
#'
#' @param revenue Revenue of the company
#' @param opmargin Margin on operating revenue
#' @param expenses Additional fixed costs
#' @param PE of the company
#' @param symbol Currency
#' @param use reporting values in "B" (billion) or "M" (millions)
#' @return - value
#' @export
#' @family Miscellaneous Functions
#' @seealso - [fin_interest()], [fin_NI()], [fin_percent()]
#' @md
#' @examples
#' fin_valuation(rev=7e9, opmargin=.1, PE=33)
#' # Market cap = $18,480,000,000
#' # (Based on PE= 33, operating Income of $0.70 B, and net income =$0.56B
#'
fin_valuation <- function(revenue=6e6*30e3, opmargin=.08, expenses=.2, PE=30, symbol = "$", use = c("B", "M")) {
use = match.arg(use)
if(use=="B"){
divisor=1e9
} else {
divisor=1e6
}
operatingIncome = revenue * opmargin
netIncome = operatingIncome *(1-expenses)
marketCap = netIncome*PE
class(marketCap) = 'money'; attr(marketCap, 'symbol') = symbol
class(netIncome) = 'money'; attr(netIncome, 'symbol') = symbol
class(operatingIncome) = 'money'; attr(operatingIncome, 'symbol') = symbol
cat("Market cap = ", print(marketCap, cat=F))
cat("\n(Based on PE= ", PE, ", operating Income of ", print(operatingIncome/divisor, cat=F), " ", use, ", and net income =", print(netIncome/divisor, cat=F), use, "\n", sep = "")
invisible(marketCap)
}
#' Compute the value of a principal & annual deposits at a compound interest over a number of years
#' @description
#' Allows you to determine the final value of an initial `principal` (with optional
#' periodic `deposits`), over a number of years (`yrs`) at a given rate of `interest`.
#' Principal and deposits are optional. You control compounding periods each year (n) and whether deposits occur at the beginning or end of the year.
#' The function outputs a nice table of annual returns, formats the total using a user-settable currency `symbol`. Can also `report` using a web table.
#'
#' *notes*: Graham valuation: fair P/E = 9 + (1.5 * growth%). e.g. $INTEL fair P/E = 9+.5*3 = 10.5 up to 9+2*10 = 29
#' Can move the weighting between a conservative .5 and an optimistic 2 (in terms of how long the growth will last and how low the hurdle rate is)
#'
#'
#' @param principal The initial investment at time 0 (default 100)
#' @param deposits Optional periodic additional investment each *year*.
#' @param interest Annual interest rate (default .05)
#' @param inflate How much to inflate deposits over time (default 0)
#' @param yrs Duration of the investment (default 10).
#' @param n Compounding intervals per year (default 12 (monthly), use 365 for daily)
#' @param when Deposits made at the "beginning" (of each year) or "end"
#' @param symbol Currency symbol to embed in the result.
#' @param report "markdown" or "html",
#' @param table Whether to print a table of annual returns (default TRUE)
#' @param largest_with_cents Default = 0
#' @param baseYear Default = current year (for table row labels)
#' @param final if set (default = NULL), returns the rate required to turn principal into final after yrs (principal defaults to 1)
#' @return - Value of balance after yrs of investment.
#' @export
#' @family Miscellaneous Functions
#' @seealso - [umx_set_dollar_symbol()], [fin_percent()], [fin_NI()], [fin_valuation()]
#' @references - <https://en.wikipedia.org/wiki/Compound_interest>
#' @md
#' @examples
#' \dontrun{
#' # 1. Value of a principal after yrs years at 5% return, compounding monthly.
#' # Report in browser as a nice table of annual returns and formatted totals.
#' fin_interest(principal = 5000, interest = 0.05, rep= "html")
#' }
#'
#' # Report as a nice markdown table
#' fin_interest(principal = 5000, interest = 0.05, yrs = 10)
#'
#' umx_set_dollar_symbol("£")
#' # 2 What rate is needed to increase principal to final value in yrs time?
#' fin_interest(final = 1.4, yrs=5)
#' fin_interest(principal = 50, final=200, yrs = 5)
#'
#' # 3. What's the value of deposits of $100/yr after 10 years at 7% return?
#' fin_interest(deposits = 100, interest = 0.07, yrs = 10, n = 12)
#'
#' # 4. What's the value of $20k + $100/yr over 10 years at 7% return?
#' fin_interest(principal= 20e3, deposits= 100, interest= .07, yrs= 10, symbol="$")
#'
#' # 5. What is $10,000 invested at the end of each year for 5 years at 6%?
#' fin_interest(deposits = 10e3, interest = 0.06, yrs = 5, n=1, when= "end")
#'
#' # 6. What will $20k be worth after 10 years at 15% annually (n=1)?
#' fin_interest(deposits=20e3, interest = 0.15, yrs = 10, n=1, baseYear=1)
#' # $466,986
#'
#' # manual equivalent
#' sum(20e3*(1.15^(10:1))) # 466985.5
#'
#' # 7. Annual (rather than monthly) compounding (n=1)
#' fin_interest(deposits = 100, interest = 0.07, yrs = 10, n=1)
#'
#' # 8 Interest needed to increase principal to final value in yrs time.
#' fin_interest(principal = 100, final=200, yrs = 5)
#'
fin_interest <- function(principal = 100, deposits = 0, inflate = 0, interest = 0.05, yrs = 10, final= NULL, n = 12, when = "beginning", symbol = NULL, largest_with_cents = 0, baseYear= as.numeric(format(Sys.time(), "%Y")), table = TRUE, report= c("markdown", "html")){
report = match.arg(report)
if(is.null(symbol)){symbol = umx_set_dollar_symbol(silent=TRUE)}
if(principal==0){
caption= paste0("Compounding ", bucks(deposits, symbol, cat=TRUE), " deposits over ", yrs, " years at ", interest*100, "% interest with ", inflate*100, "% inflation.")
} else {
caption= paste0("Compounding ", bucks(principal, symbol, cat=TRUE), " principle plus ", bucks(deposits, symbol, cat=TRUE), " annual deposits, ", interest * 100, "% interest and ", inflate*100, "% inflation.")
}
if(inflate != 0){
deposits = c(deposits, rep(deposits, times = yrs-1) *(1+inflate)^c(1:(yrs-1)))
}else{
deposits = rep(deposits, times = yrs)
}
if(!is.null(final)){
# final = prin*(1+rate)^y
if(principal==0){ principal=1 }
return((final/principal)^(1/(yrs+1))-1)
# rate is the years root of (final *prin?)
}
# 1. compute compounding rate per unit time n (allowing for zero interest so 1.0)
rate = ifelse(interest==0, 1, 1+(interest/n))
tableOut = data.frame(Year = NA, Deposits = NA, Interest = NA, Total_Deposits = NA, Total_Interest = NA, Total = scales::dollar(principal, prefix = symbol, largest_with_cents = 0))
balance = principal
totalDeposits = 0
totalInterest = 0
for (yr in 1:yrs) {
# 1. Compute compounding rate per unit time n (allowing for zero interest so 1.0)
if(when == "beginning"){
# Deposits at the beginning of each year
thisInterest = ((balance + deposits[yr]) * rate^n) - (balance + deposits[yr])
} else {
# Deposits at the end of the year
thisInterest = (balance * rate^n) - balance
}
totalDeposits = (totalDeposits + deposits[yr])
totalInterest = (totalInterest + thisInterest)
balance = (balance + deposits[yr] + thisInterest)
thisRow = c(Year=yr+baseYear, Deposit= deposits[yr], Interest = thisInterest, Total_Deposit = totalDeposits, Total_Interest = totalInterest, Total = balance)
thisRow = c(thisRow[1], scales::dollar(thisRow[-1], prefix = symbol, largest_with_cents = largest_with_cents))
tableOut = rbind(tableOut, thisRow)
}
if(table){
# principal = 0, deposits = 0, inflate = 0, interest = 0.05, yrs
umx_print(tableOut, justify = "right", caption = caption, report=report)
}
if(length(deposits)==1){
# 2. compute compounded value of the principal (initial deposit)
Compound_interest_for_principal = principal* rate^(n*yrs)
# 3. compute compounded value of the deposits
if(interest==0){
Future_value_of_a_series = deposits * yrs
} else {
# beginning: A = PMT * (((1 + r/n)^(nt) - 1) / (r/n))
# end : A = PMT * (((1 + r/n)^(nt) - 1) / (r/n)) * (1+r/n)
if(when == "beginning"){
# deposits at the beginning of each year
periods = (yrs:1)*n
Future_value_of_a_series = sum(deposits*(rate^periods))
} else {
# deposits at the end of the year
periods = ((yrs-1):1)*n
Future_value_of_a_series = sum(deposits*(rate^periods)) + (1*deposits)
}
}
Total = Compound_interest_for_principal+ Future_value_of_a_series
} else {
Total = balance
}
class(Total) = 'money'
attr(Total, 'symbol') = symbol
return(Total)
}
#' Compute NI given annual Earnings.
#'
#' @description
#' Employees pay contributions at 12%% on annual earnings between £9,568 and £50,270. Above that you pay at 2%%.
#' Employers pay at 13.8%% on all annual earnings of more than £8,840, although there are different thresholds
#' for those under the age of 21 and for apprentices under the age of 25.
#'
#' @param annualEarnings Employee annual earnings.
#' @param symbol Currency symbol to embed in the result.
#' @return - NI
#' @export
#' @family Miscellaneous Functions
#' @seealso - [fin_interest()], [fin_percent()], [fin_valuation()]
#' @references - <https://www.telegraph.co.uk/tax/tax-hacks/politicians-running-scared-long-overdue-national-insurance-overhaul/>
#' @md
#' @examples
#' fin_NI(42e3)
#' fin_NI(142000)
#'
fin_NI <- function(annualEarnings, symbol = "\u00A3") {
if(annualEarnings < 50270){
employee = .12 * max(0, (annualEarnings- 9568))
} else {
employee = (.12 * (annualEarnings- 9568)) + (.02 * (annualEarnings-50270))
}
employer = .138 * max((annualEarnings - 8840), 0)
Total = employer + employee
class(Total) = 'money'
attr(Total, 'symbol') = symbol
cat(paste0("Employer pays ", bucks(employer, symbol = symbol, cat = FALSE), ", and employee pays ", bucks(employee, symbol = symbol, cat=FALSE),
". So ", round((employer+employee)/annualEarnings*100, 2), " % total!\n")
)
return(Total)
}
#' Print a money object
#'
#' @description Print function for "money" objects, e.g. [fin_interest()].
#'
#' @aliases bucks print
#' @param x money object.
#' @param symbol Default prefix if not set.
#' @param big.mark option defaulting to ","
#' @param decimal.mark option defaulting to "."
#' @param trim option defaulting to TRUE
#' @param largest_with_cents option defaulting to 1e+05
#' @param negative_parens option defaulting to "hyphen"
#' @param ... further arguments passed to or from other methods. also cat =F to return string
#' @return - invisible
#' @seealso - [umx::fin_percent()], [umx::fin_interest()], [scales::dollar()]
#' @md
# #' @family print
#' @export
#' @examples
#' bucks(100 * 1.05^32)
#' fin_interest(deposits = 20e3, interest = 0.07, yrs = 20)
#'
bucks <- function(x, symbol = "$", big.mark = ",", decimal.mark = ".", trim = TRUE, largest_with_cents = 1e+05, negative_parens = c("hyphen", "minus", "parens"), ...) {
dot.items = list(...) # grab all the dot items cat
cat = ifelse(is.null(dot.items[["cat"]]), TRUE, dot.items[["cat"]])
if(is.null(dot.items[["cat"]])){
cat = TRUE
} else {
cat = FALSE
dot.items[["cat"]] = NULL
}
if(!is.null(attr(x, 'symbol')) ){
symbol = attr(x, 'symbol')
}
formatted = scales::dollar(as.numeric(x), prefix = symbol, big.mark = big.mark, decimal.mark = decimal.mark, trim =trim, largest_with_cents = largest_with_cents, style_negative = negative_parens, ...)
if(cat){
cat(formatted)
} else {
formatted
}
}
#' @export
#' @method print money
print.money <- bucks
#' Compute the percent change needed to return to the original value after percent off (or on).
#'
#' @description
#' Determine the percent change needed to "undo" an initial percent change. Has a plot function as well.
#' If an amount of $100 has 20% added, what percent do we need to drop it by to return to the original value?
#'
#' `fin_percent(20)` yields $100 increased by 20% = $120 (Percent to reverse = -17%)
#'
#' @param percent Change in percent (enter 10 for 10%, not 0.1)
#' @param value Principal
#' @param symbol value units (default = "$")
#' @param digits Rounding of results (default 2 places)
#' @param plot Whether to plot the result (default TRUE)
#' @param logY Whether to plot y axis as log (TRUE)
#' @return - new value and change required to return to baseline.
#' @export
#' @family Miscellaneous Functions
#' @seealso - [fin_interest()]
#' @md
#' @examples
#' # Percent needed to return to original value after 10% taken off
#' fin_percent(-10)
#'
#' # Percent needed to return to original value after 10% added on
#' fin_percent(10)
#'
#' # Percent needed to return to original value after 50% off 34.50
#' fin_percent(-50, value = 34.5)
fin_percent <- function(percent, value= 100, symbol = "$", digits = 2, plot = TRUE, logY = TRUE) {
percent = percent/100
newValue = value * (1 + percent)
percent_to_reverse = (value/newValue) - 1
class(newValue) = 'percent'
attr(newValue, 'oldValue') = value
attr(newValue, 'percent') = percent
attr(newValue, 'digits') = digits
attr(newValue, 'symbol') = symbol
attr(newValue, 'percent_to_reverse') = percent_to_reverse
if(plot){
plot(newValue, logY = logY)
}else{
return(newValue)
}
}
#' Print a percent object
#'
#' Print method for "percent" objects: e.g. [umx::fin_percent()].
#'
#' @param x percent object.
#' @param ... further arguments passed to or from other methods.
#' @return - invisible
#' @seealso - [umx::fin_percent()]
#' @md
#' @method print percent
#' @export
#' @examples
#' # Percent needed to return to original value after 10% off
#' fin_percent(-10)
#' # Percent needed to return to original value after 10% on
#' fin_percent(10)
#'
#' # Percent needed to return to original value after 50% off 34.50
#' fin_percent(-50, value = 34.5)
#'
print.percent <- function(x, ...) {
if(!is.null(attr(x, 'digits')) ){
digits = attr(x, 'digits')
}
oldValue = round(attr(x, 'oldValue'), digits)
percentChange = attr(x, 'percent')
symbol = attr(x, 'symbol')
percent_to_reverse = round(attr(x, 'percent_to_reverse'), digits)
dir = ifelse(percentChange < 0, "decreased", "increased")
cat(symbol, oldValue, " ", dir , " by ", percentChange*100, "% = ", symbol, x, " (Percent to reverse = ", percent_to_reverse*100, "%)", sep="")
}
#' Plot a percent change graph
#'
#' Plot method for "percent" objects: e.g. [umx::fin_percent()].
#'
#' @param x percent object.
#' @param ... further arguments passed to or from other methods.
#' @return - invisible
#' @seealso - [umx::fin_percent()]
#' @md
#' @method plot percent
#' @export
#' @examples
#' # Percent needed to return to original value after 10% off
#' fin_percent(-10)
#' # Percent needed to return to original value after 10% on
#' tmp = fin_percent(10)
#' plot(tmp)
#'
#' # Percent needed to return to original value after 50% off 34.50
#' fin_percent(-50, value = 34.5, logY = FALSE)
#'
plot.percent <- function(x, ...) {
tmp = list(...) # pull logY if passed in
logY = tmp$logY
symbol = attr(x, 'symbol')
digits = attr(x, 'digits')
oldValue = round(attr(x, 'oldValue'), digits)
percentChange = attr(x, 'percent')
percent_to_reverse = round(attr(x, 'percent_to_reverse'), digits)
dir = ifelse(percentChange < 0, "decreased", "increased")
# fnReversePercent(-.1)
fnReversePercent <- function(x) {
# 1/(1+.1)
percentOn = x/100
newValue = (1 + percentOn)
percent_to_reverse = 1-(1/newValue)
return(-percent_to_reverse*100)
}
# x range = -100 (%) to +500 (%)?
# y = -100 to +200?
# y range = -100 to +200?
if(percentChange>0){
p = ggplot(data.frame(x = c(0, 90)), aes(x))
lab = paste0(round(percentChange*100, 2), "% on = ", round(percent_to_reverse * 100, 2), "% off", sep = "")
labXpos = 50
labYpos = -20
logY = FALSE
} else {
p = ggplot(data.frame(x = c(-90, 0)), aes(x))
lab = paste0(round(percentChange*100, 2), "% off = ", round(percent_to_reverse * 100, 2), "% on", sep = "")
labXpos = -50
labYpos = 700
}
if(is.null(logY)||!(logY)){
p = p + ggplot2::scale_y_continuous(n.breaks = 8) + ggplot2::scale_x_continuous(n.breaks = 10)
p = p + cowplot::draw_label(lab, vjust = 1, hjust = .5, x = labXpos, y = labYpos, color= "grey")
# hor & vert
p = p + ggplot2::geom_segment(x = percentChange*100, xend=-100, y=percent_to_reverse*100, yend=percent_to_reverse*100, alpha=.5, color = "lightgrey")
p = p + ggplot2::geom_segment(x = percentChange*100, xend=percentChange*100, y=-10, yend=percent_to_reverse*100, alpha=.5, color = "lightgrey")
} else {
p = p + ggplot2::scale_y_continuous(n.breaks = 8, trans="log10") + ggplot2::scale_x_continuous(n.breaks = 10)
p = p + cowplot::draw_label(lab, vjust = 1, hjust = .5, x = labXpos, y = log10(labYpos), color= "grey")
# hor & vert
p = p + ggplot2::geom_segment(x = percentChange*100, xend=-100 , y= log10(percent_to_reverse*100), yend= log10(percent_to_reverse*100), alpha=.5, color = "lightgrey")
p = p + ggplot2::geom_segment(x = percentChange*100, xend=percentChange*100, y= -10, yend= log10(percent_to_reverse*100), alpha= .5, color = "lightgrey")
}
p = p + ggplot2::stat_function(fun = fnReversePercent, color= "lightblue")
p = p + labs(x = "Percent change", y = "Percent change to reverse", title = paste0(percentChange*100, " percent ", ifelse(percentChange>0, "on ", "off "), oldValue, " = ", (1+percentChange)*oldValue))
# p = p + ggplot2::geom_area() can't do with stat fun ...
# p = p + cowplot::draw_label("\u2B55", hjust=0, vjust=1, x = percentChange*100, y = percent_to_reverse*100, color = "lightblue")
if(umx_set_plot_use_hrbrthemes(silent = TRUE)){
# p = p + hrbrthemes::theme_ipsum()
p = p + hrbrthemes::theme_ft_rc()
} else {
# p = p + ggplot2::theme_bw()
p = p + cowplot::theme_cowplot(font_size = 11)
}
print(p)
cat(symbol, oldValue, " ", dir , " by ", percentChange*100, "% = ", symbol, x, " (Percent to reverse = ", percent_to_reverse*100, "%)", sep="")
invisible(p)
}
#' Easily plot functions in R
#'
#' @description
#' A wrapper for [ggplot2::stat_function()]
#'
#' @details Easily plot a function - like sin, using ggplot.
#'
#' @param fun Function to plot. Also takes strings like "sin(x) + sqrt(1/x)".
#' @param min x-range min.
#' @param max x-range max.
#' @param xlab = Optional x axis label.
#' @param ylab = Optional y axis label.
#' @param title Optional title for the plot.
#' @param logY Set to, e.g. "log" to set COORDINATE of y to log.
#' @param p Optional plot onto which to draw the function.
#' @return - A ggplot graph
#' @export
#' @family Plotting functions
#' @seealso - [ggplot2::stat_function()]
#' @md
#' @examples
#' \dontrun{
#' # Uses fonts not available on CRAN
#' umxPlotFun(sin, max= 2*pi)
#' umxPlotFun("sqrt(1/x)", max= 2*pi)
#' umxPlotFun(sin, max= 2*pi, ylab="Output of sin", title="My Big Graph")
#' p = umxPlotFun(function(x){x^2}, max= 100, title="Supply and demand")
#' umxPlotFun(function(x){100^2-x^2}, p = p)
#'
#' # Controlling other plot features
#' umxPlotFun(c("sin(x)", "x^3")) + ylim(c(-1,5))
#' }
#'
umxPlotFun <- function(fun= c(dnorm, "sin(x) + sqrt(1/x)"), min= -1, max= 5, xlab = NULL, ylab = NULL, title = NULL, logY = c("no", "log", "log10"), p = NULL) {
logY = xmu_match.arg(logY, c("no", "log", "log10"), check = FALSE)
if(inherits(fun, "numeric")){
stop("If you write a function symbolically, you need to put it in quotes, e.g. 'x^2'")
} else if(inherits(fun, "character")){
make_function <- function(args, body, env = parent.frame()) {
args = as.pairlist(args)
eval(call("function", args, body), env)
}
funOut = c()
for (i in fun) {
if(is.null(title)){ title = paste0("Plot of ", i) }
# failed ideas to format as latex...
# if(is.null(title)){ title = parse(text=paste0("'Plot of '", expression(i) ) ) }
# if(is.null(title)){ title = parse(text = paste0("Plot of ", i)) }
if(is.null(ylab)){ ylab = i}
thisFun = make_function(alist(x=NA), parse(text = i)[[1]] )
funOut = c(funOut, thisFun)
}
fun = funOut # 1 or more functions
}else{
# Got a bare function like sin
fun = list(fun)
}
# plot function 1
if(!is.null(p)){
if(is.na(max)){
p = p + ggplot2::stat_function(fun = fun[[1]])
} else {
p = p + ggplot2::stat_function(fun = fun[[1]], xlim= c(min, max))
}
}else{
p = ggplot(data.frame(x = c(min, max)), aes(x) )
if(logY != "no"){
p = p + ggplot2::coord_trans(y = logY)
}
p = p + ggplot2::stat_function(fun = fun[[1]])
xlab = ifelse(!is.null(xlab), xlab , "X value")
if(is.null(ylab)){
if(length(as.character(quote(fun[[1]]))) == 1){
ylab = paste0(as.character(quote(fun[[1]]), " of x"))
} else {
ylab = paste0("Function of X")
}
}
if(is.null(title)){
if(length(as.character(quote(fun[[1]]))) == 1){
pref= "Plot of function: "
}else{
pref= "Plot of Functions: "
}
result = tryCatch({
title = expression(paste0(pref,fun[[1]]))
}, error = function() {
title = paste0(pref, as.character(quote(fun[[1]]), " function"))
})
}
p = p + labs(x = xlab, y = ylab, caption = title)
}
if(length(fun) > 1){
n= 1
colorList = c("red", "green", "blue")
for (i in fun[2:length(fun)]) {
p = p + ggplot2::stat_function(fun = i, color=colorList[n])
n=n+1
}
}
if(umx_set_plot_use_hrbrthemes(silent = TRUE)){
p = p + hrbrthemes::theme_ipsum()
} else {
p = p + cowplot::theme_cowplot(font_family = "Times", font_size = 12)
}
print(p)
invisible(p)
}
#' Compute odds ratio (OR)
#'
#' @description
#' Returns the odds in each group, and the odds ratio. Takes the cases (n) and total N as a list
#' of two numbers for each of two groups.
#'
#' @details
#' Returns a list of odds1, odds2, and OR + CI. Has a pretty-printing method so displays as:
#'
#' ```
#' Group 1 odds = 0.43
#' Group 2 odds = 0.11
#' OR = 3.86 CI95[0.160, 3.64]
#' ```
#'
#' @param grp1 either odds for group 1, or cases and total N , e.g c(n=3, N=10)
#' @param grp2 either odds for group 2, or cases and total N , e.g c(n=1, N=20)
#' @param alpha for CI (default = 0.05)
#' @return - List of odds in group 1 and group2, and the resulting OR and CI
#' @export
#' @family Miscellaneous Stats Functions
#' @seealso - [umx_r_test()]
#' @references - <https://stats.oarc.ucla.edu/r/dae/logit-regression/>, <https://tbates.github.io>
#' @md
#' @examples
#' oddsratio(grp1 = c(1, 10), grp2 = c(3, 10))
#' oddsratio(grp1 = 0.111, grp2 = 0.429)
#' oddsratio(grp1 = c(3, 10), grp2 = c(1, 10))
#' oddsratio(grp1 = c(3, 10), grp2 = c(1, 10), alpha = .01)
#'
oddsratio <- function(grp1= c(n=3, N=10), grp2= c(n=1, N=10), alpha = 0.05) {
# nGrp1 = 1; nGrp2 = 2; NGrp1 = 3; NGrp2 = 4
if(length(grp1) == 2){
nGrp1 = grp1[1]
NGrp1 = grp1[2]
# odds = n/(N-n)
odds1 = nGrp1/(NGrp1-nGrp1)
} else {
odds1 = grp1
}
if(length(grp2) == 2){
nGrp2 = grp2[1]
NGrp2 = grp2[2]
odds2 = nGrp2/(NGrp2-nGrp2)
} else {
odds2 = grp2
}
OR = odds1/odds2
if(length(grp1) == 2 & length(grp1) == 2){
siglog = sqrt((1/nGrp1) + (1/NGrp1) + (1/nGrp2) + (1/NGrp2))
zalph = qnorm(1 - alpha/2)
LowerCI = exp(log(OR) - zalph * siglog)
UpperCI = exp(log(OR) + zalph * siglog)
}else{
LowerCI = NA
UpperCI = NA
}
# CI
result = list(odds1= odds1, odds2= odds2, OR= OR, LowerCI = LowerCI, UpperCI = UpperCI, alpha = alpha)
class(result) = 'oddsratio'
result
}
#' Print a scale "oddsratio" object
#'
#' Print method for the [umx::oddsratio()] function.
#'
#' @param x A [umx::oddsratio()] result.
#' @param digits The rounding precision.
#' @param ... further arguments passed to or from other methods.
#' @return - invisible oddsratio object (x).
#' @seealso - [print()], [umx::oddsratio()],
#' @md
#' @method print oddsratio
#' @export
#' @examples
#' oddsratio(grp1 = c(1, 10), grp2 = c(3, 10))
#' oddsratio(grp1 = c(3, 10), grp2 = c(1, 10))
#' oddsratio(grp1 = c(3, 10), grp2 = c(1, 10), alpha = .01)
print.oddsratio <- function(x, digits = 3, ...) {
# x = list(odds1= odds1, odds2= odds2, OR= OR, LowerCI = OR_CI_lo, UpperCI = OR_CI_hi, alpha = alpha)
charLen = nchar("Group 1 odds")
cat(sprintf(paste0("%", charLen, "s = ", format(x[["odds1"]], digits = digits)), "Group 1 odds"), fill= TRUE)
cat(sprintf(paste0("%", charLen, "s = ", format(x[["odds2"]], digits = digits)), "Group 2 odds"), fill= TRUE)
OR = round(x[[ "OR" ]], digits = digits)
if(is.na(x[[ "UpperCI" ]])){
ORstring = paste0(OR, " (input odds as c(n=, N=) to compute CI)")
} else {
LowerCI = round(x[[ "LowerCI" ]], digits=digits)
UpperCI = round(x[[ "UpperCI" ]], digits=digits)
alpha = x[[ "alpha" ]]
ORstring = paste0(OR, " CI", 100-(alpha*100), " [", LowerCI, ", ", UpperCI, "]")
}
cat(sprintf(paste0("%", charLen, "s = ", ORstring), "OR"), fill= TRUE)
cat("polite note: Remember OR is sensitive to base rate: A given odds ratio can represent very different degrees of association/correlation")
invisible(x)
}
#' Report correlations and their p-values
#'
#' For reporting correlations and their p-values in a compact table. Handles rounding, and skipping non-numeric columns.
#'
#' To compute heterochoric correlations, see [umxHetCor()].
#'
#' *note*: The Hmisc package has a more robust function called `rcorr`.
#'
#' @param X a matrix or dataframe
#' @param df the degrees of freedom for the test
#' @param use how to handle missing data (defaults to pairwise complete)
#' @param digits rounding of answers
#' @param type Unused argument for future directions
#' @return - Matrix of correlations and p-values
#' @seealso umxHetCor
#' @family Miscellaneous Stats Functions
#' @export
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' tmp = myFADataRaw[1:8,1:8]
#' umx_cor(tmp)
#' tmp$x1 = letters[1:8] # make one column non-numeric
#' umx_cor(tmp)
umx_cor <- function (X, df = nrow(X) - 2, use = c("pairwise.complete.obs", "complete.obs", "everything", "all.obs", "na.or.complete"), digits = 2, type= c("r and p-value", "smart")) {
# see also
# Hmisc::rcorr()
use = match.arg(use)
message("TODO: umx_cor assumes no missing data, n is just nrow() !!")
# nVar = dim(x)[2]
# nMatrix = diag(NA, nrow= nVar)
# for (i in 1:nVar) {
# x[,i]
# }
numericCols = rep(FALSE, ncol(X))
for (i in 1:ncol(X)) {
numericCols[i] = is.numeric(X[,i])
}
if(ncol(X) > sum(numericCols)){
message("dropped ", ncol(X) - sum(numericCols), " non-numeric column(s).")
}
R = cor(X[,numericCols], use = use)
above = upper.tri(R)
below = lower.tri(R)
r2 = R[above]^2
Fstat = r2 * df/(1 - r2)
R[row(R) == col(R)] <- NA # NA on the diagonal
R[above] = pf(Fstat, 1, df, lower.tail = FALSE)
R[below] = round(R[below], digits)
R[above] = round(R[above], digits)
# R[above] = paste("p=",round(R[above], digits))
message("lower tri = correlation; upper tri = p-value")
return(R)
}
# Return the maximum value in a row
rowMax <- function(df, na.rm = TRUE) {
tmp = apply(df, MARGIN = 1, FUN = max, na.rm = na.rm)
tmp[!is.finite(tmp)] = NA
return(tmp)
}
rowMin <- function(df, na.rm= TRUE) {
tmp = apply(df, MARGIN = 1, FUN = min, na.rm = na.rm)
tmp[!is.finite(tmp)] = NA
return(tmp)
}
#' umx_round
#'
#' A version of round() which works on dataframes that contain non-numeric data (or data that cannot be coerced to numeric)
#' Helpful for dealing with table output that mixes numeric and string types.
#'
#' @param df a dataframe to round in
#' @param digits how many digits to round to (defaults to getOption("digits"))
#' @param coerce whether to make the column numeric if it is not (default = FALSE)
#' @return - [OpenMx::mxModel()]
#' @family Miscellaneous Stats Functions
#' @export
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' head(umx_round(mtcars, coerce = FALSE))
#' head(umx_round(mtcars, coerce = TRUE))
#'
umx_round <- function(df, digits = getOption("digits"), coerce = FALSE) {
if(is.matrix(df)){
df = data.frame(df)
}
if(!is.data.frame(df)){
if(is.null(dim(df))){
if(coerce){
return(round(as.numeric(df), digits))
}else{
return(round(df, digits))
}
} else {
stop("df input for umx_round must be a dataframe")
}
}
# for each column, if numeric, round
rows = dim(df)[1]
cols = dim(df)[2]
for(c in 1:cols) {
if(coerce){
for(r in 1:rows) {
df[r, c] = round(as.numeric(df[r, c]), digits)
}
} else {
if(is.numeric(df[1, c])){
df[ , c] = round(df[ , c], digits)
}
}
}
return(df)
}
#' Compute an SE from a beta and p value
#'
#' @description
#' `SE_from_p` takes beta and p, and returns an SE.
#'
#' @param beta The effect size
#' @param p The p-value for the effect
#' @param SE Standard error
#' @param lower Lower CI
#' @param upper Upper CI
#' @return - Standard error
#' @export
#' @family Miscellaneous Stats Functions
#' @seealso - [umxAPA()]
#' @md
#' @examples
#' SE_from_p(beta = .0020, p = .780)
#' SE_from_p(beta = .0020, p = .01)
#' SE_from_p(beta = .0020, SE = 0.01)
#' umxAPA(.0020, p = .01)
SE_from_p <- function(beta = NULL, p = NULL, SE = NULL, lower = NULL, upper = NULL) {
if(is.null(beta)){
stop("beta must be given")
}
if (!is.null(p)){ # compute SE from beta and p
beta_over_SE = -log(p) * (416 * log(p) + 717)/1000
SE = abs(beta/beta_over_SE) # 3 = 5/(5/3) a/(a/b) = b
return(c(SE = SE))
# p = exp(-0.717*(beta/SE) - 0.416*(beta/SE)^2)
# p = .780
# x = log(p)
# "Of all tyrannies, a tyranny sincerely exercised for the good of its victims may be
# the most oppressive. It would be better to live under robber barons than
# under omnipotent moral busybodies." C.S. Lewis.
}else{ # compute p from beta and CI or SE
if (is.null(SE)){
if(is.null(upper) || is.null(lower)){
stop("SE_from_p: upper and lower, or SE must be provided to compute p")
}
SE = (upper - lower)/(2*1.96)
}
z = beta/SE
p_value = exp(-0.717 * z - 0.416 * z^2)
return(c(p_value = p_value))
}
}
specify_decimal <- function(x, k){
format(round(x, k), nsmall = k)
}
#' Report coefficient alpha (reliability)
#'
#' Compute and report Coefficient alpha (extracted from Rcmdr to avoid its dependencies)
#'
#' @param S A square, symmetric, numeric covariance matrix
#' @return None
#' @export
#' @family Miscellaneous Stats Functions
#' @seealso - [umx::print.reliability()],
#' @references - <https://cran.r-project.org/package=Rcmdr>
#' @examples
#' # treat car data as items of a test
#' data(mtcars)
#' reliability(cov(mtcars))
reliability <-function (S){
reliab <- function(S, R) {
k = dim(S)[1]
ones = rep(1, k)
v = as.vector(ones %*% S %*% ones)
alpha = (k/(k - 1)) * (1 - (1/v) * sum(diag(S)))
rbar = mean(R[lower.tri(R)])
std.alpha = k * rbar/(1 + (k - 1) * rbar)
c(alpha = alpha, std.alpha = std.alpha)
}
result = list()
if ((!is.numeric(S)) || !is.matrix(S) || (nrow(S) != ncol(S)) || any(abs(S - t(S)) > max(abs(S)) * 1e-10) || nrow(S) < 2)
stop("argument must be a square, symmetric, numeric covariance matrix")
k = dim(S)[1]
s = sqrt(diag(S))
R = S/(s %o% s)
rel = reliab(S, R)
result$alpha = rel[1]
result$st.alpha = rel[2]
if (k < 3) {
warning("there are fewer than 3 items in the scale")
return(invisible(NULL))
}
rel = matrix(0, k, 3)
for (i in 1:k) {
rel[i, c(1, 2)] = reliab(S[-i, -i], R[-i, -i])
a = rep(0, k)
b = rep(1, k)
a[i] = 1
b[i] = 0
cov = a %*% S %*% b
var = b %*% S %*% b
rel[i, 3] = cov/(sqrt(var * S[i, i]))
}
rownames(rel) = rownames(S)
colnames(rel) = c("Alpha", "Std.Alpha", "r(item, total)")
result$rel.matrix = rel
class(result) = "reliability"
result
}
#' Print a scale "reliability" object
#'
#' Print method for the [umx::reliability()] function.
#'
#' @param x A [umx::reliability()] result.
#' @param digits The rounding precision.
#' @param ... further arguments passed to or from other methods
#' @return - invisible reliability object (x)
#' @seealso - [print()], [umx::reliability()],
#' @md
#' @method print reliability
#' @export
#' @examples
#' # treat vehicle aspects as items of a test
#' data(mtcars)
#' reliability(cov(mtcars))
print.reliability <- function (x, digits = 4, ...){
cat(paste("Alpha reliability = ", round(x$alpha, digits), "\n"))
cat(paste("Standardized alpha = ", round(x$st.alpha, digits), "\n"))
cat("\nReliability deleting each item in turn:\n")
print(round(x$rel.matrix, digits))
invisible(x)
}
#' Convert Radians to Degrees
#'
#' @description Just a helper to multiply radians by 180 and divide by \eqn{\pi} to get degrees.
#'
#' *note*: R's trig functions, e.g. [sin()] use Radians for input!
#'
#' There are \eqn{2\pi} radians in a circle.
#' 1 Rad = \eqn{180/\pi} degrees = ~ 57.296 degrees.
#'
#' @param rad The value in Radians you wish to convert
#' @return - value in degrees
#' @export
#' @family Miscellaneous Functions
#' @seealso - [deg2rad()], [sin()]
#' @references [https://en.wikipedia.org/wiki/Radian](https://en.wikipedia.org/wiki/Radian)
#' @md
#' @examples
# Has test
#' rad2deg(pi) #180 degrees
rad2deg <- function(rad) { rad * 180/pi }
#' Convert Degrees to Degrees
#'
#' @description A helper to convert degrees (360 in a circle) to Rad (\eqn{2\pi}
#' in a circle).
#'
#' *note*: R's trig functions, e.g. [sin()] use Radians for input!
#'
#' The formula is radians = \eqn{deg x 180/\pi}.
#'
#' * 180 Degrees is equal to \eqn{\pi} radians.
#' * 1 Rad = \eqn{180/\pi} degrees = ~ 57.296 degrees.
#'
#' @param deg The value in degrees you wish to convert to radians
#' @return - value in radians
#' @export
#' @family Miscellaneous Functions
#' @seealso - [rad2deg()], [sin()]
#' @references [https://en.wikipedia.org/wiki/Radian](https://en.wikipedia.org/wiki/Radian)
#' @md
#' @examples
# Has test
#' deg2rad(180) == pi # TRUE!
deg2rad <- function(deg) { deg * pi/ 180 }
# =======================
# = Developer functions =
# =======================
#' Install OpenMx, with choice of builds
#'
#' @description
#' You can install OpenMx, including the latest NPSOL-enabled build of OpenMx. Options are:
#'
#' 1. "NPSOL": Install from our repository (default): This is where we maintain binaries supporting parallel processing and NPSOL.
#' 2. "travis": Install the latest travis built (MacOS only).
#' 3. "CRAN": Install from CRAN.
#' 4. "open travis build page": Open the list of travis builds in a browser window.
#'
#'
#' @aliases umx_update_OpenMx
#' @param loc Version to get default is "NPSOL". "travis" (latest build),CRAN, list of builds.
#' @param url Custom URL. On Mac, set this to "Finder" and the package selected in the Finder will be installed.
#' @param repos Which repository to use (ignored currently).
#' @param lib Where to install the package.
#' @return None
#' @export
#' @seealso [umxVersion()]
#' @family Miscellaneous Utility Functions
#' @references - <https://github.com/tbates/umx>, <https://tbates.github.io>
#' @md
#' @examples
#' \dontrun{
#' install.OpenMx() # gets the NPSOL version
#' install.OpenMx("NPSOL") # gets the NPSOL version explicitly
#' install.OpenMx("CRAN") # Get the latest CRAN version
#' install.OpenMx("open travis build page") # Open web page of travis builds
#' }
install.OpenMx <- function(loc = c("NPSOL", "travis", "CRAN", "open travis build page", "UVa"), url= NULL, lib, repos = getOption("repos")) {
loc = match.arg(loc)
if(loc == "UVa"){
loc = "NPSOL"
message("next time, use 'NPSOL' instead of 'UVa'")
}
oldTimeOut = getOption('timeout')
options(timeout=60*3)
if(!is.null(url)){
if(url == "Finder"){
umx_check_OS("OSX")
url = system(intern = TRUE, "osascript -e 'tell application \"Finder\" to get the POSIX path of (selection as alias)'")
message("Using file selected in front-most Finder window:", url)
} else if(url == "") {
url = file.choose(new = FALSE) # choose a file
message("Using selected file:", url)
}
install.packages(url)
} else if(loc == "NPSOL"){
if(umx_check_OS("Windows")){
detach('package:OpenMx', unload = TRUE)
}
source("https://openmx.ssri.psu.edu/getOpenMx.R")
# was source("https://openmx.ssri.psu.edu/software/getOpenMx.R")
# was https://openmx.psyc.virginia.edu/getOpenMx.R
# was source("https://openmx.ssri.psu.edu/software/getOpenMx.R")
}else if(loc == "travis"){
if(umx_check_OS("OSX")){
install.packages("https://vipbg.vcu.edu/vipbg/OpenMx2/software/bin/macosx/travis/OpenMx_latest.tgz")
# was ("https://openmx.psyc.virginia.edu/OpenMx2/bin/macosx/travis/OpenMx_latest.tgz")
# , lib = lib, repos=repos
# quit(save = "default")
} else {
stop(paste0("Sorry, travis builds are only available for MacOS :-("))
}
} else if(loc == "CRAN"){
install.packages("OpenMx", lib= lib, repos = repos)
} else if(loc == "open travis build page"){
browseURL("https://vipbg.vcu.edu/vipbg/OpenMx2/software/bin/macosx/travis/?C=M;O=D")
}
options(timeout=oldTimeOut)
}
#' @export
umx_update_OpenMx <- install.OpenMx
#' "make" the umx package using devtools: release to CRAN etc.
#'
#' @description
#' Easily run devtools "install", "release", "win", "examples" etc.
#'
#' @param what whether to "install", "release" to CRAN, "test", "check" test on "win" or "rhub", "spell", or "examples")).
#' @param pkg the local path to your package. Defaults to my path to umx.
#' @param check Whether to run check on the package before release (default = TRUE).
#' @param run If what is "examples", whether to also run examples marked don't run. (default FALSE).
#' @param start If what is "examples", which function to start from (default (NULL) = beginning).
#' @param spelling Whether to check spelling before release (default = "en_US": set NULL to not check).
#' @param which What rhub platform to use? c("mac", "linux", "win").
#' @param run_dont_test When checking.
#' @param spell for rhub, check spelling? TRUE
#' @return None
#' @export
#' @family xmu internal not for end user
#' @references - <https://devtools.r-lib.org>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' # umx_make() # Just load new code (don't rebuild help etc)
#' # umx_make(what = "quickInst") # Quick install
#' # umx_make(what = "install") # Full package rebuild and install
#' # umx_make(what = "spell") # Spellcheck Rd documents
#' # umx_make(what = "sitrep") # Are needed packages up to date?
#' # umx_make(what = "deps_install") # Update needed packages
#' # umx_make(what = "examples") # Run the examples
#' # umx_make(what = "checkCRAN") # Run R CMD check
#' # umx_make(what = "rhub") # Check on rhub
#' # umx_make(what = "win") # Check on win-builder
#' # umx_make(what = "release") # Release to CRAN
#' # tmp = umx_make(what = "lastRhub") # View rhub result
#' }
umx_make <- function(what = c("load", "quickInst", "install", "spell", "sitrep", "deps_install", "checkCRAN", "testthat", "examples", "win", "rhub", "lastRhub", "release"), pkg = "~/bin/umx", check = TRUE, run = FALSE, start = NULL, spelling = "en_US", which = c("win", "mac", "linux", "solaris"), run_dont_test = FALSE, spell=TRUE) {
what = match.arg(what)
which = match.arg(which)
if(what == "load"){
devtools::load_all(path = pkg)
changed = gert::git_status(repo = "~/bin/umx")
if(dim(changed)[1]>=1){
umx_print(gert::git_status(repo = pkg))
}
} else if(what == "quickInst"){
devtools::document(pkg = pkg);
devtools::install(pkg = pkg, quick = TRUE, dependencies= FALSE, upgrade= FALSE, build_vignettes = FALSE);
devtools::load_all(path = pkg)
changed = gert::git_status(repo = "~/bin/umx")
if(dim(changed)[1]>=1){
umx_print(gert::git_status(repo = pkg))
}
}else if(what == "install"){
devtools::document(pkg = pkg);
devtools::install(pkg = pkg);
devtools::load_all(path = pkg)
} else if (what == "spell"){
spelling::spell_check_package(pkg = pkg, vignettes = TRUE, use_wordlist = TRUE)
# }else if (what=="travisCI"){
# browseURL("https://www.travis-ci.com/tbates/umx")
}else if (what == "sitrep"){
devtools::dev_sitrep(pkg = pkg)
}else if (what == "deps_install"){
devtools::install_dev_deps(pkg=pkg)
} else if(what == "run_examples"){
devtools::run_examples(pkg = pkg, run = run, start = start)
} else if(what == "checkCRAN"){
devtools::check(pkg = pkg, run_dont_test = run_dont_test) # http://r-pkgs.had.co.nz/check.html
} else if (what =="win"){
devtools::check_win_devel(pkg = pkg)
} else if (what =="rhub"){
if(which == "mac"){
plat = "macos-highsierra-release-cran"
plat = "macos-m1-bigsur-release"
} else if(which=="linux") {
plat = "debian-gcc-patched"
# plat = "debian-clang-devel"
} else if(which=="win") {
plat = "windows-x86_64-devel" #"windows-x86_64-patched"
# plat = "windows-x86_64-devel" # broken 2021-06-12
} else if (which=="solaris"){
plat = "solaris-x86-patched-ods"
}
cat("checking ", omxQuotes(pkg), "on", omxQuotes(plat))
if(spell){
devtools::check_rhub(pkg = pkg, platforms = plat, interactive = FALSE, env_vars = c(`_R_CHECK_FORCE_SUGGESTS_` = "true", `_R_CHECK_CRAN_INCOMING_USE_ASPELL_` = "false"))
}else{
devtools::check_rhub(pkg = pkg, platforms = plat, interactive = FALSE)
}
} else if(what == "lastRhub"){
prev = rhub::list_package_checks(package = pkg, howmany = 4)
check_id = prev$id[1]
return(rhub::get_check(check_id))
}else if(what == "testthat"){
devtools::test(pkg = pkg)
} else if (what == "release"){
oldDir = getwd()
setwd(dir= pkg)
devtools::release(pkg = pkg, check = check, "--no-manual") # spelling = NULL
setwd(dir= oldDir)
}else{
stop("I don't know how to ", what)
}
}
# ================================
# = Reporting & Graphing helpers =
# ================================
#' Print the name and compact contents of variable.
#'
#' Helper function to ease debugging with console notes like: "ObjectName = \<Object Value\>".
#' This is primarily useful for inline debugging, where seeing, e.g., "nVar = 3" can be useful.
#' The ability to say \code{umx_msg(nVar)} makes this easy.
#'
#' @param x the thing you want to pretty-print
#' @return - NULL
#' @export
#' @family Miscellaneous Utility Functions
#' @references - <https://tbates.github.io>, <https://github.com/tbates/umx>
#' @md
#' @examples
#' a = "brian"
#' umx_msg(a)
#' b = c("brian", "sally", "jane")
#' umx_msg(b)
#' umx_msg(mtcars)
umx_msg <- function(x) {
nm = deparse(substitute(x) )
if(any(is.data.frame(x))){
message(nm, " = ")
str(x)
} else {
if(length(x) > 1) {
message(nm, " = ", omxQuotes(x))
} else {
message(nm, " = ", x)
}
}
}
#' Helper to make the list of vars and their shapes for a graphviz string
#'
#' @description
#' Helper to make a graphviz rank string defining the latent, manifest, and means and their shapes
#'
#' @param latents list of latent variables (including "one")
#' @param manifests list of manifest variables
#' @param preOut existing output string (pasted in front of this: "" by default).
#' @return string
#' @export
#' @family Graphviz
#' @seealso - [xmu_dot_rank()]
#' @examples
#' xmu_dot_define_shapes(c("as1"), c("E", "N"))
xmu_dot_define_shapes <- function(latents, manifests, preOut= "") {
latents = unique(latents)
manifests = unique(manifests)
preOut = paste0(preOut, "\n# Latents\n")
for(var in latents) {
if(var == "one"){
preOut = paste0(preOut, "\t", var, " [shape = triangle];\n")
} else {
preOut = paste0(preOut, "\t", var, " [shape = circle];\n")
}
}
preOut = paste0(preOut, "\n# Manifests\n")
for(thisManifest in manifests) {
preOut = paste0(preOut, "\t", thisManifest, " [shape = square];\n")
}
return(preOut)
}
#' Helper to make a graphviz rank string
#'
#' Given a list of names, this filters the list, and returns a graphviz string to force them into the given rank.
#' e.g. `"{rank=same; as1};"`
#'
#' @param vars a list of strings
#' @param pattern regular expression to filter vars
#' @param rank "same", "max", "min"
#' @return string
#' @export
#' @family Graphviz
#' @seealso - [xmu_dot_define_shapes()]
#' @md
#' @examples
#' xmu_dot_rank(c("as1"), "^[ace]s[0-9]+$", "same")
xmu_dot_rank <- function(vars, pattern, rank) {
formatted = paste(namez(vars, pattern), collapse = "; ")
ranks = paste0("{rank=", rank, "; ", formatted, "};\n")
return(ranks)
}
#' Return dot code for paths in a matrix
#'
#' @description
#' Return dot code for paths in a matrix is a function which walks the rows and cols of a matrix.
#' At each free cell, it creates a dot-string specifying the relevant path, e.g.:
#'
#' \code{ai1 -> var1 [label=".35"]}
#'
#' Its main use is to correctly generate paths (and their sources and sink objects)
#' without depending on the label of the parameter.
#'
#' It is highly customizable:
#'
#' 1. You can specify which cells to inspect, e.g. "lower".
#' 2. You can choose how to interpret path direction, from = "cols".
#' 3. You can choose the label for the from to ends of the path (by default, the matrix name is used).
#' 4. Offer up a list of from and toLabel which will be indexed into for source and sink
#' 5. You can set the number of arrows on a path (e.g. both).
#' 6. If `type` is set, then sources and sinks added manifests and/or latents output (p)
#'
#' Finally, you can pass in previous output and new paths will be concatenated to these.
#'
#' @param x a [umxMatrix()] to make paths from.
#' @param from one of "rows", "columns"
#' @param cells which cells to process: "any" (default), "diag", "lower", "upper". "left" is the left half (e.g. in a twin means matrix)
#' @param arrows "forward" "both" or "back"
#' @param fromLabel = NULL. NULL = use matrix name (default). If one, if suffixed with index, length() > 1, index into list. "one" is special.
#' @param toLabel = NULL. NULL = use matrix name (default). If one, if suffixed with index, length() > 1, index into list.
#' @param showFixed = FALSE.
#' @param digits to round values to (default = 2).
#' @param fromType one of "latent" or "manifest" NULL (default) = don't accumulate new names.
#' @param toType one of "latent" or "manifest" NULL (default) = don't accumulate new names.
#' @param model If you want to get CIs, you can pass in the model (default = NULL).
#' @param SEstyle If TRUE, CIs shown as "b(SE)" ("b \[l,h\]" if FALSE (default)). Ignored if model NULL.
#' @param p input to build on. list(str = "", latents = c(), manifests = c())
#' @return - list(str = "", latents = c(), manifests = c())
#' @export
#' @family Graphviz
#' @seealso - [plot()]
#' @md
#' @examples
#'
#' # test with a 1 * 1
#' a_cp = umxMatrix("a_cp", "Lower", 1, 1, free = TRUE, values = pi)
#' out = xmu_dot_mat2dot(a_cp, cells = "lower_inc", from = "cols", arrows = "both")
#' cat(out$str) # a_cp -> a_cp [dir = both label="2"];
#' out = xmu_dot_mat2dot(a_cp, cells = "lower_inc", from = "cols", arrows = "forward",
#' fromLabel = "fromMe", toLabel = "toYou",
#' fromType = "latent", toType = "manifest", digits = 3, SEstyle = TRUE
#' )
#' cat(out$str) # fromMe -> toYou [dir = forward label="3.142"];
#' cat(out$latent) # fromMe
#' cat(out$manifest) # toYou
#'
#' # Make a lower 3 * 3 value= 1:6 (1, 4, 6 on the diag)
#' a_cp = umxMatrix("a_cp", "Lower", 3, 3, free = TRUE, values = 1:6)
#'
#' # Get dot strings for lower triangle (default from and to based on row and column number)
#' out = xmu_dot_mat2dot(a_cp, cells = "lower", from = "cols", arrows = "both")
#' cat(out$str) # a_cp1 -> a_cp2 [dir = both label="2"];
#'
#' # one arrow (the default = "forward")
#' out = xmu_dot_mat2dot(a_cp, cells = "lower", from = "cols")
#' cat(out$str) # a_cp1 -> a_cp2 [dir = forward label="2"];
#'
#' # label to (rows) using var names
#'
#' out = xmu_dot_mat2dot(a_cp, toLabel= paste0("v", 1:3), cells = "lower", from = "cols")
#' umx_msg(out$str) # a_cp1 -> v2 [dir = forward label="2"] ...
#'
#' # First call also inits the plot struct
#' out = xmu_dot_mat2dot(a_cp, from = "rows", cells = "lower", arrows = "both", fromType = "latent")
#' out = xmu_dot_mat2dot(a_cp, from = "rows", cells = "diag",
#' toLabel= "common", toType = "manifest", p = out)
#' umx_msg(out$str); umx_msg(out$manifests); umx_msg(out$latents)
#'
#' # ================================
#' # = Add found sinks to manifests =
#' # ================================
#' out = xmu_dot_mat2dot(a_cp, from= "rows", cells= "diag",
#' toLabel= c('a','b','c'), toType= "manifest");
#' umx_msg(out$manifests)
#'
#' # ================================
#' # = Add found sources to latents =
#' # ================================
#' out = xmu_dot_mat2dot(a_cp, from= "rows", cells= "diag",
#' toLabel= c('a','b','c'), fromType= "latent");
#' umx_msg(out$latents)
#'
#'
#' # ========================
#' # = Label a means matrix =
#' # ========================
#'
#' tmp = umxMatrix("expMean", "Full", 1, 4, free = TRUE, values = 1:4)
#' out = xmu_dot_mat2dot(tmp, cells = "left", from = "rows",
#' fromLabel= "one", toLabel= c("v1", "v2")
#' )
#' cat(out$str)
#'
#' \dontrun{
#' # ==============================================
#' # = Get a string which includes CI information =
#' # ==============================================
#' data(demoOneFactor)
#' latents = c("g"); manifests = names(demoOneFactor)
#' m1 = umxRAM("xmu_dot", data = demoOneFactor, type = "cov",
#' umxPath(latents, to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = latents, fixedAt = 1.0)
#' )
#' m1 = umxCI(m1, run= "yes")
#' out = xmu_dot_mat2dot(m1$A, from = "cols", cells = "any",
#' toLabel= paste0("x", 1:5), fromType = "latent", model= m1);
#' umx_msg(out$str); umx_msg(out$latents)
#'
#' }
#'
xmu_dot_mat2dot <- function(x, cells = c("diag", "lower", "lower_inc", "upper", "upper_inc", "any", "left"), from = c("rows", "cols"), fromLabel = NULL, toLabel = NULL, showFixed = FALSE, arrows = c("forward", "both", "back"), fromType = NULL, toType = NULL, digits = 2, model = NULL, SEstyle = FALSE, p = list(str = "", latents = c(), manifests = c())) {
from = match.arg(from)
cells = match.arg(cells)
arrows = match.arg(arrows)
# Get default from and to labels if custom not set
if(is.null(fromLabel)){ fromLabel = x$name }
if(is.null(toLabel)) { toLabel = x$name }
if(inherits(x, "MxAlgebra")){
# convert to a matrix
tmp = x$result
x = umxMatrix(x$name, "Full", dim(tmp)[1], dim(tmp)[2], free = TRUE, values = tmp)
}
nRows = nrow(x)
nCols = ncol(x)
# Get parameter value and make the plot string
# Convert address to [] address and look for a CI: not perfect, as CI might be label based?
# Also fails to understand not using _std?
for (r in 1:nRows) {
for (c in 1:nCols) {
if(xmu_cell_is_on(r= r, c = c, where = cells, mat = x)){
# cell is in the target zone
if(!is.null(model)){
# Model available - look for CIs by label...
CIstr = xmu_get_CI(model, label = x$labels[r,c], SEstyle = SEstyle, digits = digits)
if(is.na(CIstr)){
value = umx_round(x$values[r,c], digits)
}else{
value = umx_round(CIstr, digits)
}
} else {
if(is.numeric(x$values[r,c])){
value = umx_round(x$values[r,c], digits)
} else {
value = x$values[r,c]
}
# tryCatch({
# value = umx_round(x$values[r,c], digits)
# }, warning = function() {
# }, error = function() {
# }, finally={
# value = x$values[r,c]
# })
# #
}
if(from == "rows"){
sourceIndex = r; sinkIndex = c; fromWidth = nRows; toWidth = nCols
} else { # from cols
sourceIndex = c; sinkIndex = r; fromWidth = nCols; toWidth = nRows
}
if(length(fromLabel) == 1){
if(fromLabel == "one"){
thisFrom = "one"
} else if(fromWidth > 1){
thisFrom = paste0(fromLabel, sourceIndex)
}else{
thisFrom = fromLabel[sourceIndex]
}
} else {
thisFrom = fromLabel[sourceIndex]
}
if(length(toLabel) == 1){
if(toLabel == "one"){
thisTo = "one"
} else if(toWidth > 1){
thisTo = paste0(toLabel, sinkIndex)
}else{
thisTo = toLabel[sinkIndex]
}
} else {
thisTo = toLabel[sinkIndex]
}
# Show fixed cells if non-0
if(x$free[r,c] || (showFixed && x$values[r,c] != 0)){
p$str = paste0(p$str, "\n", thisFrom, " -> ", thisTo, " [dir = ", arrows, " label=\"", value, "\"];")
if(!is.null(fromType)){
if(fromType == "latent"){
p$latents = c(p$latents, thisFrom)
} else if(fromType == "manifest"){
p$manifests = c(p$manifests, thisFrom)
}else{
stop("not sure what to do for fromType = ", fromType, ". Legal is latent or manifest")
}
}
if(!is.null(toType)){
if(toType == "latent"){
p$latents = c(p$latents, thisTo)
} else if(toType == "manifest"){
p$manifests = c(p$manifests, thisTo)
}else{
stop("not sure what to do for fromType = ", toType, ". Legal is latent or manifest")
}
}
}
} else {
# fixed cell
}
}
}
p$latents = unique(p$latents)
p$manifests = unique(p$manifests)
return(p)
}
#' umx_time
#'
#' A function to compactly report how long a model took to execute. Comes with some preset styles
#' User can set the format with C-style string formatting.
#'
#' The default time format is "simple", which gives only the biggest unit used. i.e., "x seconds" for times under 1 minute.
#' "std" shows time in the format adopted in OpenMx 2.0 e.g. "Wall clock time (HH:MM:SS.hh): 00:00:01.16"
#'
#' If a list of models is provided, time deltas will also be reported.
#'
#' If instead of a model the key word "start" is given in x, a start time will be recorded. "stop" gives the
#' time since "start" was called (and clears the timer)
#'
#' If a model has not been run, umx_time will run it for you.
#'
#' @param x A [OpenMx::mxModel()] or list of models for which to display elapsed time, or 'start' or 'stop'
#' @param formatStr A format string, defining how to show the time (defaults to human readable)
#' @param tz time zone in which the model was executed (defaults to "GMT")
#' @param autoRun If TRUE (default), run the model if it appears not to have been.
#' @return - invisible time string
#' @export
#' @family Reporting Functions
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
#' require(umx)
#' umx_time('stop') # alert user stop called when not yet started...
#' umx_time('stop')
#' umx_time('start')
#' data(demoOneFactor)
#' latents = c("G")
#' manifests = names(demoOneFactor)
#' myData = mxData(cov(demoOneFactor), type = "cov", numObs=500)
#' m1 = umxRAM("umx_time_example", data = myData,
#' umxPath(from = latents, to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = latents, fixedAt = 1)
#' )
#' umx_time(m1) # report time from mxModel
#' m2 = umxRun(m1)
#' umx_time(c(m1, m2)) # print comparison table
#' umx_time('stop') # report the time since timer last started, and restart
#' umx_time('stop') # report the time since timer was restarted.
#' }
#'
umx_time <- function(x = NA, formatStr = c("simple", "std", "custom %H %M %OS3"), tz = "GMT", autoRun = TRUE){
commaSep = paste0(umx_set_separator(silent = TRUE), " ")
formatStr = xmu_match.arg(formatStr, c("simple", "std", "custom %H %M %OS3"), check = FALSE)
if(is.list(x)){
# check each item is a model
if(!umx_is_MxModel(x, listOK = TRUE)){
stop("If x is a list of models, each must be a valid mxModel")
}
}else if(umx_is_MxModel(x)){
# great, we've got a model
if(!is.character(formatStr)){
stop(paste("In umx_time, pass models as a list, i.e., umx_time(c(m1, m2)), not umx_time(m1, m2)"))
}
}else if(is.character(x)){
umx_check(x %in% c('start', 'stop', "lap", "now"), "stop", "Valid time strings are 'start', 'stop', 'lap', (or a model or list of models). Leave blank for 'now'")
}else if(is.na(x)){
cat("Current time is ", format(Sys.time(), "%X, %a %d %b %Y"), "\nTry me with a list of models, or 'start', 'stop'")
invisible()
}else{
stop("You must set the first parameter to 'start', 'stop', 'now', a model, or a list of models.\nYou offered up a", class(x))
}
for(i in 1:length(x)) {
if(length(x) > 1) {
thisX = x[[i]]
} else {
if(inherits(x, "list")){
thisX = x[[i]]
} else {
thisX = x
}
}
if(inherits(thisX, "character")){
if(thisX == "start"){
options("umx_last_time" = proc.time())
thisTime = (proc.time()["elapsed"] - getOption("umx_last_time")["elapsed"])
# invisible()
} else if (thisX == "stop" || thisX == "lap" ) {
tmp = getOption("umx_last_time")
if(is.null(tmp)){
message("Timer was not yet started: I started it now...")
options("umx_last_time" = proc.time())
thisTime = (proc.time()["elapsed"] - getOption("umx_last_time")["elapsed"])
# invisible()
} else {
thisTime = (proc.time()["elapsed"] - getOption("umx_last_time")["elapsed"])
if(thisX =="stop"){
options("umx_last_time" = proc.time())
}
}
}else if(thisX == "now"){
return(format(Sys.time(), "%X, %a %d %b %Y"))
}
} else {
# Handle model
if(!umx_has_been_run(thisX) && autoRun){
thisX = mxRun(thisX)
# message("You must run the model before asking for the elapsed run time")
}
thisTime = thisX$output$wallTime
if(i == 1){
lastTime = thisTime
timeDelta = ""
percentDelta = ""
} else {
timeDelta = paste0("(\u0394: ", round(thisTime - lastTime, 3), ")")
percentDelta = paste0(round(as.numeric(thisTime) / as.numeric(lastTime)*100, 3), "%")
}
}
if(formatStr == "std"){
realFormatStr = "Wall clock time (HH:MM:SS.hh): %H:%M:%OS2"
} else if(formatStr == "simple"){
if(thisTime > (3600-1)){ # > 1 hour
realFormatStr = "%H hr, %M min, %OS2 sec"
} else if(thisTime > (60-1)){ # minutes
realFormatStr = "%M min, %OS2 sec"
} else { # < 1 minute
realFormatStr = "%OS2 sec"
}
}else{
realFormatStr = formatStr
}
if(inherits(thisX, "character")){
# Handle start-stop
timeString = format(.POSIXct(thisTime, tz), paste0("elapsed time: ", realFormatStr))
} else {
timeString = format(.POSIXct(thisTime, tz), paste0(thisX$name, ": ", realFormatStr, " ", timeDelta, commaSep, percentDelta))
}
message(timeString)
}
invisible(thisTime)
}
#' Print tables in a range of formats (markdown default, see [umx_set_table_format()] for other formats)
#' or as a web browser table.
#'
#' To aid interpretability of printed tables from OpenMx (and elsewhere)
#' you can change how NA and zero appear, and suppressing values below a certain cut-off.
#' By default, Zeros have the decimals suppressed, and NAs are suppressed altogether.
#'
#' @param x A data.frame to print (matrices will be coerced to data.frame)
#' @param digits The number of decimal places to print (getOption("digits"))
#' @param caption Optional caption.
#' @param file Whether to write to a file (defaults to NA (no file). Use "html" to open table in browser.
#' @param report How to report the results. "html" = open in browser.
#' @param kableExtra Whether to print the table using kableExtra (if report="html")
#' @param zero.print How to display 0 values (default = "0") for sparse tables, using "." can produce more readable results.
#' @param justify Parameter passed to print (defaults to "none")
#' @param na.print How to display NAs (default = "")
#' @param quote Whether or not to quote strings (FALSE)
#' @param suppress Minimum numeric value to print (NULL = print all values, no matter how small)
#' @param append If html, is this appended to file? (FALSE)
#' @param sortableDF If html, is table sortable? (TRUE)
#' @param both If html, is table also printed as markdown? (TRUE)
#' @param style The style for the table "paper","material_dark" etc.
#' @param bootstrap_options e.g. border etc.
#' @param lightable_options e.g. striped
#' @param html_font Override default font. e.g. "Times" or '"Arial Narrow", arial, helvetica, sans-s'
#' @param ... Optional parameters for print
#' @return - A dataframe of text
#' @seealso [umx_msg()], [umx_set_table_format()]
#' @family Miscellaneous Utility Functions
#' @export
#' @md
#' @examples
#' data(mtcars)
#' umx_print(mtcars[1:10,], digits = 2, zero.print = ".", justify = "left")
#' umx_print(mtcars[1,1:2], digits = 2, zero.print = "")
#' umx_print(mtcars[1,1:2], digits = 2, caption = "Hi: I'm the caption!")
#' \dontrun{
#' umx_print(mtcars[1:10,], report = "html")
#' }
umx_print <- function (x, digits = getOption("digits"), caption = NULL, report = c("markdown", "html"), file = c(NA, "tmp.html"), na.print = "", zero.print = "0", justify = "none", quote = FALSE, suppress = NULL, kableExtra = TRUE, append = FALSE, sortableDF= TRUE, html_font = NULL, style = c("paper","material_dark", "classic", "classic_2", "minimal", "material"), bootstrap_options=c("hover", "bordered", "condensed", "responsive"), lightable_options = "striped", both = TRUE, ...){
style = match.arg(style)
file = xmu_match.arg(file, c(NA, "tmp.html"), check = FALSE)
report = match.arg(report)
# Catch legacy users passing in file instead of report...
if(!is.na(file) && file == "markdown"){
cat("polite note: in future, please use report = 'markdown', not file='markdown'. file is for choosing a custom html file when writing to the browser with report= 'html'")
report = "markdown"
file = NA
}else if(!is.na(file) && file == "html"){
cat("polite note: in future, please use report = 'html', not file='html'. file is for choosing a custom html file when writing to the browser with report= 'html'")
report = "html"
file = "tmp.html"
}
!is.na(file) && file == "markdown"
if(class(x)[[1]] != "data.frame"){
if(class(x)[[1]] == "matrix" | class(x)[[1]] == "numeric"){
x = data.frame(x)
} else if(class(x)[[1]] == "tbl_df"){
x = data.frame(x)
} else {
message("Sorry, umx_print currently only prints data.frames, dreaded tibbles, matrices, and vectors.\n
File a request to print ", omxQuotes(class(x)[[1]]), " objects\n or perhaps you want umx_msg?")
return()
}
}
if(class(x)[[1]] == "character"){
print(x)
invisible(x)
}else if(is.null(dim(x)[1]) || dim(x)[1] == 0){
invisible(x)
} else {
x = umx_round(x, digits = digits, coerce = FALSE)
if(!is.null(suppress)){
# zero-out suppressed values
x[abs(x) < suppress] = 0
zero.print = "."
}
x[is.na(x)] = na.print
x[(x == 0)] = zero.print
if (is.numeric(x) || is.complex(x)){
print(x, quote = quote, right = TRUE, ...)
} else if(report == "html"){
# From report = "html"
if(both){ print(kable(x, caption= caption, format="pipe")) }
if(kableExtra){
# default html output
x = kbl(x, caption = caption, format=report)
if(zero.print != "0"){
# x = add_footnote(x, label = paste0("zero printed as ", omxQuotes(zero.print)))
x = footnote(kable_input= x, general = paste0("zero printed as ", omxQuotes(zero.print)))
}
x = xmu_style_kable(x, html_font = html_font, style = style, bootstrap_options= bootstrap_options, lightable_options = lightable_options, full_width = FALSE)
print(x)
} else {
# !kableExtra
R2HTML::HTML(x, file = file, Border = 0, append = append, sortableDF = sortableDF)
system(paste0("open ", file))
}
}else{
# markdown/latex
x = kbl(x, caption = caption, format=report)
print(x)
}
invisible(x)
}
} # end umx_print
# ===========================
# = Boolean check functions =
# ===========================
#' umx_has_been_run
#'
#' check if an mxModel has been run or not
#'
#' @param model The [OpenMx::mxModel()] you want to check has been run
#' @param stop Whether to stop if the model has not been run (defaults to FALSE)
#' @return - boolean
#' @export
#' @family Test
#' @references - <https://github.com/tbates/umx>
#' @md
#' @examples
#' \dontrun{
# has test
#' require(umx)
#' data(demoOneFactor)
#' manifests = names(demoOneFactor)
#
#' m1 = umxRAM("has_been_run_example", data = demoOneFactor, type = "cov",
#' umxPath("G", to = manifests),
#' umxPath(var = manifests),
#' umxPath(var = "G", fixedAt = 1)
#' )
#' umx_has_been_run(m1)
#' }
#'
umx_has_been_run <- function(model, stop = FALSE) {
output = model$output
if (is.null(output)){
if(stop){
stop("Provided model has no objective function, and thus no output to process further")
}else{
return(FALSE)
}
} else if (length(output) < 1){
if(stop){
stop("Provided model has no output (probably you have not yet run it?)")
} else {
return(FALSE)
}
}
return(TRUE)
}
umxCheckModel <- function(model){
# Are all the manifests in paths?
# Do the manifests have residuals?
if(any(duplicated(model@manifestVars))){
stop(paste("manifestVars contains duplicates:", duplicated(model@manifestVars)))
}
if(length(model@latentVars) == 0){
# Check none are duplicates, none in manifests
if(any(duplicated(model@latentVars))){
stop(paste("latentVars contains duplicates:", duplicated(model@latentVars)))
}
if(any(duplicated(c(model@manifestVars, model@latentVars)))){
stop(
paste("manifest and latent lists contain clashing names:", duplicated(c(model@manifestVars, model@latentVars)))
)
}
}
# Check manifests in dataframe
}
#' umx_check
#'
#' Check that a test evaluates to TRUE. If not, stop, warn, or message the user
#'
#' @param boolean.test test evaluating to TRUE or FALSE.
#' @param action One of "stop" (the default), "warning", or "message".
#' @param message what to tell the user when boolean.test is FALSE.
#' @param ... extra text will be pasted after the messages.
#' @return - boolean
#' @export
#' @family Test
#' @examples
#' umx_check(length(1:3)==3, "message", "item must have length == 3", "another comment", "and another")
#' umx_check(1==2, "message", "one must be 2", ". Another comment", "and another")
umx_check <- function(boolean.test, action = c("stop", "warning", "message"), message = "check failed", ...){
dot.items = list(...) # grab all the dot items: mxPaths, etc...
dot.items = paste(unlist(dot.items), collapse =" ") # concatenate any dot items
action = match.arg(action)
if(!boolean.test){
if(action == "stop"){
stop(paste0(message, dot.items), call. = FALSE)
} else if(action == "warning"){
warning(paste0(message, dot.items), call. = FALSE)
}else{
message(paste0("Polite message: ", message, dot.items))
}
}
return(boolean.test)
}
#' Check if a request name exists in a dataframe or related object
#'
#' Check if a list of names are in the [namez()] of a dataframe (or the [dimnames()] of a matrix), or the names of
#' the observed data of an [mzData()]
#'
#' @param namesNeeded Variable names to find (a dataframe is also allowed)
#' @param data data.frame, matrix, or mxData to search in for names (default NA)
#' @param die Whether to die if the check fails (default TRUE).
#' @param illegal Optional list of names which must NOT be present.
#' @param no_others Whether to test that the data contain no columns in addition to those in namesNeeded (default FALSE)
#' @param intersection Show the intersection of names
#' @param message Some helpful text to append when dieing.
#' @family Test
#' @export
#' @family Check or test
#' @references - <https://github.com/tbates/umx>
#' @examples
#' require(umx)
#' data(demoOneFactor) # "x1" "x2" "x3" "x4" "x5"
#' umx_check_names(c("x1", "x2"), demoOneFactor)
#' umx_check_names(c("x1", "x2"), as.matrix(demoOneFactor))
#' umx_check_names(c("x1", "x2"), cov(demoOneFactor[, c("x1","x2")]))
#' umx_check_names(c("x1", "x2"), mxData(demoOneFactor, type="raw"))
#' umx_check_names(c("z1", "x2"), data = demoOneFactor, die = FALSE)
#' umx_check_names(c("x1", "x2"), data = demoOneFactor, die = FALSE, no_others = TRUE)
#' umx_check_names(c("x1","x2","x3","x4","x5"), data = demoOneFactor, die = FALSE, no_others = TRUE)
#' # no request
#' umx_check_names(c(), data = demoOneFactor, die = FALSE, no_others = TRUE)
#'
#' \dontrun{
#' # An example error from vars that don't exist in the data
#' umx_check_names(c("bad_var_name", "x2"), data = demoOneFactor, die = TRUE)
#' }
umx_check_names <- function(namesNeeded, data = NA, die = TRUE, illegal = NULL, no_others = FALSE, intersection = FALSE, message = ""){
if(is.data.frame(data)){
namesInData = names(data)
}else if(is.matrix(data)){
namesInData = dimnames(data)[[2]]
}else if(umx_is_MxData(data)){
namesInData = namez(data)
} else if (!typeof(data) == "character"){
namesInData = data
} else {
stop("Oops: Data has to be a dataframe or matrix. You gave me a ", typeof(data))
}
if(is.data.frame(namesNeeded)){
namesNeeded = names(namesNeeded)
}else if(is.matrix(namesNeeded)){
namesNeeded = dimnames(namesNeeded)[[2]]
} else if (typeof(namesNeeded)=="character"){
namesNeeded = namesNeeded
} else if (is.null(namesNeeded) || is.na(namesNeeded)){
return(TRUE)
} else{
stop("Oops: namesNeeded has to be a list of names, a dataframe or matrix. You gave me a ", typeof(namesNeeded), "\n",
"PS: names in data were: ", omxQuotes(namesInData))
}
if(intersection){
namesFound = intersect(namesNeeded, namesInData)
message(omxQuotes(namesFound))
} else {
namesFound = namesNeeded %in% namesInData
if(any(!namesFound)){
if(die){
stop("Not all required names were found in the data. Missing were:\n", paste(namesNeeded[!namesFound], collapse = "; "), "\n", message)
} else {
return(FALSE)
}
}
if(!is.null(illegal)){
illegalFound = illegal %in% namesInData
if(any(illegalFound)){
if(die){
stop("An illegal name was found in the data:\n", omxQuotes(illegal[illegalFound]), "\n", message)
} else {
return(FALSE)
}
}
}
if(no_others & !setequal(namesInData, namesNeeded)){
if(die){
stop("Data contains columns other than those needed. Superfluous columns were:\n",
paste(namesInData[!namesInData %in% namesNeeded], collapse = "; "))
} else {
return(FALSE)
}
} else {
return(TRUE)
}
}
}
#' Get variances from a df that might contain some non-numeric columns
#'
#' Pass in any dataframe and get variances despite some non-numeric columns.
#' Cells involving these non-numeric columns are set to ordVar (default = 1).
#'
#' @param df A dataframe of raw data from which to get variances.
#' @param format to return: options are c("full", "diag", "lower"). Defaults to full, but this is not implemented yet.
#' @param use Passed to [cov()] - defaults to "complete.obs" (see param default for other options).
#' @param ordVar The value to return at any ordinal columns (defaults to 1).
#' @param digits digits to round output to (Ignored if NULL). Set for easy printing.
#' @param strict Whether to allow non-ordered factors to be processed (default = FALSE (no)).
#' @param allowCorForFactorCovs When ordinal data are present, use heterochoric correlations in affected cells, in place of covariances.
#' @return - [OpenMx::mxModel()]
#' @export
#' @family Miscellaneous Stats Functions
#' @references - <https://tbates.github.io>
#' @md
#' @examples
#' tmp = mtcars[,1:4]
#' tmp$cyl = ordered(mtcars$cyl) # ordered factor
#' tmp$hp = ordered(mtcars$hp) # binary factor
#' umx_var(tmp, format = "diag", ordVar = 1, use = "pair")
#' tmp2 = tmp[, c(1, 3)]
#' umx_var(tmp2, format = "diag")
#' umx_var(tmp2, format = "full")
#'
#' data(myFADataRaw)
#' df = myFADataRaw[,c("z1", "z2", "z3")]
#' df$z1 = mxFactor(df$z1, levels = c(0, 1))
#' df$z2 = mxFactor(df$z2, levels = c(0, 1))
#' df$z3 = mxFactor(df$z3, levels = c(0, 1, 2))
#' umx_var(df, format = "diag")
#' umx_var(df, format = "full", allowCorForFactorCovs=TRUE)
#'
#' # Ordinal/continuous mix
#' data(twinData)
#' twinData= umx_scale_wide_twin_data(data=twinData,varsToScale="wt",sep= "")
#' # Cut BMI column to form ordinal obesity variables
#' obLevels = c('normal', 'overweight', 'obese')
#' cuts = quantile(twinData[, "bmi1"], probs = c(.5, .8), na.rm = TRUE)
#' twinData$obese1=cut(twinData$bmi1,breaks=c(-Inf,cuts,Inf),labels=obLevels)
#' twinData$obese2=cut(twinData$bmi2,breaks=c(-Inf,cuts,Inf),labels=obLevels)
#' # Make the ordinal variables into mxFactors
#' ordDVs = c("obese1", "obese2")
#' twinData[, ordDVs] = umxFactor(twinData[, ordDVs])
#' varStarts = umx_var(twinData[, c(ordDVs, "wt1", "wt2")],
#' format= "diag", ordVar = 1, use = "pairwise.complete.obs")
#'
umx_var <- function(df, format = c("full", "diag", "lower"), use = c("complete.obs", "pairwise.complete.obs", "everything", "all.obs", "na.or.complete"), ordVar = 1, digits = NULL, strict = TRUE, allowCorForFactorCovs= FALSE){
format = match.arg(format)
use = match.arg(use)